Introduction

Shell structure is a distinctive feature of nuclear many-body systems and is characterized by proton and neutron magic numbers. These magic numbers are of fundamental importance for understanding nuclear stability, intrinsic nuclear structure, nucleosynthesis, and practical applications [1-5]. Magic nuclei can be identified through various observables, such as sudden increases in binding energy [6], high excitation energies, reduced E₂ transition probabilities [7], localized abrupt changes in charge radii [8, 9] and proton radii [10, 11], and reduced neutron/proton capture cross sections [12] relative to neighboring nuclei. Modern radioactive-beam facilities and detection systems have further extended the exploration of extreme nuclides, potentially revealing new magic numbers or revising the traditional shell structure [1, 13-15].

The systematic evolution of bulk properties along the oxygen-to-argon isotopic chains carries rich information on shell structure, especially the disappearance of traditional magic numbers and the emergence of new ones. In oxygen isotopes, the high excitation energy and low value of the state [16, 17], the small quadrupole transition parameter β₂ [18], the reduced proton radius [11], and the large inclusive cross section together with the broad momentum distribution observed in quasifree (p,pN) scattering [12] provide strong evidence for the N=14 subshell closure in ²²O. In addition, a large shell gap at N=16 has been clearly identified in ²⁴O from various measurements [11, 17, 19-24], indicating that ²⁴O is a doubly magic nucleus. Measurements of masses, charge radii, and Coulomb excitation in neutron-rich Ne, Na, and Mg isotopes [9, 25-27] have suggested the breakdown of the traditional magic number N=20 [28-31]. At higher mass, the observation of a low-lying state in ⁴²Si provides clear evidence for the collapse of the N=28 shell closure [32].

A variety of nuclear-structure models have been employed to elucidate the mechanisms responsible for the emergence of new magic numbers and the weakening of conventional shell closures. Shell-model studies suggest that the appearance of the new magic numbers N=14 and N=16 is driven by the strong neutron-proton tensor interaction [10, 11, 21, 33]. Ab initio calculations with modern two- and three-nucleon forces considerably improve the simultaneous description of binding energies and charge and matter radii for stable O isotopes, although deficiencies remain for the most neutron-rich systems [34]. The collapse of the N=20 shell closure is commonly attributed to occupation of the neutron pf shell together with the sd orbitals at sizable prolate deformation, i.e., the island of inversion [35, 36]. The coupled-cluster method based on nucleon-nucleon and three-nucleon potentials qualitatively reproduces the evolution of charge radii and the N=14 magicity in Ne and Mg isotopes when angular-momentum projection is included, although the isotope shifts remain challenging [37]. In mean-field theories, the ground state of ³²Mg is spherical and becomes deformed only after a 20% increase in the spin-orbit strength for the Skyrme force SLy4 [38]. A deformed ground state can also be obtained by adjusting the neutron and proton pairing gaps, but the magic number N=20 still persists in the relativistic mean-field (RMF) description of ³²Mg [39]. Angular-momentum-projection approaches based on the HFB model [40, 41] and the RMF model [42] change the spherical mean-field ground state of ³²Mg into a deformed one with β₂ close to the measured value. A similar picture is also obtained in the projected shell model [43].

As noted above, a local reduction in charge radii along an isotopic chain is one of the signatures of shell closure. However, charge radii are influenced by various mechanisms, including pairing correlations [44-48], deformation [49, 50], cluster structure [51, 52], shell evolution [53-55], and center-of_mass correlations [56-58]. Among these, pairing correlations, particularly neutron-proton pairing near the Fermi surface, are considered a key factor. To this end, a modified RMF plus BCS [RMF(BCS)*] model has been proposed, which incorporates a term for neutron-proton pairing correlations around the Fermi surface into the charge-radius formula [59, 60]. This model successfully describes the odd-even staggering and the inverted parabolic behavior observed in the charge radii of Ca isotopes and provides a reasonable description for most O, Ne, and Mg isotopes. However, it predicts odd-even staggering in the O isotopic chain together with a sudden increase at N=14, and it underestimates the charge radii of Mg nuclei with N=14 and N>18. Because the BCS approximation is not well suited to nuclei far from the β-stability line [61, 62], whereas the magic numbers 14, 16, and 20 in the sd shell appear in neutron-rich systems and even at the drip line, it is necessary to reexamine the charge radii within the same ansatz using a more appropriate treatment of pairing correlations.

The Bogoliubov transformation provides a more reliable treatment of pairing correlations than the BCS method for unstable nuclei far from the β-stability line. On this basis, we recently incorporated neutron-proton pairing correlations extracted from quasiparticle states around the Fermi surface into the multidimensionally constrained relativistic Hartree-Bogoliubov (MDC-RHB) model [63]. The charge radii of Ca and Ni isotopes, as well as the N=28, 30, 32, and 34 isotones, are then well reproduced. In the present work, we revisit the difficulties encountered in the RMF(BCS)* approach and examine the shell structure around N=14, 16, and 20 from the perspective of nuclear charge radii, thereby clarifying the relation between shell structure and neutron-proton pairing in these isotopes.

This paper is organized as follows. Section 2 briefly introduces the MDC-RHB model and the neutron-proton pairing correlation extracted from quasiparticle states around the Fermi surface. In Sect. 3, the charge radii of the Z=8, 10, 12, 14, and 18 isotopic chains are investigated within this framework, and the shell-structure phenomena around N=14, 16, and 20 are discussed. The main conclusions are summarized in Sect. 4.

Theoretical framework

In the MDC-RHB model, the RHB equation in coordinate space can be written as follows [64, 65]:(1)where λ is the Fermi energy, Δ is the pairing field, Ek is the quasiparticle energy, and (U_k(r), V_k(r))^T is the quasiparticle wave function. h is the single-particle Hamiltonian which can be written as follows:(2)where M is the nucleon mass, and S(r), V(r), and Σ_R(r) are the scalar, vector, and rearrangement potentials, respectively. For meson-exchange interactions,(3)where , , and are the coupling constants of the σ, ω₀, and ρ₀ meson fields, respectively; A₀ is the time-like component of the Coulomb field generated by photons; and e is the proton charge. The values τ₃=1 and -1 distinguish the neutron and proton components, respectively. For point-coupling interactions,(4)where α_S, α_V, α_TS, α_TV, β_S, γ_S, γ_V, δ_S, δ_V, δ_TS, and δ_TV are coupling constants for different channels. The subscripts S, V, and T denote scalar, vector, and isovector couplings, respectively.(5)which represent the isoscalar-scalar, isoscalar-vector, isovector-scalar, and isovector-vector densities, respectively. The pairing field Δ is calculated by the effective pairing interaction V and the pairing tensor κ as follows:(6)where(7)In this work, a separable pairing force of finite range is adopted. The effective pairing interaction is(8)where G is the pairing strength. and are the center-of-mass and relative coordinates, respectively. The function represents a Gaussian distribution, . The pairing strength G = 728 MeV fm³ and effective range of the pairing force a = 0.644 fm are used [66].

The modified root-mean-square (rms) charge radius r_ch is given by [67]:(9)where the first term represents the charge distributions of point-like protons, and the second term is incorporated due to the finite size of protons. The third term represents short-range correlations from the difference of the fractions of Cooper pairs between neutrons and protons [59]. A is the mass number, and is a normalization constant. The quantity is defined as(10)where is the occupation probability of the kth quasiparticle orbital for neutrons or protons, and . D_n (D_p) is used to characterize the surface diffuseness encoded in the neutron (proton) eigenfunctions [68]. In practice, the quasiparticle levels satisfying MeV are included in Eq. (10). The last term in Eq. (9) represents the neutron-proton correlation arising from a simultaneously unpaired neutron and proton. The values of a₀=0.561 and Δ=0.355(0.000) for odd-odd (even-even, odd-even, and even-odd) nuclei are the same as those in Ref. [67].

In the MDC-RHB model, multipole moments can be considered under the V₄ symmetry [69, 70]. In this work, we restrict the calculations to axial and reflection symmetries; namely, only the quadrupole deformation β₂₀ is taken into account. The effective interactions PK1 [56], NL3 [71], DD-ME2 [72], and DD-PC1 [73] are used in our calculations to avoid parameter-dependent results.

Results and discussion

As discussed in Sect. 1, shell effects at N=14 and N=16 are observed in O, N, Ne, and Mg isotopes. To investigate these characteristic magic numbers, we calculate the evolution of charge radii along the O, Ne, Mg, Si, and Ar isotopic chains within the MDC-RHB model. The results obtained without and with the neutron-proton correlations around the Fermi surface [the third and fourth terms in Eq. (9)] are denoted by RHB and RHB*, respectively. Figure 1 shows the results calculated with the PK1, DD-ME2, NL3, and DD-PC1 effective interactions. The experimental values from Refs. [8, 27, 74, 75] are also shown for comparison. The RHB calculations yield the charge radii for ^16-22O within the same order of magnitude, as shown in Fig. 1(a). Among the four interactions, PK1 provides the best description of ¹⁶O, whereas the charge radius of ¹⁸O is underestimated by all four forces. As shown in Fig. 1(b), the correction term increases the charge radius of ¹⁸O and thus improves the description obtained with the PK1 and NL3 parameter sets. The double-magic character of ¹⁶O makes it largely free of surface correlations, so its charge radius remains close to the original value. In the RHB* model, the systematic trend of the charge radii shows clear depressions at N=8 and N=16, indicating that the corresponding shell effects are enhanced after the corrections around the Fermi surface are included.

Fig. 1Full-screen View

(Color online) Charge radii of O, Ne, Mg, Si, and Ar isotopes calculated with the RHB (left) and RHB* (right) models using the DD-ME2, NL3, DD-PC1, and PK1 effective interactions. The corresponding experimental data are taken from Refs. [8, 27, 74, 75] (solid squares). The gray bands indicate the neutron numbers N=14 and 20

We attribute the kink at N=16 to the third term in Eq. (9). To clarify this point, Fig. 2 shows the single-particle levels of the ground states of ^22-26O as a function of the occupation probability v². For O isotopes, the proton magic number Z=8 implies a fully occupied proton shell and therefore D_p=0. For neutrons in ²²O to ²⁶O, pairing correlations broaden the single-particle levels around the Fermi surface and enable occupation of states above it. The values are plotted as olive solid lines together with the occupation numbers for each single-particle level; larger values correspond to partially occupied levels near the Fermi surface. The calculated D_n and D_p values derived from these quantities (see Table 1) show that D_n is largest in ²⁶O and relatively small in ²⁴O, consistent with the kink at N=16. Thus, the N=16 magic character is reinforced by enhanced neutron pairing correlations at N=18. Similarly, the sudden increase of D_n from ¹⁶O to ¹⁸O leads to the kink at N=8 along this isotopic chain. New charge-radius data for neutron-rich and proton-rich O isotopes from the planned laser-spectroscopy experiments of the COLLAPS and CRIS collaborations [13] will provide an important test of these predictions.

Fig. 2Full-screen View

(Color online) Single-particle levels of the ground states of ^22-26O as a function of the occupation probability υ², calculated with the RHB model using the NL3 effective interaction. The black dashed lines are the Fermi surfaces. The blue dashed-dotted line corresponds to the BCS formula with an average pairing gap. The olive solid line below each single-particle level denotes the corresponding value

The charge radii of Ne isotopes are shown in Fig. 1(c) and (d). The RHB results agree reasonably well with the data from Ref. [8]. In particular, the density-dependent effective interactions DD-ME2 and DD-PC1 reproduce the experimental values better, whereas the meson-exchange interactions NL3 and PK1 slightly underestimate them. Clear kinks are observed at N=8 and N=14. After the pairing correction is included, the charge radii increase, and the results obtained with the meson-exchange interactions become more consistent with the experimental data [8]. For comparison, the revised charge radii from Ref. [74] are also shown as gray squares. The RHB* predictions for ^26,28Ne obtained with DD-ME2 and DD-PC1 agree with these revised data. In both the RHB and RHB* calculations, the minimum charge radius occurs at N=14 for the DD-ME2, NL3, and PK1 interactions, whereas DD-PC1 gives the minimum at N=16. With increasing neutron number, the RHB model predicts no abrupt change in charge radius across N=20. By contrast, the RHB* model predicts that neutron-proton pairing leads to a sudden increase in the charge radius of ³²Ne, thereby signaling the magicity at N=20. The D_n and D_p values for Ne isotopes are listed in Table 1. These results suggest that the N=20 shell effect in ³²Ne is associated with enhanced pairing correlations around the neutron Fermi surface.

For the Mg isotopic chain, as shown in Fig. 1(e), the RHB calculations reproduce the charge radii of the proton-rich isotopes reasonably well, whereas those of the neutron-rich isotopes are underestimated. Kinks are observed at N=14 and N=20. As shown in Fig. 1(f), the surface correlation between protons and neutrons increases the charge radii, especially for isotopes with neutron number greater than 14. Consequently, the charge radii obtained with the RHB* model are in better agreement with the data than those from the RHB model. On the neutron-rich side of the Mg isotopic chain, this neutron-proton pairing vibration near the surface may therefore provide the microscopic origin of the observed charge-radius evolution.

In the RMF(BCS)* calculations [59], the same correction is applied to the charge radius, but the BCS approach is used in the pairing channel. As a result, the depression of the charge radius at N=14 is overestimated. By contrast, the charge radii of ²⁶Mg (N=14) obtained with the RHB* model in the present work are 3.0434 fm and 3.0415 fm for the NL3 and PK1 effective interactions, respectively, both of which are very close to the experimental value of 3.0340(26) fm [27]. This result shows that the description of the N=14 magicity is highly sensitive to the treatment of pairing correlations. Because both the pairing treatment and the form of the pairing force differ between Ref. [59] and the present work, it is necessary to examine the pairing force within the same mean-field Hamiltonian and to analyze the microscopic origin of the N=14 shell closure. We therefore calculate the charge radii of Mg isotopes by replacing the Bogoliubov quasiparticle transformation with the BCS approximation in the pairing channel. The results are shown in Fig. 3. As shown in Fig. 3(a), the charge radius obtained with the BCS treatment exhibits a pronounced kink at ²⁶Mg (N=14), both before and after the neutron-proton surface correction. This behavior is consistent with the results reported in Ref. [59] and suggests an excessively strong N=14 shell effect. In Fig. 3(b), however, the Bogoliubov quasiparticle transformation yields charge radii for ^24-29Mg that agree well with the data once the additional neutron-proton correction is included. In other words, this correction helps reproduce the charge radius at N=14.

Fig. 3Full-screen View

(Color online) Charge radii of Mg isotopes calculated with (a) the BCS approach and (b) the Bogoliubov quasiparticle transformation for the pairing treatment. The results are obtained with the NL3 effective interaction based on the same single-particle Hamiltonian. The corresponding experimental data are taken from Ref. [27] (solid squares). The gray band indicates the neutron numbers N=14 and 20

In Fig. 4, we analyze the microscopic origin of the reduction in charge radius at N=14 from the single-particle levels of ²⁶Mg. Figures 4(a) and (b) show the neutron and proton single-particle levels calculated with the BCS method, respectively. It is evident that both neutrons and protons mainly occupy the single-particle states below the Fermi surface, whereas the occupations above the Fermi surface are nearly zero. This indicates relatively small pairing gaps for both N=14 and Z=12, corresponding to enhanced nuclear stability and hence a smaller charge radius. According to Eq. (9), the neutron-proton pairing correlation around the Fermi surface for a given nucleus is determined by ΔD, which in turn depends on the noninteger occupation numbers of neutrons and protons. The BCS results show that the particles are almost integer occupied, leading to a small ΔD. The neutron single-particle levels obtained with the RHB model [Fig. 4(c)] are consistent with those from the BCS calculation. However, the RHB calculation yields fractional occupations for several proton single-particle levels around the Fermi surface [Fig. 4(d)]. Consequently, the charge radius of ²⁶Mg calculated with the RHB* model is substantially larger than that obtained with the RHB model, and the corrected result is closer to the experimental data.

Fig. 4Full-screen View

(Color online) Single-particle levels of the ground state of ²⁶Mg as a function of the occupation probability υ², calculated with the BCS approach and the Bogoliubov quasiparticle transformation for the pairing treatment. The NL3 effective interaction is used. The black dashed lines are the Fermi surfaces. The blue dashed-dotted line corresponds to the BCS formula with an average pairing gap

The absence of the N=20 shell closure in Mg isotopes is usually attributed to the island of inversion. The inversion of energy levels induces substantial deformation in the ground state of ³²Mg, thereby weakening the N=20 shell closure. However, the RHB model does not reproduce this phenomenon, because all four parameter sets yield spherical ground states. Beyond-mean-field calculations, which generate reasonably large deformation in the ground state of ³²Mg [40-42], can account for the discrepancy between the calculated and measured charge radii. In the present study, we show that the RHB* calculation also yields a more accurate charge radius for ³²Mg than the RHB calculation, especially with the density-dependent effective interactions DD-PC1 and DD-ME2. Beyond-mean-field calculations based on the present framework can be found in Refs. [76, 77], and further work in this direction is required. We also present the charge radii of ^19,20Mg in Fig. 3. The sudden change observed there indicates the persistence of the N=8 shell effect near the proton drip line, consistent with the recent invariant-mass reconstruction [78].

The charge radii of the Si isotopes are displayed in Fig. 1(g) and (h). The RHB model predicts a larger charge radius for ²⁸Si (N=14) than for ³⁰Si (N=16), contrary to the experimental trend. After the neutron-proton correlation around the Fermi surface is taken into account, the charge radii of these two nuclei agree with the data. Among these isotopes, the results obtained with the PK1 and NL3 effective interactions are comparable to the experimental values, including the new charge radius of ³²Si determined by collinear laser spectroscopy [75]. For N=20 and N=28, the kinks in the charge radii become more pronounced after the correction term is included, indicating an enhancement of the corresponding shell closures.

Figures 1(i) and (j) show the results for Ar isotopes. The charge radii calculated with the RHB model display a parabolic dependence on neutron number and deviate more strongly from the data. The RHB* results agree reasonably well with the experimental values on the neutron-rich side, except at N=22. On the proton-rich side, however, the charge radius of the N=14 isotope deviates substantially from the experimental value after the correction term is included.

To facilitate discussion of the N=14 magicity, we consider the charge-radius difference between adjacent even-even nuclei, defined as(11)for the N=14 isotones shown in Fig. 5. The results obtained with the NL3 and PK1 forces are very similar. As indicated by the green line, the empirical formula fm depends only on the mass number and therefore cannot reflect shell structure. The RHB calculations give a very small at Z=8, implying a strong N=14 shell effect in ²²O. This is consistent with the measured proton-radius data [11]. After the neutron-proton pairing correction is included, the N=14 shell closure in ²²O is weakened, whereas that in ²⁸Si is enhanced. Future measurements of weakly bound nuclei such as ²⁰C, ²²O, and ³²Ar, although experimentally challenging, will be important for clarifying the nature of the N=14 shell closure.

Fig. 5Full-screen View

(Color online) Charge-radius differences between adjacent even-even nuclei () as a function of proton number for the N=14 isotones, calculated with the (a) NL3 and (b) PK1 effective interactions. The corresponding experimental data are taken from Refs. [8, 27] (solid squares). The empirical values with (r₀=1.2 fm) are shown for reference

Finally, the average deviation and the root-mean-square deviation Δ_rms are defined as follows:(12)The corresponding values for the 25 experimental and calculated charge radii are listed in Table 2. The experimental charge radii of Mg and Ne isotopes, as well as that of ³²Si, are taken from Refs. [27], [74], and [75], respectively; all other data are taken from Ref. [8]. In Eq. (12), denotes the experimental uncertainty of the measured charge radius for the ith nucleus. For the DD-ME2, NL3, and PK1 parameter sets, the values obtained with the RHB* model are smaller than those from the RHB model, whereas the opposite trend is found for DD-PC1. The Δ_rms values for the meson-exchange interactions NL3 and PK1 decrease after the correction is included, while those for the density-dependent interactions DD-ME2 and DD-PC1 increase slightly. These results indicate that neutron-proton pairing corrections around the Fermi surface are important for mean-field calculations based on meson-exchange interactions. By contrast, density-dependent effective interactions do not appear to require this correction, because the RHB calculations already yield larger charge radii than the meson-exchange interactions.

Summary

In this work, we systematically investigate the evolution of charge radii from oxygen to argon isotopic chains within the RHB framework. A refined charge-radius formula that incorporates neutron-proton pairing correlations near the Fermi surface is adopted to describe the anomalous behavior. We analyze the emergence of new magic numbers at N=14 and N=16, together with the weakening of the traditional N=20 shell closure.

The inclusion of neutron-proton pairing correlations produces pronounced kinks in the charge radii at N=8, 20, and 28, demonstrating that these correlations enhance the corresponding shell effects. The Mg isotopic chain provides particularly useful insight into the interplay between the pairing treatment and shell structure. Whereas the conventional BCS approach artificially strengthens the N=14 shell effect, the Bogoliubov transformation predicts enhanced proton pairing correlations and yields charge radii in closer agreement with the data. Comparison with the experimental charge radii of 25 even-even nuclei further shows that the neutron-proton pairing correction significantly improves the description obtained with meson-exchange effective interactions, whereas it does not improve the results for density-dependent interactions. The interaction-dependent performance of this correction indicates that a more unified theoretical framework is still needed to reconcile meson-exchange and density-dependent interactions.