Optimized Dirac Woods-Saxon basis for covariant density functional theory

The study of exotic nuclei far from the $\beta$ stability line has been at the forefront of nuclear physics research since the discovery of the first halo nucleus ¹¹Li in the 1980s Tanihata et al. (1985). In addition to halos, novel phenomena found in exotic nuclei include the pygmy resonance Adrich et al. (2005), the disappearance of traditional magic numbers and the emergence of new ones Ozawa et al. (2000). These phenomena not only promote the worldwide development of radioactive ion beam facilities but also challenge conventional nuclear models.

The nuclear density functional theory (DFT) is able to provide a unified description for almost all nuclei in the nuclear chart and has become one of the most important microscopic methods for the study of nuclear structure Bender et al. (2003). Its relativistic version, i.e., the covariant density functional theory (CDFT), has attracted wide attention in recent years for many advantages, such as the automatic inclusion of the nucleonic spin degree of freedom and the spin-orbital interaction Ren and Zhao (2020), the explanation of the pseudospin symmetry in the nucleon spectrum Ginocchio (1997); Meng et al. (1998, 1999); Chen et al. (2003); Ginocchio (2005); Liang et al. (2015) and the spin symmetry in the antinucleon spectrum Zhou et al. (2003a); He et al. (2006); Liang et al. (2015), and the natural inclusion of the nuclear magnetism Koepf and Ring (1989), which plays an important role in nuclear magnetic moments Yao et al. (2006); Arima (2011); Li et al. (2011a, b); Li and Meng (2018) and nuclear rotations Meng et al. (2013); König and Ring (1993); Afanasjev et al. (2000); Afanasjev and Ring (2000); Afanasjev and Abusara (2010); Zhao et al. (2011a, b, 2012, 2015); Wang (2017, 2018); Ren et al. (2019, 2020).

In the nuclear DFT, the equations of motion for (quasi)nucleons are either solved in coordinate space or transformed into a matrix diagonalization problem in a complete basis, e.g., the harmonic oscillator (HO) basis. As a good approximation for the mean field of stable nuclei, the HO potential can be solved analytically, and the HO wave function in analytical form may bring convenience in the calculation of some matrix elements, e.g., the separable pairing force Tian et al. (2009). Due to the incorrect asymptotic behavior of HO wave functions, however, the expansion in a localized HO basis is incapable of describing nuclei with very diffuse spatial density distributions. To improve the asymptotic behavior of HO wave functions, a transformed HO basis has been proposed in Refs. Stoitsov et al. (1998a, b) via a local scaling transformation.

Solutions in coordinate space can describe properly the asymptotic behavior of wave functions. When dealing with systems having small separation energy, the coordinate-space calculations in large boxes were found to be more effective than the transformed HO basis Zhang et al. (2013). In the CDFT, the relativistic Hartree equation that neglects pairing correlations has been successfully solved on a three-dimensional lattice Ren et al. (2017). For the treatment of pairing correlations in weakly bound nuclei, it was shown that the Bogoliubov transformation is more suitable than the widely used Bardeen-Cooper-Schrieffer (BCS) method Dobaczewski et al. (1984); Meng (1998). However, the relativistic Hartree-Bogoliubov equation in coordinate space has only been solved by assuming spherical symmetry Meng and Ring (1996); Pöschl et al. (1997); Meng and Ring (1998). Solving the deformed relativistic Hartree-Bogoliubov equation in coordinate space is extremely difficult if not impossible Zhou et al. (2000).

A Woods-Saxon basis was proposed in Ref. Zhou et al. (2003b) as a reconciler between the HO basis and coordinate space. It is obtained by solving the Schrödinger equation or the Dirac equation containing spherical Woods-Saxon potentials with box boundary conditions, and is correspondingly referred to as the Schrödinger Woods-Saxon (SWS) or Dirac Woods-Saxon (DWS) basis. The Woods-Saxon basis has the advantage in providing appropriate asymptotic behaviors of wave functions because of the nature of the Woods-Saxon type potential Woods and Saxon (1954). It was shown in Ref. Zhou et al. (2003b) that for spherical systems the solution of the relativistic Hartree equations in the Woods-Saxon basis is almost equivalent to the solution in coordinate space. Up to the present, the SWS basis has been applied to the spherical Hartree-Fock-Bogoliubov theory Schunck and Egido (2008a, b), and the DWS basis to the spherical relativistic Hartree theory Zhou et al. (2003b), spherical relativistic Hartree-Fock-Bogoliubov theory Long et al. (2010a, b), deformed relativistic Hartree-Bogoliubov theory in continuum (DRHBc) Zhou et al. (2010); Li et al. (2012), deformed relativistic Hartree-Fock theory Geng et al. (2020), and deformed relativistic Hartree-Fock-Bogoliubov theory Geng and Long (2022).

Recently, the DRHBc theory has been successful in studying deformed halo nuclei Zhou et al. (2010); Li et al. (2012); Sun et al. (2018); Zhang et al. (2019); Sun et al. (2020); Yang et al. (2021); Sun (2021); Zhong et al. (2022), investigating the deformation effects on the location of the neutron drip line In et al. (2021), predicting stability peninsulas beyond the primary neutron drip line Zhang et al. (2021); Pan et al. (2021); He et al. (2021), exploring rotational excitations of exotic nuclei with the angular momentum projection Sun and Zhou (2021a, b), and revealing the shape coexistence from light to heavy nuclei In et al. (2020); Choi et al. (2022); Kim et al. (2022). In particular, stimulated by the success of the fist nuclear mass table with continuum effects Xia et al. (2018), many efforts have been made to construct an upgraded mass table including simultaneously deformation and continuum effects based on the DRHBc theory Zhang et al. (2020), and the DRHBc mass table for even-even nuclei has just been published Zhang et al. (2022). The next step of the DRHBc mass table is the systematical calculation of odd-mass and odd-odd nuclei in the nuclear chart Pan et al. (2022). As mentioned in Ref. Zhang et al. (2020), since the DRHBc theory with the DWS basis is numerically very complicated, large-scale DRHBc calculations are extremely time consuming, especially for odd nuclei where the blocking effects should be considered.

This paper is devoted to exploring methods to reduce the computing cost of the DRHBc theory from the perspective of the DWS basis space. On the one hand, the DWS basis space includes both bases in the Fermi sea and in the Dirac sea for completeness Zhou et al. (2003b). In order to guarantee the convergence with respect to the basis space, it was suggested to take the number of DWS bases in the Dirac sea as the same as that in the Fermi sea Zhou et al. (2010); Li et al. (2012), and later all DRHBc studies followed this way. However, in the spherical relativistic Hartree-Fock-Bogoliubov Long et al. (2010a, b) and deformed relativistic Hartree-Fock Geng et al. (2020) calculations, the number of DWS bases in the Dirac sea is only about one third of that in the Fermi sea. On the other hand, usually the Woods-Saxon potential for the basis is parameterized Koepf and Ring (1991). It is natural to consider the possibility of an optimized DWS basis whose corresponding potential is closer to the mean field of the calculated nucleus, and as a result the basis space required for convergence would be reduced. Following this idea we examine in detail the convergence behavior against the basis truncation in the Dirac sea and propose an optimized Woods-Saxon basis for the large-scale DRHBc calculations as well as other DFT calculations in the present paper.

This paper is organized as follows: In Sec. II, we introduce the DRHBc theory and the DWS basis. The numerical details are given in Sec. III. The results and discussion are presented in Sec. IV. Finally, a summary is given in Sec. V.

The details of the DRHBc theory with meson-exchange and point-coupling density functionals can be found in Refs. Li et al. (2012) and Zhang et al. (2020) respectively. In the DRHBc theory, the relativistic Hartree-Bogoliubov (RHB) equation reads Kucharek and Ring (1991)

where $\lambda$ is the Fermi energy, and $E_{k}$ and $(U_{k},V_{k})^{\rm T}$ are the quasiparticle energy and wave function, respectively. $h_{D}$ is the Dirac Hamiltonian,

$\Delta$ is the pairing potential,

with a density-dependent force of zero range,

and the pairing tensor $\kappa(\bm{r}_{1},\bm{r}_{2})$ Ring and Schuck (1980). In principle, the zero-range pairing force is not naturally converged with the pairing window, and both appropriate pairing strength and pairing window are significant in the present DRHBc calculations Zhang et al. (2020). It is expected to implement the finite-range pairing force, e.g., the Gogny Meng (1998) or separable pairing force Tian et al. (2009), in the DRHBc theory in future work.

In the DRHBc theory, the scalar potential $S(\bm{r})$ and vector potential $V(\bm{r})$ in Eq. (2) are expanded in terms of the Legendre polynomials,

so are the pairing potential and various densities. In order to properly describe the large spatial extension of weakly bound nuclei, the RHB equation, (1), is solved in the DWS basis Zhou et al. (2003b). After solving the RHB equations self-consistently, the total energy, rms radii, quadrupole deformation, and other physical quantities can be calculated.

The wave function of the DWS basis can be written as

where $n$, $\kappa$, and $m$ are its quantum numbers, which will be introduced below. $\bm{r}$, $s$, and $p$ are spatial coordinate, spin, and index ($p=1$ or $2$) for upper or lower component. $R_{n\kappa}$ is the radial wave function,

satisfying the radial Dirac equations

where $\epsilon_{n}$ is the eigenvalue (energy), and $n$ is the node number. $V_{\mathrm{WS}}+S_{\mathrm{WS}}$ and $V_{\mathrm{WS}}-S_{\mathrm{WS}}$ are parameterized Woods-Saxon potentials for the basis,

where $V^{\pm}$, $R^{\pm}$, and $a^{\pm}$ depict the depth, width, and diffuseness of the potential, respectively. A usual choice of the parametrization can be found in Ref. Koepf and Ring (1991). $\mathcal{Y}_{\kappa m}^{l}$ is the spinor spherical harmonic,

where $Y_{lm_{l}}$ is the spherical harmonic function, $\chi_{\frac{1}{2}m_{s}}$ is the spin wave function, $l$ is the orbital angular momentum, $j$ is the total angular momentum, and $m_{l}$, $m_{s}$, and $m$ are the third components of $l$, $s$, and $j$, respectively. The quantum number $\kappa$ is given by the parity $\pi$ and $j$, $\kappa=\pi(-1)^{j+1/2}(j+1/2)$, running over positive and negative integers $\pm 1,\pm 2,\cdots$. The orbital angular momenta for upper and lower components are

respectively.

Taking light ²⁰Ne, medium-heavy ¹¹²Mo, and heavy ³⁰⁰Th nuclei as examples, a detailed convergence check of the total energy against the energy cutoff in the Dirac sea $E_{\mathrm{cut}}^{D}$ will be performed. The radial Dirac equations for the DWS basis, (8), are solved with a box size $R_{\mathrm{box}}=20$ fm and a mesh size $\Delta r=0.1$ fm Zhou et al. (2003b); Zhang et al. (2020). The angular momentum cutoff for the DWS basis is $J_{\mathrm{max}}=\frac{23}{2}~{}\hbar$ Zhang et al. (2020). The Legendre expansion truncations in Eq. (5) are chosen as $\lambda_{\mathrm{max}}=6$ and $8$ for nuclei with $8\leq Z\leq 70$ and $72\leq Z\leq 100$, respectively Pan et al. (2019); Zhang et al. (2020). The above numerical details are the same as those used in the DRHBc mass table calculations for even-even nuclei Zhang et al. (2022). Since a zero-range pairing interaction (4) is adopted in the DRHBc theory, in the numerical check of the convergence with respect to the basis space pairing correlations should be neglected. In the DRHBc mass table calculations Zhang et al. (2022), the energy cutoff for the DWS basis in the Fermi sea $E_{\mathrm{cut}}^{F}=300$ MeV and the number of DWS bases in the Dirac sea is the same as that in the Fermi sea, which are able to give converged results as shown in Ref. Zhang et al. (2020) and thus are taken as the benchmark in the following.

Figure 1 shows the convergence of the total energy with respect to the energy cutoff in the Dirac sea $E_{\mathrm{cut}}^{D}$ (whose zero point is set to be the continuum threshold in the Dirac sea) for ²⁰Ne, ¹¹²Mo, and ³⁰⁰Th. From panels (a), (b), and (c) one finds that, although the total energies of the three nuclei are very different, they show a similar convergence trend with the increase of $E_{\mathrm{cut}}^{D}$. When $E_{\mathrm{cut}}^{D}=50$ MeV, the difference of the calculated total energy from the converged one (horizontal dashed line) becomes marginal in their respective scale. It is clearly seen in panel (d) that their relative differences from the converged results are on the order of $10^{-5}$. Therefore, $E_{\mathrm{cut}}^{D}=50$ MeV for the DWS basis can lead to a satisfactory accuracy in the DRHBc calculation. The $E_{\mathrm{cut}}^{D}$ value greater than zero suggests that the contribution of the bases from continuum in the Dirac sea is non-negligible with the adopted DWS basis. A similar conclusion was found in a recent study based on the axially deformed relativistic Hartree-Fock-Bogoliubov theory with the density functional PKA1 Geng and Long (2022).

To compare the number of bases in the Dirac sea $N_{D}$ with $E_{\mathrm{cut}}^{D}=50$ MeV and that in the Fermi sea $N_{F}$ with $E_{\mathrm{cut}}^{F}=300$ MeV, Fig. 2 exhibits $N_{F}$ and the ratio $N_{D}/N_{F}$ as functions of the absolute value of the quantum number $\kappa$ in the calculations of ²⁰Ne, ¹¹²Mo, and ³⁰⁰Th. The $N_{F}$ as a function of $|\kappa|$ is shown in Fig. 2(a), where the basis numbers for neutrons/protons and positive/negative parities are close to each other and summed. First, we can see that the $N_{F}$s are not very different for ²⁰Ne, ¹¹²Mo, and ³⁰⁰Th because of the large energy cutoff $E_{\mathrm{cut}}^{F}=300$ MeV, and they decrease slowly and almost linearly with $|\kappa|$. Second, as shown in Figs. 2(b) and 2(c), all the ratios of $N_{D}/N_{F}$ are smaller than 1, meaning that the basis space in the Dirac sea is not necessarily as large as that in the Fermi sea in order to obtain converged results. Third, the ratio of $N_{D}/N_{F}$ shows an odd-even staggering with the increase of $|\kappa|$. This is because, for the same $|\kappa|$, the orbital angular momentum $l$ of the basis in the Fermi sea and that in the Dirac sea differs by $1$, e.g., for positive parity the ratio is calculated in the sequence of $s_{1/2}/p_{1/2}$, $d_{3/2}/p_{3/2}$, $d_{5/2}/f_{5/2}$, $\cdots$ for $|\kappa|=1,2,3,\cdots$. In a basis space truncated by energy, the number of bases corresponding to spin-orbit partners is in general the same, and it decreases with the increasing $l$. As a result, the numerator $N_{D}$ and denominator $N_{F}$ do not decrease synchronously, leading to the odd-even staggering. Fourth, the ratio decreases with $|\kappa|$ on the whole, which reflects the influence of the different depths of the potentials in the Fermi sea and Dirac sea. Finally, the ratios of $N_{D}/N_{F}$ for the heavy nucleus ³⁰⁰Th are larger than those for lighter ones ²⁰Ne and ¹¹²Mo, which shows that the bases in the Dirac sea are more important for heavy nuclei and is consistent with the finding in Ref. Geng and Long (2022). It is also noted that the ratios for protons are obviously larger than those for neutrons in the calculation of ³⁰⁰Th, which will be further discussed below.

Is it possible to optimize the usual DWS basis to achieve a better performance? We note that Ref. Zhou et al. (2003b) has found when the Woods-Saxon potentials for the basis become more different from the converged mean-field potentials of the calculated nucleus, the contribution of the basis in the Dirac sea will be larger. To illustrate this difference in potentials, Fig. 3 compares for ³⁰⁰Th the Woods-Saxon potentials for the basis and the converged mean-field ones, and for the latter the angular dependence has been averaged. For the $V+S$ potential of neutrons, the difference mainly appears in 6 $\leq r\leq$ 12 fm and is not remarkable. For the $V-S$ potential of neutrons, distinguishable difference can also be found in 0 $<r\leq$ 6 fm. For protons, however, one finds significant differences in both $V+S$ and $V-S$ potentials. The depth of the $V+S$ potential differs by about 10 MeV, and the difference in that of the $V-S$ potential is larger than 100 MeV. This explains why the ratios $N_{D}/N_{F}$ for protons are obviously larger than those for neutrons in Fig. 2(c).

According to the above discussion, a possible way to optimize the DWS basis and reduce the basis space is using the Woods-Saxon potentials close to the converged mean-field ones to generate a new basis replacing the original DWS basis. As an attempt, for deformed nuclei, we solve a Dirac equation containing the angle averaged converged mean-field potentials to obtain the new basis, referred to as the optimized Dirac Woods-Saxon (ODWS) basis. This optimized basis is in spirit somewhat similar to the method proposed in Refs. Zhang and Onley (1988, 1991), where deformed relativistic Hartree equations are solved by expanding the Dirac spinors in terms of basis functions calculated from the self-consistent spherical Hartree potential of the same or a nearby nucleus. At the end of this section, we will see the unique advantage making the ODWS basis superior to the recipe in Refs. Zhang and Onley (1988, 1991).

To test the performance of the ODWS basis with respect to the convergence behavior, Fig. 4(a) shows the difference of the total energy from the converged result as a function of $E_{\mathrm{cut}}^{F}$ in both calculations with the DWS and ODWS bases for ³⁰⁰Th. It is easy to see that although both calculations lead to convergence with increasing cutoff energies, the results obtained with the ODWS basis are always closer to the converged one. In particular, when $E_{\mathrm{cut}}^{F}=100$ MeV, the total energy calculated with the ODWS basis is $\geq 1$ MeV lower than that with the DWS basis. The results with the ODWS basis not only exhibit better convergence with the increasing $E_{\mathrm{cut}}^{F}$, but also hardly change with the increasing $E_{\mathrm{cut}}^{D}$. This suggests that we can obtain converged results even without the bases from continuum in the Dirac sea, which would significantly reduce the basis space.

To more finely evaluate the convergence with respect to $E_{\mathrm{cut}}^{D}$, in Fig. 4(b) we reduce the range of $E_{\mathrm{cut}}^{F}$ from [100, 300] MeV in Fig. 4(a) to [200, 300] MeV, and thus we rescale the vertical axis. It is found that for the usual DWS basis, the calculated total energy changes by about 0.5 MeV from $E_{\mathrm{cut}}^{D}=0$ to 50 MeV, suggesting the importance of the bases from continuum in the Dirac sea. In contrast, for the ODWS basis, the calculated total energy is almost a constant under different cutoff energies. For example, the difference between the results from $E_{\mathrm{cut}}^{F}=200$ MeV, $E_{\mathrm{cut}}^{D}=0$ MeV and $E_{\mathrm{cut}}^{F}=300$ MeV, $E_{\mathrm{cut}}^{D}=50$ MeV is only 0.005 MeV, about $10^{-6}$ of its total energy. Therefore, the calculation using the ODWS basis can provide converged total energies for heavy nuclei such as ³⁰⁰Th, with $E_{\mathrm{cut}}^{F}=200$ MeV and $E_{\mathrm{cut}}^{D}=0$ MeV that are both smaller than those required for the DWS basis. We have further reduced $E_{\mathrm{cut}}^{D}$ to negative values to see if a smaller ODWS basis space in the Dirac sea could give converged results. It turns out that $E_{\mathrm{cut}}^{D}\lesssim-100$ MeV leads to a ground-state deformation completely different from the converged one, and the difference between the total energy from $E_{\mathrm{cut}}^{D}\approx-50$ MeV and the converged one is more than 0.3 MeV. This indicates that some bound bases within 50 MeV from the continuum threshold in the Dirac sea are important for convergence. Thus, an energy cutoff of $E_{\mathrm{cut}}^{D}=0$ MeV, i.e., the truncation right at the continuum threshold in the Dirac sea, is suggested for the ODWS basis to guarantee convergence for different nuclei.

It is also interesting to check whether the ODWS basis space truncated by $E_{\mathrm{cut}}^{F}=200$ MeV and $E_{\mathrm{cut}}^{D}=0$ MeV can guarantee the convergence of single-particle properties. Figure 5 compares the single-neutron levels of ³⁰⁰Th near the continuum threshold calculated using DWS and ODWS bases as well as those from the converged result. Distinguishable differences can be seen between the result from the DWS basis and the other two. For example, the difference for the single-neutron energy of the $5/2^{+}$ state is nearly 0.1 MeV. On the other hand, the result from the ODWS basis and the converged one are almost the same, with the largest difference less than 0.01 MeV. Therefore, according to this comparison and the convergence of the total energy in Fig. 4, $E_{\mathrm{cut}}^{F}=200$ MeV and $E_{\mathrm{cut}}^{D}=0$ MeV for the ODWS basis are enough to give converged results in the DRHBc calculation. They correspond to a certainly smaller basis space than that of the DWS basis and, thus, require less computing resource.

Finally, one may argue that using the converged mean-field potentials to construct the ODWS basis makes little sense, since one never knows a priori the converged potentials of the calculated nucleus. This is not true. A natural way out is to construct a new basis set with new potentials in each step of the iterative solution of RHB equations. This has been done and can lead to the same results as those calculated by the ODWS basis, but it is found that the reconstruction of the basis in each step is time consuming. In practice, we find that during the iteration only one reconstruction of the basis is enough to get the same results. Here the reconstruction happens when a reference parameter $\sigma\approx 1$ MeV, in which $\sigma$ means the largest difference of the mean-field potentials between the present step and the previous step during the iteration. (When $\sigma$ is less than a certain value, e.g., $10^{-4}$ MeV, the self-consistent iteration converges.) Note that there is no need to change the cutoff energies during the iteration in the above approach; namely, $E_{\mathrm{cut}}^{F}=200$ MeV and $E_{\mathrm{cut}}^{D}=0$ MeV are used for the original DWS basis from the beginning and for the ODWS basis after $\sigma\approx 1$ MeV, up to convergence. It is also worthwhile to mention that, although this approach replaces bases during the iteration, there is almost no difference in the convergence steps compared with the calculation using only the DWS basis. Therefore, the ODWS basis is practicable, and its application to large-scale DRHBc mass table calculations is expected.

In summary, based on the deformed relativistic Hartree-Bogoliubov theory in continuum, we check in detail the convergence with respect to the basis space in the Dirac sea, and propose an optimized Woods-Saxon basis for nuclear density functional theory. By checking the convergence with respect to the energy cutoff, we demonstrate that the basis space in the Dirac sea is not necessarily as large as that in the Fermi sea in order to obtain results with a satisfactory convergence accuracy of $10^{-5}$. Then we use the converged mean-field potential to generate a new basis set, referred to as the optimized Dirac Woods-Saxon basis. It is found that the new basis can lead to convergence in both total energy and single-particle levels with a basis space truncated by the energy cutoff $E_{\mathrm{cut}}^{F}=200$ MeV in the Fermi sea and $E_{\mathrm{cut}}^{D}=0$ MeV in the Dirac sea, which is substantially smaller than the original basis space required for convergence. Finally, we suggest a method to construct the optimized Dirac Woods-Saxon basis during the iterative solution of relativistic Hartree-Bogoliubov equations, which guarantees the practicability of the new basis with a suitable basis space.

The application of the optimized Dirac Woods-Saxon basis to large-scale mass-table-type calculations based on the covariant density functional theory is expected to significantly reduce computing resource. Similarly, its counterpart, the optimized Schrödinger Woods-Saxon basis, is also available for nonrelativistic density functional theory.