Pairing Properties In Relativistic Mean Field Models Obtained From Effective Field Theory

We apply recently developed effective field theory nuclear models in mean field approximation (parameter sets G1 and G2) to describe ground-state properties of nuclei from the valley of $\beta$-stability up to the drip lines. For faster calculations of open-shell nuclei we employ a modified BCS approach which takes into account quasi-bound levels owing to their centrifugal barrier, with a constant pairing strength. We test this simple prescription by comparing with available Hartree-plus-Bogoliubov results. Using the new effective parameter sets we then compute separation energies, density distributions and spin--orbit potentials in isotopic (isotonic) chains of nuclei with magic neutron (proton) numbers. The new forces describe the experimental systematics similarly to conventional non-linear $\sigma-\omega$ relativistic force parameters like NL3.


Introduction
The relativistic field theory of hadrons known as quantum hadrodynamics (QHD) has become a very useful tool for describing bulk and single-particle properties of nuclear matter and finite nuclei in the mean field approximation [1,2,3,4]. Compared with the nonrelativistic approach to the nuclear many-body problem, the relativistic model explicitly includes the mesonic degrees of freedom and treats the nucleons as Dirac particles. At the mean field (Hartree) level, nucleons interact in a relativistic covariant way by exchanging virtual mesons: an isoscalar-vector ω meson, an isoscalar-scalar σ meson and an isovectorvector ρ meson. With these ingredients the mean field treatment of QHD automatically takes into account the spin-orbit force, the finite range and the density dependence of the nuclear force. Adjusting some coupling constants and meson masses from the properties of a small number of finite nuclei, the relativistic mean field (RMF) model produces excellent results for binding energies, root-mean-square radii, quadrupole and hexadecapole deformations and other properties of spherical and deformed nuclei [5,6].
The original linear σ − ω model of Walecka [7] was complemented with cubic and quartic non-linearities of the σ meson [8] (non-linear σ − ω model) to improve the results for the incompressibility and for finite nuclei. Since these models were proposed to be renormalizable, the scalar self-interactions were limited to a quartic polynomial and scalar-vector or vectorvector interactions were not allowed. Recently, and inspired by effective field theory (EFT), Furnstahl, Serot and Tang [9,10] abandoned the idea of renormalizability and extended the RMF theory by including other non-linear scalar-vector and vector-vector self-interactions as well as tensor couplings [4,9,10,11,12,13].
The EFT Lagrangian has an infinite number of terms since it contains all the nonrenormalizable couplings consistent with the underlying QCD symmetries. Therefore, it is mandatory to develop a suitable scheme of expansion and truncation. At normal nuclear densities the scalar (Φ) and vector (W ) meson fields are small compared with the nucleon mass (M), and they vary slowly with position in finite nuclei. This indicates that the ratios Φ/M, W/M, |∇Φ|/M 2 and |∇W |/M 2 can be used as the expansion parameters. With the help of the concept of naturalness, it is then possible to compute the contributions of the different terms in the expansion and to truncate the effective Lagrangian at a given level of accuracy [4,10,12,13]. None of the couplings should be arbitrarily dropped out to the given order without a symmetry argument.
References [10,12,13] have shown that it suffices to go to fourth order in the expansion.
At this level one recovers the standard non-linear σ−ω model plus a few additional couplings, with thirteen free parameters in all. These parameters have been fitted (parameter sets G1 and G2) to reproduce some observables of magic nuclei [10]. The fits display naturalness (i.e., all coupling constants are of the order of unity when written in appropriate dimensionless form), and the results are not dominated by the last terms retained. This evidence confirms the utility of the EFT concepts and justifies the truncation of the effective Lagrangian at the first lower orders.
Recent applications of the models based on EFT include studies of pion-nucleus scattering [14] and of the nuclear spin-orbit force [15], as well as calculations of asymmetric nuclear matter at finite temperature with the G1 and G2 sets [16]. In a previous work [17] we have analyzed the impact of each one of the new couplings introduced in the EFT models on the nuclear matter saturation properties and on the nuclear surface properties. In Ref. [18] we have looked for constraints on the new parameters by demanding consistency with DBHF calculations and the properties of finite nuclei. During the last years a large amount of work has been devoted to measuring masses of nuclei far from stability [19]. This body of experimental data has been used as a benchmark to test the predictions of the currently existent (relativistic and non-relativistic) nuclear effective forces [20]. This fact motivates us to investigate in the present work the behaviour of the parameter sets G1 and G2 derived from EFT in regions far from the stability line.
To study ground-state properties of spherical open-shell nuclei one has to take into account the pairing correlations. Relativistic mean field calculations near the β-stability line have usually included pairing in a constant gap BCS approximation [5,21,22], with the gaps fitted to empirical odd-even mass differences. This approach works properly when the main effect of the pairing correlations is a smearing of the Fermi surface. Since the BCS pairing energy diverges for large momenta, a cut-off has to be introduced in the pairing 3 channel to simulate phenomenologically the finite range of the particle-particle force. The limitations of this simple BCS method appear when one deals with nuclei far from the βstability line. Close to the drip lines the Fermi level falls near the particle continuum and it is known that the BCS model does not provide a correct description of the coupling between bound and continuum states [23,24]. In the non-relativistic framework this difficulty was overcome by the unified description of the mean field and the pairing correlations provided by the Hartree-Fock-Bogoliubov (HFB) theory [25,26], with Skyrme [23,24] or Gogny forces [27].
The same unified treatment was developed by Kucharek and Ring [28] in the relativistic framework. However, a quantitative description of the pairing correlations in nuclei cannot be achieved with relativistic mean field parametrizations because the meson exchange forces are not properly adapted to large momentum transfer [28,29]. Later, Ring and coworkers [29,30,31,32] have used the RMF interaction for the particle-hole channel plus the pairing part of the Gogny force [27] (with the D1S parameters [33]) for the particle-particle channel, in relativistic Hartree-plus-Bogoliubov (RHB) calculations. Other authors have employed a density-dependent zero-range pairing force [34] instead of the Gogny pairing force [35,36].
Recent calculations with non-relativistic Skyrme forces and a zero-range force in the particle-particle channel have shown that a BCS approach is able to provide a good qualitative estimate of the drip lines if some quasibound states due to their centrifugal barrier (plus the Coulomb barrier for protons) are included in the calculation [37,38,39]. In this work we will use a similar BCS approach with quasibound states to approximately take into account the effects of the continuum contributions near the drip lines. We will employ a constant pairing strength which can be considered as a simplification of the zero-range pairing force and which gives similar results to those obtained with a delta force for spherical nuclei [40].
The paper is organized as follows. We summarize the mean field approximation to the EFT nuclear model in the second section. In the third section we describe our modified BCS approach with quasibound states, and perform some calculations to test its possibilities and limitations by comparing with Bogoliubov results available from the literature. The fourth section is devoted to the detailed study with the EFT parametrizations G1 and G2 of properties such as separation energies, particle densities and spin-orbit potentials of nuclei belonging to chains of isotopes (isotones) with magic proton (neutron) number. Our conclusions are laid in the last section.
2 Relativistic mean field approach from effective field theory The effective field theory approach to QHD has been developed in the recent years. The theory and the equations for nuclear matter and finite nuclei can be found in the literature [4,9,10] and here we shall only outline the formalism. We start from Ref. [9] where the field equations were derived from an energy density functional containing Dirac baryons and classical scalar and vector mesons. This functional can be obtained from the effective Lagrangian in the Hartree approximation, but it can also be considered as an expansion in terms of the ratios of the meson fields and their gradients to the nucleon mass of a general energy density functional that contains the contributions of correlations within the spirit of density functional theory [4,10].
According to Refs. [4,10] the energy density for finite nuclei can be written as (r) and A ≡ eA 0 (r). Variation of the energy density (1) with respect to ϕ † α and the meson fields gives the Dirac equation fulfilled by the nucleons and the meson field equations, which are solved self-consistently by numerical iteration. We refer the reader to Ref. [10] for the expressions of the variational equations.
The terms with g γ , λ, β s and β v take care of effects related with the electromagnetic structure of the pion and the nucleon (see Ref. [10]). Specifically, the constant g γ concerns the coupling of the photon to the pions and the nucleons through the exchange of neutral vector mesons. The experimental value is g 2 γ /4π = 2.0. The constant λ is needed to reproduce the magnetic moments of the nucleons. It is defined by with λ p = 1.793 and λ n = −1.913 the anomalous magnetic moments of the proton and the neutron, respectively. The terms with β s and β v contribute to the charge radii of the nucleon [10].
In this work we will employ the EFT parameter sets G1 and G2 of Refs. [4,10]. The , ζ 0 , f v , α 1 and α 2 were fitted by a least-squares optimization procedure to twenty-nine observables (binding energies, charge form factors and spin-orbit splittings near the Fermi surface) of the nuclei 16 O, 40 Ca, 48 Ca, 88 Sr and 208 Pb, as described in Ref. [10].
The constants β s , β v and f ρ were then chosen to reproduce the experimental charge radii of the nucleon. The fits yielded two best, distinct parameter sets (G1 and G2) with essentially the same χ 2 value [10].
We report in Table 1 the values of the parameters and the saturation properties of G1 and G2. One observes that the fitted parameters differ significantly between both interactions.
For example, G2 presents a positive value of κ 4 , as opposite to G1 and to many of the most successful RMF parametrizations, such as the NL3 parameter set [41]. Formally a negative value of κ 4 is not acceptable because the energy spectrum then has no lower bound [42].
Furthermore, the wrong sign in the Φ 4 coupling constant may cause troubles in obtaining stable solutions in light nuclei like 12 C. We note that the value of the effective mass at saturation M * ∞ /M in the EFT sets (∼ 0.65) is somewhat larger than the usual value in the RMF parameter sets (∼ 0.60). This fact is related with the presence of the tensor coupling f v of the ω meson to the nucleon, which has an important bearing on the spin-orbit force [10,15,17].
One should mention that the EFT perspective also has been helpful to elucidate the empirical success of the usual non-linear σ − ω models that incorporate less couplings (just up to cubic and quartic self-interactions of the scalar field): the EFT approach accounts for the success of these RMF models and provides an expansion scheme at the mean field level and for going beyond it [4,10,12]. In practice it has been seen that the mean field phenomenology of bulk and single-particle nuclear observables does not constrain all of the new parameters of the EFT model unambiguosly. That is, the constants of the EFT model are underdetermined by the observables currently included in the fits and different parameter sets with low χ 2 (comparable to G1 and G2) can be found [10,12,13,14]. However, the extra couplings could prove to be very useful for the description of further observables. Indeed, for densities above the normal saturation density, and owing to the additional non-linear couplings, the EFT models are able [18] to give an equation of state and nuclear matter scalar and vector self-energies in much better agreement with the microscopic Dirac-Brueckner-Hartree-Fock predictions than the standard non-linear σ − ω parametrizations (the latter completely fail in following the DBHF trends as the nuclear density grows [18,22]).
The sets G1 and G2 were fitted including centre-of-mass corrections in both the binding energy and the charge radius. Therefore, we will utilize the same prescription of Ref. [10] in our calculations with G1 and G2. Namely, a correction to the binding energy and a correction to the mean-square charge radius. 7

The pairing calculation
It is well known that pairing correlations have to be included in any realistic calculation of medium and heavy nuclei. In principle the microscopic HFB theory should be used for this purpose. However, for pairing calculations of a broad range of nuclei not too far from the β-stability line, a simpler procedure is usually considered in which a seniority potential acts between time-reversed orbitals. In this section we want to discuss and test a straightforward improvement of this simple approximation to be able to describe in addition nuclei near the drip lines, at least on a qualitative level. Without the complications intrinsic to a full Bogoliubov calculation, our faster approximation will allow us later on to perform extensive calculations of chains of isotopes and isotones with the relativistic parameter sets.
The pairing correlation will be considered in the BCS approach [25,26]. One assumes that the pairing interaction v pair has non-zero matrix elements only between pairs of nucleons invariant under time reversal: where |α = |nljm and |α = |nlj − m (with G > 0 and m > 0). Most often the BCS calculations in the RMF model have been performed using a constant gap approach [5,21,22]. Instead, here we choose a seniority-type interaction with a constant value of G for pairs belonging to the active pairing shells.
The contribution of the pairing interaction to the total energy, for each kind of nucleon (neutrons or protons), is where n α is the occupation probability of a state with quantum numbers α ≡ {nljm} and the sum is restricted to positive values of m. One has The Lagrange multiplier µ is called the chemical potential and the gap ∆ is defined by 8 As usual the last term in Eq. (6) will be neglected. It is not a very important contribution and its only effect is a renormalization of the pairing energies [25,26].
Assuming constant pairing matrix elements (5) in the vicinity of the Fermi level one gets [25,26] where A is the number of neutrons or protons involved in the pairing correlation. The solution of these two coupled equations allows one to find µ and ∆. Using Eqs. (7) and (8) the pairing energy for each kind of nucleon can be written as This simple approach breaks down for nuclei far from the stability line. The reason is that in this case the number of neutrons (for isotopes) or protons (for isotones) increases, the corresponding Fermi level approaches zero and the number of available levels above it is clearly reduced. Moreover, in this situation the particle-hole and pair excitations reach the continuum. Ref. [23] showed that if one performs a BCS calculation using the same quasiparticle states as in a HFB calculation, then the BCS binding energies are close to the HFB ones but the r.m.s. radii (i.e., the single-particle wave functions) dramatically depend on the size of the box where the calculation is performed. This is due to the fact that there are neutrons (protons) that occupy continuum states for which the wave functions are not localized in a region, thus giving rise to an unphysical neutron (proton) gas surrounding the nucleus.
Recent non-relativistic calculations near the drip lines with Skyrme forces [38,39] have shown that the above problem of the BCS approach can be corrected, in an approximate manner, by taking into account continuum effects by means of the so-called quasibound states, namely, states bound because of their own centrifugal barrier (centrifugal-plus-Coulomb barrier for protons). When the quasibound states are included in the BCS calculation (from now on a qb-BCS calculation), it is necessary to prevent the unrealistic pairing of highly excited states and to confine the region of influence of the pairing potential to the vicinity of the Fermi level. Instead of using a cutoff factor as in Ref. [38], in our calculations we will restrict the available space to one harmonic oscillator shell above and below the Fermi level.
In order to check this approach we have performed with the G1 parameter set (G n = 21/A MeV, see next section) calculations of the binding energy and r.m.s. radius of the 120 Sn and 160 Sn nuclei in boxes of sizes between 15 and 25 fm (as in the non-relativistic calculations of Ref. [23]). The results taking into account the quasibound levels 1h 9/2 , 2f 5/2 and 1i 13/2 for 120 Sn, and 1i 11/2 and 1j 13/2 for 160 Sn, are compared in Figure 1 with the output of a standard BCS calculation with only bound levels. It turns out that in the qb-BCS case the results are essentially independent of the size of the box where the calculations are carried out. When the quasibound levels are included the binding energies are larger than when only the bound levels are taken into account, due to the damping of the pairing correlation caused by disregarding the continuum states in the standard BCS calculation [23]. We also show in Figure 1 the results of a BCS calculation using all bound and unbound levels (i.e., without restricting ourselves to quasibound levels) in the considered range. It is obvious that in this case the results are box dependent, as the binding energy and neutron r.m.s. radius of 160 Sn evidence.
Another test of the qb-BCS approach concerns the asymptotic behaviour of the particle densities [24]. In Figure 2 we display the radial dependence of the neutron density of 150 Sn (as in Ref. [24]) calculated with the G1 parameter set in boxes of radii between 15 and 25 fm. For large enough distances the density decreases smoothly when the size of the box increases (except very near of the edge, where the density suddenly drops to zero because of the ϕ α = 0 boundary condition). This means that no neutron gas surrounding the nucleus has appeared. In a Bogoliubov calculation the asymptotic behaviour of the particle density is governed by the square of the lower component of the single-quasiparticle wave function corresponding to the lowest quasiparticle energy [24]. This asymptotic behaviour is displayed by the (almost straight) dotted line in Figure 2. It can be seen that the density obtained with our approach decreases more slowly than the RHB density, i.e., asymptotically the qb-BCS density is not able to follow the RHB behaviour. This coincides with the conclusion of Ref. [24] (see Figure 19 of that work) where non-relativistic HFB densities are compared for large distances with the densities obtained in the qb-BCS approach with a state-dependent pairing [37].
Although the qb-BCS densities do not display the right asymptotic behaviour, it was conjectured in Ref. [24] that such an approach could allow one to compute properties of nuclei much closer to the drip lines than in a standard BCS calculation. Very recently, RHB calculations up to the drip lines of the two-neutron separation energy S 2n for nickel isotopes [35] and of the charge and neutron r.m.s. radii for tin isotopes [36] have been carried out using the NL-SH parameter set [43] plus a density-dependent zero-range pairing force. We have repeated these calculations with our qb-BCS method for both isotopic chains (with a pairing interaction strength G n = 22.5/A in the case of NL-SH).
We display the values of the S 2n separation energies for the Ni chain in Figure 3a. The can be compared with the RHB values displayed in Figure 2 of Ref. [35]. The tendencies are the same, though the qb-BCS pairing energies are slightly larger than in the RHB calculation. In Figure 4 we draw our results for the radii of the Sn isotopes, and compare them with the RHB values. In the case of the charge radii the agreement is excellent. The neutron radii obtained in our method closely follow the behaviour of the RHB neutron radii and the kink at N = 132 is qualitatively reproduced.
We have furthermore computed the binding energies of nuclei of the N = 20 isotonic chain for which RHB results exist with the NL3 parameter set [31]. We present the extracted two-proton separation energies S 2p in Table 2. The agreement between the qb-BCS and RHB approaches again is very good. In both models the last stable nucleus is 46 Fe, as in experiment. Notice that in the present case the first levels with positive energy correspond to those of the pf shell. Due to the Coulomb barrier all these levels become quasibound in our approach, and it is expected that they will lie close to the canonical levels. This explains the goodness of the qb-BCS energies for this isotonic chain.
¿From the previous comparisons we see that the simple qb-BCS calculation is able to reasonably follow the main trends of the more fundamental RHB pairing calculation. One can also conclude that the consideration of quasibound states in the BCS approach is, actually, a key ingredient to eliminate the spurious nucleon gas arising near the drip lines.

Results for EFT parameter sets
We want to analyze the ability of the G1 and G2 parameter sets based on effective field theory [4,10] to describe nuclear properties far from the stability line, i.e., far from the region where the parameters were fitted. To our knowledge such calculations have not been explored so far. We will contrast the results with experiment and with those predicted by the NL3 set, that we take as one of the best representatives of the usual RMF model with only scalar self-interactions.
As indicated, we shall use a schematic pairing with a state-independent matrix element MeV and C p = 25 MeV for G2, and finally C n = 20.5 MeV and C p = 23 MeV for NL3. Figure 5 shows that the neutron and proton state-independent gaps (∆ n and ∆ p ) predicted by our calculation with G1 are scattered around the empirical average curve 12/ √ A [44].
A similar picture is found with the parameter sets G2 and NL3.

Two-particle separation energies
In Figure 6a we present the two-neutron separation energies S 2n for the chain of Ni isotopes.
Clear shell effects arise at N = 28 and 50. The three relativistic interactions (G1, G2 and NL3) slightly overestimate the shell effect at N = 28 as compared with the experimental value, which also happens in more sophisticated RHB calculations with NL3 [30,32]. In our qb-BCS approach some disagreement with experiment is found for the N = 38 and N = 40 isotopes. Again, this also occurs in the RHB calculations of Refs. [30,32] with NL3.
However, if we compare Figure 6a with the results that we have shown in Figure 3a for the NL-SH parameter set, we see that NL-SH achieves a better agreement with experiment for these N = 38 and N = 40 isotopes.
We stop our calculation towards the neutron drip line when the two-neutron separation energy vanishes or when the neutron chemical potential becomes positive. The fact that S 2n is not always zero at the drip line is connected with the quenching of the shell structure with N, which is a force-dependent property [24]. This effect is illustrated in Figure 25 of Ref.
[24] for HFB calculations with different non-relativistic forces. We find similar situations with the considered relativistic sets in our qb-BCS calculations of separation energies. In the case of the Ni isotopes we reach the drip line at N = 66 with the G1 and NL3 sets and at N = 68 with the G2 set. This agrees nicely with the value N = 66 obtained in HFB calculations with the Skyrme forces SIII [34,45] and SkP [34]. For NL-SH our qb-BCS scheme predicts the drip line at N = 72 (see Figure 3a), the same value found in the RHB calculations of Ref. [35].
In Figure 7a we display our qb-BCS results for the two-neutron separation energies of the Sn isotopic chain. In Ref. [32] it was claimed that pure BCS calculations in the constant gap approach (with NL3) are not suitable for the Sn isotopes. We observe in Figure 7a that below N = 60, as one moves towards N = 50, some discrepancies with the experimental values appear, which also arise in the RHB calculations [32]. The three forces slightly overestimate the shell effect at N = 82 (as the RHB results of Refs. [30,32] for NL3). We have computed Sn isotopes up to A = 176, when S 2n vanishes for NL3 (in good agreement with RHB results for NL-SH [36] and HFB results for the Skyrme force SkP [23]). For G1 and G2 we find that S 2n does not yet vanish at N = 126, and it is not possible to increase the neutron number due to the shell closure at N = 126 (the neutron chemical potential becomes positive for the N = 128 isotope). This means that the quenching of the shell effect at N = 126 for NL3 (and NL-SH) is larger than for the G1 and G2 parameter sets.
Our calculated S 2n energies for Pb isotopes are shown in Figure 8. Ref. [38]. It is slightly larger for G2 than for G1 and NL3. Experimental information for this shell effect is not available. NL3 would predict another shell effect at Z = 58, which does not appear experimentally. The effect is less pronounced in G1 and it does not show up in G2. The three forces indicate that the proton drip line is reached after the 156 W isotope, in agreement with experimental information [46]. Figures 11a and 11b show, respectively, the calculated S 2p separation energies for the N = 50 and N = 126 isotone chains. Note that we did not use any information about these nuclei in our fit of the G p pairing strength. For N = 50 the set G2 follows the experimental data very well, specially for the larger Z. The trend of G1 and NL3 is only a little worse.
The proton drip line is located at 100 Sn in the three parametrizations, in good accordance with experiment. The quenching of the shell effect at Z = 50 is larger for G2 than for G1 and NL3. The available data for two-proton separation energies of N = 126 isotones are reasonably well estimated by the relativistic sets. However, the trend of NL3 is worse than that of G1 and G2. It would then be very interesting to perform RHB calculations of this chain to confirm the behaviour of NL3. The last nucleus of the chain stable against two-proton emission is 218 U according to G1 and NL3, and 220 Pu according to G2. The three sets predict a shell effect at Z = 92, though it is relatively quenched for G2.

One-particle separation energies
We have computed one-neutron (one-proton) separation energies for Ni and Sn isotopes  [25,26]. In the spherical approximation one replaces the blocked single-particle state by an average over the degenerate states in its j shell. This way the rotational and time-reversal invariance of the many-body system is restored in the intrinsic frame [47]. In this approach the contribution of the j shell that contains the blocked state to the number of active particles and the pairing energy is respectively. The remaining active shells contribute in the usual manner [Eqs. (6) and (10)].
Due to rearrangement effects, blocking the single-particle state with smallest quasiparticle energy E α = (ε α − µ) 2 + ∆ 2 in the even A − 1 nucleus, does not necessarily lead to the largest binding energy of the odd A nucleus. Therefore, in some cases one has to repeat the calculation blocking in turn the different single-particle states that lie around the Fermi level to find the configuration of largest binding energy [23,27,47].  (Figure 10a), NL3 predicts a shell effect at Z = 58 which is not found experimentally, whereas for G1 this effect is clearly smaller and it does not appear for G2. The last stable nucleus against one-proton emission is 151 Tm according to the three parameter sets.

One-body densities and potentials
The nuclear densities for chains of isotopes of light and medium size nuclei have recently been studied in the RHB approximation [30,31,35,36]. As N grows the neutron and mass densities extend outwards and the r.m.s. radii and the surface thickness increase. Special attention has been payed to the isospin dependence of the spin-orbit interaction. The magnitude of the spin-orbit potential is reduced when one approaches the neutron drip line and, as a consequence, there is a reduction of the energy splittings between spin-orbit partner levels [30,31,36]. To our knowledge, for isotones such an study has only been carried out in the N = 20 chain [31]. It is to be remarked that the EFT parametrizations G1 and G2 contain a tensor coupling of the ω meson to the nucleon which plays a very important role in the spin-orbit force because there exists a trade-off between the size of this coupling and the size of the scalar field [15,17].
In Figures 12a and 12b we display, respectively, the neutron and proton densities of some N = 28 isotones from Z = 16 to Z = 32 as predicted by the G2 set in our qb-BCS approach. Figures 13a and 13b show the results for some N = 82 isotones from Z = 40 to Z = 70.
Since N is fixed in each isotonic chain, the spatial extension of the neutron densities is very similar for the different nuclei of the chain. In any case, as one goes from the lightest to the heaviest isotone of the chain, the neutron densities tend to be depressed in the interior region and their surface thickness (90%-10% fall-off distance) shows a decreasing tendency.
The proton densities of the isotones exhibit a strong dependence on Z: by adding more This behaviour may be related with the shell effect for protons at Z = 50.
The spin-orbit interaction is automatically included in the RMF approximation. It appears explicitly when the lower spinor of the relativistic wave function is eliminated in favour of the upper spinor. This way one obtains a Schrödinger-like equation with a term V SO (r) that has the structure of the single-particle spin-orbit potential. Including the contribution of the tensor coupling of the ω meson, the spin-orbit term reads [15,30] where M = M − 1 2 (Φ + W ). We have checked numerically that the contribution to the spin-orbit potential of the f ρ tensor coupling of the ρ meson is very small, even when one approaches the drip lines. Hence we have not written this contribution in Eq. (15).
The spin-orbit potential (15) for some lead isotopes computed with G2 and NL3 is displayed in Figures 15a and 15b, respectively. As a general trend, for both G2 and NL3, when the number of neutrons is increased the depth of the spin-orbit potential decreases gradually and the position of the bottom of the well is shifted outwards, which implies a significant weakening of the spin-orbit interaction. The same effect arises in other isotopic chains in RHB calculations [30,35,36]. Comparing the spin-orbit potentials obtained with the G2 and NL3 sets, one sees that they have a similar strength for all the isotopes analyzed and that the minima of the wells are located at similar positions (slightly shifted to larger values of r in G2). The higher effective mass of G2 at saturation (M * ∞ /M = 0.664) with respect to NL3 (M * ∞ /M = 0.595) is compensated by the tensor coupling included in G2 (f v = 0.692). To ascertain the relative importance of the tensor coupling we have drawn in the insert of Figure 15a, for 228 Pb, the full potential (15) and the contribution resulting from setting f v = 0 in Eq. (15). We see that the full V SO (r) is much deeper and wider. The maximum depth of V SO (r) changes from −68 MeV fm −2 (right scale of the insert) to −44 MeV fm −2 when f v = 0. That is, the tensor coupling accounts for roughly one third of the total spin-orbit strength in the G2 parameter set.
One expects that the weakening of the spin-orbit potential in going to the neutron drip line will bring about a reduction of the spin-orbit splittings ∆ε = ε nl,j=l−1/2 − ε nl,j=l+1/2 (16) of the neutron levels [30]. Figure 16 displays the energy splittings of some spin-orbit partner levels of neutrons for lead isotopes, obtained with the G2 and NL3 parameter sets. The splittings predicted by G2 and NL3 are very close as a consequence of the similarity of the corresponding spin-orbit potentials. Partner levels with high angular momentum undergo some reduction in the splitting along the Pb isotopic chain, but partners with small angular momentum show an almost constant splitting. By comparison of their RHB results for Ni and Sn, the authors of Ref. [30] pointed out that the weakening of the spin-orbit interaction should be less important for heavier isotopic chains. Our calculations for Pb would confirm this statement. All the single-particle levels involved in Figure 16 are bound. Of course, one should not expect the results for ∆ε to be so reliable in our qb-BCS approach if one, or both, of the partner levels lies at positive energy. The reason is that the single-particle energies of the quasibound levels do not exactly reproduce the energies of the corresponding canonical states of a RHB calculation.
In Figures 17a and 17b we show the spin-orbit potential for isotones of N = 82 from Z = 40 to Z = 70, for the G2 and NL3 parametrizations. Similarly to what is found for isotopes, the results obtained from G2 and from NL3 are comparable and the spin-orbit potential well V SO (r) moves outwards with the addition of protons, following the tendency of the proton density. However, for isotones we find that the behaviour of the depth of the spin-orbit potential well is not so monotonous: it increases when one goes from the neutron drip line up to the β-stable region, while it decreases afterwards as more protons are added.

Summary and conclusion
We have analyzed the pairing properties of some chains of isotopes and isotones with magic Z and N numbers in the relativistic mean field approach. The study has been performed for the G1 and G2 parametrizations that were obtained in Ref. [10] from the modern effective field theory approach to relativistic nuclear phenomenology. We have compared the results with those obtained with the NL3 parameter set which is considered to be very successful for dealing with nuclei beyond the stability line.
For accurate calculations of pairing far from the valley of β-stability in the relativistic models, the relativistic Hartree-Bogoliubov approach should be applied. However, we have presented a simpler modified BCS approach which allows one to obtain pairing properties near the drip lines fast and confidently. The method has been used previously in nonrelativistic calculations with Skyrme forces [38,39]. The key ingredient is to take into account the continuum contributions through quasibound levels due to their centrifugal barrier.
To further simplify the calculations we have assumed pairing matrix elements of the type for Sn isotopes is larger in NL3 than in G1 and G2, while for Pb isotopes none of the three sets exhibits a quenching of the shell effect at N = 184 in our qb-BCS calculation.                        (b) 28