-
The neutron- or proton-rich nuclei far away from the
β− 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]
∫dr′(h−λΔΔ−h+λ)(ψUψV)=E(ψUψV),
(1) where h is the s.p. Hamiltonian, Δ is the pairing potential, λ is the Fermi energy, E is the q.p. energy, and
ψU andψV are the upper and lower components of the corresponding q.p. wave functions. Note that, for a bound system (λ<0 ), according to the asymptotic behavior of the wave functionsψU andψV in Eq. (1), the q.p. states with energy0<E<|λ| are defined as discrete q.p. states, and those with energyE>|λ| are continuum q.p. states [4, 9].In the RHB theory, the Dirac Hamiltonian can be expressed as
hD(r,r′)=[α⋅p+V(r)+β(M+S(r))]δ(r−r′),
(2) where
V(r) andS(r) are the vector and scalar potentials, respectively. For simplicity, we also consider a δ-function-type pairing potential Δ asΔ(r,r′)=Δ0(r)δ(r−r′).
(3) In the following, we assume spherical symmetry for the nuclear system. Then, the Dirac spinor wave functions
ψU(r) andψV(r) can be expressed asψU(r)=1r(GU(r)ΩljmiFU(r)Ωl′jm),ψV(r)=1r(GV(r)ΩljmiFV(r)Ωl′jm),
(4) where
l+l′=2j andl=j±12 , andΩljm is the spherical spinor. Then, the RHB equation can be reduced to four coupled equations for the radial wave functions, namelyGU ,FU andGV ,FV [17], as{−dFUdr+κrFU+(V+S−λ)GU+Δ0GV=EGUdGUdr+κrGU−(−V+S+2M+λ)FU+Δ0FV=EFUdFVdr−κrFV−(V+S−λ)GV+Δ0GU=EGV−dGVdr−κrGV+(−V+S+2M+λ)FV+Δ0FU=EFV,
(5) where the mean-field potentials are usually shifted by M for comparison with non-relativistic results, and
κ=(−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 forV(r) andS(r) ,V(r)±S(r)=V±01+er−r±0a±,
(6) where the parameters
V+0=−78.43 MeV,r+0=4.049 fm,a+=0.9254 fm, andV−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 56Ca [42]. The Fermi energy in Eq. (5) takes the valueλ=−4.00 MeV, read from the self-consistent result for neutrons in 56Ca. The pairing strengthΔ0(r) is studied according to the two following cases:Case1,constantpairingstrength:Δ0(r)=d0,
(7) Case2,Gaussianpairingstrength:Δ0(r)=d0e−x0(r−r0r0)2,
(8) where
d0=2 MeV,x0=5 , andr0=1.2A1/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,df(r)dr=f(r+h)−f(r−h)2h+O(h2),
(9) df(r)dr=−f(r+2h)+8f(r+h)−8f(r−h)+f(r−2h)12h+O(h4),
(10) 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),
df(r)dr=−3f(r)+4f(r+h)−f(r+2h)2h+O(h2),
(11) df(r)dr=3f(r)−4f(r−h)+f(r−2h)2h+O(h2),
(12) and those of the 5-point ADF (5PADF),
df(r)dr=−25f(r)+48f(r+h)−36f(r+2h)+16f(r+3h)−3f(r+4h)12h+O(h4), (13) df(r)dr=25f(r)−48f(r−h)+36f(r−2h)−16f(r−3h)+3f(r−4h)12h+O(h4).
(14) With these differential formulas, the RHB equation (Eq. (5)) can be expressed as a matrix in the coordinate space,
(AB1D0B2C0DD0−A−B10D−B2−C)(GUFUGVFV)=Hqp(GUFUGVFV)=E(GUFUGVFV),
(15) where
GU (FU ), andGV (FV ) are vectors of large (small) components of q.p. wave functionsψU andψV , respectively. The matrix elements in A,B1 ,B2 , and C are almost the same as those resulting from solving the Dirac equation in Ref. [50], except that the Fermi energy−λ 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Δ0 . Once the matrixHqp in Eq. (15) is diagonalized, one can obtain the eigenvalues of the q.p. energies E and the corresponding wave functionsGU ,FU ,GV , andFV in the coordinate space. In the following calculations, we set the box sizes in the coordinate space asRbox=20 fm andRbox=40 fm with lattice intervaldr=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 forr=0 and outside the boxr>Rbox . Here,f(r) denotes the radial wave functionsGU(FU) andGV(FV) in Eqs. (4) and (5).It is evident that by using the different differential formula in Eqs. (9)−(14), the matrix form of
Hqp 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 pointf(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 largeG(r) and smallF(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 forG(r) andF(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
Hqp 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
Δ0 , it can be proved that the solution of the q.p. energy for the RHB equation (Eq. (5)) can be expressed as [4]E=√(ε−λ)2+Δ20,
(16) with the corresponding wave functions
(GUFU)=c1(GF),
(17) (GVFV)=c2(GF),
(18) 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,
(V+Sκr−ddrκr+ddrV−S−2M)(GF)=ε(GF).
(19) The coefficients
c1 andc2 can be calculated as{c21=12(1+ε−λ√(ε−λ)2+Δ20),c22=12(1−ε−λ√(ε−λ)2+Δ20).
(20) 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]
v2=∫4πr2(G2V+F2V)dr
(21) of the neutrons for
κ=−1 (s1/2 ) andκ=+1 (p1/2 ) in 56Ca, calculated by the FDM using the 3PCDF, 5PCDF, 3PADF, and 5PADF to establish the q.p. HamiltonianHqp in the boxRbox=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 κ=−1 κ=1 κ=−1 κ=1 κ=−1 κ=1 κ=−1 κ=1 κ=−1 κ=1 49.481 (0.9996) \begin{array}{c}49.481\ ( 0.9996 )\end{array} 49.466(0.9996) \begin{array}{c}44.761 \( 0.9995 )\end{array} 49.436(0.9996) 49.466(0.9996) 49.466(0.9996) \begin{array}{c}30.446\ (0.9989)\end{array} 30.446(0.9989) 30.401(0.9989) 30.309(0.9989) 30.400(0.9989) 30.400(0.9989) 14.968(0.9955) \begin{array}{c}14.968 \( 0.9955)\end{array} 14.896(0.9955) 14.732(0.9954) 14.892(0.9955) 14.890(0.9954) \begin{array}{c}2.343 \( 0.7606 )\end{array} 2.343(0.7606) \begin{array}{c}14.645 \( 0.9953)\end{array} 2.302(0.7475) 2.226(0.7197) 2.302(0.7476) 2.301(0.7473) Table 1. Values of q.p. energy E and corresponding occupation probability
v2 in the brackets below for the neutrons in 56Ca withκ=−1 (s1/2 ) andκ=+1 (p1/2 ) calculated with constant pairing strength by using the FDM with different differential formulas and the shooting method in the boxRbox=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
κ=−1 andκ=+1 are exactly the same. To find out the reason, we present their corresponding wave functions in Fig. 1. Given that theGU(FU) andGV(FV) components of the wave functions differ by factorsc1 andc2 , respectively, as expected from Eqs. (17) and (18), here we only show theGV(FV) components of the wave functions. Figs. 1 (a) and (b) show the wave functions of the first q.p. state solutions forκ=−1 andκ=+1 , respectively, with the same q.p. energyE=49.481 MeV. Note that the wave functions in Fig. 1 (a) withκ=−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 the1s1/2 s.p. state owing to the Bogoliubov transformation. However, in Fig. 1 (b), forκ=+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. states1p1/2 ,2s1/2 , and2p1/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)
GV andFV components of q.p. wave functions of neutrons in 56Ca for the states withκ=−1 in (a) (c) (e) (g) andκ=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
κ=−1 andκ=+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 atf(r−h) andf(r+h) , but not at the middle pointf(r) . This leads to a unitary transformation matrix with alternative±1 diagonal elements between the Dirac HamiltonianHκ andH−κ , which produces energy-degenerate physical and spurious states. It can be demonstrated that such unitary transformation still exists between the q.p. HamiltonianHqp,κ andHqp,−κ . 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
Hqp . 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
v2 in brackets for neutrons withκ=−1 andκ=+1 in 56Ca obtained by the FDM with the 5PADF and the shooting method in the boxesRbox=20 fm andRbox=40 fm. First, for the states withκ=−1 , one can find that all the results for the largest q.p. energy of49.435 MeV provided by the different methods in box sizesRbox=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. state1s1/2 .κ=−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) κ=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
εcan with the corresponding occupation probabilitiesv2 in the brackets for neutrons in 56Ca withκ=−1 andκ=+1 calculated with Gaussian pairing strength by using the FDM with the 5PADF and the shooting method in the boxesRbox=20 fm andRbox=40 fm. The unit of the energy values is MeV.Taking the FDM results for
κ=−1 calculated in the boxRbox=20 fm as examples, one can find two nearly degenerate states with q.p. energies of14.945 MeV (v2=0.6441 ) and14.539 MeV (v2=0.3517 ) and comparable occupation probabilitiesv2 . Similar nearly degenerate states at14.972 MeV (v2=0.5618 ) and14.579 MeV (v2=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 probabilitiesv2 are clearly different. Therefore, we plot the wave functions ofGU ,FU ,GV , andFV in Figs. 2 (a) and (c) for the state of14.945 MeV (v2=0.6441 ) obtained by the FDM and14.972 MeV (v2=0.5618 ) obtained by the shooting method in the boxRbox=20 fm. The wave functions of the close states with q.p. energies of14.539 MeV obtained by the FDM and14.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 ofE≈15 MeV, which is larger than the Fermi energy|λ|=4.00 MeV, it should be a continuum q.p. state. Their wave functionsGU andFU must be oscillating, andGV andFV must have exponential decaying asymptotic behaviors, as clearly shown in Figs. 2 (a) and (c). Besides, the large component of the wave functionsGV with one node shows that this q.p. state must correspond to the2s1/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 boxRbox=20 fm, a smaller lattice interval, such asdr=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 fordr=0.1 fm as an example. Additionally, in the inset of Fig. 2 (e), we enlarge the component ofGV obtained by both methods within the same region,r=15 to20 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 componentsGU , 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
GU ,FU in (a) (b), andGV ,FV in (c) (d) of neutrons in 56Ca for the2s1/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)Rbox=20 fm and (b) (d)Rbox=40 fm. The insets (e) and (f) show the asymptotic wave functions ofGV up tor=20 fm. The unit of the wave functions is fm-1/2.When we calculate with a larger box size
Rbox=40 fm, the nearly degenerate q.p. states with comparable occupation probabilities seem to disappear. Explicitly, one can obtain the states with14.791 MeV (v2=0.9699 ) and15.797 MeV (v2=0.0188 ) from the FDM, and those with14.790 MeV (v2=0.9708) and15.829 MeV (v2=0.0177 ) from the shooting method. In particular, the states with relatively larger occupation probabilities, namely14.791 MeV (v2=0.9699 ) obtained from the FDM and14.790 MeV (v2=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 ofGV are also shown in a logarithmic scale in the inset of Fig. 2 (f). Their asymptotic behaviors are reasonable atr=20 fm. However, nearr=40 fm, their asymptotic behaviors are similar to those nearr=20 fm calculated in the boxRbox=20 fm in Fig. 2 (e), with much smaller wave function values. By comparing the results of the wave functionsGU obtained withinRbox=20 fm andRbox=40 fm in Figs. 2 (a) and (b), one can clearly see that they significantly differ, in particular atr≈5 fm andr≈20 fm. This demonstrates that a box size ofRbox=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ψU components.Using the model potentials in Eqs. (6)−(8), we can only find continuum q.p. states with
κ=−1 . In Table 2, we also list the results for the states withκ=+1 . Here, one can find a discrete q.p. state with an energy ofE≈1.65 MeV<|λ|=4.00 MeV, and occupation probability ofv2≈0.837 . For this state, all the calculations from the FDM and shooting method with boxesRbox=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 componentsGU ,FU ,GV , andFV 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 ofRbox=20 fm is large enough to describe these wave functions. From the nodes of theGU andGV components, we know that this q.p state should correspond to the2p1/2 state near the Fermi energy. In the inset of Figs. 3 (b) and (d), we enlarge the asymptotic wave functions fromr=15 to20 fm for the componentsGU andGV , respectively. It is interesting to note that, with a smaller box sizeRbox=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 the2s1/2 state in Fig. 2 (e), where the asymptotic wave functions ofGV components given by both the FDM and shooting method quickly drop to zero nearRbox=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 thes1/2 state has odd parity, so the forward asymmetric difference formula at the box boundaryRbox makes the wave functions quickly drop to zero, given that the boundary condition requires the wave functions atr≥Rbox to be zero. However, for the G component of the2p1/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 isRbox=40 fm, the shooting method also makes the wave functionsGU andGV 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
GU ,FU in (a), andGV ,FV in (c) of neutrons in 56Ca for the2p1/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 boxesRbox=20 fm andRbox=40 fm. The insets (b) and (d) show the asymptotic wave functions ofGU andGV up tor=20 fm. The unit of the wave functions is fm-1/2.In Table 2, for
κ=+1 , in addition to the discrete q.p. state with energyE≈1.65 MeV, there are also several states with energyE≈30 MeV and a large occupation probabilityv2>0.9 resulting from both the FDM and shooting method withRbox=20,40 fm. According to the aforementioned analysis, we know that these continuum q.p. states should correspond to the most deeply bound1p1/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 withE≈15 MeV and non-negligible occupation probabilitiesv2≈0.02 forκ=−1 calculated withRbox=40 fm. The wave functions of these 'split' q.p. states also have oscillatingψU and localψV components similar to those nearby states withv2>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
εcan and corresponding occupation probabilitiesv2 for the two most deeply bound states are listed in Table 2. For the states withκ=−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 boxRbox=20 fm as examples, it seems from the occupation probabilities that the canonical s.p. state2s1/2 'splits' into two nearly degenerate q.p. states withE≈15 MeV and comparablev2 . The similar 'splitting' can be also found in the results calculated with a larger boxRbox=40 fm, although one q.p. stateE≈14.79 MeV has a much largerv2≈0.97 , whereas the other one nearby has a much smallerv2≈0.02 . Forκ=+1 , one can also find this similar 'splitting' in the q.p. states atE≈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
2s1/2 and2p1/2 obtained from the FDM and shooting method in the boxesRbox=20 fm andRbox=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. For2s1/2 , the results obtained by both the FDM and the shooting method withRbox=20 fm are forced to be zero near the box boundary. The results calculated withRbox=40 fm have better asymptotic behaviors nearr=20 fm. However, nearr=40 fm, their asymptotic behaviors are similar to those nearr=20 fm calculated in the boxRbox=20 fm, with much smaller wave function values. For2p1/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) andF(r) for (a)2s1/2 and (b)2p1/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 boxesRbox=20 fm andRbox=40 fm. The insets (b) and (d) show the asymptotic wave functions ofG(r) up tor=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 56Ca. 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
Rbox=20 fm. The asymptotic wave functions of the large componentsGU andGV of the2p1/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 boxRbox=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 boxRbox=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 the2s1/2 state, but they are better for the2p1/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.
Solving the relativistic Hartree-Bogoliubov equation with the finite-difference method
- Received Date: 2024-08-10
- Available Online: 2025-01-15
Abstract: The relativistic Hartree-Bogoliubov (RHB) theory is a powerful tool for describing exotic nuclei near drip lines. The key technique is to solve the RHB equation in the coordinate space to obtain the quasi-particle states. In this paper, we solve the RHB equation with the Woods-Saxon-type mean-field and Delta-type pairing-field potentials by using the finite-difference method (FDM). We inevitably obtain spurious states when using the common symmetric central difference formula (CDF) to construct the Hamiltonian matrix, which is similar to the problem resulting from solving the Dirac equation with the same method. This problem is solved by using the asymmetric difference formula (ADF). In addition, we show that a large enough box is necessary to describe the continuum quasi-particle states. The canonical states obtained by diagonalizing the density matrix constructed by the quasi-particle states are not particularly sensitive to the box size. Part of the asymptotic wave functions can be improved by applying the ADF in the FDM compared to the shooting method with the same box boundary condition.