Introduction

The shell structure of nuclei is a fundamental framework in nuclear physics. However, for neutron-rich nuclei far from the stability valley, the traditional magic numbers N = 8, 20, 28, and 40 vanish [1-10] and new magic numbers such as N = 14, 16, 32, and 34 emerged [11-14]. The evolution of shell structures in exotic nuclei has significantly deepened our understanding of nuclear quantum many-body systems and their underlying nuclear forces. Light neutron-rich nuclei are particularly intriguing, as the development of radioactive beam facilities has facilitated unprecedented exploration of isotopic chains extending to the dripline region [15].

The neutron-rich oxygen isotopes ²²O and ²⁴O are doubly magic nuclei with a high excitation energy of the first 2⁺ excited state, providing clear evidence for the emergence of the N=14 and 16 subshell closures [14, 16-18]. However, a comparison of the systematic behavior of the 2⁺ energy levels in the oxygen and carbon isotopic chains suggests that the N=14 shell gap disappears in the carbon isotopes [19, 20]. In ¹⁵C, the ground state and first excited state are 1/2⁺ and 5/2⁺, respectively, indicating an inversion of the neutron ν1s_1/2 and ν0d_5/2 orbits in ¹⁵C compared to ¹⁷O [21, 22]. A similar inversion is anticipated when neutrons are added to ²⁰C [19, 20]. As the transition region between oxygen and carbon isotopes, nitrogen isotopes also signal erosion of the N=14 shell gap [23-26]. Moreover, the evolution of the N=16 shell gap in neutron-rich boron, carbon, nitrogen, and oxygen isotopic chains remains an intriguing question. Further experimental and theoretical studies are required to fully understand the origin and evolution of N=14 and N=16 shell gaps.

Neutron-rich boron, carbon, and nitrogen isotopes exhibit many exotic phenomena. For instance, two-neutron halo structures were observed in ¹⁷B [27], ¹⁹B [28], and ²²C [29], while one-neutron halo structures were found in ¹⁵C [30] and ¹⁹C [31]. In addition, the observation of the most neutron-rich unbound nuclei, ²⁰B and ²¹B, was recently reported in Ref. [32]. Furthermore, a thick neutron skin has been detected in ^17-22N [26]. These light neutron-rich systems, which straddle the neutron dripline, have accumulated a wealth of experimental data, thereby providing an ideal testing ground for theoretical models [33]. These studies have inspired extensive theoretical studies, including the shell model (SM) [21, 34-36], Gamow shell model (GSM) [37-39], antisymmetrized molecular dynamics (AMD) [40, 41], and ab initio no-core shell model (NCSM) calculations [42-44]. However, traditional SM calculations, which employ phenomenological interactions, have limited predictive power as the nuclei approach the dripline. In contrast, ab initio calculations based on high-resolution interactions derived from the chiral effective field theory (χEFT) of quantum chromodynamics have made great progress over the past decades [45-52], offering a more robust framework for understanding and predicting the properties of neutron-rich nuclei.

Systematic calculations of neutron-rich boron, carbon, and nitrogen isotopes enable us to gain a more comprehensive understanding of their exotic structures and shell evolution while simultaneously testing our ab initio methods and nuclear forces. Therefore, in this study, we performed systematic calculations of neutron-rich boron, carbon, and nitrogen isotopes using the ab initio valence-space in-medium similarity renormalization group (VS-IMSRG) [53-58] based on nucleon-nucleon (NN) and three-nucleon (3N) interactions derived from χEFT.

The remainder of this paper is organized as follows: First, we provide a brief introduction to the VS-IMSRG framework. Next, we present systematic calculations of the spectra of neutron-rich boron, carbon, and nitrogen isotopes. Subsequently, we discuss the shell evolution of the N=14 and N=16 shell gaps in these isotopes utilizing the calculated effective single-particle energies (ESPE). Finally, we conclude with a summary of this work.

Method

The intrinsic Hamiltonian of an A-nucleon system can be written as $H={\displaystyle \sum _{i<j}^A}\left(\frac{{\left(p_i-p_j\right)}^2}{2mA}+v_{ij}^{\text{NN}}\right)+{\displaystyle \sum _{i<j<k}^A}{v}_{ijk}^{\text{3N}}\mathrm{,}$ (1) where p is the nucleon momentum in the laboratory, m is the nucleon mass, $v^{\text{NN}}$ and $v^{\text{3N}}$ are the NN and 3N interactions, respectively. In this study, we used three different χEFT NN + 3N interactions: EM1.8/2.0 [59, 60], N³LO+3N(lnl) [9, 61] and N⁴LO+3N(lnl) [62]. The EM1.8/2.0 interaction, consisting of a next-to-next-to-next-to-leading order (N³LO) NN interaction softened by the similarity renormalization group (SRG) evolution with momentum resolution scale λ=1.8 fm^-1 and a next-to-next-to-leading order (N²LO) 3N interaction with momentum cutoff Λ=2.0 fm^-1, can reproduce well the ground-state energies up to ¹³²Sn [59, 60]. The N³LO+3N(lnl) interaction, which presents good description of ground-state energies up to nickel isotopes is composed of an N³LO NN interaction and an N²LO 3N interaction, adopting a large SRG scale of λ=2.6 fm^-1 for the NN interaction without including the induced 3N force [9, 61]. The N⁴LO+3N(lnl) interaction, composed of a next-to-next-to-next-to-next-to-next-to-leading order (N⁴LO) NN interaction and an N²LO 3N interaction, also adopts a large SRG scale of λ=2.6 fm^-1 for the NN interaction, without including the induced 3N force [62]. In this study, we take a harmonic-oscillator basis at $\hslash \omega =16\text{MeV}$ with e_max=2n+l=14 and E_3max=14, which is sufficiently large for convergence.

In practical calculations, Hamiltonian (1) is rewritten as normal-order operators with respect to the Hartree-Fock reference state, $\begin{array}{l}H=E_0+{\displaystyle \sum _{ij}}{f}_{ij}:a_i^{†}{a}_j:+\frac{1}{4}{\displaystyle \sum _{ijkl}}{\text{Γ}}_{ijkl}:a_i^{†}{a}_j^{†}{a}_l{a}_k:\\ +\frac{1}{36}{\displaystyle \sum _{ijklmn}}{W}_{ijklmn}:a_i^{†}{a}_j^{†}{a}_k^{†}{a}_n{a}_m{a}_l:\mathrm{,}\end{array}$ (2) where E₀, f, Γ, W refer to the normal-ordered zero-, one-, two- and three-body terms respectively. Because the contribution of 3N force can be captured well at the normal-ordered two-body level [63], the residual normal-ordered three-body term W is neglected.

Next, using the continuous unitary transformation $U(s)$, the Hamiltonian (2) is decoupled from the large Hilbert space to a small valence space using VS-IMSRG. This is achieved by solving the flow equation, $\frac{\text{d}H(s)}{\text{d}s}=[\eta (s)\mathrm{,}H(s)]\mathrm{,}$ (3) with an anti-Hermitian generator, $\eta (s)\equiv \frac{dU(s)}{ds}{U}^{†}\left(s\right)=-{\eta }^{†}\left(s\right)\mathrm{.}$ (4) In the present study, we adopted protons in the p shell and neutrons in the sd shell as the valence space. In practical calculations, the Magnus formalism is employed with all operators truncated at the two-body level [64]. The effective Hamiltonians for this model space are consistently derived by VS-IMSRG with ensemble normal ordering (ENO) [55]. Then, we use the large-scale shell model code KSHELL [65] to diagonalize the effective Hamiltonians.

Results

3.1

Energy Spectra

Nuclear spectra provide invaluable insights into nuclear structure. In this study, the spectra of neutron-rich boron, carbon, and nitrogen isotopes were systematically calculated using ab initio VS-IMSRG with three sets of χEFT NN + 3N interactions, that is, EM1.8/2.0 [59, 60], N³LO+3N(lnl) [9, 61], and N⁴LO+3N(lnl) [62]. For comparison, we also performed phenomenological SM calculations using the YSOX interaction within the full psd valence space, which accurately reproduced the properties of boron, carbon, nitrogen, and oxygen isotopes [34]. The results were compared with the available experimental data, as illustrated in Fig. 1, 2 and 3.

Fig. 1Full-screen View

Energy spectra of ^16-24N calculated by ab initio VS-IMSRG using the EM1.8/2.0, N³LO+3N(lnl), and N⁴LO+3N(lnl) interactions, compared with the SM calculations using the YSOX interaction and the available experimental data [66]

Fig. 2Full-screen View

Energy spectra of ^15-22C calculated by ab initio VS-IMSRG using the EM1.8/2.0, N³LO+3N(lnl), and N⁴LO+3N(lnl) interactions, compared with the SM calculations using the YSOX interaction and the available experimental data [66]

Fig. 3Full-screen View

Energy spectra of ^14-21B calculated by ab initio VS-IMSRG using the EM1.8/2.0, N³LO+3N(lnl), and N⁴LO+3N(lnl) interactions, compared with the SM calculations using the YSOX interaction and the available experimental data [66]

Nitrogen isotopes have abundant experimental data, which provide a valuable test ground for our ab initio calculations. In addition, nitrogen nuclei lie between oxygen and carbon isotopes, offering essential insights into the shell evolution. For ¹⁶N, the ordering of the first four states is experimentally determined as 2^-, 0^-, 3^-, and 1^- [66]. Our VS-IMSRG calculations with the EM1.8/2.0 and N³LO+3N(lnl) interactions showed good agreement with the experimental results, whereas the SM calculations using the YSOX interaction failed to reproduce the spectra of ¹⁶N. For ¹⁷N, the VS-IMSRG calculations with the three NN + 3N interactions as well as the SM calculations using the YSOX interaction reproduced the lowest three states well. In the case of ¹⁸N, the ground state was correctly reproduced by the VS-IMSRG calculations with the EM1.8/2.0 interaction and the SM calculations using the YSOX interaction. However, the VS-IMSRG calculations with N³LO+3N(lnl) and N⁴LO+3N(lnl) interactions give 2^- as the ground state. For ^19,21N, all calculations correctly gave the ground state as 1/2^- and predicted the first excited state as 3/2^-, with the excited states yet to be experimentally confirmed. For ²⁰N, the ground state is 2^- given by all calculations. In the case of ²²N, the VS-IMSRG calculations with the N³LO+3N(lnl) and N⁴LO+3N(lnl) interactions predicted a 0^- ground state, whereas the VS-IMSRG calculations with the EM1.8/2.0 interaction and the SM calculations using the YSOX interaction predicted the ground state to be 2^- and 1^-, respectively. The size of the N=16 shell gap in ²³N remains unclear [67]. All calculations predicted the ground state of ²³N to be 1/2^-. However, the SM calculations with the YSOX interaction yielded higher excitation spectra than the VS-IMSRG results, indicating an overly large N=16 shell gap in the SM calculations using the YSOX interaction. This behavior was also observed in ²¹B, ²²C, and ²⁴O, as discussed above. For ²⁴N, all the calculations give the ground state as 2^- and predict the first excited state as 1^-.

Neutron-rich carbon isotopes have attracted considerable experimental and theoretical attention. The ground state and first excited state of ¹⁵C are 1/2⁺ and 5/2⁺, respectively, indicating an inversion of the neutron ν1s_1/2 and ν0d_5/2 orbits in carbon isotopes [21, 22]. Our VS-IMSRG calculations, along with the SM calculations, accurately reproduced the spectra of ¹⁵C. The VS-IMSRG calculations with the N³LO+3N(lnl) interaction for ^16,17C did not converge; therefore, we only present the results for ^16,17C using the EM1.8/2.0 and N⁴LO+3N(lnl) interactions in Fig. 2. For ^16,18,20C, the first 2⁺ excited states obtained from our VS-IMSRG calculations are in good agreement with the experimental data, whereas SM calculations using the YSOX interaction generally provide higher energies for the first 2⁺ excited states of these nuclei. In particular, for ²²C, the energy of the first 2⁺ excited state calculated by SM using the YSOX interaction is significantly higher than that predicted by VS-IMSRG. Notably, the first 2⁺ excited state of ²⁴O, calculated by SM with YSOX interaction, is 5.41 MeV, which is higher than the experimental value of 4.76 MeV. In contrast, our VS-IMSRG calculations with the EM1.8/2.0, N³LO+3N(lnl), and N⁴LO+3N(lnl) interactions give the first 2⁺ excited state of ²⁴O to be 4.87, 5.26, and 4.39 MeV, respectively, which are in better agreement with the experimental value. Moreover, the first 2⁺ excited state of ²²C calculated by VS-IMSRG is significantly higher than that of the neighboring carbon isotopes, implying the presence of the N=16 shell gap in ²²C. However, the calculated first 2⁺ excited state of ²²C is lower than that of ²⁴O, suggesting a reduction in the N=16 shell gap in ²²C compared to ²⁴O. In addition, the low energy of the first 2⁺ excited state in ²⁰C indicates the disappearance of the N=14 shell gap in ²⁰C [19]. For ¹⁷C and ¹⁹C, our VS-IMSRG calculations failed to reproduce the experimental spectra because of the absence of continuum coupling in the current method, as discussed in Ref. [37, 68].

In the case of neutron-rich boron isotopes, experimental data are limited. For ¹⁴B, the ordering of the first five states is experimentally assigned as 2^-, 1^-, 3^-, 2^-, and 4^- [66]. Our ab initio VS-IMSRG calculations with the three NN + 3N interactions as well as the SM calculations using the YSOX interaction are all in good agreement with the experimental results. The spin and parity of the ground state for isotopes heavier than ¹⁴B have not yet been confirmed experimentally. For 15,17,19,21B, our VS-IMSRG calculations with the three NN + 3N interactions predict that the ground states of these isotopes are all 3/2^-, which is consistent with SM calculations using the YSOX interaction. For ¹⁶B, the VS-IMSRG calculations with the N³LO+3N(lnl) and N⁴LO+3N(lnl) interactions predict the ground state to be 2^-, whereas the VS-IMSRG calculations with the EM1.8/2.0 interaction predict a ground state of 0^- with a 2^- excited state that is nearly degenerate. In addition, the ground state of ¹⁶B obtained from the SM calculations using the YSOX interaction is also 0^-. For the unbound nucleus ¹⁸B, SM calculations using YSOX and WBP [69] interactions predicted the ground state to be 2^-. Our VS-IMSRG calculations with the EM1.8/2.0 and N⁴LO+3N(lnl) interactions also predicted the ground state of ¹⁸B to be 2^-, but the VS-IMSRG calculations with the N³LO+3N(lnl) interaction gave the ground state of ¹⁸B to be 3^-. Both ²⁰B and ²¹B are unbound nuclei that exist as resonances and decay via one or two neutron emissions [32]. For ²⁰B, our VS-IMSRG calculations with the three NN + 3N interactions predict the ground state to be 1^-, whereas SM calculations using the YSOX interaction yield the ground state 2^-. For ²¹B, our VS-IMSRG calculations with the three NN + 3N interactions predicted a significantly lower energy for the first excited state 1/2^- compared to the SM calculations using the YSOX interaction, a situation similar to that observed for ²²C and ²⁴O above. The 1/2^- state of ²¹B calculated with VS-IMSRG has larger neutron excitation components than those in the SM calculations. However, the subsequent analysis of the ESPE calculated by VS-IMSRG suggests that the N=16 subshell still exists in ²¹B.

Overall, our ab initio VS-IMSRG calculations based on the three χEFT NN + 3N interactions, particularly the EM1.8/2.0 interaction, were in good agreement with the available experimental data. These results provide critical assistance for future research.

3.2

Shell evolution

In neutron-rich boron, carbon, nitrogen, and oxygen isotopes, the evolution of N = 14 and N = 16 shell gaps has attracted significant experimental and theoretical attention. The ESPE provides a good reflection of the size of shell gaps, which defined as [70, 71] ${\epsilon }_i={\epsilon }_i^{\text{core}}+{\displaystyle \sum _j}{n}_j{V}_{ij}^m\mathrm{,}$ (5) where ${\epsilon }_i^{\text{core}}$ is the valence space single-particle energy, Vm is the monople interaction, and nj is the occupation number calculated consistently using the SM.

To study the evolution of the N=14 and N=16 shell gaps, we calculated the ESPE of the ν1s_1/2, ν0d_5/2, and ν0d_3/2 orbits in these nuclei employing the VS-IMSRG method, based on the EM1.8/2.0, N³LO+3N(lnl), and N⁴LO+3N(lnl) interactions, compared with the SM calculations using the YSOX interaction. As shown in Fig. 4, for the VS-IMSRG calculations, when protons are removed from the oxygen isotopes, the ESPE of the ν1s_1/2 orbital is gradually lowered relative to the ν0d_5/2 orbital, and a reversal occurs in the boron and carbon isotopes. This indicates that the N=14 neutron subshell is present in the oxygen isotopes but disappears in the boron, carbon, and nitrogen isotopic chains. The SM calculations using the YSOX interaction show the same trend as the VS-IMSRG results, whereas the N=14 shell gaps in ²¹N and ²²O calculated by SM are significantly larger than those obtained from VS-IMSRG. Notably, the N=16 sub-shell is present in the boron, carbon, nitrogen, and oxygen isotopic chains and gradually decreases from ²⁴O to ²¹B. Moreover, the N=16 shell gap calculated by SM using the YSOX interaction was significantly larger than that calculated by VS-IMSRG. Additionally, the N=16 shell gap calculated by VS-IMSRG using the N³LO+3N(lnl) interaction is generally larger than that calculated using the EM1.8/2.0 and N⁴LO+3N(lnl) interactions. Correspondingly, the energy spectra exhibited higher excitation energies, particularly in ²⁴O, as discussed above.

Fig. 4Full-screen View

ESPEs relative to the neutron ν0d_5/2 orbital in boron, carbon, nitrogen, and oxygen isotopes, calculated by VS-IMSRG using the EM1.8/2.0, N³LO+3N(lnl) and N⁴LO+3N(lnl) interactions, respectively, compared with the SM calculations using the YSOX interaction

To further understand the contribution of the 3N force to the shell evolution, Fig. 5 shows the VS-IMSRG calculations using the full NN + 3N EM1.8/2.0 interaction (solid lines) compared with the VS-IMSRG calculations without the 3N force (dashed lines). As shown in Fig. 5, the 3N force enhances the N=14 shell gap in the oxygen isotopes while reducing the N=16 shell gap in the boron, carbon, nitrogen, and oxygen isotopes. A similar trend was observed for the results using N³LO+3N(lnl) and N⁴LO+3N(lnl) interactions.

Fig. 5Full-screen View

ESPEs relative to the neutron ν0d_5/2 orbital in boron, carbon, nitrogen, and oxygen isotopes, calculated by VS-IMSRG with the full NN + 3N EM1.8/2.0 interation (solid lines), compared to the results obtained without the 3N force included (dashed lines)

To analyze the roles of different components of the interaction in shell evolution, we employed the spin-tensor decomposition method [72-74] to decompose the effective interaction derived from VS-IMSRG into central, tensor, and spin-orbit (LS) parts. First, we conducted systematic calculations for boron isotopes and compared the contributions of different components of the interaction derived from VS-IMSRG across three sets of χEFT NN + 3N interactions: N³LO+3N(lnl), N⁴LO+3N(lnl), and EM1.8/2.0. We also compared the results with the SM calculations using the YSOX interaction. As shown in Fig. 6, the solid lines represent the results with the full interaction, the dashed lines represent the results without the contribution of the tensor force, and the dot-dashed lines represent the results without the contribution of the two-body LS part. In our VS-IMSRG calculations, we can see that the tensor force slightly lowers the ESPE of the ν0d_3/2 orbital relative to the ν0d_5/2 orbital, whereas the two-body LS component increases the ESPE of both the ν0d_3/2 and ν1s_1/2 orbits relative to the ν0d_5/2 orbital. Moreover, the contributions of the tensor force and the two-body LS part were similar across the three sets of NN + 3N forces. Additionally, SM calculations using the YSOX interaction demonstrated the same trend as the VS-IMSRG results, although the strength of the tensor force and the two-body LS force showed minor differences. Next, we present the contributions of the different components of the interactions in various isotopes. For clarity, we only show the results calculated by VS-IMSRG using the EM1.8/2.0 interaction. In Fig. 7, the solid and dashed lines represent the ESPE with and without the contribution of the tensor force, respectively. It can be seen that the tensor force only slightly reduces the N=16 shell gap for all the isotopes. Fig. 8 compares the results obtained with and without the two-body LS part. The two-body LS part significantly lowers the ESPE of the ν0d_5/2 orbital and causes inversion of the ν1s_1/2 and ν0d_5/2 orbits in oxygen isotopes compared to boron and carbon isotopes. This inversion is the origin of the N=14 shell gap in the oxygen isotopes. In addition, the two-body LS part significantly increases the ESPE of the ν0d_3/2 orbital, leading to the appearance of the N=16 sub-shell.

Fig. 6Full-screen View

ESPEs relative to the neutron ν0d_5/2 orbital in boron isotopes, calculated by VS-IMSRG using the N³LO+3N(lnl), N⁴LO+3N(lnl) and EM1.8/2.0 interactions with (solid lines) and without the contribution of the tensor force (dashed lines)/LS force (dot-dashed lines), as well as the SM calculations using the YSOX interaction

Fig. 7Full-screen View

ESPEs of the neutron ν1s_1/2, ν0d_5/2, and ν0d_3/2 orbits in boron, carbon, nitrogen, and oxygen isotopes, calculated by VS-IMSRG using the EM1.8/2.0 interaction with (solid lines) and without (dashed lines) the contribution of the tensor force

Fig. 8Full-screen View

ESPEs of the neutron ν1s_1/2, ν0d_5/2, and ν0d_3/2 orbits in boron, carbon, nitrogen, and oxygen isotopes, calculated by VS-IMSRG using the EM1.8/2.0 interaction with (solid lines) and without (dashed lines) the contribution of the spin-orbit (LS) force

In Ref. [4], the inversion between the neutron ν1s_1/2 and ν0d_5/2 orbits in ¹⁵C and ¹⁷O is depicted by analyzing the monopole matrix elements of the YSOX interaction, assuming that, from ¹⁵C to ¹⁷O, the proton $\pi 0p_{1/2}$ orbit is fully occupied, and the last neutron is in either the ν1s_1/2 or ν0d_5/2 orbital. Similar to the discussion in Refs. [4], we extracted the contributions of the monopole matrix elements from various components of the effective interaction derived by VS-IMSRG, including the central, tensor, and LS parts, as shown in Fig. 9. Here, the rest is due to the Adependence of the Hamiltonian. The two-body LS part significantly contributes to the shift of the ν1s_1/2 orbital relative to the ν0d_5/2 orbital, even when determining the inversion of the ν1s_1/2 and ν0d_5/2 orbits. In contrast, the contribution from the tensor force was much smaller than that observed in the YSOX interaction in Ref. [4].

Fig. 9Full-screen View

ESPEs of the neutron ν1s_1/2 orbit relative to the ν0d_5/2 orbit in ¹⁵C and ¹⁷O, calculated by VS-IMSRG with the EM1.8/2.0 interaction. The contributions of the tensor, LS, and central forces to the ESPE were decomposed. The rest part is due to our Hamiltonian is A-dependent

Conclusion

Utilizing three sets of chiral NN + 3N interactions, EM1.8/2.0, N³LO+3N(lnl), and N⁴LO+3N(lnl), we present ab initio VS-IMSRG calculations for neutron-rich boron, carbon, and nitrogen isotopes. We systematically calculated and predicted the low-lying spectra of these nuclei. For comparison, we also performed SM calculations using the YSOX interaction. Based on the calculated spectra and effective single-particle energies, we studied the evolution of N=14 and N=16 shell gaps. Our results suggest that the N=14 neutron subshell is present in the oxygen isotopes but collapses in the boron, carbon, and nitrogen isotopic chains. Furthermore, the N=16 sub-shell is predicted to be present in boron, carbon, nitrogen, and oxygen isotopes, but its size gradually decreases from ²⁴O to ²¹B. Additionally, we investigated the roles of the different components of the interaction in shell evolution by employing the spin-tensor decomposition method. We find that the two-body spin-orbit force plays a significant role in the formation of the N=14 and N=16 shell gaps. In general, our ab initio VS-IMSRG calculations agree well with the available experimental data, and theoretical predictions for these nuclei will be helpful for future experiments.