{'fraction_non_alphanumeric': 0.030906697170483242, 'fraction_numerical': 0.016467744372261672, 'mean_word_length': 4.510048757283863, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analysis of bulk and surface contributions in the neutron skin of nuclei

The neutron skin thickness of nuclei is a sensitive probe of the nuclear symmetry energy having multiple implications for nuclear and astrophysical studies. However, precision measurements of this observable are difficult. The analysis of the experimental data may imply some assumptions about the bulk or surface nature of the formation of the neutron skin. Here, we study the bulk or surface character of neutron skins of nuclei following from calculations with Gogny, Skyrme, and covariant nuclear mean-field interactions. These interactions are successful in describing nuclear charge radii and binding energies but predict different values for neutron skins. We perform the study by fitting two-parameter Fermi distributions to the calculated self-consistent neutron and proton densities. We note that the equivalent sharp radius is a more suitable reference quantity than the half-density radius parameter of the Fermi distributions to discern between the bulk and surface contributions in neutron skins. We present calculations for nuclei in the stability valley and for the isotopic chains of Sn and Pb.


I. INTRODUCTION
The description of the sizes and shapes of atomic nuclei is among the oldest problems in nuclear physics. The rms radius of the charge distribution in nuclei surely is the most prominent example of this kind of observable. Owing to the high degree of accuracy achieved by the elastic electron-nucleus and muon-nucleus scattering experiments, the nuclear charge radius is nowadays known with uncertainties that are for many nuclei smaller than 1% [1,2]. In contrast, our knowledge about the neutron distribution and its rms radius in nuclei is not so precise. Actually, accurate determinations of the rms radius of the neutron density are still lacking. This implies that the so-called neutron skin thickness, generally defined as the neutron-proton rms radius difference in the atomic nucleus, is not precisely known either. The neutron skin thickness observable, Eq. (1), is a very sensitive probe of the pressure difference that exists between neutrons and protons in the atomic nucleus. As such, ∆R np is intimately correlated with the density dependence of the nuclear symmetry energy and with the equation of state of pure neutron matter [3][4][5][6][7][8][9][10][11][12]. Owing to this fact, the accurate calibration of ∆R np is a problem of significant implications for studies that embrace diverse facets of both nuclear physics and nuclear astrophysics (see for example Refs. [13][14][15][16][17][18][19][20][21][22][23]). It should be mentioned that Eq. (1) is not the only useful prescription to characterize the different spatial extension of the neutron and proton densities in a nucleus. The neutron-proton radius difference has also been computed using Helm radii instead of rms radii for the discussion of nuclear halos in the literature [24,25]. Nevertheless, in the present work we use the conventional and more frequent definition (1) for ∆R np , Eq. (1).
Because neutrons are uncharged particles, the measurement of their spatial distribution in the nucleus is more difficult than for the positively charged protons (see, e.g., Ref. [26]). The experimental access to neutron densities and neutron rms radii usually involves strongly interacting hadronic probes, for example, in the case of experiments that perform proton-nucleus elastic scattering [27][28][29][30][31], α-particle elastic scattering [32], and techniques based on the inelastic scattering excitation of the giant dipole and spin-dipole resonances [33,34]. The neutron skin thickness of nuclei has also been investigated through radiochemical and x-ray techniques in antiprotonic atoms [35][36][37][38][39][40], taking advantage of the fact that the nuclear periphery is very sensitive to antiprotons in the normally electronic shell. Though the scope of our analysis of neutron skins is of theoretical, in this article we refer to some extent to the experimental investigations in antiprotonic atoms to set the stage for our calculations.
A few years ago, Trzcińska et al. [38,39] extracted the neutron skin thickness of a large set of nuclei in experiments with antiprotons conducted at the former LEAR facility of CERN. The measurements were made for 26 stable isotopes distributed across the mass table, from the light and symmetric nucleus 40 Ca to the heavy and asymmetric nucleus 238 U. Because of the fact that the antiproton-nucleon interaction is very strong, antiprotons are able to interact with the atomic nucleus at distances where the nuclear density is much smaller than its central value. Slow enough antiprotons can form a hydrogen-like atom with the nucleus. When the antiproton annihilates with a nucleon producing pions that may miss the nucleus, it leaves a residue that is one neutron or proton fewer. From the analysis of these yields, information about the neutron distribution in the nucleus can be obtained [35][36][37][38][39]. A second experimental method measures antiprotonic x-rays from where the atomic level widths and shifts owing to the strong interaction are determined [38][39][40]. By combining the results obtained with these two experimental techniques, the neutron-proton rms difference ∆R np can be deduced provided that the charge density of the nucleus is known [38][39][40].
The extraction of ∆R np values from antiprotonic atoms assumes nucleon densities in the form of twoparameter Fermi (2pF) distributions. The procedure involves interpreting whether the difference between the peripheral neutron and proton densities arises from an increase of the mean location of the surface of the neutron density (i.e., from an increase of the bulk radius of neutrons) or, rather, from an increase of the surface diffuseness of the neutron density. This question is also instrumental in studies of properties of nuclei by parity-violating electron scattering [41]. The radiochemical data in antiprotonic atoms were shown to be in favor of interpreting ∆R np as an increase of the neutron surface diffuseness [38,39] but with the caveat that some room existed within assigned errors for an intermediate situation [38]. Indeed, a comparison with the droplet model [42] suggests that the neutron skin sizes of the antiprotonic measurements can be described similarly well in the droplet model theory by a difference in the diffuseness as by a difference in the bulk radius of the neutron and proton densities.
In the present work we theoretically investigate the bulk and surface components in the neutron skin thickness of nuclei by parametrizing self-consistently calculated nucleon densities by 2pF distributions. The 2pF form is also common in experimental analyses. Our calculations are performed with some representative effective nuclear forces. These are the finite-range Gogny D1S interaction and the zero-range Skyrme SLy4 force from the nonrelativistic framework, and the NL3 and FSUGold parameter sets from the relativistic mean-field (RMF) framework. We discuss the predictions that these meanfield models make for the decomposition of the neutron skin thickness in bulk and surface contributions for the set of stable nuclei that were analyzed in the experiments in antiprotonic atoms. We also study the theoretical predictions along the isotopic chains of Sn and Pb and the variation of the results as one moves from the proton to the neutron drip line in these isotopic chains.
The structure of this article is the following. In the Sec. II we present a brief summary of the experimental methodology and results in antiprotonic atoms to high-light some interesting aspects for our study. In Sec. III we discuss the common definitions of nuclear radii and characterize the bulk and surface contributions in the neutron skin thickness of nuclei. In Sec. IV we analyze the theoretical mean-field results in stable isotopes and the predictions for nuclei across the Sn and Pb isotopic chains. We present the summary and conclusions in Sec. V. Some relations for 2pF functions are collected in the Appendix.

II. SOME ASPECTS OF THE METHODS AND RESULTS IN ANTIPROTONIC ATOMS
The radiochemical study of antiprotonic atoms [35][36][37][38][39][43][44][45] consists of the analysis of the nuclei with a mass number one unit smaller than the target mass number A t where the antiproton annihilation takes place. These products that have one less neutron or proton are short lived and their decay is followed by emission of γ rays. Standard nuclear spectroscopy methods allow one to determine the absolute number of these residual nuclei with mass number A t − 1. The probability of producing a cold product with mass number A t − 1 is calculated by using the antiproton-nucleus optical potential fitted to reproduce the x-ray experimental data (atomic level shifts and level widths) [38,40,46]. The probability distribution for obtaining a nucleus with mass number A t −1 has a maximum located about 2-3 fm outside the halfdensity radius R 1/2 of the target nucleus.
To compare the experimental data for any target it is convenient to introduce the peripheral halo factor [47]: The first term in brackets on the right-hand side of this equation, gives the ratio of thep annihilations on peripheral neutrons to thep annihilations on peripheral protons, normalized with the Z t /N t value of the target nucleus. The quantities δp n and δp p are thep-n and p-p scattering amplitudes, respectively, and the factor F = Im(δp n )/Im(δp p ) gives the ratio of thep annihilation probabilities on a neutron and on a proton. Consequently, halo factor (2) essentially measures the change of the ratio of the neutron-to-proton concentration in the peripheral region with respect to the bulk value represented by the N t /Z t ratio in the target nucleus. A halo factor larger (smaller) than 1 means an increase of the relative neutron (proton) concentration in the nuclear periphery. Assuming the ratio F is known, the measurement of the N t −1 and Z t −1 product yields allows one to obtain the experimental value of the halo factor. It was concluded that the comparison of the results from the antiprotonic x-ray method with the radiochemical method favors a value of F ≈ 1 [38]. The number of nuclei withone fewer neutron or proton is proportional to the quantities Im(δp n )ρ n or Im(δp p )ρ p , respectively, integrated over the region where the annihilation process takes place. Thus, approximately, one can estimate theoretically the halo factor as [37,48] f theor halo (r) ≈ This factor reaches the experimental value f expt halo at a distance r ≈ R 1/2 + 2.5 fm beyond the half-density radius of the charge distribution, where the antiproton annihilation takes place. One can use this fact to reproduce the neutron density distribution from some experimental observables.
The charge distribution in many stable nuclei is known very precisely. Usually the experimental charge density is given in some analytical form that is or may be converted to a 2pF distribution [49,50]. The simple 2pF formula has only two free parameters with clear physical meaning: on the one hand, a describes the diffuseness of the surface of the density profile: on the other hand, the half density or central radius C describes the mean location of this surface (i.e., C is indicative of the extension of the bulk part of the density distribution). The other, more sophisticated, parametrizations of the density distributions [49,50] describe better the central part of the nucleus or modify the surface part with higher-order terms. We do not use them because the interpretation of the multiple parameters is harder and less direct.
To study the differences at the nuclear periphery between the neutron and proton densities in the antiprotonic atoms experiments, the authors of Ref. [38] used 2pF functions. They used the notation "neutron skintype" distribution and "neutron halo-type" distribution to describe the two extreme cases of 2pF shapes having either C n > C p and a n = a p , or C n = C p and a n > a p , respectively. As we mentioned in the Introduction, the radiochemical data gave support to understanding the difference between the neutron and proton rms radii in the antiprotonic atoms from an increase of a n rather than from an increase of C n , compared to the a p and C p values. Thus, the neutron halo-type distribution (C n − C p = 0) was assumed [38]. (Note that here "halo type" is a useful notation, but the nuclei in the antiprotonic experiments are stable isotopes and it is not meant that they have halos like very light or exotic nuclei, such as 11 Li [24,25]). It was noted in the same work [38] that some room was left, nevertheless, within the error bars for intermediate cases with C n > C p and a n > a p .
The 2pF shape can be applied for the description of charge, proton, or neutron densities. Indeed, if the charge density is known in the 2pF form, the corresponding point proton density can be easily found and it also takes a 2pF form, with the parameters obtained through the deconvolution procedure [51,52]. The more relevant expressions are given in the Appendix. From these proton D1S SLy4 C n -C p = -0.20 C n -C p = -0.10 C n -C p = 0.00 C n -C p = 0.10 C n -C p = 0.20 FIG. 1: (Color online) The halo factor in 208 Pb as a function of the distance from the center of the nucleus, calculated using the 2pF formulas from the experimental charge density [53] assuming a neutron skin thickness ∆Rnp = 0.16 fm and different values of Cn − Cp (thin lines). The results of the theoretical predictions with the Gogny D1S and Skyrme SLy4 nuclear forces are also shown. The values of the halo factor deduced from experiment [38,40] are marked by crosses. The value of the proton half-density radius Cp is indicated by an arrow.
density profiles one can try to deduce the neutron density profiles with some additional information. For instance, if one knows the C p and a p parameters of the proton density and the value of the neutron skin thickness ∆R np defined in Eq. (1) (e.g., from independent experimental measurements), a relationship between the half-density radius C n and the diffuseness a n of the neutron 2pF distribution can be found (see the Appendix). Such a procedure creates a family of neutron density profiles that depend on one free parameter, which can be chosen to be either C n − C p or a n − a p . Now it is possible to check which set of values of C n and a n reproduces the experimental halo factor (2). This allows one to determine the 2pF neutron distribution. Next we apply this idea to the 208 Pb nucleus as an illustrative test case.
In Fig. 1 we display the halo factor, Eq. (3), for the 208 Pb nucleus computed with 2pF density distributions. The 2pF neutron densities are obtained from the experimental charge density taken from Ref. [53], which we convert to a 2pF proton density using the formulas given in the Appendix. We consider a family of 2pF neutron densities that differ one from the other in the half-density radius C n and diffuseness a n , but they all have the same value for the neutron skin thickness: ∆R np = 0.16 fm (which corresponds to the value extracted from the results of the x-ray method in the 208 Pb antiprotonic atom [40]). We allow the difference C n − C p to vary in the range from −0.20 to 0.20 fm. The value of the diffuseness a n of the 2pF neutron density is then obtained using Eq. (A10) from the Appendix.
The values of the 2pF parameters in our present calculation are given in Table I. Agreement with the results from previous experiments [38,40] (indicated by crosses in Fig. 1) is found for the parameters of the 2pF distributions in the range from the values C n − C p = 0 fm and a n − a p = 0.13 fm to the values C n − C p = 0.10 fm and a n − a p = 0.07 fm. In the first limit, the neutron skin thickness ∆R np is basically due to an enhancement of the diffuseness of the neutron density, whereas the second limit corresponds to a mixed configuration where both the neutron diffuseness and the neutron half-density radius are larger than the values of the proton density. We have to stress that even small differences between C n and C p may have a meaningful influence on the value of the neutron skin thickness. The discussed results indicate an intermediate character of the neutron skin thickness in 208 Pb, but with a preference for the density pattern where C n ≈ C p and a n > a p . A more precise conclusion is difficult because of the large experimental error bars. For comparison, in Fig. 1 we also plot the theoretical halo factor (3) in 208 Pb calculated directly from the neutron and proton densities of the Gogny D1S [54] and Skyrme SLy4 [55] effective nuclear forces. At the distances relevant for antiproton annihilation, the halo factor obtained with these theoretical mean-field densities agrees considerably well with the experimental data. The numerical values of the half-density radius C and of the diffuseness parameter a of the D1S and SLy4 equivalent 2pF neutron and proton densities for 208 Pb are reported in Table I. For either force, one observes that the halfdensity radii C n and C p of the 2pF distributions are different and that the values of a n and a p are also different. Thus, the results of these models correspond to some intermediate situation between the "neutron halo-type" 2pF distribution (where C n = C p and a n > a p ) and the "neutron skin-type" 2pF distribution (where C n > C p and a n = a p ) discussed in Ref. [38]. In any case, the two forces, D1S and SLy4, show some preference, as also do the calculations considered earlier, for the "neutron halo-type" distribution in 208 Pb, especially in the case of the D1S force, where C n − C p = 0.04 fm and a n − a p = 0.08 fm.
From the discussions of this section it is easily seen that the characterization of the "bulk" or "surface" formation of the neutron skin thickness in nuclei is a nontrivial problem. In the following, we wish to analyze this question according to the predictions derived from mean-field models of nuclear structure that are tested to be successful for charge radii, binding energies, and a wealth of phenomena in nuclei. First, we need to establish a prescription to separate bulk and surface contributions in neutron skins calculated with mean-field nuclear densities.

III. DISCERNING BULK FROM SURFACE IN NEUTRON SKINS
Nuclear radii can be defined in several independent ways. The most popular formulas and their relations are discussed in the book by Hasse and Myers [49] and we recall and analyze them in this section. One of the simplest ways to describe the size of nuclei is to define a central radius C in terms of the integral of the nuclear density profile ρ(r) (for either neutrons or protons) as Another option is the half-density radius R 1/2 which is defined from the local condition For density profiles that have a symmetric surface, such as the 2pF distribution defined in Eq. (4), the central radius C and the half-density radius R 1/2 coincide. The equivalent sharp radius R is the radius of a uniform sharp distribution that has a constant density equal to the bulk value of the actual density and that contains the same number of nucleons as the considered nucleus: For its importance in experimental techniques, the equivalent rms radius Q is commonly used. It describes a uniform sharp distribution with the same rms radius as the given density profile: Because the neutron skin thickness (1) is defined through rms radii, it can be expressed easily with Q: In uniform, sharp-edge nuclear distributions all of the above-mentioned definitions coincide, but in realistic leptodermous density profiles they give distinct values. The 132 Sn obtained in the RMF theory (the NL3 interaction [56] has been used).
relation between the C, Q, and R radii can be expressed by expansion formulas in powers of b/R, where b is the so-called surface width of the density profile. The latter quantity is defined by To leading order, one has the relationships [49] and These expansions are useful as long as b/R is small. This condition is well fulfilled in many nuclei because b ∼ 1 fm and R ∼ r 0 A 1/3 ; therefore, b 2 /R 2 is typically no larger than A −2/3 . Moreover, to consider no further corrective terms in the above relations of the C and Q radii with R may be quite accurate for specific shapes. For example, in the case of 2pF distributions, the first nonvanishing corrections to the terms within brackets in Eqs. (11) and (12) are of order (b/R) 6 and (b/R) 4 , respectively. The multiple definitions of nuclear radii may sometimes cause misleading conclusions, especially if nuclear properties sensitive to the nuclear matter distribution in nuclei are concerned or two different models are compared. Therefore one has to be careful about the suitable choice of the radius definition. To compare the foregoing definitions of radii for heavy nuclei, in Fig. 2 we plot the neutron densities obtained in a self-consistent meanfield calculation and the fitted 2pF distributions for the neutron-deficient nucleus 100 Sn and the neutron-rich nucleus 132 Sn. These profiles are compared with sharp-edge density distributions having radii C, R, and Q, calculated from the above expressions with the 2pF function. The central densities of the sharp surface spheres are fixed so to fulfill the particle number normalization. Figure 2 illustrates the fact that the central radius C does not allow the bulk density to be reproduced. A sharp sphere of radius C overestimates the self-consistent density in the whole nuclear interior. The equivalent rms radius Q also fails, because it clearly underestimates the original density in the bulk. Only the equivalent sharp radius R is able to properly reproduce the bulk part of the self-consistent and 2pF density profiles of the nucleus. As discussed in Ref. [49], the equivalent sharp radius R is the quantity of basic geometric importance of the three radii C, R, and Q. A sharp distribution of radius R has the same volume integral as the actual density of the finite nucleus and differs from it only in the surface region. Therefore, the radius R appears to be the suitable quantity to be used to measure the size of the bulk part of the nucleus.
On account of Eq. (9) for the neutron skin thickness and relationship (12) between Q and R, one obtains the expression up to terms of order O(b 4 /R 3 ). Thus, one can make a meaningful distinction between a bulk contribution and a surface (diffuseness) contribution to the neutron skin thickness of nuclei as follows: with independent of surface properties and The nucleus may develop a neutron skin by separation of the bulk radii R of neutrons and protons or by modification of the diffuseness b of the neutron and proton surfaces. In the general case a combination of both effects can be found. To which degree the different patterns arise in mean-field calculations of finite nuclei is a question we address in the following sections.
Experimental [50] and theoretical mean-field density distributions present oscillations in their inner bulk region, which implies that some suitable average is needed to determine the bulk density value ρ(bulk) of Eq. (7) to compute the equivalent sharp radius R. The difficulty may be easily solved by fitting a Fermi function to the original density, as we have illustrated in Fig. 2. In this case, the bulk density value ,that is, ρ 0 /[1 + exp (−C/a)], is, to excellent accuracy the ρ 0 parameter of the Fermi function. An effect not described by functions like the 2pF distribution is the occurrence of a nonsymmetric surface shape, which is possible in real nuclear densities. However, one plausibly expects that neutron skins of nuclei are dominated by the difference existing between neutrons and protons in the location of the surface (i.e., R) and/or in the spatial extent of this surface (i.e., b); the difference in the degree of asymmetry between the shapes of the neutron and proton surfaces is expected to be a less important, higher-order correction.
It may be practical to rewrite expressions (15) and (16) for the bulk and surface contributions to the neutron skin thickness directly in terms of the parameters of the 2pF function of Eq. (4). Using the fact that the diffuseness parameter a in a 2pF function is related to the surface width b by the formula and inverting relation (11) between C and R, one easily finds to the given order of approximation that Eqs. (15) and (16) become, respectively, and We note from these results that that is, the quantities C n − C p (easily obtained from the 2pF distributions themselves) and R n − R p [obtained through Eq. (7)] differ by surface diffuseness terms and should not be mixed. We have seen in Fig. 2 that a sharp sphere having a radius C significantly distorts the appearance of the actual nuclear density by overshooting it in the whole bulk region. It is thus preferable to use 3/5(R n −R p ) rather than 3/5(C n −C p ) as a measure of the bulk contribution to the neutron skin thickness.

A. Nuclei of the experiments in antiprotonic atoms
To get information about the "bulk" or "surface" character of the thickness of neutron skins in theoretical mean-field calculations, we parametrize the selfconsistently calculated proton and neutron densities with 2pF distributions. This procedure can be applied very suitably for many heavy nuclei and gives a clear distinction between bulk and surface properties of nuclei. However there is no universal method to do this parametrization. A popular prescription is to use a χ 2 minimization of the differences between the density to be reproduced and the 2pF profile, or of the differences between their logarithms. These methods may somewhat depend on conditions given during minimization (number of mesh points, limits, etc.). We prefer to extract the parameters of the 2pF profiles by imposing that they reproduce the same quadratic r 2 and quartic r 4 moments of the selfconsistent mean-field densities. These two conditions, together with the normalization to the proton and neutron numbers, allow us to determine in a unique way the equivalent 2pF densities. This method can be applied to any density distribution. Its focus is on a good reproduction of the surface region of the original density because the local distributions of the quantities r 2 ρ(r) and r 4 ρ(r) are peaked at the periphery of the nucleus. As an example of the results of the present determination of the 2pF profiles, in Fig. 3 we display in logarithmic scale the self-consistent neutron and proton densities for 208 Pb together with their equivalent 2pF distributions calculated with the NL3 interaction. It can be seen that there is overall good agreement in the central and surface regions of the nucleus. In particular, the 2pF densities reproduce the mean-field densities well at the distances that are relevant for antiproton annihilation.
In this section, we would like to get some insight about the "bulk" or "surface" character of the neutron skin predicted by the mean-field approach in the nuclei for which the neutron skin thickness values were extracted in Refs. [38,39] from the measurements in antiprotonic atoms. These nuclei range from 40 Ca to 238 U and all lie along the valley of stability. We carry out the calculations with the aforementioned nonrelativistic interactions D1S [54] and SLy4 [55] plus the relativistic interactions NL3 [56] and FSUGold [57], as representative examples of successful nuclear mean-field models. We impose spherical symmetry in all nuclei described in this article. The effect of deformation on the neutron skin thickness was discussed elsewhere [58]. The two nonrelativistic forces have a soft symmetry energy [15,18]. On the contrary, the covariant NL3 parameter set has a stiff symmetry energy, which is usual in the relativistic models [15,18]. Note that the notation soft or stiff refers to whether the symmetry energy of the model increases slowly or rapidly as a func- tion of the nuclear density. The covariant parameter set FSUGold was devised to have a softer density dependence of the symmetry energy [57] than the typical relativistic models. Thus, FSUGold is, in this aspect, closer to the nonrelativistic models than NL3. In Fig. 4 we display the experimental data with error bars determined from the antiprotonic atoms. There is a relatively clear correlation between the experimental value of the neutron skin thickness of these 26 stable nuclei and the overall relative neutron excess I = (N −Z)/A of the nucleus. This trend has been fitted by the linear relationship ∆R np = (0.90 ± 0.15)I + (−0.03 ± 0.02) fm with a χ 2 factor of 0.5 [38,39]. In the same figure we plot the theoretical ∆R np value calculated according to Eq. (1) using the mean-field densities of the indicated effective nuclear models for the 23 even-even nuclei that exist in the experimental data set. It is obvious that the theoretical models make largely different predictions for the neutron skin thickness. This is especially visible at large values of I. It is seen that the models that have a softer symmetry energy give smaller ∆R np values, whereas the models with a stiffer symmetry energy give larger ∆R np values [11,12]. Note that as soon as I = 0, even when it is small, discrepancies arise among the ∆R np values of the models. In contrast, all of the considered interactions make an almost identical prediction of ∆R np ≈ −0.05 fm for the 40 Ca nucleus.
One observes in Fig. 4 that, similarly to the situation in the experimental data set, the theoretical neutron skin values computed with each nuclear interaction show an average linear behavior as a function of the neutron excess I. Actually, the linear correlation factor of ∆R np with I in the present models is considerably high (between 0.95 and 0.97). The fit of the neutron skin values calculated with the Skyrme SLy4 force yields ∆R np = 1.01I − 0.035 fm, whereas ∆R np = 0.88I − 0.04 fm is obtained in the case of the Gogny D1S interaction. We find a linear fit of ∆R np = 1.20I − 0.03 fm using the values of the neutron skin thickness computed with the relativistic FSUGold model. If we consider the relativistic NL3 parameter set, which has a stiffer symmetry energy than the other models, we find the linear fit ∆R np = 1.50I − 0.03 fm. Compared with the slope of 0.90 for the fit of the experimental data [38,39], the slopes in forces like D1S and SLy4 is much closer to it, whereas the agreement deteriorates as the model has a stiffer symmetry energy. Thus, one can conclude that the comparison of the slopes of ∆R np with I between theory and the antiprotonic measurements for nuclei across the mass table favors the models that have a soft symmetry energy.
In Fig. 5 we display the bulk and surface contributions to the theoretical neutron skin thickness that are computed by applying Eqs. (15) and (16), respectively, against the neutron excess I. These values are obtained from the 2pF distributions associated with the mean-field neutron and proton densities calculated with the D1S force as well as with the Skyrme SLy4 force and the two RMF parameter sets FSUGold and NL3. The numerical calculations show that the value of ∆R bulk np + ∆R surf np obtained through Eqs. (18) and (19) can be slightly less accurate in some nuclei [compared to the exact value of Eq. (1)] than the result obtained through Eqs. (15) and (16). The small differences, whenever they arise, are due almost entirely to the replacement of a 2 q /R q by a 2 q /C q (q = n, p) in Eqs. (18) and (19). Of course, Eqs. (18) and (19) are quite practical because they do not require the additional calculation of the R q values and can be applied straightforwardly using the C q and a q parameters  15) and (16), respectively] for the nuclei presented in Fig. 4, calculated in various models.
of the 2pF distributions themselves.
The four panels presented in Fig. 5, as a consequence of the differences between models, predict a slightly different splitting of the neutron skin thickness into their bulk and surface contributions for each considered nucleus. Therefore, the values of the bulk and surface contributions to ∆R np are to some extent model dependent. The total neutron skin thickness roughly follows an increasing tendency with the relative neutron excess I, as discussed previously. This tendency is also seen in the bulk contribution ∆R bulk np , especially in the nonrelativistic interactions. However the surface contribution does not clearly follow the same trend and shows a less definite behavior as a function of I.
From Fig. 5 we see that for relatively neutron-rich nuclei (I 0.15), the surface part generally contributes 50% or more to the total neutron skin thickness. If the nuclear model has a stiff symmetry energy and I is large, the bulk part of ∆R np may become larger than the surface part, as can be seen in the NL3 panel of Fig. 5 at I > 0.2. More symmetric nuclei (I 0.15) do not present a definite tendency. The theoretical calculations show that in these nuclei the sharp radius is larger for protons than for neutrons (R p > R n ) while the surface width is larger for neutrons than for protons (b n > b p ). Consequently, the bulk contribution to the neutron skin thickness becomes negative, as it can be seen in Fig. 5, and the relatively small value of the neutron skin thickness is basically because of a strong cancellation between the bulk and surface parts. In the lightest nuclei, both contributions are negative and produce a "proton skin" rather than a neutron skin.

B. Medium and heavy mass isotopic chains
The set of nuclei chosen in Sec. IV A have given us some insight into the bulk and surface contributions to the neutron skin thickness in stable isotopes. To investigate in a systematic way how the bulk and surface contributions evolve with the neutron number we study the chains of even-even Sn and Pb isotopes from the proton to the neutron drip line.
First, Fig. 6 displays the neutron skin thickness ∆R np along the Sn and Pb isotopic chains computed in the mean-field approximation with the Gogny D1S and Skyrme SLy4 forces and using the relativistic mean-field models NL3 and FSUGold. For small and moderate values of the relative neutron excess (I 0.2), the neutron skin thickness grows almost linearly with I in each isotopic chain, as in the case of the stable nuclei analyzed in the previous section. However, ∆R np shows a rather pronounced kink at I ≈ 0.20 − 0.25. Beyond this value, it increases again almost linearly as a function of the relative neutron excess, but with a larger slope. Finally, a new departure from the linear behavior of ∆R np as a function of I can be observed when the isotopes are on the edge of the neutron drip line. The kinks and changes of slope in the neutron skin thickness as a function of the relative neutron excess are clearly connected with the doubly magic nuclei 132 Sn (I ≈ 0.24) and 208 Pb (I ≈ 0.21) in the stable nuclei region and with another two doubly magic nuclei, namely 176 Sn (I ≈ 0.43) and 266 Pb (I ≈ 0.38), near the neutron drip line. Therefore, it is obvious that these changes of slope are produced by quantal effects, which modify the linear trend of ∆R np as a function of I in a non-negligible way. The average slope of ∆R np with I for various forces is clearly different. As has been shown in Refs. [11,12], for each nuclear model, ∆R np is strongly correlated with the slope of the symmetry energy with respect to the density computed at saturation. As in Sec. IV A, we fit the mean-field proton and neutron densities by 2pF distributions to investigate the bulk and surface contributions to the neutron skin. To analyze the changes that occur in ∆R np first we look at the parameters that characterize the 2pF distributions. In Fig. 7 we display the central density ρ 0 , the half-density radius C and the diffuseness a, which are obtained from the D1S densities along the Sn and Pb isotopic chains. The 2pF parameters show an overall smooth behavior as a function of the relative neutron excess I but with local modulations near the shell closures. The central density ρ 0 grows (for neutrons) or declines (for protons) almost linearly as a function of I in both isotopic chains. It is just a simple consequence of the increasing asymmetry along the isotopic chains. The central radii C n and C p both show a global increasing tendency with the neutron number. One can notice that C n grows faster with I than C p . However, the a n and a p diffuseness parameters behave in a completely different way. Whereas a n grows on average with increasing neutron excess, a p remains roughly constant with I in the two isotopic chains. On top of these general trends we see that the 2pF parameters associated with the neutron densities show kinks around the shell closures at N = 82 in Sn and N = 126 in Pb (that is, in the region I ≃ 0.20 − 0.25), as well as changes of slope near the neutron drip lines. One observes that some signature of the neutron shell effect also appears in the 2pF parameters of the proton densities around magic neutron numbers. Figure 8 shows the evolution with I of the 2pF parameters corresponding to the quantal densities computed with the SLy4 Skyrme force and with the relativistic mean-field NL3 parameter set. In the two models the global trends are similar to the ones discussed before for 2pF parameters of D1S distribution. However, they show some differences that can be attributed in part to the different properties of the isovector channel of the interaction.
In Figs. 9(a) and 9(b) we display the bulk and the surface contributions to the nuclear skin thickness calculated with the D1S force along the Sn and Pb isotopic chains. These contributions show for Sn isotopes two well-defined regions as a function of I. One of them covers the neutron major shell between 100 Sn and 132 Sn, i.e. 0 ≤ I 0.25, and the other region corresponds to the next major shell between 132 Sn and 176 Sn in the range 0.25 I 0.43. For nearly symmetric Sn isotopes close to 100 Sn, which are neutron deficient, C p is larger than C n (see Fig. 7) and, consequently, the bulk part of the neutron skin is negative or at most it takes very small positive values.
For these values of I, the surface contribution is positive and relatively small, reducing the opposite effect of "proton skin" due to negative values of the bulk part. Looking at more asymmetric isotopes, one can see that, in the magic nuclei 132 Sn and 176 Sn, the bulk and the surface contributions are roughly equal. This implies a rather compact neutron density distribution with a relatively stiff surface as a consequence of the kinks exhibited by the neutron diffuseness parameter a n around the neutron magic numbers, which can be seen in Fig. 7. In the regions between magic numbers, the surface contribution to the neutron skin thickness is larger than the bulk contribution. The splitting between the surface and bulk contributions reaches its maximum value roughly at midshell. This behavior of the bulk and surface contributions to the neutron skin thickness points out that Sn isotopes in the middle of major shells develop a larger surface region in the density distribution than the magic ones. In other words, these isotopes with neutron number in between magic values are more of "halo" type than the limiting magic nuclei which show a mixed character between "halo" and "neutron skin." The bulk and the surface contributions to the neutron skin of the Pb isotopes calculated with the D1S force show a similar behavior to the case of the Sn isotopes analyzed before. However, for this heavy isotopic chain, the bulk part gives a more important contribution to the total neutron skin for nuclei in between the magic 208 Pb and 266 Pb nuclei.
The discussed differences in the behavior of the neutron skin and its bulk and surface contributions along the neutron-rich Sn and Pb isotopes can be qualitatively understood as follows. Within a major shell, the rms radii of the different single-particle orbits are spread around their average value. The rms radius of orbits with low angular momentum l are larger than the average, while the rms radius of orbits with high l are slightly smaller but close to the average in the shell. Consequently, the outermost region of the density is basically provided by the orbits with low l in the last populated major shell. On the contrary, orbits with high l in this major shell have their most important contribution at shorter distances from the center of the nucleus, increasing more the bulk than the surface part of the nuclear density. Hence, the filling order of the last single-particle orbits is crucial to determine if the growing of the neutron radius, and consequently the neutron skin, is due to an increase of the surface or the bulk of the density. In the case of neutron-rich isotopes of Sn above N = 82, the first filled orbits are 2f 7/2 , 3p 3/2 and 3p 1/2 . This ordering produces an important enhancement of the nuclear surface, which can be appreciated from Figs. 7 and 9. Once midshell is filled, the 1h 9/2 and 1i 13/2 orbits start to be appreciably populated, increasing the bulk more than the surface of the densities and hence of the neutron skin, as can be seen in the aforementioned figures. The situation in Pb isotopes is just the contrary; above N = 126, the first occupied levels are 2g 9/2 and 1j 15/2 , which increase the bulk more than the surface. Only near the drip line, N = 184, are the low-momentum orbits 3d 5/2 , 4s 1/2 , and 3d 3/2 relevant, increasing the surface and quenching the bulk contributions to the density and consequently of the neutron skin. This behavior in Pb isotopes can also be seen in Figs. 7 and 9.
In Figs. 9(c)-9(f) we display the results calculated with the Skyrme SLy4 force and the NL3 parameter set. The bulk and surface contributions qualitatively behave as the ones computed with the D1S force discussed previously. Of course, some differences are found in comparing the results obtained with the different models because of their different isovector properties. In models with a stiff density dependence of the symmetry energy (for instance, NL3), the bulk contribution to the neutron skin is more important than in models with a soft symmetry energy, as in the case of the SLy4 and D1S forces. This tendency can be especially noted in the Pb isotonic chain.

V. CONCLUSIONS
Ground-state properties of stable nuclei, such as charge radii and binding energies, can be reproduced fairly well by using mean-field models with effective interactions such as the Skyrme or Gogny forces and with relativistic Lagrangians. Although the isoscalar part of these models is well constrained, their isovector properties are much less determined and the predictions in this sector can differ considerably, even for stable nuclei. A typical example is the neutron skin thickness of nuclei. For this observable theoretical predictions for 208 Pb using nonrelativistic forces give a value around 0.15 fm while relativistic mean-field parametrizations predict values that are almost twice as large.
The neutron skin thickness of a set of 26 nuclei, distributed over the whole periodic table, has been obtained from the analysis of experiments with antiprotonic atoms [38,39] combined with the charge radii obtained from electron scattering data. One important result of these experiments with antiprotonic atoms is that there exists a rather clear linear correlation between the neutron skin thickness and the overall neutron excess I. Theoretical mean-field calculations of the neutron skin also show this tendency even more clearly. It is found that the relativistic parametrizations systematically predict larger neutron skins than the ones computed with nonrelativis-tic interactions. This is because the symmetry energy is stiffer in the relativistic models than in the nonrelativistic models.
To analyze the experimental data of antiprotonic atoms an ansatz of the nuclear densities is needed. The two-parameter Fermi distributions have been used often to this end. It is found that the experimental data can be reproduced by a variety of these Fermi distributions with different values of C n −C p and a n −a p . The experimental values of the halo factor in 208 Pb are well reproduced by distributions with C n ≈ C p (halo model) and by Fermi densities with both halo and neutron skin (a n ≈ a p ) contributions. This latter scenario is also well predicted by the nonrelativistic mean-field densities obtained with the D1S and SLy4 forces.
We have also parametrized the mean-field densities via two-parameter Fermi distributions. We do this by imposing that both mean-field and parametrized densities give the same quadratic and quatric moments. This parametrization of the mean-field densities also allows the neutron skin thickness to be split easily into two contributions, namely the bulk part and the surface part. It is found that the mean-field neutron skins computed in nuclei with I > 0.1 can be shared between non-negligible surface and bulk parts. This applies both for stable nuclei investigated in antiprotonic experiments and for drip line isotopes in all the theoretical models considered in this work.
To analyze the neutron skin in neutron-rich nuclei, we have theoretically studied its variation along the Sn and Pb isotopic chains up to the neutron drip line using selected mean-field models. As expected, ∆R np shows generally linear growing trend with I. However, shell effects, which are always present in mean-field calculations, produce noticeable departures from this linear dependence in nuclei with large neutron excesses.
Regarding the bulk and surface contributions to the neutron skin thickness in Sn isotopes, it can be seen that the considered mean-field models point toward more of a surface character in stable nuclei. This effect is reinforced in the neutron-rich region. In the case of the Pb isotopic chain, bulk and surface contributions have similar values in stable isotopes, whereas the bulk part is larger than the surface part in the more neutron-rich region of this isotopic chain.