Influence of the symmetry energy on the giant monopole resonance of neutron-rich nuclei

We analyze the influence of the density dependence of the symmetry energy on the average excitation energy of the isoscalar giant monopole resonance (GMR) in stable and exotic neutron-rich nuclei by applying the relativistic extended Thomas-Fermi method in scaling and constrained calculations. For the effective nuclear interaction, we employ the relativistic mean field model supplemented by an isoscalar-isovector meson coupling that allows one to modify the density dependence of the symmetry energy without compromising the success of the model for binding energies and charge radii. The semiclassical estimates of the average energy of the GMR are known to be in good agreement with the results obtained in full RPA calculations. The present analysis is performed along the Pb and Zr isotopic chains. In the scaling calculations, the excitation energy is larger when the symmetry energy is softer. The same happens in the constrained calculations for nuclei with small and moderate neutron excess. However, for nuclei of large isospin the constrained excitation energy becomes smaller in models having a soft symmetry energy. This effect is mainly due to the presence of loosely-bound outer neutrons in these isotopes. A sharp increase of the estimated width of the resonance is found in largely neutron-rich isotopes, even for heavy nuclei, which is enhanced when the symmetry energy of the model is soft. The results indicate that at large neutron numbers the structure of the low-energy region of the GMR strength distribution changes considerably with the density dependence of the nuclear symmetry energy, which may be worthy of further characterization in RPA calculations of the response function.


I. INTRODUCTION
The departure of the equation of state (EOS) of neutron-rich matter (N > Z) from the symmetric limit (N = Z) is governed by the nuclear symmetry energy. This makes the knowledge of the symmetry energy important for a variety of areas of nuclear physics (masses, neutron densities, heavy ion collisions, etc.) and of astrophysics (supernovae, neutron stars, nucleosynthesis, etc.) [1][2][3][4][5][6]. However, the available experimental data on nuclei, which correspond mostly to isotopes with small and moderate neutron-proton asymmetries, do not constrain very precisely the value of the symmetry energy at the saturation density. Rather, it seems to be better constrained at a lower density around 0.10 fm −3 [7][8][9][10], corresponding to an average between the bulk and surface symmetry contributions in the semi-empirical mass formula [11][12][13]. All in all, the density dependence of the symmetry energy is a basic nuclear property that still is not sufficiently well known [1,4]. This situation has prompted a recent ongoing effort to characterize more closely the symmetry energy in the regime of subsaturation densities, relevant for finite nuclei, from the phenomenology of different nuclear observables [1,4,[11][12][13][14][15][16][17][18][19]. A proper understanding of the symmetry energy is particularly desirable in order to predict ground-state and excitation properties of atomic nuclei far from the naturally-occurring region of stability [20][21][22][23][24].
The study of collective excitations of nuclei near and away from the stability line is a salient source of information for nuclear structure [25][26][27][28][29]. For instance, the isoscalar giant monopole resonance (GMR) constrains the compression modulus K 0 of the EOS of symmetric nuclear matter. This resonance has been accurately measured in heavy and medium mass nuclei through recently improved α-particle scattering experiments [27,28]. By comparing the 1p-1h RPA results predicted by Gogny [30], Skyrme [31][32][33] and relativistic mean field (RMF) interactions [34,35] with the excitation energy measured in 208 Pb (14.17± 0.28 MeV [27]), it is concluded that the value of the nuclear incompressibility K 0 must lie in the range of 200-280 MeV. Performing measurements of the GMR in exotic nuclei may deliver other valuable informations and enhance our knowledge of the physics of isospin-rich and loosely-bound nuclear systems [26,29,32,36,37]. However, it represents an important experimental challenge due to e.g. the low intensities in the radioactive ion beams. With the advent of new technologies it is hoped that measurements of the breathing mode in short-lived unstable nuclei will become feasible in the future [26]. Note that promising first successes have already been demonstrated [29]. Although the value of the excitation energy of the GMR is driven to a large extent by the compression modulus of symmetric matter, it is known to be influenced also by the density dependence of the symmetry energy in a neutron-rich nucleus such as 208 Pb [33,[38][39][40][41][42]. This influence is expected to be amplified in isotopes with large neutron numbers, but how and to which extent has not actually been assessed.
Our purpose in this paper is to investigate the effects of the density dependence of the symmetry energy on the average excitation energy of the GMR in isospin-rich nuclei, of heavy and medium mass, when moving away from the stability valley. We ought to mention, however, that our aim here is not pursuing to constrain the nuclear EOS from the present calculations, but rather to analyze the effects that arise from the existing uncertainty in the knowledge of the density dependence of the symmetry energy. We advance that the calculations will show that the impact on the excitation energy of the GMR is somehow moderate for isotopes not far from stability, while the more notable effects are found for nuclei with large neutron numbers, thereby (unfortunately) rendering them more difficult to see in experiment. To describe the effective nuclear interaction, we found convenient to use the model of [8][9][10] that consists of a standard RMF Lagrangian plus an isoscalar-isovector nonlinear interaction between the ω and ρ meson fields. This coupling provides an efficient way to change in a controlled manner the density dependence of the symmetry energy without modifying the EOS of symmetric matter. In particular, the compression modulus K 0 is not modified. When the model is applied to finite nuclei, the variation of the strength of the isoscalar-isovector coupling changes the properties of the neutron density distribution (for example, the neutron rms radius), which have large experimental uncertainties, but the predicted ground-state binding energies and charge rms radii, which are accurately known in experiment, remain basically unaltered. Thus, this RMF model allows one to modify the density content of the symmetry energy while preserving the agreement with the accurate ground-state information on which the model has been calibrated previously (see [8-10, 41, 43]).
The consistent microscopic framework for the calculation of the properties of the GMR is furnished by the Hartree-Fock plus RPA theory, which corresponds to the small-amplitude limit of the time-dependent Hartree-Fock method. The generalization of this theory to deal with pairing in open-shell nuclei is known as quasiparticle RPA (QRPA). The QRPA provides the full information on the response function and centroid energies of the GMR. These microscopic calculations are manageable but rather demanding if they are fully selfconsistent, which is accentuated when they are performed using the canonical basis or in isotopes with large neutron numbers as it has been recently discussed in [37]. We will carry out our analysis by studying the estimates of the average excitation energy of the GMR that can be obtained without performing explicit RPA calculations but which carry meaningful physical information on the RPA strength. The typical example of this procedure is the sum-rule approach to the RPA developed in the nonrelativistic framework using Skyrme forces [44][45][46]. The sum-rule techniques provide average properties of the GMR such as the centroid energy and the resonance width and have been successfully applied to study the GMR in different nuclei and isotopic chains [30,33,37,[44][45][46][47][48]. In the sum-rule approach one obtains the cubic and inverse energy-weighted moments of the RPA strength function by means of, respectively, scaling and constrained calculations of the uncorrelated Hartree-Fock nuclear ground state. Note that these ground-state scaling and constrained values are the exact RPA moments [44][45][46]. Because highly-accurate ground-state calculations are simpler, the sum-rule techniques also have been used in the literature as a means to check the self-consistency of explicit RPA calculations of the strength function [32,37,47,48].
The validity of the sum-rule theorems when pairing is active (QRPA) has been proven in recent publications [37,49]. Capelli et al. [37] have made a careful study of constrained Hartree-Fock (CHF), constrained Hartree-Fock-Bogoliubov (CHFB) and QRPA calculations of the excitation energy of the GMR in the light isotopic chains of Ca and Ni. The calculations have reached up to the theoretical neutron drip line of the Skyrme SkM* force, i.e., up to the isotopes 76 Ca and 84 Ni. Admitting that these isotopic chains are known to terminate earlier in experiment than in the predictions of SkM* and of most of the nuclear effective forces, and knowing also that it will be difficult to measure the GMR in even less exotic isotopes, the calculations spanning the whole theoretical isotopic chains are very useful to test the predictions of nuclear forces under extreme isospin conditions (which can take place in astrophysical scenarios [50]), as well as to test the different many-body frameworks employed to deal with the problem. It may be noted also that Khan [51] has recently concluded that knowing the GMR in a single isotope such as 208 Pb may not suffice to extract K 0 of nuclear matter and that the whole isotopic chain should be measured in order to have a grasp on the various effects acting on the GMR [51].
From the results of [37] for Ca and Ni, it is seen that the inclusion of pairing (CHFB and QRPA calculations) does not modify significantly the behaviour shown by the GMR excitation energies obtained in the CHF calculations along the isotopic chains. In the case of the Ni chain, the atomic number not being as small as in Ca, the CHFB and CHF values agree fairly well for most isotopes from 56 Ni up to the neutron drip line included, as noticed in figure 5 of [37]. Of course, superfluidity is to be taken into account if one wants to describe in detail the excitation spectra or if precise values of the excitation energy within a few hundred keV in specific nuclei are needed for comparison with accurate measurements [37,51,52]. It has been pointed out in [37] that the QRPA calculations for isotopes having large neutron numbers become much time consuming and that ensuring proper numerics in such isospin-rich nuclei becomes somewhat difficult due to the weakly bound neutrons. Indeed, the authors of [37] state that in the case of large neutron numbers very accurate QRPA calculations may become almost prohibitive in keeping with the size of the basis.
The GMR is a collective excitation whose average properties in general vary smoothly with the mass number A and are rather insensitive to shell effects, at least for not too light isotopes. Thus, a description of the average properties of this resonance through semiclassical approaches like the extended Thomas-Fermi (ETF) method can be, for some purposes, a practical and interesting alternative. It has been applied in different studies in the nonrelativistic Skyrme-Hartree-Fock framework [46,47,[53][54][55]. The results have shown that the average excitation energies of the GMR obtained from scaling and constrained quantal HF calculations are nicely described by the corresponding semiclassical ETF counterparts. The ETF predictions of GMR excitation energies become increasingly closer to the quantal predictions with increasing atomic number, as it can be expected in semiclassical methods. Indeed, in the calculations along the Pb isotopic chain with the Skyrme force shown in figure 12 of [47], the ETF scaling and constrained excitation energies differ from the corresponding quantal values by less than 0.4 MeV and they reproduce quite accurately the change of the quantal results in going from the proton to the neutron drip line. Thus, from the experience in nonrelativistic calculations, we believe that the semiclassical scaling and constrained estimates of the average excitation energy will be a useful and realistic procedure to have a quick representation of the general trends of the GMR along isotopic chains also in the relativistic mean field framework. And in heavy systems such as the Pb isotopes, we expect that the semiclassical calculations of the GMR average excitation energy will yield reliable predictions practically on the quantitative level. Certainly, for a precise microscopic description, the more demanding QRPA calculations would have to be performed.
The semiclassical relativistic extended Thomas-Fermi (RETF) approach [56][57][58][59][60] has recently been applied in scaling and constrained calculations of the GMR energy in stable isotopes in [61,62]. The results turn out to be in good agreement with the values extracted from full relativistic RPA calculations [61,62]. Let us mention as a side comment that the RETF scaling procedure has been helpful in obtaining the virial theorem for the relativistic mean field model [61] and that it has allowed self-consistent calculations of the surface contribution K surf to the nuclear incompressibility in the relativistic theory [63]. In this paper we will use the relativistic nuclear mean field model with an isoscalar-isovector coupling mentioned previously together with the RETF approach for describing the average excitation energy of the GMR in the Pb and Zr isotopic chains, through scaling and constrained calculations. The present results are expected to provide a realistic indication on how the density dependence of the symmetry energy impacts on the behaviour of the excitation en-ergies and widths of the GMR along isotopic chains of heavy and medium mass. In the next section the basic theory is presented, in the third section the results are discussed and the summary and conclusions are laid in the fourth section. Some details about the semiclassical RETF calculations are given in the appendices.

II. FORMALISM
A detailed description of the scaling and constrained calculations of the excitation energy of the GMR in the relativistic framework using the RETF approach has been presented elsewhere [61,62]. Here, we shall restrict ourselves to outline the more relevant expressions, incorporating the new contributions that arise from the mixed ω-ρ meson interaction.
A. ETF energy density of the relativistic mean field model The basic model for the present study is the relativistic mean field theory of nucleons and mesons. As in the nonrelativistic problem, the theory can be formulated in the Thomas-Fermi approach. When corrections of order 2 (in gradients of the nucleon densities and of the Dirac effective mass) to the pure Thomas-Fermi contribution are included in the energy density, it is known as the relativistic extended Thomas-Fermi approach [56][57][58][59][60].
The local energy density that we employ in the present work can be written as follows: where ρ = ρ p +ρ n is the baryon density, ρ 3 = 1 2 (ρ p −ρ n ) is the isovector density, and the term E will be described immediately below. The energy density (1) contains scalar Φ = g s φ 0 (r), vector W = g v V 0 (r) and vector-isovector B = g ρ b 0 (r) meson fields, as well as the Coulomb field A = eA 0 (r). These fields represent, respectively, σ-meson, ω-meson, ρ-meson and photon exchange. The cubic κΦ 3 and quartic λΦ 4 self-interactions of the scalar field soften the EOS of symmetric nuclear matter around the saturation point. In addition, the model is supplemented by a nonlinear mixed coupling between the ω-meson and the ρ-meson fields. This Λ V B 2 W 2 interaction was introduced in [8] to soften the density dependence of the symmetry energy and of the EOS of neutron matter, which are two quantities that are predicted to be stiff by the majority of relativistic mean field models.
In the quantum approach, the contribution E to the RMF energy density in equation (1) is given by Here, the ϕ ν are single-particle Dirac wave functions, α and γ 0 are the usual Dirac matrices and m denotes the rest mass of the nucleons. The RETF representation of E consists of a pure Thomas-Fermi part E 0 plus a part E 2 , which is of order 2 compared to E 0 . The nucleon variables of the RETF energy density are the neutron and proton densities (ρ n and ρ p ) instead of the single-particle wave functions ϕ ν . Thereby, E is a functional of the nucleon densities and of the Dirac effective mass m * = m − Φ in the RETF approach. The gradient corrections contained in the E 2 term provide an improved description of the nuclear surface with respect to the Thomas-Fermi treatment with only the E 0 contribution. The reader can find the full RETF expression of the functional E in appendix A. We also provide there the variational Euler-Lagrange equations derived from equation (1) in the RETF approach for the nucleon densities and for the meson and photon fields.
The energy density H is to be integrated over space to compute the total energy of the nucleus. By means of partial integrations, and on account of the variational meson field equations given in appendix A, it is possible to write the relativistic energy density of a finite nucleus in the following form: Here, we have introduced the effective scalar density ρ eff s defined as The form given in equation (2) for H turns out to be more practical to perform the scaling of the equations in order to compute the excitation energy of the GMR in the scaling formalism.

B. Energy of the GMR in the scaling approach
The excitation energy of the isoscalar giant monopole resonance of a finite nucleus of mass number A is customarily written as where K A is called the compression modulus or incompressibility of the finite nucleus, and B M is called the mass or inertia parameter of the monopole oscillation. In the scaling approach, the incompressibility K A is obtained from the restoring force of the monopole vibration modeled through a scaling transformation of the nucleon densities: Here, α denotes the scaling parameter and q = n, p refers to neutrons or protons. The meson fields do not scale as simple powers of α because of the finite-range character of the meson interactions; we collect in appendix B the variational field equations for the scaled meson fields of the RMF model with the mixed ω-ρ interaction. Following [61,62], it is helpful to write the scaled Dirac effective mass m * because the energy density E and the scalar density ρ s of the RETF model scale, respectively, as E α (r) = α 4 E[ρ q (αr),m * (αr)] ≡ α 4Ẽ (αr) and as ρ sα (r) = α 3 ρ s [ρ q (αr),m * (αr)] ≡ α 3ρ s (αr). The tilded quantitiesẼ andρ s have the same expressions as E and ρ s given in appendix A if one replaces m * bym * in them. Taking into account the above, the scaled form of the energy density of the present model becomes with the definitionρ eff All densities and fields on the r.h.s. of (7) depend on the variable αr.
The calculation of the first and second derivatives of the scaled energy E(α) = dr H α (r) with respect to α proceeds along the same lines described in [61,62]. The final result for the scaling incompressibility is It is interesting to note that all of the nucleonic contributions to K scal A are summarized in the term with the scalar density and that the massless photon field and the quartic meson couplings such as λΦ 4 and Λ V B 2 W 2 do not provide explicit contributions. The expression (8) requires the calculation of the derivatives of the meson fields with respect to the scaling parameter α at equilibrium (α = 1). These derivatives can be computed by differentiation of the scaled meson field equations with respect to α and by solving them at α = 1, as we describe in appendix B. We will use the notation E scal M to refer to the excitation energy of the GMR calculated in the scaling approach, namely, The mass parameter of the monopole resonance is given by the expression [62,64,65] We note that in the nonrelativistic limit it reads where m is the nucleon rest mass. Because of the large contribution of the nucleon rest mass term mρ to the energy density H in the integrand of equation (10), the value of the excitation energy of the resonance is little modified by using either B M or B nr M in the calculations.

C. Energy of the GMR in the constrained approach
One also can study the average excitation energy of the giant monopole resonance through constrained calculations [44,62,[66][67][68][69]. For that purpose, one solves the problem associated with the constrained functional The densities, fields and energy obtained from the solution of the variational equations deduced from equation (12) depend on the value of the parameter η at which one performs the calculation. By expanding E(η) in a harmonic approximation about its minimum, located at the ground-state rms radius R 0 (corresponding to η = 0), one obtains the compression modulus of the finite nucleus in the constrained approach: Here, R 2 η = r 2 η is the square mass radius of the nucleus computed with the density ρ of equation (12). One can show that equation (13) can be rewritten in the following equivalent forms: We have used (14) to check the numerical accuracy of our calculations of K cons A . Finally, the excitation energy of the constrained isoscalar monopole vibration is computed as

D. Comparison with the expressions in the nonrelativistic sum-rule approach
We would like to point out that there is a parallelism between the scaling and constrained expressions for the GMR energy derived above for the relativistic model and the known results in the nonrelativistic formalism of sum rules. That is, one can define suitable average excitation energies of the resonance such as in terms of ratios of the integral moments m k = ∞ 0 dω ω k S(ω) of the RPA strength function S(ω). In the nonrelativistic RPA theory of the GMR [44][45][46], it has been demonstrated that the cubic energy-weighted moment m 3 = ∞ 0 dω ω 3 S(ω) of the RPA strength function is exactly equivalent to the second derivative of the energy obtained from a monopole scaling transformation. It is also known that the RPA energy-weighted moment ("Thouless theorem"), and that the exact RPA inverse energy-weighted moment m −1 = ∞ 0 dω S(ω)/ω can be obtained from a constrained calculation of the Hartree-Fock ground state: ("dielectric theorem"). Taking into account these results and the fact that m 1 = (2/m 2 ) B nr M , one sees an interesting analogy between E 3 and E 1 of the nonrelativistic RPA formalism of sum rules and E scal M of equations (9) and (8) and E cons M of equations (15) and (14) given previously for the GMR average excitation energy in the relativistic model, though to our knowledge the sum-rule approach has not been derived in the relativistic RPA theory. It is worth mentioning that the integral moments of the nonrelativistic RPA strength function obey certain inequalities [44][45][46]. In particular, one has where E x is the centroid energy of the GMR evaluated from the ratio of the m 1 moment to the m 0 moment. That is, E 1 is a lower limit for the mean energy (centroid) of the resonance, whereas E 3 is an upper limit. It is clear from the definition of the m 3 and m −1 moments that E 3 and E 1 explore different energy domains of the collective excitation, i.e., E 3 is more influenced by the high-energy region of the strength function S(ω), whereas E 1 is more sensitive to the low-energy components of S(ω). It can be shown [44] that the quantity provides an upper bound for the width of the giant monopole resonance in the RPA approximation. The calculations carried out for stable nuclei show that E 3 and E 1 become increasingly close with increasing nuclear mass. Therefore, in a stable heavy nucleus the value of the width σ is small and the monopole strength becomes concentrated around a single peak [44], which is the situation found in experiment. Following equation (20) and the above comments, in the discussion of our results the spread between the calculated E scal First, in this section we test with some examples the results of the relativistic extended Thomas-Fermi (RETF) approach vis-á-vis relativistic Hartree and RPA calculations. To this end, we report in table 1 semiclassical RETF results for the binding energy and the excitation energy of the GMR in the nuclei 208 Pb and 90 Zr. We also include in the comparison the value of the neutron skin thickness ∆R np = R n − R p (the difference between the neutron and proton rms radii of the nucleus) because this quantity is very sensitive to the density dependence of the symmetry energy [7,12,70,71]. The results have been computed with the NL3 mean field interaction [72] supplemented by the nonlinear Λ V coupling between the ω-meson and ρ-meson fields [8]. As pointed out in the Introduction, this model, through the isoscalar-isovector meson coupling, allows one to change the symmetry energy leaving unchanged the properties of the EOS of symmetric nuclear matter. By modifying the value of Λ V one can investigate the changes in finite nuclei due to the variation of the density dependence of the symmetry energy almost without affecting the previous quality of the model for binding energies and charge radii [8-10, 41, 43].
For the coupling constant Λ V we use the values 0.00, 0.01, 0.02, and 0.03 (models denoted as NL3.00, NL3.01, NL3.02 and NL3.03 in table 1). It is known that the binding energy of nuclei constrains an average between the symmetry energy at saturation and the surface symmetry energy, rather than the individual value of the symmetry energy at saturation  [8][9][10][11][12]. Thus, as done in earlier literature [8-10, 41, 43], we constrain the symmetry energy S(ρ) of the models to take a value of 25.67 MeV at a Fermi momentum k F = 1.15 fm −1 (i.e., at a subsaturation density ρ ≈ 0.10 fm −3 ). This implies that the coupling constant g ρ of the NL3 interaction has to be suitably adjusted for each value of the isoscalar-isovector Λ V coupling. As a consequence, the value S(ρ 0 ) of the symmetry energy at the saturation density ρ 0 changes with the Λ V parameter. We find the results S(ρ 0 ) = 37.4, 35.0, 33.2 and 31.7 MeV when we vary Λ V from 0.00 to 0.03. That is, the symmetry energy S(ρ) of the model around saturation becomes softer (it increases more slowly with the nuclear density) when the value of the Λ V parameter is larger. Further details on the behaviour of S(ρ) in the NL3 model supplemented by the Λ V coupling can be found in [10]. Table I shows that by increasing the Λ V coupling the neutron skin thickness of a neutronrich nucleus is significantly reduced, i.e., the models with a softer symmetry energy give rise to smaller neutron skins. Here, this is achieved while maintaining the proton rms radius of the nucleus essentially intact. For example, in the quantal Hartree calculation with the NL3.00 model we find R p = 5.466 fm for 208 Pb and 4.199 fm for 90 Zr, while the respective values using NL3.03 are 5.478 fm and 4.210 fm. The changes in R p are thus less than 0.5%. Moreover, one sees in Table I that also the binding energies are little altered by Λ V . Both 208 Pb and 90 Zr are slightly more bound at large Λ V but the changes in B/A are again no larger than 0.5%. Thus, the predictions for charge radii and binding energies of the original NL3 model are not seriously compromised by the additional Λ V coupling, which otherwise largely softens the symmetry energy.
From table 1 we observe that in the semiclassical RETF approach the binding energies of 208 Pb and 90 Zr are larger and the neutron skins smaller, than in the corresponding quantal Hartree calculations. In principle, this fact is not to be considered a lack of accuracy of the Thomas-Fermi-like methods because semiclassical and quantal calculations must differ by shell effects. It is well known that the semiclassical ground-state calculations by construction smooth out the quantal shell effects and provide the average part of the quantal binding energy [73][74][75]. Actually, in the semiclassical calculations one works with the particle densities and completely avoids to calculate single-particle wave functions. By this same reason, the semiclassical approach yields smooth average densities and energies that do not show the quantum shell oscillations. The nuclear ETF functionals with gradient corrections up to order 2 always produce some overbinding (between a ∼2% in heavy nuclei and a ∼10% in light nuclei) and yield smaller rms radii (between ∼1% and ∼4%) than the quantal mean field calculations [53,56,57,74]. The quantal corrections, which are more prominent in magic nuclei (our situation in table 1), can eventually be added perturbatively to the semiclassical binding energies using techniques such as the Strutinsky shell-correction method [75] or the expectation value method [56,57,74].
Our focus here is somehow more about the change of the results with the softness of the symmetry energy than on the absolute values themselves. In this respect, one has to note that the quantal and semiclassical results reported in table 1 show both practically the same variation as a function of the Λ V coupling. The binding energies remain similarly constant and the neutron skins decrease by almost the same amount of 0.08 fm in 208 Pb and of 0.03 fm in 90 Zr in both the Hartree and RETF calculations when Λ V is changed from 0.00 to 0.03. Therefore, the RETF approach predicts practically the same variation with Λ V as the quantal Hartree calculations. This implies that the rate of change of the results with the softness of the symmetry energy in the present model is well described by the RETF method. As a further test, which is motivated also because later we will study RETF results along isotopic chains, in figure 1 we plot for Pb isotopes the difference between the value of the neutron skin thickness of the nucleus calculated with Λ V = 0.00 and with Λ V = 0.03, as a function of the parameter I = (N − Z)/A. We observe by comparison with the corresponding quantal Hartree results, that the RETF approach performs quite satisfactorily along the Pb isotopic chain in predicting the right size of the modification induced by the change of the symmetry energy of the model.
With regard to the average excitation energies of the GMR, of Section II C) for 208 Pb and 90 Zr calculated in the RETF approach. We compare our RETF results with the values of the centroid energy m 1 /m 0 of the GMR extracted from the relativistic RPA (RRPA) strength, which have been reported in [41]. We see that for the 208 Pb nucleus, the RRPA centroid energies very nicely lie in between the RETF constrained and scaled estimates, as it is required in the nonrelativistic sum-rule approach to the RPA [44][45][46] (recall equation (19) of section II D). In the case of a lower-mass nucleus such as 90 Zr, the semiclassical excitation energies are not expected to be as successful as for the heavier 208 Pb. Indeed, we see in table 1 that the RETF constrained and scaled estimates in 90 Zr overestimate a little the RRPA centroid energies. The discrepancies are nevertheless moderate, with differences no larger than 2% for E cons M and 5% for E scal M with respect to the RRPA m 1 /m 0 value of 90 Zr. And, as it happens in the case of 208 Pb, the semiclassical average excitation energies of 90 Zr as a function of the Λ V coupling display the tendency of the RRPA excitation energies. In summary, we can conclude from the results discussed in the present section that the RETF approach is a reasonable starting point to study the influence of the density dependence of the symmetry energy on the GMR of neutron-rich nuclei.

B. GMR in heavy neutron-rich nuclei
As mentioned in the Introduction, we found in [47] within the Skyrme-Hartree-Fock framework that the semiclassical calculations are able to describe realistically the tendencies of the excitation energy of the GMR in regions away from the stability valley, especially for heavy nuclei such as the lead isotopes where a semiclassical approach performs better. In the present section we analyze with the RETF approach to the relativistic mean field model the effects of the density dependence of the symmetry energy on the excitation energy of the GMR in heavy neutron-rich nuclei by contrasting the results for the nucleus 266 Pb, which we take as an example of a heavy system with a large neutron excess away from stability, with the results for the stable isotope 208 Pb.
We display in figure 2 for 208 Pb and in figure 3 for 266 Pb the calculated values of the energy per particle E/A, the neutron skin thickness ∆R np = R n − R p , the compression modulus of the finite nucleus K A and the average excitation energy E x of the GMR as a function of the coupling constant Λ V in the relativistic model. We recall that the coupling Λ V induces a change of the density dependence of the symmetry energy. The model has a stiff symmetry energy at Λ V = 0.00, while it has a soft symmetry energy at Λ V = 0.03. The RMF Lagrangian supplemented by the nonlinear interaction between the ω and ρ meson fields leaves the binding energy of the nucleus almost independent of Λ V , at least for moderate values of this coupling. We note from the results for 266 Pb in figure 3 that this happens not only in the case of stable nuclei, which we have seen previously for 208 Pb and 90 Zr in table 1, but also in the case of much more neutron-rich nuclei. The neutron thickness ∆R np decreases almost linearly with an increase of Λ V in both 208 Pb and 266 Pb. The reduction of the ∆R np value at Λ V = 0.03 amounts to about a 30% of the original neutron skin thickness at Λ V = 0.00. The proton rms radius R p of 208 Pb and 266 Pb has been found to remain nearly constant with Λ V , and therefore the decrease of ∆R np is due almost exclusively to the decrease of R n . It is worth mentioning that the alluded variation of the neutron skin thickness of 208 Pb in the present calculations is close to roughly a 1% change in the neutron rms radius R n of this nucleus. That is, the range of variation we have allowed in the density dependence of the symmetry energy by changing Λ V between 0.00 and 0.03 induces a change in R n of 208 Pb around 1%. In this sense, it is consistent with the uncertainty that will remain in this observable after the PREX experiment at the Jefferson Laboratory [76], which aims to measure R n of 208 Pb to 1% accuracy via parity-violating electron scattering.
The average excitation energy of the GMR calculated in the RETF approach involves the knowledge of the finite nucleus compression modulus K A , or restoring force of the resonance. We display K A in figures 2(c) and 3(c) for 208 Pb and 266 Pb, respectively. Let us discuss first the change of K A with the mass number A in isotopes. We see from the results for 208 Pb and 266 Pb that the finite nucleus compression modulus decreases when A increases. This behaviour may be easily understood from the so-called leptodermous expansion of K A [77]: where K 0 is the compression modulus of symmetric infinite nuclear matter and K surf , K τ and K Coul are the leading surface, volume-symmetry and Coulomb corrections. The coefficients K surf , K τ , and K Coul in the scaling model are negative, and in particular K τ may be large [33,63,[77][78][79][80]. Within an isotopic chain, the nuclear charge Z is fixed and therefore the neutron-proton asymmetry I = (N − Z)/A increases when A increases. Then, although the surface and Coulomb contributions in (21) become a little less negative with increasing A, the increase of the negative volume-symmetry contribution K τ I 2 dominates. Therefore, K A in isotopes decreases with larger A. We stress that the expansion (21) is provided here for indicative purposes only and that in all results we have computed K A self-consistently.
The scaling value K scal A of the finite nucleus incompressibility, given by equation (8) of section II B, is displayed in figures 2(c) and 3(c) by the upper curves with solid symbols. For both isotopes 208 Pb and 266 Pb, K scal A shows an upward trend with larger Λ V . This variation of K scal A with Λ V is due to the softening of the density dependence of the symmetry energy. It modifies the value of K τ and consequently K A (of course, K τ does not describe alone the whole effect as there are surface-symmetry and other corrections neglected in (21) [63]). The present models have all K 0 = 271.5 MeV, while the K τ coefficient changes from −699 MeV in NL3.00 to −381 MeV in NL3.03. The recent empirical extractions of K τ suggest K τ values in this range. For example, isospin diffusion in nuclear reactions favours the constraint K τ = −500 ± 50 MeV [4], the breathing mode of Sn isotopes suggests K τ = −550 ± 100 MeV [28] (or −500 ± 50 MeV [79]) and neutron skins of nuclei point to K τ = −500 +125 −100 MeV [12].
The excitation energy E scal M of the monopole oscillation evaluated in the scaling approach is determined by the ratio of the restoring force K scal A to the mass parameter B M (section II). As indicated, to replace B M by its nonrelativistic limit B nr M is a very good approximation. In practice, this implies that B M varies with the parameters of the nuclear interaction essentially in the same way as r 2 of the nucleus does, see equation (11). When Λ V is increased, r 2 of the matter distribution of the nucleus decreases because of the reduction of the neutron rms radius originated by the softer symmetry energy. As a consequence, the mass parameter of the GMR (analogously, the m 1 sum rule (17) in the nonrelativistic language) also decreases. This effect combined with the growing tendency of K scal A with Λ V implies that the scaling excitation energy of the GMR will be enhanced with a softer symmetry energy, which is confirmed in figures 2(d) and 3(d) where E scal M is plotted. We recall that this enhancement of E scal M is due to the different density dependence of the symmetry energy alone because the incompressibility K 0 of symmetric matter is forced to be the same in all the interactions considered here.
Let us now analyze the results obtained through the constrained calculations described in section II C. For the stable nucleus 208 Pb, we see in figure 2 (20)), and therefore this quantity is basically constant with the symmetry energy in 208 Pb. In contrast to this situation in the stable 208 Pb nucleus, in the more neutron-rich isotope 266 Pb we see in figure 3 that the change of K cons  (15) by approximating B M with its nonrelativistic form B nr M = m A r 2 , and to write K cons A using the first expression given in equation (14). This results in which in terms of the defined quantities m 1 and m −1 has the same form that E cons M takes in the nonrelativistic RPA formalism of sum rules (i.e., E 1 of equation (16)) [44][45][46]. It is to be noted that the denominator of (22) is very sensitive to the presence of loosely bound nucleons in the nucleus [47], which is precisely the situation that occurs in 266 Pb. When one is approaching the neutron drip line, the negative chemical potential µ n of neutrons becomes close to vanishing, implying that the system contains weakly bound neutrons. These neutrons are very soft against compression and in the constrained calculations they give rise to a large increase of the absolute value of (∂R 2 η /∂η) η=0 , i.e., of m −1 in (22). This effect has been found in [47] in constrained HF (and ETF) calculations for nonrelativistic Skyrme forces, where it can be related to the enhancement of the RPA strength in the lowenergy region produced by the increase of the transitions to the continuum at large neutron numbers. The same effect is found in the recently published QRPA and CHFB Skyrme calculations, which take account of the pairing correlations, performed in the light isotopic chains of Ca and Ni [37]. In our present calculations, the increase of m −1 at large neutron numbers is accentuated the softer is the symmetry energy (note that parameter sets having a softer symmetry energy have also weaker binding energy of the last bound neutrons [10] and therefore reinforce the effect). This explains the decreasing trend of E cons M with Λ V in 266 Pb noticed in figure 3. The decrease of K cons A with Λ V seen in the same figure is similarly understood by recalling equation (14) that relates K cons A with (∂R 2 η /∂η) −1 η=0 .

C. GMR along the Pb and Zr isotopic chains
In this section we would like to investigate at which values of the neutron excess of the nucleus the effects coming from the symmetry energy become more relevant in the excitation energies of the GMR in isotopic chains of heavy (Pb) and medium (Zr) mass. Within the Skyrme-Hartree-Fock framework, in [47] we analyzed the excitation energy and the resonance width of the GMR in some isotopic chains with the sum-rule approach (i.e., through scaling and constrained calculations). In the case of a small atomic number such as the Ca isotopes, the semiclassical ETF results were able to describe the global variation of the quantal HF excitation energies but could not reproduce them quantitatively. However, the ETF results for the GMR excitation energies become progressively more accurate with larger atomic number. In the case of the Pb chain, the ETF scaling and constrained excitation energies were found to be in very good agreement with the quantal HF values on the quantitative level from the proton to the neutron drip line (the differences between ETF and quantal values not exceeding 0.4 MeV) [47]. Thus, from the experience in the nonrelativistic calculations, as well as suggested by the comparison of RETF and RRPA results made in table 1, the present semiclassical predictions for the GMR are expected to be sufficiently realistic, and rather accurate quantitatively (compared to quantal calculations) in the case of the Pb isotopes.
We depict in figure 4 the scaling and constrained average excitation energies of the GMR for the Pb isotopic chain. Here we only use the values 0.00 and 0.03 for the nonlinear ω-ρ interaction Λ V , representing stiff and soft density dependences of the symmetry energy. In the semiclassical calculations, the vanishing of the neutron chemical potential µ n defines the location of the neutron drip line. We see in figure 4 that the latter is reached at a smaller mass number when the model has a soft symmetry energy, in agreement with the earlier observations of [10]. Note that shell effects are absent in these semiclassical predictions of the drip line (a quantal Hartree calculation with NL3 predicts the Pb neutron drip line at A = 266). The proton drip line occurs at neutron number N not much larger than Z. Therefore, because one deals with almost symmetric matter, the position of the proton drip line is unaffected by the density dependence of the symmetry energy.
The calculated GMR excitation energies show in figure 4 a downward tendency with increasing mass number A when one proceeds along the isotopic chain, by similar reasons to those pointed out in our previous discussion about the 208 Pb and 266 Pb isotopes. The downward trend is more noticeable for E cons M than for E scal M in the isotopes with larger neutron numbers. The scaling excitation energies E scal M obtained using a parameter set with a soft symmetry energy (Λ V = 0.03) are systematically higher in the whole isotopic chain than if the symmetry energy is stiff (Λ V = 0.00). Similarly to the scaling results discussed for 208 Pb and 266 Pb, this behaviour of E scal M is a consequence of the fact that the matter radius of the nucleus is smaller and the scaling incompressibility K scal A larger for a soft symmetry energy. At variance with the uniformity displayed by the E scal M results, the constrained estimate E cons M of the GMR excitation energy has a distinct feature with the symmetry energy below and above A ≃ 254. If the mass number is below A ≃ 254, E cons M is larger when the parameter set has a soft symmetry energy. However, this tendency is reversed from A ≃ 254 to the neutron drip line when the isotopes are more neutron rich. Actually, in a given nucleus the Λ V parameter reduces its neutron rms radius, which means that the system becomes on average more compact. As long as the nucleus does not be very neutron rich, if it is more compact, the effect of the constraining parameter η on the nuclear rms radius is less, and therefore the absolute value of (∂R 2 η /∂η) η=0 in the denominator of equation (22) for  As a result of their specific contribution to R 2 η when one applies an external constraint, the absolute value of (∂R 2 η /∂η) η=0 increases sharply instead of decreasing. This originates the large reduction of E cons M seen in figure 4 for the more neutron-rich isotopes. The reduction of E cons M at large neutron numbers is enhanced when the symmetry energy is soft compared to the case of a stiff symmetry energy because the outer neutrons are more weakly bound in the former case. A comparable decrease of the excitation energy of the GMR, however, is not present in the scaling calculations because the discussed effect does not arise in a self-similar scaling of the density of the nucleus. As we mentioned in section II D, the scaling and constrained average energies pertain to different energy regimes of the GMR. While E scal M informs more about the relatively high-energy contributions to the RPA strength, E cons M carries more information about the low-energy sector which is the region that shows the more noticeable changes in largely neutron-rich nuclei because of the presence of weakly bound neutrons.
We may take σ should display a sharp increase for the Pb isotopes with large neutron numbers. We see from figure 4 that this effect is predicted to be stronger when the nuclear interaction has a soft symmetry energy. As discussed in [47] in nonrelativistic calculations with Skyrme forces, the value of the estimate σ of the GMR width provides a qualitative idea about the distribution of the RPA strength. A small width σ indicates that the RPA strength of the GMR is basically concentrated in a narrow single peak, while a large width σ suggests that the peak is broad, or even that the RPA strength may be fragmented into several peaks. The large value of m −1 ∝ −(∂R 2 η /∂η) η=0 (see equation (22)) that one finds in the constrained calculations for the more neutron-rich Pb isotopes, points towards an enhancement of the RPA strength in the low-energy region in these nuclei, due to transitions from the last weakly bound single-particle levels to the continuum [47]. (This effect, which has a quantal origin, is accounted at least on average by the semiclassical calculations [47].) The discussed scenario is actually confirmed by the QRPA calculations of [37] where the fragmentation as well as the enhancement of the low-energy part of the RPA strength in Ca and Ni isotopes approaching the neutron drip line can be seen in figures 3 and 6.
Lastly, we present in figure 5 the calculated results for the excitation energies of the GMR in the Zr isotopic chain. The variation of the excitation energies obtained in the scaling scheme as a function of the mass number A and of the Λ V coupling is qualitatively similar to the pattern exhibited by the Pb isotopes. That is, the excitation energies E scal M show a monotonic decrease with increasing mass number, and they are larger when calculated with a soft symmetry energy. The relative differences between the values of E scal M at Λ V = 0.03 and at Λ V = 0.00 are smaller compared to the chain of the Pb isotopes because of the smaller isospin content of the Zr isotopes. Concerning the results for the constrained excitation energies in Zr, the value of E cons M shows a stronger dependence on the mass number A at large neutron numbers than the scaling value E scal M . As long as A 120, the excitation energy E cons M for zirconium is to a large extent independent of the symmetry energy. But afterwards (A 120), E cons M becomes visibly lower for the soft symmetry energy. In consonance with the behaviour of E scal M and E cons M found in figure 5, it is predicted that the width σ of the GMR for the Zr nuclei should remain practically constant till A ≃ 100. The width σ would start increasing moderately between A ≃ 100 and A ≃ 120, because of the slight gradual separation between E scal M and E cons M in these isotopes. Finally, the resonance width would be much enhanced in the region of the most neutron-rich Zr isotopes, where the effect is largely reinforced by having a soft symmetry energy in the nuclear interaction.

IV. SUMMARY AND CONCLUSIONS
Though the excitation energy of the isoscalar giant monopole resonance is mainly ruled by the compression modulus of symmetric nuclear matter, it is also affected by the symmetry energy in neutron-rich nuclei. In this paper we have analyzed the influence of the density dependence of the symmetry energy on the average excitation energy of the GMR. In order to do this study, we have used a relativistic mean field model supplemented by an additional ω-ρ nonlinear interaction driven by a coupling constant Λ V . This additional coupling allows one to generate different parameter sets that preserve the properties of the symmetric matter and differ among them, only, in the density dependence of the symmetry energy, which is softened by increasing Λ V . As a function of the coupling constant Λ V , the model predicts in finite nuclei almost the same binding energy and proton rms radius as in the case where Λ V = 0, while the neutron rms radius gets significantly reduced.
The excitation energy of the isoscalar giant monopole resonance has been computed using the semiclassical relativistic extended Thomas-Fermi theory. We obtain the energy of the breathing mode in two different ways: through the scaling model and performing constrained calculations. It is known that these calculations with nonrelativistic Skyrme forces nicely reproduce the average energies E 3 and E 1 obtained from the more microscopic RPA calculations, not only for stable nuclei but also for largely neutron-rich nuclei [47]. In our study, we have calculated the scaling and constrained semiclassical average energies with the relativistic model along the isotopic chains of Pb and Zr. We can also obtain some qualitative information about the distribution of the RPA strength from the resonance width computed using these two semiclassical estimates of the excitation energy.
As it is known, the excitation energy of the breathing mode decreases, at least for stable nuclei, following roughly an A −1/3 law [77] due to the geometrical enhancement of the mass denominator. However, when the neutron-proton asymmetry increases, departures from this behaviour appear. The reason is that the isospin-dependent part of the finite nucleus incompressibility K A plays a non-negligible role by reducing the restoring force and making nuclei with large neutron numbers softer than the stable isotopes. This effect in the average excitation energy is present in the scaling estimate but it is particularly strong in the constrained calculation. When the neutron-proton asymmetry increases, the outer neutrons become more and more weakly bound starting to behave almost as extremely asymmetric nuclear matter, which is very soft compared to symmetric matter [81]. This effect makes the compression modulus K cons A smaller and, therefore, the average excitation energy of the GMR decreases. As discussed in the case of constrained HF calculations with Skyrme forces [47], this scenario reflects the important contribution to the low-energy part of the RPA strength distribution coming from the particle-hole transitions to the continuum from the weakly bound neutron energy levels in nuclei having large neutron numbers. A similar situation has been found recently by Capelli, Colò and Li [37] in QRPA and constrained HFB Skyrme calculations performed along the light isotopic chains of Ca and Ni.
For a given nucleus the excitation energy of the GMR also depends on the density content of the symmetry energy of the model used to compute it. In stable nuclei the scaled and constrained estimates of the excitation energy are larger for parameter sets that have a symmetry energy with a soft density dependence. In the case of the constrained estimate, this tendency is reversed in the nuclei with large neutron-proton asymmetry. The reason lies in the fact that in these nuclei the outermost and weakly bound neutrons are very compressible and consequently the constrained compression modulus K cons A becomes lower when the density dependence of the symmetry energy is soft (the outer neutrons are then less bound) than when the symmetry energy is stiff. Thus, the average excitation energies based on the inverse energy-weighted sum rule are predicted to be significantly reduced for large neutron numbers when computed with a parameter set that has a softer density dependence of the symmetry energy.
In conclusion, the density dependence of the symmetry energy has a moderate influence on the value of the excitation energy of the isoscalar giant monopole resonance in stable nuclei and in nuclei that are not far from stability. It has a larger impact on the GMR of isotopes with large isospins, where the calculations reveal that the constrained energies become much reduced in models having a soft symmetry energy. The results point towards a considerable enhancement of the low-energy part of the RPA strength and to likely different structures in the excitation spectra of the GMR in unstable neutron-rich nuclei if the nuclear symmetry energy is soft. It may be worth pursuing a detailed characterization of these predictions related with the softness of the symmetry energy by computing the strength distribution in a few selected isotopes at large neutron numbers with the more demanding but more microscopic QRPA theory of the GMR response. Unfortunately, it is difficult that the indicated effects be tested in the laboratory as they become prominent at neutron numbers beyond present reach. Such nuclear systems are otherwise expected to exist in astrophysical conditions, such as those prevailing in the inner crust of neutron stars, where exotic nuclear clusters can be stabilized by an enveloping gas of neutrons and electrons [50]. It is thought that low-energy collective excitations play an important role in astrophysical processes that involve nuclei far from stability [82], and in that context the modification of the low-energy region of the strength distribution of the monopole resonance by the softness of the symmetry energy could have implications.

Appendix A: Equations of the relativistic extended Thomas-Fermi model
We have introduced the energy density of the relativistic nuclear mean field model supplemented by an isoscalar-isovector mixed ω-ρ interaction in equation (1) of the text. Following [56][57][58][59], in the relativistic extended Thomas-Fermi approach the expression of the nucleonic part E is written as E = E 0 + E 2 , where and E 2 = q=n,p C 1q (k Fq , m * )(∇ρ q ) 2 + C 2q (k Fq , m * ) (∇ρ q · ∇m * ) + C 3q (k Fq , m * )(∇m * ) 2 . (A2) The leading term E 0 is the contribution of the pure Thomas-Fermi part, while the term E 2 contains the gradient corrections of order 2 . The variables k Fq , m * , and ǫ Fq are, respectively, the local Fermi momentum k Fq = (3π 2 ρ q ) 1/3 , the Dirac effective mass m * = m − Φ, and the dispersion relation ǫ Fq = (k 2 Fq + m * 2 ) 1/2 of the effective particle, with q = n for neutrons and q = p for protons. The coefficients C iq stand for the following analytic functions of k Fq and m * = m − Φ: The RETF ground-state densities and meson fields are calculated by solving by numerical iteration the Euler-Lagrange equations associated to the energy density (1), with the constraint of baryon number conservation (which is imposed through Lagrange multipliers µ q , the chemical potentials of neutrons and protons): We have τ 3q = +1 for protons and τ 3q = −1 for neutrons. The baryon density is ρ = ρ p + ρ n , while ρ 3 = 1 2 (ρ p − ρ n ) is the isovector density of the nucleus. Finally, the RETF expression of the scalar density ρ s entering equation (A5) is given by In the case of the scalar field one has to take account of the terms that arise due to the fact that the scalar density ρ s is itself a function of the scalar field. Following the same steps as above, from the Klein-Gordon equation for the scaled scalar field Φ α , one readily arrives at the equation whose solution at α = 1 provides the value of ∂Φ α /∂α| α=1 .