Neutron skin thickness in droplet model with surface width dependence: indications of softness of the nuclear symmetry energy

We analyze the neutron skin thickness in finite nuclei with the droplet model and effective nuclear interactions. The ratio of the bulk symmetry energy J to the so-called surface stiffness coefficient Q has in the droplet model a prominent role in driving the size of neutron skins. We present a correlation between the density derivative of the nuclear symmetry energy at saturation and the J/Q ratio. We emphasize the role of the surface widths of the neutron and proton density profiles in the calculation of the neutron skin thickness when one uses realistic mean-field effective interactions. Next, taking as experimental baseline the neutron skin sizes measured in 26 antiprotonic atoms along the mass table, we explore constraints arising from neutron skins on the value of the J/Q ratio. The results favor a relatively soft symmetry energy at subsaturation densities. Our predictions are compared with the recent constraints derived from other experimental observables. Though the various extractions predict different ranges of values, one finds a narrow window L ~ 45-75 MeV for the coefficient L that characterizes the density derivative of the symmetry energy which is compatible with all the different empirical indications.


I. INTRODUCTION
Neutron skin thickness is the name in common usage to refer to the difference between the root-mean-square (rms) radii of the neutron and proton density distributions of atomic nuclei: ∆R np = r 2 1/2 n − r 2 1/2 p .
Experimentally, the value of the proton rms radius r 2 1/2 p is obtained from the charge radius. The latter has been measured by electron-nucleus elastic scattering with high accuracy (often the accuracy in charge radii is better than 1% [1]). In contrast, our knowledge of the neutron distribution in nuclei and of its rms radius r 2 1/2 n , as well as our knowledge of ∆R np , is till date less precise. This situation looks inadequate in anticipation of the next generation of rare ion accelerator facilities that are planned for the synthesis and study of exotic nuclei, as in RIKEN centers (Japan), in FAIR at GSI (Germany), in the CSR at HIRFL (China), or in FRIB at MSU (USA). Without securing sufficient knowledge of the neutron distribution in stable nuclei, the prospects of nuclear structure theory in this thriving domain may be compromised. It is expected that parity-violating electron scattering will provide in the nearby future a leap forward in the quest for high precision determinations of the neutron radius in heavy nuclei [2].
The calibration of the neutron skin thickness of nuclei also is one of the problems at the forefront of nuclear structure by reason of being intimately correlated with the nuclear symmetry energy. Indeed, the symmetry energy is a fundamental quantity in nuclear physics and astrophysics, because it governs at the same time important properties of very small entities like atomic nuclei and of very large objects like neutron stars [3]. One of the crucial properties of the symmetry energy, which still is not sufficiently well constrained, is its dependence on the nuclear density. It is relevant in the astrophysical context in order to understand a wealth of phenomena [3][4][5], including supernova explosions, neutrino emission, and the cooling mechanism of protoneutron stars, as well as mass-radius relations in neutron stars. Moreover, the density content of the symmetry energy eventually relates to basic issues in physics. This is the case of precision tests of the Standard Model through atomic parity non-conservation observables [6], and even of studies on constraining a possible time variation of the gravitational constant [7]. In terrestrial laboratories, the available tools to delineate the density dependence of the symmetry energy at saturation and subsaturation densities, include the interaction potential between neutron-rich nuclei [8], observables like isospin diffusion and isoscaling in heavy-ion reactions at intermediate energies [9][10][11][12][13][14][15][16][17][18][19][20][21], different modes of collective excitations of nuclei [22][23][24][25][26], and, of course, data on the binding and structure of neutron-rich nuclei and on their neutron skin thickness.
In the literature there exist several theoretical formulations to investigate the neutron skin thickness of neutronrich nuclei and its connections with the symmetry energy. This is the case, for instance, of methods based on the droplet model [27,28], on the concept of surface symmetry energy [5,29,30], thermodynamical arguments [31], nucleonic density form factors [32], mean-field analyses [33][34][35], or studies in the spirit of the Landau-Migdal approximation [36]. It has been shown that the neutron skin thickness in heavy nuclei, like 208 Pb, calculated in mean-field models with either non-relativistic or relativistic effective nuclear interactions, displays a linear correlation with the slope of the neutron equation of state (EOS) obtained with the same interactions at a neutron density ρ ≈ 0.10 fm −3 [37][38][39]. A similar correlation exists between ∆R np and the density derivative of the bulk symmetry energy [13-16, 33, 40, 41], as the latter is a measure of the pressure difference between neutrons and protons. These correlations have been exploited in recent years to gain a deeper understanding of the isospin properties of the effective nuclear interaction and to relate them with nuclear and astrophysical observations.
The rms radius of neutron densities in nuclei has been measured with hadronic probes such as proton-nucleus elastic scattering [42][43][44][45] or inelastic scattering excitation of the giant dipole and spin-dipole resonances [46,47]. On the other hand, antiprotonic atoms are helpful to probe the size of the neutron skin of nuclei from the fact that the nuclear periphery is very sensitive to antiprotons in the normally electronic shell. Experimentalists combine two different techniques in this case [48][49][50], namely, the measurement of the antiprotonic X-rays which determine the atomic level shifts and widths due to the strong interaction, and the radiochemical analysis of the yields after the antiproton annihilation. The values of the neutron skin thickness of 26 stable nuclei from 40 Ca to 238 U deduced from antiprotonic atoms data by Trzcińska et al. [48,49] follow a roughly linear trend with the overall relative neutron excess I = (N − Z)/A of these nuclei. This trend can be fitted by the relationship ∆R np = (0.90 ± 0.15)I + (−0.03 ± 0.02) fm as discussed in Refs. [48,49]. As mentioned, all neutron skin thickness measurements have relatively large uncertainties in comparison with charge radii, and sometimes the results from different experimental techniques are not totally consistent among them [47,49]. The neutron skin sizes determined in Refs. [48,49] from the analysis of antiprotonic atoms are till date the largest set of uniformly measured values of ∆R np all over the periodic table (40 ≤ A ≤ 238). Due to this reason, we shall use hereinafter these data as the experimental benchmark for our calculations.
The droplet model (DM) describes in a physically transparent way nuclear radii and relates them directly with basic properties of the nuclear interactions. In the present paper we study the neutron skin thickness of atomic nuclei with the DM using various effective nuclear interactions of the Skyrme, Gogny, and relativistic mean-field (RMF) type. The present work extends with a new analysis and perspective a first presentation of our study made in Ref. [51]. Here, we will show that the ratio of the DM parameters J and Q, which drives the value of the neutron skin thickness in heavy nuclei, is correlated with the slopes in density of the nuclear symmetry energy and of the EOS of neutron matter. We compare the DM values for the neutron skin thickness with the results obtained in self-consistent extended Thomas-Fermi (ETF) calculations of finite nuclei [54,58,62,63], since both methods are free of shell effects. A non-negligible role of the contribution of the difference in the surface widths of the neutron and proton density profiles is noticed. Next, we use the experimental neutron skin thickness measured in antiprotonic atoms to explore the range of possible values of the ratio J/Q that are favored by neutron skins. With these values we can predict some properties of the density dependence of the nuclear symmetry energy. Our results are compared with the constraints recently obtained in the literature using other observables and methods.
The present paper is arranged as follows. In the second section we study the neutron skin thickness of heavy nuclei on the basis of the DM [27,52,53] and show a correlation that links the value of the ratio J/Q, which governs the neutron skin thickness of nuclei, with the slope of the symmetry energy in bulk matter at saturation. In the third section, the contribution of the surface widths of the neutron and proton density distributions to the neutron skin thickness is analyzed with the DM using nonrelativistic and covariant mean-field nuclear interactions. In the fourth section we estimate possible constraints on the density dependence of the nuclear symmetry energy on the basis of the DM and the experimental data on neutron skin sizes derived from antiprotonic atoms. We discuss the present results in comparison with the recent constraints obtained from various observables and methods. Finally, the summary and our conclusions are laid in the fifth section. We outline the procedure for the calculation of the Q coefficient in the Appendix.

II. THE FRAMEWORK
A. Neutron skin thickness in the droplet model In the DM of average nuclear properties [52][53][54] the neutron skin thickness of a finite nucleus is computed from the expression [27,52] where e 2 Z/70J is a correction due to the Coulomb interaction, R = r 0 A 1/3 is the nuclear radius, and b n and b p are the surface widths of the neutron and proton density profiles. In the "standard" version of the DM it is assumed that b n = b p = 1 fm [27,28,52], which implies a vanishing surface width correction to the neutron skin thickness. The quantity t in Eq. (2) represents the distance between the neutron and proton mean surface locations.  This distance is computed as [27,52] where I = (N − Z)/A, J is the symmetry energy coefficient at saturation, Q is the surface stiffness coefficient, and c 1 = 3e 2 /5r 0 . The coefficient J represents with a very good accuracy the energy cost per nucleon to convert all protons into neutrons in symmetric infinite nuclear matter at saturation density ρ 0 . The surface stiffness coefficient Q measures the resistance of the system against separation of neutrons from protons to form a neutron skin. To extract Q from an effective nuclear interaction requires to perform calculations of asymmetric semi-infinite nuclear matter (ASINM). Therefore, the calculated value of Q may depend somewhat on the type of approach, such as the Hartree-Fock or Thomas-Fermi methods, employed to describe the nuclear surface [30,[54][55][56][57][58]. From Eq. (3), one sees that the leading contribution to t in large nuclei is the term 3 2 r 0 (J/Q)I. Thus, the DM suggests that one can expect a correlation between ∆R np and J/Q in heavy nuclei. We illustrate this fact in the left panel of Fig. 1. We depict there the neutron skin thickness of 208 Pb obtained from self-consistent quantal calculations with the Skyrme and Gogny Hartree-Fock methods as well as with the RMF Hartree approach. The results are shown as a function of the value of the J/Q ratio for various mean-field effective interactions. The J values of the selected interactions (about 27-32 MeV in the non-relativistic forces and about 32-45 MeV in the covariant forces, see Table I) cover widely the plausible physical range of the bulk symmetry energy. The values of Q used in this work have been extracted from ASINM calculations performed in the extended Thomas-Fermi (ETF) approach as described in the Appendix (see also Ref. [58]). Even if the shell effects, present in the mean-field calculations of ∆R np , are not built in the DM [27], one observes a considerably linear correlation between the values of ∆R np and J/Q. It should be pointed out that while all the effective interactions have been accurately calibrated to data on binding energies and charge radii, and describe these properties very successfully, they predict widely different values of ∆R np , as we see in the present case of 208 Pb. This underlines the fact that the isospin sector of the effective interactions is little constrained.

B. Properties of the nuclear symmetry energy
Let us consider the energy per particle e(ρ, δ) in asymmetric infinite nuclear matter of total density ρ = ρ n +ρ p and relative neutron excess δ = (ρ n − ρ p )/ρ, where ρ n and ρ p stand for the neutron and proton densities, respectively. The general expression e(ρ, δ) = e(ρ, δ = 0) + c sym (ρ) defines the symmetry energy coefficient c sym (ρ) of a nuclear EOS at the density ρ. This expression is particularly useful because c sym (ρ) dominates the corrections to the symmetric limit for all values of δ, especially at the subsaturation densities of relevance for finite nuclei [59]. Actually, c sym (ρ) provides with excellent accuracy the difference between the binding energies of pure neutron matter (δ = 1) and symmetric matter (δ = 0). It is customary, and insightful, to characterize the behavior of an EOS around the saturation density ρ 0 by means of a few bulk parameters calculated at the saturation point, as in the formula [10, 13-16, 52, 59-61] where ǫ = (ρ 0 − ρ)/3ρ 0 expresses the relative density displacement from ρ 0 . Here, the quantities a v and K v denote the energy per particle and the incompressibility modulus of symmetric nuclear matter. One has c sym (ρ 0 ) = J. The DM coefficients L and K sym are, respectively, proportional to the slope and the curvature of the symmetry energy coefficient c sym (ρ) at saturation density: . (6) The quadratic expansion c sym (ρ) ≈ J − Lǫ + 1 2 K sym ǫ 2 in Eq. (5) is often a reliable representation of the actual value of the c sym (ρ) coefficient at densities roughly between ρ 0 /2 and 2ρ 0 [59]. For instance, in the case of a typical subsaturation density value ρ = 0.10 fm −3 , one finds that the above quadratic expansion of c sym (ρ) differs from the exact c sym (ρ) by less than 1% in many different nuclear mean field forces [51]. These facts point to the usefulness of investigating parameters such as L and K sym for the characterization of the density dependence of the symmetry energy.
The values of the DM coefficients J, L, and K sym for the non-relativistic forces and the RMF parameter sets considered in this work are given in Table I. In the Skyrme and Gogny effective interactions the symmetry energy coefficient at saturation J takes values around 30 MeV. The RMF parameterizations have larger values of J, also with a larger spread. The slope (L) and the curvature (K sym ) of the symmetry energy at saturation take even more widely scattered values among the different interactions. The consequence is that all mentioned nuclear models predict a different behavior of the symmetry energy at subsaturation densities, what can be seen e.g. in Fig. 1 of Ref. [13].
As is known, the density dependence of the symmetry energy near saturation tends to be much softer in the non-relativistic forces than in the covariant mesonexchange models of nuclear structure (see the values of L in Table I). In particular, taking into account Eq. (5), the slope of the neutron EOS and of the symmetry energy calculated at densities ρ close to the saturation value ρ 0 differ considerably between models. We also know that ∆R np in 208 Pb shows a linear dependence with these slopes at some subsaturation density ρ ≃ 0.10 fm −3 [37][38][39][40]. Therefore it is reasonable that the neutron skin thickness in 208 Pb and the leading term L/3ρ 0 of Eqs. (7) and (8) are related. As far as ρ 0 does not change much in the different effective forces, a correlation between ∆R np in 208 Pb and the DM coefficient L is thus expected [13-15, 33, 41] and we display it in the middle panel of Fig. 1.
From the discussed results of ∆R np versus J/Q and of ∆R np versus L in Fig. 1, a correlation between the DM coefficient L and the J/Q ratio is to be expected too. Note that L rules the density dependence of the symmetry energy of the nuclear equation of state [Eq. (8)], and that Q governs the thickness of the neutron skin of finite nuclei [Eqs. (2) and (3)]. The correlation between L and J/Q can be seen in the right panel of Fig. 1. This correlation shows that the value of the slope L of the symmetry energy at the saturation density increases with the value of the ratio between the bulk symmetry energy coefficient J and the surface stiffness coefficient Q. The trend is considerably linear among the various nuclear effective interactions.
Previous literature [56,58] has shown that the systematics of experimental binding energies relates increasing values of J with decreasing values of Q in nuclear effective interactions whose parameters have been adjusted to describe experimental data. We note this same trend in Table I, where the RMF sets that in general have larger J values also tend to have smaller Q values than their nonrelativistic counterparts. A smaller Q coefficient means that it is easier to develop a neutron skin in finite nuclei. Consistently, the neutron skin thickness in 208 Pb (or any other heavy nucleus) is usually larger when computed with a RMF parameter set than when computed with a non-relativistic force. These facts, and the noticed correlation between L and J/Q, allow one to interpret in a qualitative way within the DM, the correlation pointed out in Refs. [37][38][39][40] between the slope of the symmetry energy at some subsaturation density ρ ≃ 0.10 fm −3 and the neutron skin thickness in a heavy nucleus. In a recent work [51] we have investigated further the relations of the neutron skin thickness with the parameters L and K sym that characterize the density dependence of the symmetry energy around saturation.

III. SURFACE WIDTH CONTRIBUTION TO THE NEUTRON SKIN THICKNESS
The neutron skin thickness values derived from measurements performed in antiprotonic atoms have been ob-tained in Refs. [48,49]. It is assumed that the neutron skin is due to an enhancement of the neutron surface width with respect to the proton surface width, and that the mean location of the proton and neutron surfaces in these nuclei are the same. This situation corresponds to the so-called "neutron halo-type" distribution [48]. It has been shown that the same set of experimental values of neutron skin thickness can be explained with similar quality as in [48] by means of the "standard" version of the DM (where b n = b p ) [28]. The latter case assumes that the peripheral neutrons are concentrated at the neutron surface, which is shifted with respect to the proton surface, and that both the neutron and proton density distributions have the same surface width. This is rather the pattern of the so-called "neutron skin-type" distribution according to Ref. [48].
The analysis of neutron and proton densities calculated with nuclear mean field interactions carried out in Ref. [32] by means of the Helm model, points out that self-consistent mean field densities show a mixed character between the "neutron halo" and "neutron skin" patterns. This means that, actually, the self-consistent neutron and proton density profiles obtained with nuclear effective interactions differ not only in the mean location of their surfaces but also in their surface widths. In the following we shall see that similar conclusions are found from the calculations of the neutron skin thickness performed in the DM with formula (2). It will turn out that the surface width contribution ∝ (b 2 n − b 2 p ) in the DM expression (2) for the neutron skin thickness, which arises from b n = b p , is necessary to reproduce the neutron skin thickness values calculated from the definition (1) using self-consistent densities of finite nuclei obtained with the ETF approach in mean-field theory, for non-relativistic forces as well as for relativistic parameterizations.
In Fig. 2 we display by empty symbols, as a function of the overall relative neutron excess I = (N − Z)/A, the neutron skin thickness predicted by the "standard" version of the DM (namely, Eq. (2) with b n = b p ) using some well-known effective forces. The nuclei are those from 40 Ca to 238 U measured in the experiments with antiprotonic atoms [48,49], and which were studied with the DM in Ref. [28]. The values shown in Fig. 2 have been computed using the SIII and SkM* Skyrme forces and the NL-SH and NL3 RMF parameter sets, as suitable examples. We have chosen these four parameter sets for display in Fig. 2 because they span the whole range of values of the ratio J/Q of nuclear interactions that describe reasonably well the ground-state properties of finite nuclei, having a bulk symmetry energy coefficient J between 28 MeV and 37 MeV (see Fig. 1 and Table I). The DM results for ∆R np obtained with the other mean field interactions considered in this work that have a J coefficient between the values of SIII and NL3, also lie within the window of results delimited by the SIII and NL3 interactions in Fig. 2.
The values of ∆R np predicted by the DM are compared in Fig. 2 with the values that we obtain from self-consistent ETF calculations in finite nuclei (filled symbols). Both models do not incorporate shell effects. In Fig. 2 we have used the ETF approach in the nonrelativistic [62] and relativistic [63] frameworks to compute ∆R np . We have calculated these ETF values of ∆R np by application of Eq. (1) with the rms radii of the self-consistent neutron and proton densities obtained in each isotope. From Fig. 2 two significant points stem. First, the predictions of the DM in the "standard" form (b n = b p ) systematically undershoot the ETF neutron skin thickness computed in finite nuclei with the selected effective nuclear interactions. In particular, this trend is reinforced with growing neutron excess I. Second, it can be observed that for a given nucleus the difference between the ETF value of ∆R np computed with (1) and the value provided by the "standard" DM prescription is slightly larger in the RMF parameter sets than in the Skyrme forces. Altogether, these facts suggest that in the mean-field interactions the surface width contribution to the DM formula for ∆R np does not vanish, and that this contribution has some dependence on the ratio J/Q of the force.
To apply the full Eq. (2) to compute the neutron skin thickness including the surface width correction, one needs to evaluate the neutron and proton surface widths in finite nuclei. In practice, there is not conclusive experimental evidence on the difference of the surface widths b n and b p of the nuclear density distributions [27,28,32,50,54]. To estimate b 2 n − b 2 p we will therefore rely on theoretical guidance as a surrogate. Within the context related to the DM and the leptodermous expansion of a finite nucleus [52][53][54], the surface properties in finite nuclei can be extracted from ASINM calculations. The semi-infinite geometry does not include shell, Coulomb or finite-size effects. We will use here the ETF method including 2 corrections for describing ASINM, since the ETF method is free of the Friedel oscillations of the quantal densities [30,57]. We summarize in the Appendix the basic aspects of the procedure. More details can be found in Ref. [58]. From Eqs. (14)-(17) of the Appendix, we can obtain the values of the surface widths b n and b p in ASINM with a given relative neutron excess in the bulk δ 0 . In the DM, these surface widths correspond to the values b n and b p in finite nuclei if δ 0 is calculated from the overall relative neutron excess I of the nucleus through the following relation [52,54,58]: which takes into account the Coulomb correction.
Once the neutron and proton surface widths in finite nuclei are known, we can compute their contribution to the neutron skin thickness, which reads [cf. Eq. (2)] The corresponding values of ∆R sw np for the nuclei considered in Fig. 2 are displayed in the bottom panel of  [62,63]. One now observes an improved and remarkable agreement between the DM predictions and the self-consistent ETF values computed with the same interaction, stemming from the inclusion of the calculated ∆R sw np contribution. It is interesting to note that the neutron skin thickness obtained with Eq. (1) from the rms radii of the self-consistent ETF calculations in finite nuclei shows a well defined increasing linear tendency with the relative neutron excess I, similarly to the case of the results of the DM and in consonance with the trend of the experimental values derived from antiprotonic atoms [48].
The lower panel of Fig. 3 also suggests that, for a given nucleus, ∆R sw np grows when the J/Q ratio of the nuclear interaction increases (see Table I). To analyze this behavior in more detail, we fit ∆R sw np by means of a law σ sw I, which defines the slope σ sw of ∆R sw np with respect to the relative neutron excess I. This slope is displayed in Fig. 4 as a function of the J/Q ratio for different interactions. The slopes σ sw lie inside a band limited by two straight lines, corresponding to the equations σ sw = 0.3J/Q+0.07 fm (left) and σ sw = 0.3J/Q − 0.05 fm (right). It is worth to notice that all considered Skyrme forces have slopes σ sw below 0.25 fm, whereas the analyzed RMF models have slopes σ sw always above this value.
One relevant conclusion is that bulk and finite nuclei properties described through successful theoretical meanfield models constrain the possible values of the surface width contribution to the neutron skin thickness between the limits portrayed in Fig. 4. From this figure it can be deduced that the surface width contribution to the neutron skin thickness ∆R sw np has, on top of a global increasing trend with J/Q, a more involved dependence on the parameters of the effective nuclear interactions. For instance, on the one hand, the RMF force FSUGold and the Skyrme force SkM* have almost the same J/Q ratio (see Table I). However, the predicted values of σ sw are clearly different for both interactions as can be seen in Fig 4. On the other hand, it is possible to find different interactions which have almost the same slope σ sw , and therefore the same ∆R sw np , but with different values of the J/Q ratio. Some examples of this fact in Fig. 4 are the Skyrme forces SIII and SGII (with slopes of 0.16 fm and 0.15 fm, respectively), and the RMF parameterizations FSUGold and NL3 (with slopes 0.30 fm and 0.31 fm). These facts suggest that it is possible to find the same total neutron skin thickness by combining a small value of the J/Q ratio in t [see Eq. (3)] with a large ∆R sw np contribution or, vice versa, by combining a large J/Q ratio in t with a small ∆R sw np contribution.

IV. ESTIMATES OF THE DENSITY CONTENT OF THE SYMMETRY ENERGY
As we have pointed out, the nuclear symmetry energy and the properties of the EOS of neutron-rich matter are of increasing importance in both nuclear physics and astrophysics. There is a significant recent effort in the community towards constraining the values of the parameters that characterize the density dependence of the symmetry energy in the subsaturation regime of the EOS. The new developments come in particular from the investigation of isospin-sensitive observables in intermediateenergy heavy-ion collisions [9][10][11][12][13][14][15][16][17][18][19][20][21] and in nuclear resonances [22][23][24][25]. Obviously, the different studies do not deal with exactly the same regimes of density and energy. One also has to keep in mind that the connection of the experiments with the EOS is not at all trivial; it often requires of extrapolations of the measured data, which imply a model dependence. Therefore, it is important to further investigate indications from those and other experimental probes of the symmetry energy, with different methodologies, as well as to study the interplay between the constraints derived from the different analyses.
In the present section we want to apply the experience gained in the DM study of the neutron skin thickness performed in the previous sections, to estimate possible constraints on the density dependence of the nuclear symmetry energy as suggested by neutron skin data. To this end, we shall first obtain the range of values of the J/Q ratio which are compatible with the neutron skin thickness derived from the experimental data in antiprotonic atoms. We recall that this set of data for 26 stable nuclei is till date the largest set of uniformly measured neutron skins spanning the mass table. Once the estimated range of values for the ratio J/Q will be found, the constraints derived from the neutron skin data on the DM parameter L, related to the density derivative of the symmetry energy, will be determined from the existing linear correlation between the values of L and J/Q. This correlation has been displayed in Fig. 1 and is a general feature of mean field interactions that have been adjusted to reproduce with good accuracy binding energies and charge radii (often among other properties) of nuclei across the periodic table.
A. Constraints on the J/Q ratio From the previous section we know that the surface width part in Eq. (2) gives a non-negligible contribution to neutron skins in effective nuclear interactions. This contribution is needed to reproduce the neutron skin thickness values computed self-consistently in ETF calculations of finite nuclei. We have also seen that to leading order, both the mean location of the neutron and proton surfaces (3) and the surface width correction (10) are basically driven by the value of the J/Q ratio. The discussions in the previous sections suggest to fit the experimental ∆R exp np data by means of the following DM inspired ansatz: where t is given by Eq. (3). The second term is the surface width contribution. It is parameterized to reproduce the dashed lines on Fig. 4, with c = 0.07 fm or c = −0.05 fm. With the ansatz (11), we will use J/Q as an open parameter. It will be constrained by a least-squares minimization from the experimental values ∆R exp np derived from the analysis of antiprotonic atoms [48,49]. We note that Eq. (11), as well as t given by Eq. (3), depends on the particular values of the symmetry energy at saturation J and of the nuclear matter radius r 0 . We fix these quantities to the empirical values J = 31.6 MeV and r 0 = 1.143 fm (the latter corresponds to a saturation density ρ 0 = 0.16 fm −3 ). We consider the values c = 0.07 fm and c = −0.05 fm in Eq. (11), discussed in connection with Fig. 4, to simulate the upper and lower bounds of the window of the theoretical predictions for σ sw obtained according to mean-field models of nuclear structure. In the χ 2 -minimization we have weighted each ∆R exp np − ∆R np difference by the inverse of the associated experimental uncertainties. That is, in practice we have minimized the quantity where ∆R np is calculated with Eq. value of one standard deviation associated with the fit made through Eq. (12). To check our method of minimization and error estimation, we have applied the same procedure to make a linear fit mI + n of the experimental data ∆R exp np . In this case we have obtained ∆R np = (0.901 ± 0.147)I + (−0.034 ± 0.023) fm, which fully agrees with the result quoted by the experimentalists [48,49].
The results for ∆R np from our fits compared with the experimental data ∆R exp np are displayed as a function of I in Fig. 5. Both extractions of J/Q, for c = 0.07 and c = −0.05 fm, predict basically the same total neutron skin thickness with a similar quality and they are close to the average ∆R np = (0.90 ± 0.15)I + (−0.03 ± 0.02) fm [48,49] of the experimental data. However, the splitting of the neutron skin thickness into a part coming from the distance t and another part coming from the surface width ∆R sw np is different in both cases, as we have discussed in the previous section. Therefore, it becomes clear that the experimental neutron skin thickness data, by themselves, may be able to constrain the total value but not its partition into a bulk and a surface width contribution.
It is generally acknowledged that the value of the bulk nuclear symmetry energy coefficient J is about 30-32 MeV, but there is some uncertainty. To assess the dependence of the extraction of J/Q on the assumed value for the J coefficient (which in the above we have taken as 31. We recall that the DM of nuclei does not incorporate shell effects and averages the corresponding quantal magnitudes. With the method described we have fitted ∆R exp np that in general contains shell effects, and possible correlation and deformation contributions. But, as mentioned, the neutron skin data analyzed show a well defined linear trend with the relative neutron excess I (namely, ∆R np = (0.90 ± 0.15)I + (−0.03 ± 0.02) fm [48,49]). This trend, and the agreement of the DM values for ∆R np with the self-consistent ETF calculations in finite nuclei, which also are free of shell effects and have a linear trend with I (Figs. 2 and 3), gives more reliability to the predictions obtained with the DM formula from the experimental data.
B. Constraints on the L parameter As we have discussed previously, the parameter L has a direct relation with the slope of the symmetry energy of the nuclear EOS, see Eq. (8). Having determined the values of the J/Q ratio compatible with the experimental neutron skins measured in antiprotonic atoms, we can use the linear correlation between L and J/Q found in meanfield effective nuclear interactions, to obtain an insight on the values of the parameter L favored by neutron skins.
We have displayed the linear correlation L = mJ/Q+n in the rightmost panel of Fig. 1 (the linear correlation coefficient for the shown interactions is r = 0.978). The values of the m and n coefficients have some dependence on the set of interactions chosen to make the linear regression. We have checked that this dependence is rather weak. Namely, we have tested the correlation of L with J/Q by taking into account successively 10, 14, 18, and 24 interactions and we have found the linear regressions to be comprised between L = 139J/Q − 52 MeV and L = 150J/Q − 57 MeV. Considering these two limiting cases and the constraint 0.6 J/Q 0.9 found in the previous section, leads to a variation of L between 31 MeV and 78 MeV. Thus, our estimate for the L coefficient, which takes into account the surface width correction in ∆R np obtained in the calculations with meanfield interactions, basically lies in the range 30 L 80 MeV. Had we kept the value of J fixed at 31.6 MeV, the extracted range for L would be a little narrower: In a previous work [51] we have investigated the correlations between the symmetry energy coefficient in finite nuclei and in the EOS at subsaturation densities. These correlations allow one to derive an aproximate formula for the neutron skin thickness with explicit dependence on the L coefficient. By comparison of that result for the neutron skin thickness with the experimental data set of Refs. [48,49], in Ref. [51] we found a range of values L = 55 ± 25 MeV (displayed in Fig. 3 of Ref. [51]) when one includes the surface width contribution ∆R sw np in the calculations. That prediction is consistent with the val-ues obtained here by the present procedure. A somewhat higher range L = 75 ± 25 MeV was obtained [51] when one neglects ∆R sw np . Although a vanishing ∆R sw np value, corresponding to b n = b p in the nucleon density distributions, is not favored by the mean field interactions (see Section III), it cannot be discarded without having more experimental evidence on the value of b n as we have noted in Section III (see also Refs. [28,48]).
In recent years, considerable advances in probing experimentally the density dependence of the symmetry energy at subsaturation have been achieved in heavy ion collisions (HIC) at intermediate energies. It has been found that the symmetry energy can be modelized around the saturation density with reasonable good approximation by [11][12][13][14][15][16][17][18][19][20]  Our estimates for the stiffness γ can be compared with alternative predictions derived in the recent literature. For instance, our range 0.32 γ 0.84 overlaps with the values γ ∼ 0.55-0.77 that Danielewicz [29] obtains from the study of binding energies, neutron skins, and isospin analog states of selected nuclei. It also contains the value γ ∼ 0.55 that is inferred from the analysis of neutron-proton emission ratios in HIC carried out by Famiano et al. [20], as well as with the value γ ∼ 0.69 obtained by Shetty et al. [17][18][19] from isotopic scaling in intermediate-energy nuclear reactions.
The stiffness of the symmetry energy at subsaturation densities also has been investigated from isospin diffusion data in HIC, by means of simulations with an isospin-and momentum-dependent transport model with in-medium nucleon-nucleon cross sections [12][13][14][15][16]. In this case the prediction is 0.69 γ 1.05. It corresponds to a behavior of c sym (ρ) that nearly ranges between the familiar ρ 2/3 dependence of the purely kinetic symmetry energy of a free Fermi gas, in the lower limit, and linearity in the density ρ, in the upper limit. The constraints on the stiffness of the symmetry energy derived from isospin diffusion, combined with an analysis of the properties of Skyrme interactions, are found to lead to a constraint Ref. [29] n-p emission ratios [20] isoscaling [17][18][19] isospin diffusion [12][13][14][15][16] Thomas-Fermi model [64,72] GDR [25] PDR [24] FIG. 6: Comparison of the estimated values of the parameter L from different observables and methods. Some of the estimates have been analyzed through Eq. (13) for csym(ρ).
63 L 113 MeV [12][13][14][15][16]. The predictions from isospin diffusion are thus a little stiffer, though the lower limits of γ and L obtained using this method are in agreement with the upper limits of γ and L obtained from our study of neutron skins with inclusion of the surface width contribution.
Another valuable reference comes from the celebrated Thomas-Fermi model of Myers andŚwiatecki [64,72]. This model was fitted very precisely to the binding energies of a comprehensive set of 1654 nuclei. It predicts an EOS that leads to a coefficient L = 49.9 MeV. Note that if we compare c sym (ρ) calculated from the EOS of the Thomas-Fermi model with Eq. (13), an exponent γ = 0.51 is obtained. Additional information on the density content of the symmetry energy arises from the constraints on the symmetry pressure P sym = ρ 0 L/3 extracted by Klimkiewicz et al. [24] from the properties of pygmy dipole resonances in nuclei. These are indicative of a value γ ∼ 0.35-0.65 if one assumes P sym = ρ 0 γJ following from Eq. (13) given above. On the other hand, Trippa et al. [25] have obtained the constraint 23.3 < c sym (ρ = 0.1 fm −3 ) < 24.9 MeV from consideration of the giant dipole resonance in 208 Pb, which implies a range ∼ 0.5-0.65 for the γ exponent. We depict in Fig. 6 the estimated ranges of values for the L parameter from the discussed analyses.
In summary, in spite of the discrepancies in the details, the various findings from experimental isospin-sensitive signals, including ours, agree all on a rather soft nuclear symmetry energy at subsaturation densities. Recent studies of pure neutron matter at low densities based on universal properties of dilute Fermi gases lead to a similar conclusion [73,74]. One may mention that there exists recent circumstantial evidence [75], derived from π − /π + ratios in central HIC collisions at SIS/GSI ener-gies, hinting at that the nuclear symmetry energy is soft also in the regime of suprasaturation densities (ρ ≥ 2ρ 0 ). However, further experimental and theoretical confirmations of this fact need to be awaited [75].

V. SUMMARY AND CONCLUSIONS
The droplet model predicts that the neutron skin thickness of atomic nuclei is correlated with the ratio J/Q, where J is the symmetry energy in bulk matter and Q is the surface stiffness coefficient. We have shown that the J/Q ratio displays a linear relationship with the DM parameter L in nuclear mean field models that are calibrated to experimental ground-state properties such as binding energies, charge radii, and single-particle data. In this way, the known correlation between the neutron skin thickness in a heavy nucleus and the density derivative of the symmetry energy (or of the neutron equation of state) evaluated at a subsaturation density, can be interpreted in the context of the DM.
According to the droplet model, the neutron skin thickness is correlated with the overall relative neutron excess I = (N − Z)/A of nuclei. This fact is in agreement with the experimental findings using information from antiprotonic atoms. The DM expression for the neutron skin thickness contains "bulk" and "surface" parts. The bulk part corresponds to the contribution proportional to the distance t between the neutron and proton mean surface locations. This part is quite dependent on the Skyrme force or RMF parameterization used to compute it (see Fig. 2). In finite nuclei, this DM bulk contribution systematically underestimates the neutron skin thickness extracted directly as the difference of the neutron and proton rms radii of the nucleus from self-consistent ETF calculations with effective forces. This evidence indicates that the surface part, due to the b n = b p contribution, is necessary in the DM formula to properly estimate the neutron skin thickness in finite nuclei in mean-field models using effective nuclear interactions.
The DM surface contribution ∆R sw np to the neutron skin thickness is smaller than the bulk part and shows a well defined linear increasing tendency with the overall relative neutron excess I. We have investigated the dependence of the slope σ sw with respect to I of this surface contribution using various Skyrme and RMF forces. We have found that the slopes σ sw lie in a region of the σ sw -J/Q plane that can be roughly limited by two straight lines as a function of the J/Q value. It implies that nuclear properties, which are included in the calibration of the free parameters of the Skyrme and RMF interactions, constrain the possible values of the surface width contribution to the neutron skin thickness in the DM. If the same nuclear interaction is used, a good agreement between both the DM formula and the selfconsistent ETF calculations of ∆R np in finite nuclei is found along the whole periodic table when the contribution ∆R sw np is included in the DM.
To analyze possible bounds suggested by experimental neutron skin data on the value of the J/Q ratio, we have adjusted the DM neutron skin thickness formula to the neutron skin sizes measured in antiprotonic atoms [48,49]. We have determined a window 0.6 J/Q 0.9 for the J/Q ratio by using the largest and smallest surface contributions ∆R sw np obtained from successful Skyrme forces and RMF parameterizations. These two fits reproduce the experimental data with almost the same quality. In other words, the experimental data of the neutron skin thickness in finite nuclei constrain the total theoretical estimate but not its partition into a bulk and a surface contribution. Once the window of J/Q values is known, the compatible range of values of the parameter L can be estimated from the linear correlation between L and J/Q shown in Fig. 1. From our analysis we find the constraints 30 L 80 MeV.
If a model symmetry energy c sym (ρ) = J(ρ/ρ 0 ) γ is assumed, a prediction for the value of the stiffness γ of the symmetry energy can be obtained, with use of the empirical values of J and ρ 0 . In this way we find the estimate 0.32 γ 0.84. Thus, our analysis of the experimental neutron skins deduced from antiprotonic atoms suggests a relatively soft symmetry energy, in good accord with the recent indications from pygmy [24] and giant [25] dipole resonances. Our prediction for the stiffness γ of the symmetry energy also is in reasonable agreement with the constraints derived by Danielewicz [29], by Famiano et al. [20], and by Shetty et al. [19] using different observables, as well as with the value γ = 0.51 of the EOS of the Thomas-Fermi model of Myers andŚwiatecki [64,72]. In the upper limit, our prediction overlaps with the lower limit provided by the analysis of isospin diffusion data in intermediate-energy heavy ion collisions [12][13][14][15][16].
In summary, different techniques of extracting the parameters that describe the density dependence of the symmetry energy predict different values. However, taking the average of the central values of the predictions displayed in Fig. 6, with account of their uncertainties when available, there exist narrow windows of the parameters L ∼ 45-75 MeV and γ ∼ 0.5-0.8 which are, actually, compatible with all the different methods mentioned to obtain them. These ranges of values of the indicated parameters for describing the leading density dependence of the symmetry energy in bulk matter seem to be the "optimal" ones according to present experimental evidence from nuclear data.
To conclude, we are aware that the neutron skin thickness data derived from antiprotonic atoms are to some extent model dependent, and have for some nuclei large error-bars. Also, our theoretical method represents just an average approximation. In spite of these limitations, we hope to have shown that from neutron skin data it is possible to make a reasonable estimate of the density dependence of the symmetry energy with uncertainties that are not significantly much larger than those currently obtained from other experimental observables. One expects that future data from the planned parity-violating elec-tron scattering experiment for measuring the neutron radius in 208 Pb [2] will contribute to narrow down the constraints derived here from the thickness of neutron skins of nuclei. by grant AP2005-4751 from MEC (Spain). M.W. gratefully acknowledges the warm hospitality extended to him by the Nuclear Theory Group and the Departament d'Estructura i Constituents de la Matèria at the University of Barcelona during his stay in Barcelona.

VI. APPENDIX
To compute the surface stiffness coefficient Q, and also the neutron and proton surface widths b n and b p that appear in the contribution ∆R sw np of Eq. (10) to the neutron skin thickness, we obtain the self-consistent neutron and proton density profiles in asymmetric semi-infinite nuclear matter (ASINM). To do that we consider a semiinfinite slab with a plane interface separating a mixture of protons and neutrons at the left, whose densities decrease smoothly to zero at the right as empty space is reached. The axis perpendicular to the interface is taken to be the z axis. Thus, the relative neutron excess δ = (ρ n − ρ p )/(ρ n + ρ p ) depends locally on the z coordinate. When z goes to minus infinity, the neutron and proton densities approach the values of asymmetric uniform nuclear matter in equilibrium, corresponding to an interior bulk neutron excess δ 0 .
In order to obtain the proton and neutron densities in ASINM one has to minimize the total energy per unit area with respect to arbitrary variations of the densities, with the constraint of conservation of the number of protons and neutrons. When δ 0 is not very large, so that occurrence of drip nucleons does not take place (which is the situation in all cases considered in the present work), the constrained energy per unit area reads [56,58,68] In this equation, ε(z) represents the energy density functional of the nuclear effective interaction under investigation, and µ n and µ p are the neutron and proton chemical potentials. The explicit expressions of ε(z) in the extended Thomas-Fermi (ETF) method for Skyrme forces and relativistic mean field interactions can be found in e.g. Appendix A of Ref. [58]. The proton and neutron densities obey the coupled local Euler-Lagrange equations δε(z) δρ n − µ n = 0, δε(z) δρ p − µ p = 0.
We solve them fully self-consistently by numerical iteration [62]. We note in this respect that we have not used any parameterized form of the densities such as e.g. Fermi shapes. In the relativistic model the variational equations (15) are supplemented with additional field equations for the meson fields; the calculational details for the relativistic problem can be found in Refs. [58,62,63]. From the calculated density profiles, one obtains the mean locations of the surfaces, z 0q (q = n, p), as where the primes indicate a derivative with respect to the z coordinate. The proton and neutron surface widths in ASINM are obtained as the second moment of ρ ′ (z): The distance t = z 0n − z 0p between the mean surface locations of the neutron and proton density profiles allows one to extract the surface stiffness coefficient Q from the relation [52] t = z 0n − z 0p = 3r 0 2 which is valid in the limit of small asymmetries. For each given nuclear interaction we solve the ASINM problem for 5 different values of δ 0 , between 0.005 and 0.025, and we then evaluate Q from the slope of t.
There exists a second way of computing Q. It is based on the fact that in dividing the energy in bulk and surface parts, as soon as δ 0 = 0, there are two possibilities to define the bulk reference energy [58,68]. One definition is based on the chemical potentials µ n and µ p of each nucleon species; this definition is thermodynamically consistent and it is the one that we have given in Eq. (14). The second definition of the reference energy is based on taking the value of the energy per particle in bulk asymmetric nuclear matter. Accordingly, upon the bulk reference energy chosen, there exist two forms of the surface energy in ASINM, which are called E surf,µ and E surf,e [58,68]. In the small asymmetry limit, it can be shown that the difference between these two quantities behaves as Thus, the slope of E surf,e − E surf,µ with respect to δ 2 0 provides another means to extract the value of the surface