Introduction

The study of exotic nuclei far from the β stability line is a challenging frontier in experimental [1-4] and theoretical [5-10] nuclear physics. The unstable nuclei with extreme N/Z ratios, which are weakly bound systems, have exhibited many exotic properties that differ from those of the stable nuclei, such as halo structures [11-17], changes in the traditional magic numbers [18-23], and new nuclear excitation modes [24, 25], which may herald new physics. The study of exotic nuclei is crucial for not only comprehensively understanding the rich nuclear structures and properties but also investigating element synthesis and nuclear astrophysics [26]. However, because of the short lifespan, the cross sections for synthesizing exotic nuclei are small, which makes it difficult to create them experimentally. Thus, advanced large-scale radioactive beam facilities and updated detector techniques have been developed, upgraded, or planned worldwide [27-34]. Meanwhile, abundant theoretical studies on exotic nuclei provide valuable guidance for the design of experiments and analysis of experimental results [35-37].

In the exotic nuclei—particularly the drip line nuclei—the neutron or proton Fermi surfaces are typically close to the continuum threshold. With the effects of the pairing correlation, the valence nucleons have a certain probability to be scattered into the continuum and occupy the resonant states therein, making the nuclear density distributions very diffuse and extended. It is therefore essential to treat the pairing correlations and the couplings to the continuum properly in the theoretical descriptions of exotic nuclei [38-42]. Additionally, for one-neutron halo nuclei, e.g., ³¹Ne [43] and ³⁷Mg [44], the blocking effect [45] should be considered to treat the unpaired odd nucleon.

The Hartree–Fock–Bogliubov (HFB) theory has achieved great success in describing exotic nuclei with a unified description of the mean field and pairing correlation via Bogoliubov transformation [45]. Different models based on the HFB theory have been used to study exotic nuclei, such as the Gogny–HFB theory [46], the Skyrme–HFB theory [47], the relativistic continuum Hartree–Bogoliubov (RCHB) theory [39], and the density-dependent relativistic Hartree–Fock–Bogoliubov (RHFB) theory [48]. To explore the halo phenomena in deformed nuclei, these models have been extended to the deformed framework, e.g., the deformed relativistic Hartree–Bogoliubov (DRHB) theory [16, 49, 50] and the coordinate-space Skyrme–HFB approach [51-54].

Traditionally, these H(F)B equations are often solved in configuration spaces, i.e., via the basis expansion method [55]. However, the calculations are closely related to the space size and the shape of the expanded basis. Although the harmonic oscillator basis, which has a tail in the shape of a Gaussian function, can efficiently describe stable nuclei, significant difficulties were encountered when it was applied to exotic nuclei. Bases with proper shapes, such as the Woods–Saxon basis [56] and the transformed harmonic oscillator basis [57, 58], are often used for the exotic nuclei. For example, to explore deformed halos [59, 60], the DRHB theory based on a Woods–Saxon basis [16, 49] was developed. In contrast to the basis expansion method, solving the HFB equation in the coordinate space is believed to be more effective. In the coordinate space, the discretized method with the box boundary condition has been widely used, whereby a series of discrete quasiparticle levels can be easily obtained. However, flaws in this method have been identified, such as the nonphysical drops in the nuclear densities at the box boundary, the discretization of continua and resonant states, and the inclusion of unphysical states. In contrast, the Green's function (GF) method [61-63] in the coordinate space can avoid these problems and has significant advantages, i.e., it can describe the asymptotic behaviors of wave functions properly, provide the energies and widths of the resonant states directly, and treat the bound states and the continua on the same footing.

Owing to these advantages, the GF method has been applied extensively in nuclear physics to study the contribution of the continuum to the nuclear structures and excitations. In 1987, Belyaev et al. formulated the GF for the HFB equation [64]. Subsequently, this HFB GF was applied to the quasiparticle random-phase approximation [65], which was further used to describe the collective excitations coupled to the continuum [66-71]. In 2009, the continuum HFB theory in a coupled channel representation was developed to explore the effects of the continuum and pairing correlation in deformed neutron-rich Mg isotopes [72]. In 2011, Zhang et al. developed the fully self-consistent continuum Skyrme–HFB theory with the GF method [73], which was applied to investigate the giant halos [74] and the effects of the pairing correlation on the quasiparticle resonances [75, 76]. In 2019, to explore the halo phenomena in neutron-rich odd-A nuclei, the self-consistent continuum Skyrme–HFB theory was extended by including the blocking effect [77]. In recent years, the GF method has also been adopted for the covariant density functional theory [78-86] in studies on nuclear structure. For example, by introducing the GF method to the relativistic mean field theory (GF-RMF), the single-particle level structures, including the bound states and resonant states and the pseudospin symmetries therein, were investigated for neutrons [87, 88], protons [89], and Λ hyperons [90]. Additionally, it was confirmed that exact values of the energies and widths could be obtained by searching for the poles of the GF or the extremes of the density of states in disregard of the widths of resonant states [91, 92]. By combining the GF method with the RCHB theory, the pairing correlation and continuum are well described in the giant halos of the Zr isotopes [93]. By extending the GF-RMF model to the coupled channel representation, the halo candidate nucleus ³⁷Mg reported experimentally was analyzed and confirmed to be a p-wave one-neutron halo according to the Nilsson levels [94]. In addition, the complex-scaled GF method [95] has been established as a powerful tool for the exploration of resonant states, which was extended to the framework of the relativistic mean field [96] and deformed nuclei [97]. Additionally, the RMF-CMR-GF approach was developed by combining the complex momentum representation method with the GF method in the relativistic mean-field framework to study the halo structures in neutron-rich nuclei [98]. These studies proved the effectiveness of the GF method for the description of the continuum.

In this study, the neutron-rich Ca, Ni, Zr, and Sn isotopes are investigated systematically using the continuum Skyrme–HFB theory formulated with the GF method in the coordinate space, in which the pairing correlations, the couplings with the continuum, and the blocking effect for the odd unpaired nucleon are treated properly. The remainder of the paper is organized as follows. In Sec. 2, we briefly introduce the continuum Skyrme–HFB theory. Numerical details are presented in Sec. 3. The results and discussions are presented in Sec. 4, and conclusions are drawn in Sec. 5.

Theoretical framework

In the Hartree–Fock–Bogoliubov (HFB) theory [45], the pair correlated nuclear system is described in terms of independent quasiparticles by the Bogoliubov transformation. The HFB equation in the coordinate space [38] is $\left(\begin{array}{cc}h-\lambda & \tilde{h}\\ {\tilde{h}}^{*}& -h^{*}+\lambda \end{array}\right){\varphi }_i\left(r\sigma \right)=E_i{\varphi }_i\left(r\sigma \right)\mathrm{,}$ (1) where Ei represents the quasiparticle energy, ${\varphi }_i(r\sigma )$ represents the quasiparticle wave function, and λ represents the Fermi energy, which is determined by constraining the expectation value of the nucleon number. The HF Hamiltonian $h(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime })$ and the pair Hamiltonian $\tilde{h}(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime })$ are obtained by the variation of the total energy functional with respect to the particle density $\rho (r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime })$ and the pair density $\tilde{\rho }(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime })$, respectively. The solutions of the HFB equations have two symmetric branches: one is positive (E_i>0) with the quasiparticle wave function ${\varphi }_i(r\sigma )$, and the other is negative (-E_i<0) with the conjugate wave function ${\overline{\varphi }}_{\tilde{i}}\left(r\sigma \right)$. The wave functions ${\varphi }_i\left(r\sigma \right)$ and ${\overline{\varphi }}_{\tilde{i}}\left(r\sigma \right)$ are in two components and are related as follows: ${\varphi }_i\left(r\sigma \right)=\left(\begin{array}{c}{\phi }_{\mathrm{1,}i}\left(r\sigma \right)\\ {\phi }_{\mathrm{2,}i}\left(r\sigma \right)\end{array}\right),{\overline{\varphi }}_{\tilde{i}}(r\sigma )=\left(\begin{array}{c}-{\phi }_{\mathrm{2,}i}^{*}\left(r\tilde{\sigma }\right)\\ {\phi }_{\mathrm{1,}i}^{*}\left(r\tilde{\sigma }\right)\end{array}\right),$ (2) with $\phi \left(r\tilde{\sigma }\right)=-2\left(r\mathrm{,}-\sigma \right)$. In this paper, we follow the notation used in Ref. [65] for convenience.

For an even-even nucleus, the ground state $|{\Phi }_0\rangle$ is a quasiparticle vacuum with the particle density $\rho \left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)$ and pairing density $\tilde{\rho }\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)$ determined as follows: $\rho \left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)\equiv \langle {\Phi }_0\left|c_{r\text{'}{\sigma }^{\prime }}^{†}{c}_{r\sigma }\right|{\Phi }_0\rangle \mathrm{,}$ (3a) $\tilde{\rho }\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)\equiv \langle {\Phi }_0\left|c_{r\text{'}\tilde{\sigma }\text{'}}{c}_{r\sigma }\right|{\Phi }_0\rangle \mathrm{,}$ (3b) where $c_{r\sigma }^{†}$ and $c_{r\sigma }$ are the particle creation and annihilation operators, respectively. The densities can be unified as a generalized density matrix $R\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)$, with $\rho \left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)$ and $\tilde{\rho }\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)$ being the “11” and “22” elements, respectively. With the quasiparticle wave functions, $R\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)$ can be written in a simple form: $R\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)={\displaystyle \sum _i{\overline{\varphi }}_{\tilde{i}}}\left(r\sigma \right){\overline{\varphi }}_{\tilde{i}}^{†}\left(r^{\prime }{\sigma }^{\prime }\right)\mathrm{.}$ (4) For an odd-A nucleus, the last odd nucleon is unpaired, for which the blocking effect should be considered. The nuclear ground state in this case is a one-quasiparticle state $|{\Phi }_1\rangle$, which can be constructed on the basis of a HFB vacuum $|{\Phi }_0\rangle$ as $|{\Phi }_1\rangle ={\beta }_{i_{\text{b}}}^{†}|{\Phi }_0\rangle \mathrm{,}$ (5) where ${\beta }_{i_{\text{b}}}^{†}$ is the quasiparticle creation operator and i_b denotes the blocked quasiparticle level occupied by the odd nucleon. Accordingly, the particle density $\rho (r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime })$ and the pairing density $\tilde{\rho }\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)$ are $\rho \left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)\equiv \in \langle {\Phi }_1\left|c_{r\text{'}{\sigma }^{\prime }}^{†}{c}_{r\sigma }\right|{\Phi }_1\rangle \mathrm{,}$ (6a) $\tilde{\rho }\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)\equiv \langle {\Phi }_1\left|c_{r\text{'}\tilde{\sigma }\text{'}}{c}_{r\sigma }\right|{\Phi }_1\rangle \mathrm{,}$ (6b) and the generalized density matrix $R\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)$ becomes $\begin{array}{c}R\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)={\displaystyle \sum _{i:\text{all}}{\overline{\varphi }}_{\tilde{i}}}\left(r\sigma \right){\overline{\varphi }}_{\tilde{i}}^{†}\left(r^{\prime }{\sigma }^{\prime }\right)\\ -{\overline{\varphi }}_{{\tilde{i}}_{\text{b}}}\left(r\sigma \right){\overline{\varphi }}_{{\tilde{i}}_{\text{b}}}^{†}\left(r^{\prime }{\sigma }^{\prime }\right)+{\varphi }_{i_{\text{b}}}\left(r\sigma \right){\varphi }_{i_{\text{b}}}^{†}\left(r^{\prime }{\sigma }^{\prime }\right)\mathrm{,}\end{array}$ (7) where two more terms are introduced compared with those for the even-even nuclei.

In the conventional Skyrme–HFB theory, the HFB equation (1) in the coordinate space is often solved with the box boundary condition, and a series of discretized eigensolutions including the quasiparticle energy Ei and the corresponding wave functions ${\varphi }_i(r\sigma )$ can be obtained. Then, the generalized density matrix $R\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)$ can be calculated by summing these discretized quasiparticle states, in accordance with Eqs. (5) and (9). We call this method the box-discretized approach. However, the applicability of the box boundary condition in the description of exotic nuclei—particularly those close to the drip line—is poor. A sufficiently large coordinate space (or box size) should be used to describe the extended density distribution.

The GF method can avoid these problems of the box-discretized approach, as it imposes the correct asymptotic behaviors on the wave functions—particularly for the weakly bound sates and the continuum. The GF $G\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime };E\right)$ with an arbitrary quasiparticle energy E defined for the coordinate-space HFB equation obeys $\begin{array}{l}\left[E-\left(\begin{array}{cc}h-\lambda & \tilde{h}\\ {\tilde{h}}^{*}& -h^{*}+\lambda \end{array}\right)\right]G\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime };E\right)\\ =\delta \left(r-r^{\prime }\right){\delta }_{\sigma {\sigma }^{\prime }}\mathrm{,}\end{array}$ (8) which is a 2× 2 matrix. The generalized density matrix $R\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)$ in Eq. (9) can be calculated by taking the integrals of the GFs on the complex quasiparticle energy plane, as follows: $\begin{array}{c}R\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)=\frac{1}{2\pi i}[{\displaystyle {\oint }_{C_{E<0}}}\text{d}EG\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime };E\right)\\ -{\displaystyle {\oint }_{C_{\text{b}}^{-}}}\text{d}EG\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime };E\right)\\ +{\displaystyle {\oint }_{C_{\text{b}}^{+}}}\text{d}EG\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime };E\right)]\mathrm{,}\end{array}$ (9) where the contour path C_E<0 encloses all the negative quasiparticle energies -E_i<0, $C_{\text{b}}^{-}$ encloses only the pole of $-E_{i_{\text{b}}}$, and $C_{\text{b}}^{+}$ encloses only the pole of $E_{i_{\text{b}}}$.

In the spherical case, the quasiparticle wave functions ${\varphi }_i\left(r\sigma \right)$ and ${\overline{\varphi }}_{\tilde{i}}\left(r\sigma \right)$ are only dependent on the radial parts, and they can be expanded as follows: ${\varphi }_i\left(r\sigma \right)=\frac{1}{r}{\varphi }_{nlj}\left(r\right)Y_{jm}^l\left(\widehat{r}\sigma \right)\mathrm{,}$ (10a) ${\overline{\varphi }}_{\tilde{i}}\left(r\sigma \right)=\frac{1}{r}{\overline{\varphi }}_{nlj}\left(r\right)Y_{jm}^{l*}\left(\widehat{r}\tilde{\sigma }\right)\mathrm{,}$ (10b) ${\varphi }_{nlj}\left(r\right)=\left(\begin{array}{c}{\phi }_{\mathrm{1,}nlj}\left(r\right)\\ {\phi }_{\mathrm{2,}nlj}\left(r\right)\end{array}\right)\mathrm{,}$ (10c) ${\overline{\varphi }}_{nlj}\left(r\right)=\left(\begin{array}{c}-{\phi }_{\mathrm{2,}nlj}^{*}\left(r\right)\\ {\phi }_{\mathrm{1,}nlj}^{*}\left(r\right)\end{array}\right)\mathrm{,}$ (10d) where $Y_{jm}^l\left(\widehat{r}\sigma \right)$ is the spin spherical harmonic, and $Y_{jm}^l\left(\widehat{r}\tilde{\sigma }\right)=-2\sigma Y_{jm}^l\left(\widehat{r}-\sigma \right)$. Similarly, the generalized density matrix $R\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)$ and the GF $G\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime };E\right)$ can be expanded as $R\left(r\sigma \mathrm{,}{r}^{\prime }{\sigma }^{\prime }\right)={\displaystyle \sum _{ljm}{Y}_{jm}^l}\left(\widehat{r}\sigma \right)R_{lj}\left(r\mathrm{,}{r}^{\prime }\right)Y_{jm}^{l*}\left(\widehat{r}\text{'}{\sigma }^{\prime }\right)\mathrm{,}$ (11a) $G(r\sigma \mathrm{,}{r}^{\prime }{\sigma }^{\prime };E)={\displaystyle \sum _{ljm}}{Y}_{jm}^l(\widehat{r}\sigma )\frac{{\mathcal{G}}_{lj}(r\mathrm{,}{r}^{\prime };E)}{r{r}^{\prime }}{Y}_{jm}^{l*}(\widehat{r}\text{'}{\sigma }^{\prime })\mathrm{,}$ (11b) where $R_{lj}\left(r\mathrm{,}{r}^{\prime }\right)$ and ${\mathcal{G}}_{lj}\left(r\mathrm{,}{r}^{\prime };E\right)$ are the radial parts of the generalized density matrix and GF, respectively.

Thus, the radial local generalized density matrix R(r)=R(r,r) can be expressed by the radial HFB GF ${\mathcal{G}}_{lj}\left(r\mathrm{,}{r}^{\prime };E\right)$ as follows: $\begin{array}{c}R\left(r\right)={\displaystyle \sum _{lj}{R}_{lj}\left(r\mathrm{,}r\right)}\\ =\frac{1}{4\pi r^2}[{\displaystyle \sum _{lj:\text{all}}\left(2j+1\right)}{\displaystyle \sum _{n:\text{all}}{\overline{\varphi }}_{nlj}^2}\left(r\right)-{\overline{\varphi }}_{n_{\text{b}}{l}_{\text{b}}{j}_{\text{b}}}^2\left(r\right)\\ +{\varphi }_{n_{\text{b}}{l}_{\text{b}}{j}_{\text{b}}}^2\left(r\right)]\\ =\frac{1}{4\pi r^2}\frac{1}{2\pi i}[{\displaystyle \sum _{lj:\text{all}}\left(2j+1\right)}{\displaystyle {\oint }_{C_{E<0}}}dE{\mathcal{G}}_{lj}\left(r\mathrm{,}r;E\right)\\ -{\displaystyle {\oint }_{C_{\text{b}}^{-}}}dE{\mathcal{G}}_{l_{\text{b}}{j}_{\text{b}}}\left(r\mathrm{,}r;E\right)+{\displaystyle {\oint }_{C_{\text{b}}^{+}}}dE{\mathcal{G}}_{l_{\text{b}}{j}_{\text{b}}}\left(r\mathrm{,}r;E\right)]\mathrm{.}\end{array}$ (12) From the radial generalized matrix R(r), one can easily obtain the radial local particle density ρ(r) and pair density $\tilde{\rho }(r)$, which are the “11” and “12” components of R(r), respectively. In the same way, one can express other radial local densities needed in the functional of the Skyrme interaction, such as the kinetic-energy density τ(r) and the spin–orbit density J(r), in terms of the radial GF. For the construction of the GF, see Refs. [73, 77].

Numerical details

For the Skyrme interaction in the ph channel, the SLy4 parameter set [100] is adopted. For the pairing interaction in the pp channel, a density-dependent δ interaction (DDDI) is used: $v_{\text{pair}}\left(r\mathrm{,}{r}^{\prime }\right)=\frac{1}{2}\left(1-P_{\sigma }\right)V_0\left[1-\eta {\left(\frac{\rho \left(r\right)}{{\rho }_0}\right)}^{\alpha }\right]\delta \left(r-r^{\prime }\right)\mathrm{.}$ (13) The pair Hamiltonian $\tilde{h}\left(r\sigma \mathrm{,}r\text{'}{\sigma }^{\prime }\right)$ is then reduced to a local pair potential [47]: $\Delta \left(r\right)=\frac{1}{2}{V}_0\left[1-\eta {\left(\frac{\rho \left(r\right)}{{\rho }_0}\right)}^{\alpha }\right]\tilde{\rho }\left(r\right)\mathrm{.}$ (14) The strength of the pairing force V₀=-458.4 MeV⋅fm³, density ρ₀=0.08 fm^-3, and parameters η=0.71 and α=0.59 are constrained by reproducing the experimental neutron pairing gaps for the Sn isotopes [101, 102, 68]. With these parameters, the DDDI can reproduce the scattering length a=-18.5 fm in the ¹S channel of the bare nuclear force in the low density limit [101]. The cutoff of the quasiparticle states is set to a maximal angular momentum of $j_{\text{max}}=25/2$ and a maximal quasiparticle energy of E_cut=60 MeV.

The HFB equation is solved in the coordinate space with the space size $R_{}=20$ fm and mesh size dr=0.1 fm. To calculate the densities with the GF, the integrals of the GFs are performed along a contour path C_E<0, which is set as a rectangle with height γ = 0.1 MeV and length E_cut = 60 MeV to enclose all the quasiparticle states with negative energies. For the odd-A nuclei, two more contour paths $C_{\text{b}}^{+}$ and $C_{\text{b}}^{-}$, which only enclose the blocked quasiparticle states at energies $E_{i_{\text{b}}}$ and $-E_{i_{\text{b}}}$, are introduced owing to the blocking effect of the odd unpaired nucleon. Details are presented in Ref. [77]. To perform the contour integration, an energy step of Δ E = 0.01 MeV on the contour path is adopted.

Results and discussion

In Fig. 1, the two-neutron separation energies S_2n(N,Z) =E(N-2,Z) - E(N, Z) are plotted for the even-even and odd-even Ca, Ni, Zr, and Sn isotopes. Red circles indicate those calculated according to the continuum Skyrme–HFB theory with the SLy4 parameter set, blue triangles indicate the results of the discretized method, and black squares indicate the available experimental data [99]. The differences between the S_2n values obtained via the GF method and discretized method are small. Good agreement with the experimental data is observed, indicating the reliability of the continuum Skyrme–HFB theory for the prediction of neutron drip line. The traditional shell closures, i.e., N=28 in Ca isotopes, N=50 in Ni isotopes, N=50, 82 in Zr isotopes, and N=82 in Sn isotopes, can be identified, where the S_2n decreases sharply. For example, in the Sn chain, S_2n decreases from 13.25 MeV at ¹³²Sn to 4.94 MeV at ¹³⁴Sn with the neutron number exceeding the magic number N=82. In the Ca, Ni, and Zr chains, the two-neutron separation energies quickly reach 0 as mass increasing, resulting in relatively short neutron drip lines, which are ⁶⁷Ca, ⁸⁹Ni, and ¹²³Zr, respectively. In contrast, in the Sn chain, S_2n remains below 1.0 MeV in a wide mass region after the gap of N=82 and finally becomes negative until A=178, suggesting that ¹⁷⁷Sn is a neutron drip line nucleus. These weakly bound nuclei are interesting owing to the possible appearance of neutron halos, although this is experimentally difficult to achieve. In addition, the exploration of the neutron drip line and the determination of the limit of the nuclear landscape are important in nuclear physics. However, various theoretical studies indicate that the predicted neutron drip line is very model-dependent [103]. Moreover, different physical quantities and criteria yield different neutron drip line predictions.

Fig. 1Full-screen View

(Color online) Two-neutron separation energies S_2n in the Ca, Ni, Zr, and Sn isotopes as a function of the mass number A calculated according to the continuum Skyrme–HFB theory with the SLy4 parameter set (filled red circles), in comparison with the results of the discretized method (open blue triangles) and experimental data (filled black squares) [99].

To explore the neutron drip lines in the Ca, Ni, Zr, and Sn isotopes, in Fig. 2 the single-neutron separation energies $S_{\text{n}}\left(N\mathrm{,}Z\right)=E\left(N-\mathrm{1,}Z\right)-E\left(N\mathrm{,}Z\right)$ are plotted. The results obtained using the continuum Skyrme–HFB theory with the SLy4 parameter set are indicated by red circles, which are consistent with the experimental data [99] indicated by the black squares. Strong odd-even staggering is observed in all isotopes. In general, the S_n in the even-even nucleus is approximately 2~3 MeV larger than those in the neighboring odd-A nuclei, which is attributed to the unpaired odd neutron with vanishing pairing energy. Consequently, compared with those in Fig. 1, the neutron drip lines determined via the one-neutron separation energy are significantly shortened. In the Ca, Ni, Zr, and Sn isotopes, the drip line nuclei are ⁶⁰Ca, ⁸⁶Ni, ¹²²Zr, and ¹⁴⁸Sn, respectively, whose positions are indicated by the black arrows. Outside the neutron drip line determined by S_n, the bound even-even nuclei behave as interesting Borromean systems. For example, considering the bound nucleus ⁶⁰Ca, ⁶⁰Ca+n is unbound, while ⁶⁰Ca+n+n is bound.

Fig. 2Full-screen View

(Color online) Single-neutron separation energies S_n as a function of the mass number A in the (a) Ca, (b) Ni, (c) Zr, and (d) Sn isotopes calculated according to the continuum Skyrme–HFB theory with the SLy4 parameter set (filled red circles), in comparison with the results of the discretized method (open blue triangles) and experimental data (filled black squares) [99].

In Fig. 3, we plot the neutron pairing energy E_pair, which is expressed as $E_{\text{pair}}=\frac{1}{2}{\displaystyle \int }\text{d}r\Delta \left(r\right)\tilde{\rho }\left(r\right)\mathrm{.}$ (15) The red solid symbols correspond to the even-even nuclei, and the open symbols correspond to the odd-A nuclei. The neutron pairing energies E_pair of the odd-A nuclei are obviously lower than those of the neighboring even-even nuclei, owing to the absent contribution of pairing energy from the unpaired neutron. This also explains why the drip line determined via the single-neutron separation energy S_n is much shorter than that obtained via the two-neutron separation energy S_2n. In addition, obvious shell effects are observed in the pairing energy. For example, in Sn isotopes, the pairing energy is 0 at N=82 and N=126, and it is maximized at the half-shell N=102. As a result, the traditional shell closures, i.e., N=28, 40 in Ca isotopes, N=40, 50 in Ni isotopes, N=50, 82 in Zr isotopes, and N=82,126 in Sn isotopes, can also be clearly observed, which are consistent with those shown in Fig. 1. In addition, a sub-shell N=32 is observed in Ca isotopes.

Fig. 3Full-screen View

(Color online) Neutron pairing energy -E_pair as a function of the mass number A in the (a) Ca, (b) Ni, (c) Zr, and (d) Sn isotopes calculated according to the continuum Skyrme–HFB theory with the SLy4 parameter set. The filled and open squares indicate the results for the even-even and odd-even nuclei, respectively.

In the following, the possible neutron halos in the Ca, Ni, Zr, and Sn isotopes are examined—particularly those in the weakly bound nuclei close to the neutron drip line, where the Fermi surfaces are very close to the continuum threshold and the valence neutrons can be easily scattered to the continuum by the pairing correlation.

In Fig. 4, the neutron canonical single-particle structure as a function of the mass number is plotted for the (a) Ca, (b) Ni, (c) Zr, and (d) Sn isotopes. Details on obtaining the canonical single-particle levels are presented in Refs. [104, 105]. The neutron Fermi energy λ_n as well as the canonical single-particle energies ε are shown. As the neutron number increases, the Fermi energy λ_n in each chain increases, finally reaching the continuum threshold, while all the HF single-particle levels decrease. The traditional shell closures, i.e., N=28, 40 in Ca isotopes, N=40, 50 in Ni isotopes, N=82 in Zr isotopes, and N=82 in Sn isotopes, can be observed, exhibiting large gaps. Different single-particle structures are revealed in the Ca, Ni, Zr, and Sn chains, whereby halos may be formed. In the Ni chain, above the shell closure of N=50, there are several weakly bound states and low-lying positive canonical states in the continuum with small angular momenta, which favor the formation of halos. For example, in ⁸⁶Ni, around the Fermi surface, there are two weakly bound states, i.e., 2d_5/2 and 3s_1/2, and one low-lying positive canonical state 2d_3/2. The Sn chain is similar to the Ni chain but has more advantages for halo formation, where weakly bound states and the low-lying positive states in the continuum exist above the N=82 shell closure. Additionally, the Fermi surface λ_n gradually approaches 0, and the Sn isotopes in a large mass region are weakly bound. In the Ca chain, above the shell closure of N=40, the main state is 1g_9/2, which evolves from a canonical positive state in the continuum ($A\le 62$) to a weakly bound level ($A\ge 64$). Although there is a possibility of valence neutrons occupying the 1g_9/2 orbital, the contributed density is very localized owing to the large central barrier. In the neutron-rich Ca isotopes, the low-lying positive canonical states in the continuum 3s_1/2 and 2d_5/2 also play important roles. The formation of halos is most unlikely for the Zr chain, where the shell closure of N=82 is located around the threshold of the continuum and it is difficult for the valence neutrons to overcome the large gap and occupy the continuum.

Fig. 4Full-screen View

(Color online) Single-particle levels for neutrons around the Fermi surface in the canonical basis for the (a) Ca, (b) Ni, (c) Zr, and (d) Sn isotopes calculated according to the continuum Skyrme–HFB theory. The short dashed lines represent the Fermi surface.

In Fig. 5, the neutron root-mean-square (rms) radii $r_{\text{rms}}=\sqrt{\frac{{\displaystyle \int }\text{d}r4\pi r^4\rho \left(r\right)}{{\displaystyle \int }\text{d}r4\pi r^2\rho \left(r\right)}}\mathrm{,}$ (16) are plotted for the Ca, Ni, Zr, and Sn isotopes, which are based on Skyrme–HFB calculations with the GF method (filled circles and solid lines) and box-discretized method (open triangles and dashed lines) employing the SLy4 parameter set. The radii r=b₀ A^1/3 in the traditional liquid-drop model (black lines) are shown as well, with the coefficient b₀ determined via the radii of deeply bound nuclei. In the Ca, Ni, Zr, and Sn chains, they are $r\approx 0.991A^{1/3}$, 0.984A^1/3, 0.957A^1/3, and 0.961 A^1/3, respectively. In each chain, with the addition of neutrons, the nuclear rms radii r_rms increase steeply and deviate from the radii A^1/3 rule. For example, in the Ca chain, compared with the isotopic trend in $N⩽20$ with $r\approx 0.991A^{1/3}$, the neutron rms radii in ⁵⁰Ca and the heavier isotopes exhibit steep increases with an increase in N. In this mass range, possible neutron halos may occur. Additionally, odd-even staggering of rms radii can be clearly observed in the Sn chain, where the odd-A nuclei ^151-165Sn have larger rms radii than the neighboring two even-even nuclei. Details are presented in Ref. [77]. The odd-even staggering phenomena in nuclear radii and nuclear mass have attracted considerable research interest in recent years [106-112].

Fig. 5Full-screen View

(Color online) Neutron rms radii r_rms for the (a) Ca, (b) Ni, (c) Zr, and (d) Sn isotopes calculated according to the Skyrme–HFB theory with the GF method (filled circles and solid lines), in comparison with the results of the discretized method (open triangles and dashed lines).

In exotic nuclei, diffuse density distributions in the coordinate space are often observed. Thus, in Fig. 6, to explore the exotic structures in the (a) Ca, (b) Ni, (c) Zr, and (d) Sn isotopes, we also plot the neutron density distributions ρ(r), where the solid and dashed lines indicate those obtained via the GF method and the box-discretized method, respectively. As a global trend, the neutron density distributions are extended with an increase in the neutron number. The shell structures significantly affect the density distribution; i.e., compared with the bound nuclei, the density distributions of the neutron-rich nuclei in the (a) Ca, (b) Ni, (c) Zr, and (d) Sn isotopes with the neutron number exceeding the neutron closure N=28, 50, 50, 82 are far more extended, which is consistent with the behaviors of the rms radii plotted in Fig. 5. In addition, compared with the Ca and Zr chains, the Ni and Sn chains exhibit more diffuse density distributions, which can be explained by their small two-neutron separation energies S_2n in a large mass range, as shown in Fig. 1. For the Zr isotopes, the density distributions are relatively localized. The large S_2n shown in Fig. 1 suggests the absence of halos in Zr isotopes. However, according to the RCHB theory with the NLSH parameter set [113] and the continuum Skyrme–HFB theory with the SKI4 parameter set [114, 74], giant halos in Zr isotopes have been predicted. In all isotopes, compared with the box-discretized method predicting nonphysical sharp reductions in density at the space boundary, the GF method can better describe the extended density distributions—particularly for very neutron-rich isotopes. Additionally, the densities obtained via the GF method can be independent of the space sizes, as discussed in Ref. [77], which are mainly determined by the proper boundary conditions of the bound states, weakly bound sates, and the continuum employed when constructing the GFs.

Fig. 6Full-screen View

(Color online) Neutron density ρ(r) for the (a) Ca, (b) Ni, (c) Zr, and (d) Sn isotopes calculated according to the Skyrme–HFB theory with the GF method (solid lines) and the discretized method (dashed lines).

To explore the contributions of different partial waves to the extended density distributions in Fig. 6, taking the neutron-rich (a) ⁶⁴Ca, (b) ⁸⁶Ni, (c) ¹²⁰Zr, and (d) ¹⁷⁴Sn as examples, we plot in Fig. 7 the compositions ${\rho }_{lj}\left(r\right)/\rho \left(r\right)$ as functions of the radial coordinate r. As shown, outside the nuclear surface referring to the right boundary of the shallow regions of the total nuclear density distributions, the orbitals located around the Fermi surface have the most significant effect on the density distributions. For example, in neutron-rich ⁶⁴Ca, the partial waves p_1/2, f_5/2, g_9/2, s_1/2, and d_5/2 contribute significantly to the total density in the area of $5\text{fm}<r<15$. These levels are located within ~5 MeV around the Fermi surface, as shown in Fig. 8, where the neutron canonical single-particle levels as well as the occupation probabilities are presented. As we move further in the coordinate space with r>15 fm, the contributions of the partial waves s_1/2 and d_5/2 with small angular momenta increase significantly, while the contributions of other partial waves decrease. In the case of ⁸⁶Ni, the single-particle levels 2d_5/2, 3s_1/2, and 2d_3/2 are located above the neutron shell of N=50 and close to the Fermi surface, playing the main role in the neutron density distribution. Although the positive state 1g_7/2 is also very close to the Fermi surface, it contributes little to the density in the large coordinate space owing to the large centrifugal barrier. For the nucleus ¹²⁰Zr with the neutron number very close to the closure of N=82, the single-particle levels between the closures N=50 and N=82, including 2d_5/2, 1g_7/2, 3s_1/2, 2d_3/2, and 1h_11/2, contribute significantly to the densities in the large coordinate space. The insignificant occupation of the positive state 2f_7/2 leads to a small contribution to the density. Regarding the neutron-rich ¹⁷⁴Sn with the neutron number exceeding the closure of N=82, the weakly bound single-particle levels 2f_7/2, 2f_5/2, 3p_3/2, and 3p_1/2 with small angular momenta play the key role in the extended density distributions in the large coordinate space. From this analysis, we can conclude that the single-particle levels around the Fermi surface—particularly the waves with small angular momenta—are the main cause of the extended nuclear density distributions.

Fig. 7Full-screen View

(Color online) Contributions of different partial waves to the total neutron density ${\rho }_{lj}\left(r\right)/\rho \left(r\right)$ as a function of the radial coordinate r for the neutron-rich nuclei (a) ⁶⁴Ca, (b) ⁸⁶Ni, (c) ¹²⁰Zr, and (d) ¹⁷⁴Sn. The shallow regions are for the nuclear density distributions ρ(r), which are rescaled by multiplying by a factor of 5.

Fig. 8Full-screen View

(Color online) Neutron canonical single-particle levels in (a) ⁶⁴Ca, (b) ⁸⁶Ni, (c) ¹²⁰Zr, and (d) ¹⁷⁴Sn calculated according to the Skyrme–HFB theory with the GF method. The blue dashed lines represent the Fermi surface. The occupation probabilities (v²) for the canonical orbitals are proportional to the lengths of the lines (logarithmic scale).

In Fig. 8, the particle occupation probabilities v² on different canonical levels ${\epsilon }_{\text{N}}^{\text{can}}$ are indicated by the lengths of the lines. Without pairing, the values of the occupation probabilities v² should be either 1 or 0, separated by the Fermi surface. With the effect of pairing, the nucleons occupying the levels below the Fermi surface can be scattered to higher levels, which results in the occupation of the weakly bound states above the Fermi surface and even levels in the continuum. In the neutron-rich nuclei ⁶⁴Ca, ⁸⁶Ni, ¹²⁰Zr, and ¹⁷⁴Sn, the numbers of neutrons $N_{\lambda }={\displaystyle {\sum }_{{\epsilon }_k>{\lambda }_{\text{n}}}\left(2j+1\right)}{v}_k^2$ scattered above the Fermi surface λ_n are 4.33, 2.51, 0.159, and 0.173, respectively. In the case of ⁶⁴Ca, the weakly bound single-particle 1g_9/2 contributes approximately 3.7 neutrons.

To explore the effects of pairing, we plot in Fig. 9 the numbers of neutrons N_λ scattered above the Fermi surface for the (a) Ca, (b) Ni, (c) Zr, and (d) Sn chains obtained via the Skyrme–HFB theory with the GF method. As shown, large numbers of neutrons are scattered from the single-particle levels below the Fermi surface to the weakly bound states above the Fermi surface and even levels in the continuum owing to the pairing—particularly in the nuclei with the neutron number filling the half-full shells. Additionally, an obvious shell structure is observed. When the number of neutrons reaches a magic number, i.e., N=28, 40 in Ca isotopes, N=40, 50 in Ni isotopes, N=50, 82 in Zr isotopes, and N=82,126 in Sn isotopes, N_λ is almost 0 owing to the absence of the pairing for the closed-shell nuclei. Additionally, at the points of N=32 in the Ca isotopes, N=54,68 in the Zr isotopes, and N=88 in the Sn isotopes, very small numbers of neutrons (N_λ) are obtained, indicating weak pairing in these nuclei, along with the possible existence of sub-shells and new magic numbers. Furthermore, the evolution of N_λ is consistent with the trend of the pair energy -E_pair in Fig. 3.

Fig. 9Full-screen View

(Color online) Number of neutrons N_λ occupying the single-particle levels above the Fermi surface λn as a function of the mass number A for the (a) Ca, (b) Ni, (c) Zr, and (d) Sn isotopes calculated according to the Skyrme–HFB theory with the GF method.

In Fig. 10, we further investigate the number of neutrons occupying the continuum with single-particle energies of ε>0 MeV, i.e., $N_0={\displaystyle {\sum }_{{\epsilon }_k>0}\left(2j+1\right)}{v}_k^2$. Compared with the number of neutrons N_λ occupying the levels above the Fermi surface, the number of neutrons occupying the continuum N₀ is significantly smaller. For example, in the Ca chain, N₀ is less than 1 in all isotopes except ⁶²Ca. In ⁶²Ca, the single-particle level 1g_g/2 appears as a low-lying canonical positive state in the continuum with energy of ε=0.198 MeV and a high occupation probability of v² = 0.197, resulting in almost 1.979 neutrons occupying it. However, in the neighboring ⁶⁰Ca, a very low occupation probability v²=0.015 of 1g_9/2 is obtained, and in ⁶⁴Ca, the 1g_9/2 state drops to a weakly bound level with energy of ε=-0.072 MeV. After the neutron contribution is removed from 1g_g/2 in ⁶²Ca, only 0.38 neutrons are in the continuum, which is denoted by an empty circle in panel (a). Except for the Sn chain, the shape of the evolution of N₀ is very close to those of the pairing energy and pairing gap. Although the case of the Sn chain becomes very complex, we can observe the shell structure at N=82 and N=126, where N₀ is almost 0.

Fig. 10Full-screen View

(Color online) Number of neutrons N₀ occupying in the continuum (above the threshold ε=0 MeV) as a function of the mass number A for the (a) Ca, (b) Ni, (c) Zr, and (d) Sn isotopes calculated according to the Skyrme–HFB theory with the GF method.

Summary

In this study, the exotic nuclear properties of neutron-rich Ca, Ni, Zr, and Sn isotopes were examined systematically according to the continuum Skyrme–HFB theory in the coordinate space formulated with the GF method, in which the pairing correlations, the couplings to the continuum, and the blocking effects for the unpaired nucleon in odd-A nuclei are treated properly.

First, the two-neutron separation energies S_2n and one-neutron separation energies S_n were calculated, which were consistent with experimental data. Significant differences exist for the drip lines determined by S_2n and S_n. In the Ca, Ni, Zr, and Sn isotopes, the drip line nuclei are ⁶⁷Ca, ⁸⁹Ni, ¹²³Zr, and ¹⁷⁷Sn according to S_2n and ⁶⁰Ca, ⁸⁶Ni, ¹²²Zr, and ¹⁴⁸Sn according to S_n. Owing to the absent contribution of pairing energy from the single unpaired odd neutron, the neutron pairing energies (E_pair) of the odd-A nuclei are approximately 2 MeV lower than those of the neighboring even-even nuclei. This explains why the drip lines determined via S_n are much shorter than those determined via S_2n. In addition, from the fluctuation trends of the pairing energy, the traditional neutron magic numbers are clearly displayed, i.e., N=28, 40 in Ca isotopes, N=40, 50 in Ni isotopes, N=50, 82 in Zr isotopes, and N=82,126 in Sn isotopes.

Second, to explore the possible halo structures in the neutron-rich Ca, Ni, Zr, and Sn isotopes, the neutron single-particle structures, the rms radii, and the density distributions were investigated. In the neutron-rich Ca, Ni, Sn nuclei—particularly the weakly bound nuclei close to neutron drip line—the rms radii increase sharply, with significant deviations from the traditional $r\propto A^{1/3}$ rule. Additionally, very diffuse spatial density distributions are observed in these nuclei, possibly indicating a halo phenomenon therein. By analyzing the contributions of different partial waves to the total density, we found that the orbitals located around the Fermi surface—particularly those with small angular momenta—are the main cause of the extended nuclear density and large rms radii.

Finally, the numbers of halo nucleons that can reflect the effects of pairing were examined. Two different numbers of neutrons were defined: N_λ neutrons occupying the single-particle levels above the Fermi surface λn and N₀ neutrons occupying the continuum. We found that the evolutions of N_λ and N₀ with respect to the mass number A are consistent with the trend of the pairing energy -E_pair, which supports the key role of the pairing correlations in the halo phenomena.

The theoretical calculations were supported by the nuclear data storage system at Zhengzhou University.