Solving the relativistic Hartree-Bogoliubov equation with the finite-difference method

The neutron- or proton-rich nuclei far away from the $ \beta- $stability line found by the newly developed radioactive ion beam facilities and detection techniques show exotic characteristics, challenging the available nuclear many-body theories [1]. Given that the Fermi energies of these nuclei are close to the continuum threshold, the valence neutrons or protons could be scattered into continuum states owing to the pairing correlation [2, 3]. A powerful tool to describe such exotic nuclei is the Hartree-Fock-Bogoliubov (HFB) theory, which transforms pair-correlated neutrons and protons to independent quasi-particles (q.p.) [4]. Using the energy-density functional theory, the Hartree-Fock mean field and pairing field can be solved self-consistently by the HFB equation. The two-component q.p. wave functions solved by the HFB equation can properly describe the nucleon densities including the continuum states. One can also transform the q.p. states into canonical states, which helps to understand the single-particle (s.p.) levels with influence of the pairing correlation. Within both the non-relativistic and relativistic frameworks, this theory has been widely applied and succeeded in describing exotic nuclei [5–8].

The key technique in the HFB theory is to solve the HFB equation. It is preferable to solve this equation in the coordinate space because this directly provides the q.p. wave functions and thus the nucleon densities in the coordinate space. To this end, one can usually employ the shooting or Numerov [9–17], finite-element [17, 18], finite-difference [19], Green's function [20–23], Jost function [24], and Lagrange-mesh [25, 26] methods to solve the HFB equation in the coordinate space. Given that exotic nuclei usually have an extensive density distribution, one has to work with a large coordinate space. Some proper basis can also be used to solve the HFB equation [27–34]. In the relativistic Hartree(-Fock)-Bogoliubov (RH(F)B) theory, the Dirac Woods-Saxon basis is also helpful to improve the description of the extensive density distribution in the coordinate space [35–41]. In particular, the Dirac Woods-Saxon basis has been applied in the deformed relativistic Hartree-Bogoliubov theory in the continuum (DRHBc). With this theory, the large-scale calculation of the nuclear mass table is currently proceeding and has already obtained fruitful results [42–45].

In the coordinate space, the finite-difference method (FDM) is a simple and efficient method to solve differential equations. However, when this method is applied to solve the Dirac equation, one will inevitably encounter the problem of spurious states with the commonly used symmetric central difference formula (CDF) [46]. This is also known as the fermion-doubling problem in the lattice quantum chromodynamics (QCD) theory [47, 48]. Recently, this problem was resolved by using the asymmetric difference formula (ADF) [49, 50]. In this paper, we use the FDM to solve the relativistic Hartree-Bogoliubov (RHB) equation with given mean field and pairing field. The problem of spurious states will be described and resolved by the same technique suggested by Ref. [50].

This paper is organized as follows. In Sec. II, we briefly review the RHB equation, and the FDM to solve this equation with the CDF and ADF. The details of the mean field and pairing field are also provided in this section. In Sec. III, we discuss the calculation results correspondingly. In particular, we show the results calculated by assuming a constant pairing strength in Subsec. III.A, and those calculated by assuming a Gaussian pairing strength in Subsec. III.B. A summary is provided in Sec. IV.

In the Bogoliubov theory, a pair-correlated nuclear system is described in terms of independent quasi-particles. The RHB equation in the coordinate space is [9, 10]

where h is the s.p. Hamiltonian, Δ is the pairing potential, λ is the Fermi energy, E is the q.p. energy, and $ \psi_U $ and $ \psi_V $ are the upper and lower components of the corresponding q.p. wave functions. Note that, for a bound system ($ \lambda<0 $), according to the asymptotic behavior of the wave functions $ \psi_U $ and $ \psi_V $ in Eq. (1), the q.p. states with energy $ 0<E<|\lambda| $ are defined as discrete q.p. states, and those with energy $ E>|\lambda| $ are continuum q.p. states [4, 9].

In the RHB theory, the Dirac Hamiltonian can be expressed as

where $ V({\bf{r}}) $ and $ S({\bf{r}}) $ are the vector and scalar potentials, respectively. For simplicity, we also consider a δ-function-type pairing potential Δ as

In the following, we assume spherical symmetry for the nuclear system. Then, the Dirac spinor wave functions $ \psi_U({\bf{r}}) $ and $ \psi_V({\bf{r}}) $ can be expressed as

where $ l+l'=2j $ and $ l=j\pm \frac{1}{2} $, and $ \Omega_{ljm} $ is the spherical spinor. Then, the RHB equation can be reduced to four coupled equations for the radial wave functions, namely $ G_U $, $ F_U $ and $ G_V $, $ F_V $ [17], as

where the mean-field potentials are usually shifted by M for comparison with non-relativistic results, and $\kappa=(-1)^{j+l+{1}/{2}}(j+{1}/{2})$. In this paper, to test the validity of the FDM to solve the RHB equation, we consider the Woods-Saxon-type potential for $ V(r) $ and $ S(r) $,

where the parameters $ V^+_0=-78.43 $ MeV, $ r^+_0=4.049 $ fm, $ a^+=0.9254 $ fm, and $ V^-_0=787.70 $ MeV, $ r^-_0=4.048 $ fm, $ a^-=0.6523 $ fm are fitted to the neutron mean-field potential obtained by the self-consistent DRHBc calculation for ⁵⁶Ca [42]. The Fermi energy in Eq. (5) takes the value $ \lambda=-4.00 $ MeV, read from the self-consistent result for neutrons in ⁵⁶Ca. The pairing strength $ \Delta_0(r) $ is studied according to the two following cases:

where $ d_0=2 $ MeV, $ x_0=5 $, and $ r_0=1.2A^{1/3} $ fm, and A is the total mass number of the nucleus. In the following, we proceed with a one-step calculation to solve the RHB equation (Eq. (5)) with the above model potentials by using the FDM and shooting method [17] for comparison. Once the validity of the FDM is demonstrated, it can be used in the self-consistent RHB calculation and extended to study realistic drip-line nucleus in the near future.

In the FDM, the first-order derivative of a given function $ f(r) $ can be approximated by a numerical differential formula on the equal interval lattices. The commonly used symmetric central difference formulas (CDFs) are the 3-point (3PCDF) and 5-point (5PCDF) formulas,

where h is the lattice interval. Besides, one can also use an asymmetric difference formula (ADF), such as the forward or backward 3-point ADF (3PADF),

and those of the 5-point ADF (5PADF),

With these differential formulas, the RHB equation (Eq. (5)) can be expressed as a matrix in the coordinate space,

where $ G_U $ ($ F_U $), and $ G_V $ ($ F_V $) are vectors of large (small) components of q.p. wave functions $ \psi_U $ and $ \psi_V $, respectively. The matrix elements in A, $ B_1 $, $ B_2 $, and C are almost the same as those resulting from solving the Dirac equation in Ref. [50], except that the Fermi energy $ -\lambda $ is added in the diagonal elements of A and C. Using the δ-function-type pairing potential in Eq. (3), the matrix D is also diagonal with the elements $ \Delta_0 $. Once the matrix $ H_{\rm{qp}} $ in Eq. (15) is diagonalized, one can obtain the eigenvalues of the q.p. energies E and the corresponding wave functions $ G_U $, $ F_U $, $ G_V $, and $ F_V $ in the coordinate space. In the following calculations, we set the box sizes in the coordinate space as $ R_{\rm{box}}=20 $ fm and $ R_{\rm{box}}=40 $ fm with lattice interval ${\rm{d}}r=h=0.1$ fm. This lattice interval is checked to provide precise enough results for the following discussion. Furthermore, we assume the same boundary condition for the wave functions as that in Ref. [50], i.e., $ f(r)=0 $ for $ r=0 $ and outside the box $ r>R_{\rm{box}} $. Here, $ f(r) $ denotes the radial wave functions $G_U\; (F_U)$ and $G_V \; (F_V)$ in Eqs. (4) and (5).

It is evident that by using the different differential formula in Eqs. (9)−(14), the matrix form of $ H_{\rm{qp}} $ will be different. According to Refs. [46, 49, 50], by using CDFs such as Eqs. (9) and (10), one will inevitably encounter the problem of spurious states when solving the Dirac equation, given that the information of the wave functions at the middle point $ f(r) $ is missing when calculating the first-order derivatives. This problem can be completely resolved by using ADFs such as Eqs. (11)−(14). However, note that the same ADF for both the large $ G(r) $ and small $ F(r) $ components of Dirac wave functions may lead to a non-Hermitian matrix of the Dirac Hamiltonian. To solve this problem, Ref. [50] pointed out that one should use the forward or backward ADF alternatively for $ G(r) $ and $ F(r) $ according to their different parities. This prescription can not only guarantee the Hermiticity of the Dirac Hamiltonian, but also include the full wave function information while doing the first-order derivatives, thereby eliminating the spurious states completely.

In the following, we use the CDFs in Eqs. (9)−(10) and ADFs in Eqs. (11)−(14) to establish the q.p. Hamiltonian $ H_{\rm{qp}} $ and obtain the q.p. eigenstates. First, we check whether spurious states appear in the q.p. states and whether the prescription pointed in Ref. [50] is helpful to resolve the spurious q.p. state problem.

For the pairing field potential with constant pairing strength $ \Delta_0 $, it can be proved that the solution of the q.p. energy for the RHB equation (Eq. (5)) can be expressed as [4]

with the corresponding wave functions

where the s.p. energy ε, and the corresponding wave functions G and F, are the solutions of the Dirac equation with the same mean-field potentials V and S,

The coefficients $ c_1 $ and $ c_2 $ can be calculated as

Therefore, we can use the above solutions as a benchmark to check the validity of the FDM to solve the RHB equation.

The results for the q.p. energy E, and the corresponding occupation probabilities [17]

of the neutrons for $ \kappa=-1 $ ($ s_{1/2} $) and $ \kappa=+1 $ ($ p_{1/2} $) in ⁵⁶Ca, calculated by the FDM using the 3PCDF, 5PCDF, 3PADF, and 5PADF to establish the q.p. Hamiltonian $ H_{\rm{qp}} $ in the box $ R_{\rm{box}}=20 $ fm with constant pairing strength are listed in Table 1. For comparison, the results calculated by the shooting method [17] are also listed in this table.

3PCDF			5PCDF			3PADF			5PADF			Shooting
$ \kappa=-1 $	$ \kappa=1 $		$ \kappa=-1 $	$ \kappa=1 $		$ \kappa=-1 $	$ \kappa=1 $		$ \kappa=-1 $	$ \kappa=1 $		$ \kappa=-1 $	$ \kappa=1 $
49.481 (0.9996)	$\fbox{\begin{array}{c}49.481\\ ( 0.9996 )\end{array}}$		${\begin{array}{c}49.466 \\(0.9996)\end{array}}$	$\fbox{\begin{array}{c}44.761 \\( 0.9995 )\end{array}}$		$\begin{array}{c}49.436\\(0.9996)\end{array}$			$\begin{array}{c}49.466\\ (0.9996)\end{array}$			$\begin{array}{c}49.466 \\(0.9996)\end{array}$
$\fbox{\begin{array}{c}30.446\\ (0.9989)\end{array}}$	$\begin{array}{c}30.446 \\(0.9989)\end{array}$			$\begin{array}{c}30.401 \\(0.9989)\end{array}$			$\begin{array}{c}30.309\\ (0.9989)\end{array}$			$\begin{array}{c}30.400\\ (0.9989)\end{array}$			$\begin{array}{c}30.400\\ (0.9989)\end{array}$
$\begin{array}{c}14.968 \\(0.9955)\end{array}$	$\fbox{\begin{array}{c}14.968 \\( 0.9955)\end{array}}$		$\begin{array}{c}14.896\\ (0.9955)\end{array}$			$\begin{array}{c}14.732\\ (0.9954)\end{array}$			$\begin{array}{c}14.892 \\(0.9955)\end{array}$			$\begin{array}{c}14.890 \\(0.9954)\end{array}$
$\fbox{\begin{array}{c}2.343 \\( 0.7606 )\end{array}}$	$\begin{array}{c}2.343 \\(0.7606)\end{array}$		$\fbox{\begin{array}{c}14.645 \\( 0.9953)\end{array}}$	$\begin{array}{c}2.302 \\(0.7475)\end{array}$			$\begin{array}{c}2.226 \\(0.7197)\end{array}$			$\begin{array}{c}2.302 \\(0.7476)\end{array}$			$\begin{array}{c}2.301 \\(0.7473)\end{array}$

Table 1.  Values of q.p. energy E and corresponding occupation probability $ v^2 $ in the brackets below for the neutrons in ⁵⁶Ca with $ \kappa=-1 $ ($ s_{1/2} $) and $ \kappa=+1 $ ($ p_{1/2} $) calculated with constant pairing strength by using the FDM with different differential formulas and the shooting method in the box $ R_{\rm{box}}=20 $ fm. The spurious q.p. states are marked by boxes. The unit of the energy values is MeV.

Taking the results obtained by 3PCDF in Table 1 as examples, one can clearly see that the energies and occupation probabilities of the q.p. states for $ \kappa=-1 $ and $ \kappa=+1 $ are exactly the same. To find out the reason, we present their corresponding wave functions in Fig. 1. Given that the $G_U\; (F_U)$ and $G_V\; (F_V)$ components of the wave functions differ by factors $ c_1 $ and $ c_2 $, respectively, as expected from Eqs. (17) and (18), here we only show the $G_V \;(F_V)$ components of the wave functions. Figs. 1 (a) and (b) show the wave functions of the first q.p. state solutions for $ \kappa=-1 $ and $ \kappa=+1 $, respectively, with the same q.p. energy $ E=49.481 $ MeV. Note that the wave functions in Fig. 1 (a) with $ \kappa=-1 $ are normal, which proves that this state is the physical q.p. state. The one peak structure shows that this q.p. state corresponds to the $ 1s_{1/2} $ s.p. state owing to the Bogoliubov transformation. However, in Fig. 1 (b), for $\kappa=+1$, the rapidly oscillating wave functions on neighboring lattices definitely show that this is a spurious q.p. state. Note also that half of the envelope of these oscillating wave functions is identical to that of the physical state. Similarly, we can find the wave functions of physical q.p. states corresponding to the s.p. states $ 1p_{1/2} $, $ 2s_{1/2} $, and $ 2p_{1/2} $ in Figs. 1 (d), (e) and (h), and the spurious states in Fig. 1 (c), (f), and (g), respectively. All the spurious states are marked in boxes in Table 1.

Figure 1.  (color online) $ G_V $ and $ F_V $ components of q.p. wave functions of neutrons in ⁵⁶Ca for the states with $ \kappa=-1 $ in (a) (c) (e) (g) and $ \kappa=1 $ in (b) (d) (f) (h), obtained by solving the RHB equation with Woods-Saxon-type mean-field and constant-pairing-strength potentials and using the FDM with the 3PCDF. The unit of the wave functions is fm^-1/2.

Figure 1 shows that the physical and spurious states with degenerate energies appear alternatively in $ \kappa=-1 $ and $ \kappa=+1 $, solved by the FDM with the 3PCDF for the RHB equation. This is similar to the the solutions of the s.p. states in the Dirac equation solved with the same methods [50]. Note also that the 3PCDF in Eq. (9) includes the information of the wave functions at $ f(r-h) $ and $ f(r+h) $, but not at the middle point $ f(r) $. This leads to a unitary transformation matrix with alternative $ \pm 1 $ diagonal elements between the Dirac Hamiltonian $ H_{\kappa} $ and $ H_{-\kappa} $, which produces energy-degenerate physical and spurious states. It can be demonstrated that such unitary transformation still exists between the q.p. Hamiltonian $ H_{{\rm{qp}}, \kappa} $ and $ H_{{\rm{qp}}, -\kappa} $. According to the benchmark established by Eqs. (16)-(18), the spurious s.p. states in the Dirac equation are transformed into spurious q.p. states solved by the same methods for the RHB equation. The number of spurious q.p. states can be reduced by using the 5PCDF, as shown in Table 1. The energy-degenerate physical and spurious q.p. states cannot be found because the unitary transformation matrix mentioned above disappears when the Hamiltonian matrix is established using the 5PCDF. However, the spurious states still exist, given that the information of the wave functions at the middle point to calculate the first-order derivative is still missing.

To remove the problem of spurious states completely in the Dirac equation, Ref. [50] proposed to use the ADF by including critical information of the wave functions at the middle point to construct the Dirac Hamiltonian. Following this prescription, we used the 3PADF and 5PADF to construct the q.p. Hamiltonian $ H_{\rm{qp}} $. The results are listed in Table 1. Note that all the spurious q.p. states disappear. In particular, the 5PADF produces results closer to those given by the shooting method. This demonstrates that the 5PADF can provide results with higher precision than those from the 3PADF. Furthermore, we checked that all the physical q.p. states obtained by the FDM are identical to those from the benchmarks established by Eqs. (16)−(18), transformed from the s.p. states obtained by the same FDM.

In the previous subsection, the case of constant pairing strength is simple because the RHB equation can be decoupled into the Dirac equation. In this section, we discuss the case of Gaussian pairing strength (Eq. (8)) with the peak near the nucleus surface. In this case, we checked that the FDM with the CDF also lead to spurious q.p. states. In the following, to remove these spurious states, we use the 5PADF in the FDM with higher precision and compare with the results obtained by the shooting method.

In Table 2, we list the values of q.p. energies E with the corresponding occupation probabilities $ v^2 $ in brackets for neutrons with $ \kappa=-1 $ and $ \kappa=+1 $ in ⁵⁶Ca obtained by the FDM with the 5PADF and the shooting method in the boxes $ R_{\rm{box}}=20 $ fm and $ R_{\rm{box}}=40 $ fm. First, for the states with $ \kappa=-1 $, one can find that all the results for the largest q.p. energy of $ 49.435 $ MeV provided by the different methods in box sizes $R_{\rm{box}}=20,\; 40$ fm are exactly the same. Besides, the occupation probability of this state is almost one. It can be inferred from the approximate relation between q.p. and s.p. energies (Eq. (16)) that this state with the largest q.p. energy should correspond to the most deeply bound s.p. state $ 1s_{1/2} $.

$ \kappa=-1 $
5PADF (20 fm)			Shooting (20 fm)			5PADF (40 fm)			Shooting (40 fm)
E	εcan		E	εcan		E	εcan		E	εcan
49.435 (0.9998)	−52.086 (1.0000)		49.435 (0.9998)	−52.023 (1.0000)		49.435 (0.9997)	−52.012 (1.0000)		49.435 (0.9997)	−51.944 (1.0000)
14.945 (0.6441)	−20.097 (0.9980)		14.972 (0.5618)	−20.048 (0.9980)		14.791 (0.9699)	−20.169 (0.9980)		14.790 (0.9708)	−20.110 (0.9980)
14.539 (0.3517)			14.579 (0.4339)			15.797 (0.0188)			15.829 (0.0177)
$ \kappa=1 $
5PADF (20 fm)			Shooting (20 fm)			5PADF (40 fm)			Shooting (40 fm)
E	εcan		E	εcan		E	εcan		E	εcan
30.382 (0.9856)	−34.309 (0.9997)		30.363 (0.9961)	−34.254 (0.9997)		30.370 (0.9920)	−34.309 (0.9997)		30.348 (0.9610)	−34.246 (0.9997)
29.054 (0.0131)			33.328 (0.0021)			28.887 (0.0049)			30.868 (0.0363)
1.650 (0.8374)	−5.162 (0.8406)		1.649 (0.8372)	−5.106 (0.8403)		1.650 (0.8375)	−5.161 (0.8406)		1.649 (0.8372)	−5.097 (0.8403)

Table 2.  Values of q.p. energies E and canonical energies $ \varepsilon_{\rm{can}} $ with the corresponding occupation probabilities $ v^2 $ in the brackets for neutrons in ⁵⁶Ca with $ \kappa=-1 $ and $ \kappa=+1 $ calculated with Gaussian pairing strength by using the FDM with the 5PADF and the shooting method in the boxes $ R_{\rm{box}}=20 $ fm and $ R_{\rm{box}}=40 $ fm. The unit of the energy values is MeV.

Taking the FDM results for $ \kappa=-1 $ calculated in the box $ R_{\rm{box}}=20 $ fm as examples, one can find two nearly degenerate states with q.p. energies of $ 14.945 $ MeV ($ v^2=0.6441 $) and $ 14.539 $ MeV ($ v^2=0.3517 $) and comparable occupation probabilities $ v^2 $. Similar nearly degenerate states at $ 14.972 $ MeV ($ v^2=0.5618 $) and $ 14.579 $ MeV ($ v^2=0.4339 $) are also obtained by the shooting method. Although the q.p. energies obtained by both methods are close to each other, their occupation probabilities $ v^2 $ are clearly different. Therefore, we plot the wave functions of $ G_U $, $ F_U $, $ G_V $, and $ F_V $ in Figs. 2 (a) and (c) for the state of $ 14.945 $ MeV ($ v^2=0.6441 $) obtained by the FDM and $ 14.972 $ MeV ($ v^2=0.5618 $) obtained by the shooting method in the box $ R_{\rm{box}}=20 $ fm. The wave functions of the close states with q.p. energies of $ 14.539 $ MeV obtained by the FDM and $ 14.579 $ MeV by the shooting method are similar to those in Figs. 2 (a) and (c), respectively; they are not shown for simplicity. As mentioned before, given that this state has a q.p. energy of $ E\approx 15 $ MeV, which is larger than the Fermi energy $ |\lambda|=4.00 $ MeV, it should be a continuum q.p. state. Their wave functions $ G_U $ and $ F_U $ must be oscillating, and $ G_V $ and $ F_V $ must have exponential decaying asymptotic behaviors, as clearly shown in Figs. 2 (a) and (c). Besides, the large component of the wave functions $ G_V $ with one node shows that this q.p. state must correspond to the $ 2s_{1/2} $ s.p. state. Note that although the q.p. energies of this state obtained by the FDM and shooting method are close, their wave functions are clearly different. Note also that within the box $ R_{\rm{box}}=20 $ fm, a smaller lattice interval, such as ${\rm{d}}r=0.05$ fm, can provide more consistent wave functions for similar q.p. energy states resulting from the 5PADF and the shooting method. Here, we show the results for ${\rm{d}}r=0.1$ fm as an example. Additionally, in the inset of Fig. 2 (e), we enlarge the component of $ G_V $ obtained by both methods within the same region, $ r=15 $ to $ 20 $ fm, in a logarithmic scale. Their asymptotic wave functions are exactly the same, forced to be zero at the boundary of the box. Despite this, other parts of their wave functions notably differ, as shown in Fig. 2 (c). This is owing to the different components $ G_U $, given that the four components together should be normalized to one.

Figure 2.  (color online) Representation of the q.p. wave functions for the components $ G_U $, $ F_U $ in (a) (b), and $ G_V $, $ F_V $ in (c) (d) of neutrons in ⁵⁶Ca for the $ 2s_{1/2} $ state calculated with Gaussian pairing strength by the FDM with the 5PADF (solid and dash-dotted lines) and the shooting method (dashed and dotted lines) in the boxes (a) (c) $ R_{\rm{box}}=20 $ fm and (b) (d) $ R_{\rm{box}}=40 $ fm. The insets (e) and (f) show the asymptotic wave functions of $ G_V $ up to $ r=20 $ fm. The unit of the wave functions is fm^-1/2.

When we calculate with a larger box size $ R_{\rm{box}}=40 $ fm, the nearly degenerate q.p. states with comparable occupation probabilities seem to disappear. Explicitly, one can obtain the states with $ 14.791 $ MeV ($ v^2=0.9699 $) and $ 15.797 $ MeV ($ v^2=0.0188 $) from the FDM, and those with $ 14.790 $ MeV ($ v^2=0.9708) $ and $ 15.829 $ MeV ($ v^2=0.0177 $) from the shooting method. In particular, the states with relatively larger occupation probabilities, namely $ 14.791 $ MeV ($ v^2=0.9699 $) obtained from the FDM and $ 14.790 $ MeV ($ v^2=0.9708) $ obtained from the shooting method, are almost the same for both q.p. energies and occupation probabilities. In Figs. 2 (b) and (d), we plot their wave functions. Note that all the components of their wave functions are almost the same. The asymptotic wave functions of $ G_V $ are also shown in a logarithmic scale in the inset of Fig. 2 (f). Their asymptotic behaviors are reasonable at $ r=20 $ fm. However, near $ r=40 $ fm, their asymptotic behaviors are similar to those near $ r=20 $ fm calculated in the box $ R_{\rm{box}}=20 $ fm in Fig. 2 (e), with much smaller wave function values. By comparing the results of the wave functions $ G_U $ obtained within $ R_{\rm{box}}=20 $ fm and $ R_{\rm{box}}=40 $ fm in Figs. 2 (a) and (b), one can clearly see that they significantly differ, in particular at $ r\approx 5 $ fm and $ r\approx 20 $ fm. This demonstrates that a box size of $ R_{\rm{box}}=20 $ fm is not large enough to describe these oscillating wave functions. Therefore, a larger box size is necessary to provide reliable wave functions for continuum q.p. states, especially for the $ \psi_U $ components.

Using the model potentials in Eqs. (6)−(8), we can only find continuum q.p. states with $ \kappa=-1 $. In Table 2, we also list the results for the states with $ \kappa=+1 $. Here, one can find a discrete q.p. state with an energy of $ E\approx 1.65 $ MeV $ <|\lambda|=4.00 $ MeV, and occupation probability of $ v^2\approx 0.837 $. For this state, all the calculations from the FDM and shooting method with boxes $ R_{\rm{box}}=20,\; 40 $ fm produce almost the same results. We plot their wave functions in Fig. 3. Note from Figs. 3 (a) and (c) that the four components $ G_U $, $ F_U $, $ G_V $, and $ F_V $ calculated by different methods and different box sizes are exactly the same. This is because all the components have exponential decaying asymptotic wave functions, and a box size of $ R_{\rm{box}}=20 $ fm is large enough to describe these wave functions. From the nodes of the $ G_U $ and $ G_V $ components, we know that this q.p state should correspond to the $ 2p_{1/2} $ state near the Fermi energy. In the inset of Figs. 3 (b) and (d), we enlarge the asymptotic wave functions from $ r=15 $ to $ 20 $ fm for the components $ G_U $ and $ G_V $, respectively. It is interesting to note that, with a smaller box size $ R_{\rm{box}}=20 $ fm, only the shooting method will force these wave functions to be zero, while the FDM will not. This seems to be different from the results for the $ 2s_{1/2} $ state in Fig. 2 (e), where the asymptotic wave functions of $ G_V $ components given by both the FDM and shooting method quickly drop to zero near $ R_{\rm{box}}=20 $ fm. Note that the shooting method uses the same recursion formula for both the even- and odd-parity wave functions. Meanwhile, the FDM uses the forward (backward) differential formula in Eq. (13) (Eq. (14)) for the odd(even)-parity wave functions to simultaneously guarantee the boundary condition at the origin and the Hermiticity of the Hamiltonian matrix [50]. The G component of the $ s_{1/2} $ state has odd parity, so the forward asymmetric difference formula at the box boundary $ R_{\rm{box}} $ makes the wave functions quickly drop to zero, given that the boundary condition requires the wave functions at $ r\geq R_{\rm{box}} $ to be zero. However, for the G component of the $ 2p_{1/2} $ state, the backward differential formula can provide the smoothly changing wave functions before the box boundary. This is also reported in Ref. [50]. When a box size is $ R_{\rm{box}}=40 $ fm, the shooting method also makes the wave functions $ G_U $ and $ G_V $ quickly drop to zero near r=40 fm; however, for this state, the FDM does not.

Figure 3.  (color online) Representation of the q.p. wave functions for the components $ G_U $, $ F_U $ in (a), and $ G_V $, $ F_V $ in (c) of neutrons in ⁵⁶Ca for the $ 2p_{1/2} $ state calculated with Gaussian pairing strength by the FDM with the 5PADF (solid and dash-dotted lines) and the shooting method (dashed and dotted lines) in the boxes $ R_{\rm{box}}=20 $ fm and $ R_{\rm{box}}=40 $ fm. The insets (b) and (d) show the asymptotic wave functions of $ G_U $ and $ G_V $ up to $ r=20 $ fm. The unit of the wave functions is fm^-1/2.

In Table 2, for $ \kappa=+1 $, in addition to the discrete q.p. state with energy $ E\approx 1.65 $ MeV, there are also several states with energy $ E\approx 30 $ MeV and a large occupation probability $ v^2>0.9 $ resulting from both the FDM and shooting method with $R_{\rm{box}}=20,\; 40$ fm. According to the aforementioned analysis, we know that these continuum q.p. states should correspond to the most deeply bound $ 1p_{1/2} $ state. Besides, there are some q.p. states nearby with similar energies and small but non-negligible occupation probabilities. This is similar to the results with $ E\approx 15 $ MeV and non-negligible occupation probabilities $v^2\approx 0.02$ for $ \kappa=-1 $ calculated with $ R_{\rm{box}}=40 $ fm. The wave functions of these 'split' q.p. states also have oscillating $ \psi_U $ and local $ \psi_V $ components similar to those nearby states with $v^2 > 0.9$. To understand these states, we next analyze the canonical states.

The canonical states can be obtained by diagonalizing the density matrix constructed by the q.p. states with the same κ [17]. The canonical s.p. energies $ \varepsilon_{\rm{can}} $ and corresponding occupation probabilities $ v^2 $ for the two most deeply bound states are listed in Table 2. For the states with $ \kappa=-1 $, one can see that the results of both the canonical s.p. energies and occupation probabilities obtained by different methods within different boxes are almost the same. Note that the really small differences come from the different q.p. states used to diagonalize the density matrix. Taking the calculations with the box $ R_{\rm{box}}=20 $ fm as examples, it seems from the occupation probabilities that the canonical s.p. state $ 2s_{1/2} $ 'splits' into two nearly degenerate q.p. states with $ E\approx 15 $ MeV and comparable $ v^2 $. The similar 'splitting' can be also found in the results calculated with a larger box $ R_{\rm{box}}=40 $ fm, although one q.p. state $ E\approx 14.79 $ MeV has a much larger $ v^2\approx 0.97 $, whereas the other one nearby has a much smaller $v^2\approx 0.02$. For $ \kappa=+1 $, one can also find this similar 'splitting' in the q.p. states at $ E\approx 30 $ MeV. Therefore, we know that this 'splitting' is common for the continuum q.p. states. Although these q.p. states correspond to deeply bound canonical states, they are transformed to be q.p. states in the continuum, which are resonant q.p. states with certain width owing to the pairing correlation. With the box boundary condition considered in this study, the continuum q.p. states are discretized. Then, one can obtain the 'split' q.p. states around the resonant q.p. states. To obtain the exact resonant q.p. energy and width, other techniques should be used to properly consider the continuum asymptotic boundary condition, such as the Green's function method [20–23].

Figures 4 (a) and (c) show the wave functions G and F for the canonical states $ 2s_{1/2} $ and $ 2p_{1/2} $ obtained from the FDM and shooting method in the boxes $ R_{\rm{box}}=20 $ fm and $ R_{\rm{box}}=40 $ fm. One can see that the results calculated by different methods with different boxes are almost the same. In the insets of Figs. 4 (b) and (d), the asymptotic behaviors of the large component G obtained by different calculations are shown in a logarithmic scale. For $ 2s_{1/2} $, the results obtained by both the FDM and the shooting method with $ R_{\rm{box}}=20 $ fm are forced to be zero near the box boundary. The results calculated with $ R_{\rm{box}}=40 $ fm have better asymptotic behaviors near $ r=20 $ fm. However, near $ r=40 $ fm, their asymptotic behaviors are similar to those near $ r=20 $ fm calculated in the box $ R_{\rm{box}}=20 $ fm, with much smaller wave function values. For $ 2p_{1/2} $, the asymptotic wave functions provided by the FDM seem better than those provided by the shooting method. The reason is the same as that analyzed before for Fig. 3. By comparing with the q.p. wave functions obtained by different methods and different boxes in Figs. 2 and 3, it can be concluded that, although with a not-large-enough box, the results of continuum q.p. states may not be reliable, the canonical states obtained by superposition of these q.p. states are robust except for the asymptotic wave functions.

Figure 4.  (color online) Canonical s.p. wave functions of the components $ G(r) $ and $ F(r) $ for (a) $ 2s_{1/2} $ and (b) $ 2p_{1/2} $ states calculated with Gaussian pairing strength by the FDM with the 5PADF (solid and dash-dotted lines) and the shooting method (dashed and dotted lines) in the boxes $ R_{\rm{box}}=20 $ fm and $ R_{\rm{box}}=40 $ fm. The insets (b) and (d) show the asymptotic wave functions of $ G(r) $ up to $ r=20 $ fm. The unit of the wave functions is fm^-1/2.

In this paper, we solved the RHB equation with the FDM. To check the feasibility, we compared with the results obtained by the shooting method. The mean-field potential was set as Woods-Saxon type fitted to the self-consistent calculations for ⁵⁶Ca. The pairing potential was taken as a delta-type function with constant and Gaussian pairing strengths. For the constant pairing strength case, the RHB equation decouples into the Dirac equation. The q.p. state solutions of the RHB equation exhibit a simple relation with the s.p. state solutions of the Dirac equation. We show that spurious q.p. states are obtained by the FDM using the symmetric CDF to construct the q.p. Hamiltonian matrix. This problem is solved by using the ADF in the FDM. The obtained q.p. state solutions are checked by the benchmarks related with the s.p. state solutions. Then, we solved the RHB equation in the case of Gaussian pairing strength with the shooting method and the FDM using the ADF. We show that both methods provide a consistent description for the discrete q.p. state in a box $ R_{\rm{box}}=20 $ fm. The asymptotic wave functions of the large components $ G_U $ and $ G_V $ of the $ 2p_{1/2} $ state obtained by the FDM are better than those provided by the shooting method owing to the ADF. For the continuum q.p. states, both the FDM and shooting method require a much larger box $ R_{\rm{box}}=40 $ fm to obtain reliable results. When continuum q.p. states are discretized by the box boundary condition, there are inevitably 'split' q.p. states around the resonant q.p. states with similar energies and fractional occupation probabilities. This is more evident for those q.p. states near the continuum threshold. However, these 'split' q.p. states superpose to be one definite canonical s.p. state by diagonalizing the density matrix, which is not particularly sensitive to the box size except for its asymptotic wave functions. With a relatively smaller box $ R_{\rm{box}}=20 $ fm, the asymptotic wave functions of the large components for the canonical states given by the FDM are the same as those obtained by the shooting method for the $ 2s_{1/2} $ state, but they are better for the $ 2p_{1/2} $ state, owing to the ADF used in the FDM. These investigations demonstrate the feasibility of the FDM to solve the RHB equation, which can be further used to develop future self-consistent RHB theories.