Introduction

The Dirac equation is a fundamental relativistic wave equation in quantum mechanics that describes the behavior of spin-1/2 massive particles. It is a generalization of the Schrödinger equation to account for relativistic effects. Solving Dirac equations with various potentials has played an important role in many areas of fundamental physics, including atomic, nuclear, hadron, and molecular systems. The Dirac equation has been studied in the context of nuclear physics, especially within the framework of relativistic mean-field theories [1, 2] and in explaining various single-particle effects [3-5]. One of the most important aspects of contemporary nuclear physics is the study of the shell evolution of exotic nuclei, which refers to the potentially dramatic changes in the shell structure as one approaches driplines with excessive protons or neutrons [6]. The evolution of the shell structure, which is crucial for our understanding of nuclear stability as well as the origin of heavy elements, can be induced by the isospin dependence of the spin and tensor potentials in various mean-field model approaches. The Dirac equation is unique for such studies because it provides not only insight into the origin of spin-orbital potentials as the competition between scalar and vector potentials [1] but also the possibility of adding a tensor potential at the mean field level. This study focuses on numerically solving the Dirac equation with scalar, vector, and tensor potentials. We hope that it can be a useful tool for studying the structure of exotic nuclei, as well as other quantum systems where spin-orbit and tensor effects can be important.

Dirac equation solvers are not available in the nuclear physics community, except for the Fortran subroutine embedded in relativistic mean field codes [7, 8]. Previous works in other areas on solving the Dirac equation focused particularly on the vector potential and utilized various schemes. In Ref. [9] the finite element method is applied to solve the Dirac equation, which is validated to converge to within floating point solutions. The problems inherent in the numerical solutions of the Dirac equation owing to spurious solutions are discussed in Ref. [10], where a stabilized version of the finite element method is also applied. The mapped Fourier method was used to numerically address the Dirac equation in. [11]. In Ref. [12] evolutionary algorithms are used and demonstrated in a case study of a muon orbiting the ²⁰⁸Pb nucleus. The work of Ref. [13] also focuses on studying muonic atoms, and applies the shooting method to integrate the Dirac equation. Moreover, a power series expansion method was employed in Ref.[14]. The numerical aspects of solving the Dirac equation are discussed in Ref. [15]. Other important methods include the Woods-Saxon basis [16], complex momentum representation [17], analytic continuation in the coupling constant (ACCC) [18] and Green’s function methods [19, 20].

Our work solves the Dirac equation in spherical coordinate space by utilizing the shooting method with a Runge-Kutta 4 integration scheme. We have prepared the code, named DiracSVT, in three different popular programming languages: Python, Matlab, and C++. These are gradually becoming the staple in physics, replacing Fortran as a more modern, versatile, and modular alternative. We include scalar, vector, and tensor potentials of an isospin-dependent Wood-Saxon form in the Dirac equation [21]. The code is used to solve the available single-particle states around doubly magic nuclei, including ¹⁶O, ⁴⁰Ca, ⁴⁸Ca, ⁵⁶Ni, ¹⁰⁰Sn, ¹³²Sn, and ²⁰⁸Pb. The same potential can also reproduce the shell structures of many open nuclei well over the entire nuclear chart.

The remainder of this paper is organized as follows. Section 2 describes the formulation of the Dirac equation with corresponding potentials and boundary conditions. Section 3 describes our method and the algorithm for solving the equation. Section 4 presents our results, and Sect. 5 discusses these results and presents an outlook on the possibilities for future extensions to this work.

Theory

We limit ourselves to spherical symmetry. The Dirac Hamiltonian H with a scalar potential $S$, a vector potential $V$ and a tensor potential $U$ can be given by (see, e.g., [22, 23]) $H=\overrightarrow{\alpha }\cdot \overrightarrow{p}+\beta (m+S)+V-i\beta \overrightarrow{\alpha }\cdot \widehat{r}U\mathrm{.}$ (1) The corresponding Dirac equation for the radial wave function can be expressed as $\left(\frac{\text{d}}{\text{d}r}+\frac{\kappa }{r}-U(r)\right)g_{\kappa }\left(r\right)=(E+m-\text{Δ}(r))f_{\kappa }\left(r\right),$ (2) $\left(\frac{\text{d}}{\text{d}r}-\frac{\kappa }{r}+U(r)\right)f_{\kappa }\left(r\right)=-(E-m-\text{Σ}(r))g_{\kappa }\left(r\right),$ (3) where f and g are the two components of the radial wave function and $\kappa =-(l+1)$ for $j=l+\frac{1}{2}$ and $\kappa =l$ for $j=l-\frac{1}{2}$. For convenience, we define our potentials as $\sum =V+S$ and $\text{Δ}=V-S$ as the sum of the vector and scalar potentials. The above equations can also be rewritten as ${g^{\prime }}_{\kappa }=(2m+B-\text{Δ})f_{\kappa }+\left(U-\frac{\kappa }{r}\right)g_{\kappa },$ (4) ${f^{\prime }}_{\kappa }=(-B+\text{Σ})g_{\kappa }+\left(\frac{\kappa }{r}-U\right)f_{\kappa },$ (5) where the binding energy B was obtained from the total energy $E=B+m$.

The potentials can be parameterized in a standard Woods-Saxon shape (see, e.g., [21, 24]) $\text{Σ}=\frac{{\text{Σ}}_{\mathrm{0,}n\left(p\right)}}{1+e^{\frac{r-R_{\sigma }}{a_{\sigma }}}}+V_{\text{Coulomb}},$ (6) $\text{Δ}=\frac{{\text{Δ}}_{\mathrm{0,}n\left(p\right)}}{1+e^{\frac{r-R_{\delta }}{a_{\delta }}}}+V_{\text{Coulomb}}\mathrm{,}$ (7) where $R_{\sigma \mathrm{,}\delta }=r_{\sigma \mathrm{,}\delta }{A}^{1/3}$ and a are radius and diffuseness parameters, respectively. Naturally, the Coulomb potential term exists only for protons.

Quantities ${\sum }_{\mathrm{0,}n\left(\text{p}\right)}\mathrm{,}$ and ${\text{Δ}}_{\mathrm{0,}n\left(\text{p}\right)}$ are defined for the three possible scenarios. In Scenario 1, ${\text{Δ}}_0$ is assumed to be directly proportional to ${\text{Σ}}_0$ with a scaling factor. Scenario 2 does not explicitly consider the dependence of the potential on the neutron excess. Scenario 3, on the other hand, describes dependence via a spin-orbit-dependent term ${\delta }_{\text{so}}$. For further details, please refer to Ref. [21]). For each of the three different scenarios, where n and p stand for neutron and proton states respectively, ${\text{Σ}}_{\mathrm{0,}\text{p}}=V_0\left(1+\delta \frac{N-Z}{A}\right),$ (8) ${\text{Σ}}_{\mathrm{0,}\text{n}}=V_0\left(1-\delta \frac{N-Z}{A}\right)\mathrm{,}$ (9) and ${\text{Δ}}_{\mathrm{0,}n\left(\text{p}\right)}$ is defined for Scenario 1 as ${\text{Δ}}_0=-\lambda {\text{Σ}}_0\mathrm{,}$ (10) for Scenario 2 as ${\text{Δ}}_0=-\lambda V_0\mathrm{,}$ (11) which assumes no isospin-dependence, and for Scenario 3 as ${\text{Δ}}_{\mathrm{0,}\text{p}}=-\lambda V_0\left(1-{\delta }_{\text{so}}\frac{N-Z}{A}\right),$ (12) ${\text{Δ}}_{\mathrm{0,}\text{n}}=-\lambda V_0\left(1+{\delta }_{\text{so}}\frac{N-Z}{A}\right)\mathrm{.}$ (13) The parameters $V_0\mathrm{,}\delta \mathrm{,}{\delta }_{\text{so}}\mathrm{,}\lambda \mathrm{,}{a}_{\sigma }=a_{\delta }\mathrm{,}{r}_0\mathrm{,}$ and $r_{\mathrm{0,}ls}$ are fitted to data for the three scenarios separately. In the simplest case, the diffuseness parameters $a_{\sigma }$ and $a_{\delta }$ are assumed identical.

The Coulomb potential is defined as $V_{\text{Coulomb}}=\{\begin{array}{ll}c\frac{Z}{r}\hfill & r>r_{\sigma }\hfill \\ c\frac{Z\left(3r_{\sigma }^2-r^2\right)}{2r_{\sigma }^3}\hfill & r\le r_{\sigma }\mathrm{,}\hfill \end{array}$ (14) where $c=\frac{e^2}{4\pi {ϵ}_0\hslash c}\equiv \alpha \approx \mathrm{0.007297.}$ (15) In addition, we have a tensor potential defined similarly in a Woods-Saxon shape: $U(r)=\frac{U_0}{1+e^{\frac{{\text{Δ}}^{\sigma }\left(r\right)}{a_{\sigma }}}}$ (16) The boundary condition for the bound states is defined separately for limiting cases of small ($r\to 0$) and large ($r\to \infty$) radii. For small r, we have $f=\{\begin{array}{ll}-a_0\left(\frac{-B+\text{Σ}}{k}\right){ϵ}^{l+2}\hfill & k<0\hfill \\ a_0{ϵ}^l\hfill & k>0\hfill \end{array},$ (17) $g=\{\begin{array}{ll}{a}_0{ϵ}^{l+1}\hfill & k<0\hfill \\ a_0\left(\frac{2m+B-\text{Δ}}{k}\right){ϵ}^{l+1}\hfill & k>0\hfill \end{array},$ (18) where $ϵ$ denotes the radial variable in the small radius regime. For the boundary conditions in the case where $r\to \infty$, it is assumed that both wave functions approach zero in this limit.

To handle unbound neutron resonance states with positive energy, we implement for large r similar boundary conditions to Eq. (26), as in [25]: $\psi =\left(\begin{array}{c}\left[C{j}_l\left(r\right)+D{y}_l\left(r\right)\right]Y_{j\mathrm{,}l\mathrm{,}{j}_z}\\ \frac{\sqrt{B^2+2Bm}}{B+2m}\left[C{j}_{l^{\prime }}\left(r\right)+D{y}_{l^{\prime }}\left(r\right)\right]Y_{j\mathrm{,}{l}^{\prime }\mathrm{,}{j}_z}\end{array}\right),$ where j and y indicate the spherical Bessel functions $l^{\prime }=j\pm \frac{1}{2}$ and $j_z=-j\mathrm{,}-j+\mathrm{1,}\dots j$. Here, $Y_{j\mathrm{,}l\mathrm{,}jz}$ are the spinor spherical harmonics, as detailed in [26]. In the 1D case, one of the two components of the spinor spherical harmonics vanishes, and the solution is reduced to two components f and g. As such the boundary conditions for $r\to \infty$ becomes $f=C{j}_l\left(r\right)+D{y}_l\left(r\right)\mathrm{,}$ (19) $g=\frac{\sqrt{B^2+2Bm}}{B+2m}\left[C{j}_{l^{\prime }}\left(r\right)+D{y}_{l^{\prime }}\left(r\right)\right]\mathrm{.}$ (20) The low-lying proton resonance states can be well approximated as bound states, as most of these states are pretty long-lived with negligibly small imaginary parts. The boundary condition should be refined as a Coulomb Hankel function for higher-lying proton resonance states.

Method

We attempt to combine an inner solution and an outer solution in what is known as the shooting method. The shooting method is a numerical scheme suitable for boundary-value problems. More specifically, the shooting method is a technique for solving a boundary value problem by transforming it into an Initial Value Problem (IVP). The method guesses the values of the boundary conditions at one end of the interval and then integrates the differential equation from that end using an ordinary differential equation solver. The estimated boundary values were adjusted until the solution satisfied the boundary conditions at the other end of the interval. No verification was performed for the number of nodes.

In our case, we employed shooting in two directions: outward and inward. We used the Runge-Kutta 4 integration scheme as the ordinary differential equation solver. This is akin to the methods employed in earlier studies, such as [27].

The outward shooting solution originated from the origin. However, setting r=0 causes a singularity in the model. To circumvent this, we set the inner point at a small distance $r=ϵ$ from the origin. As for initial values, the boundary conditions from Eqs. (17–20) are used. It should be noted that an initial guess for the binding energy is required. For simplicity, we used the experimental values with some variations for the example cases studied below.

Inward shooting originates from an asymptotic distance from the origin $r=\infty$. As an approximation, we consider a point that we set to a relatively large distance from the origin. We used the boundary conditions for a large r in this shooting case.

Both the outward and inward solutions can be used independently, as is. However, there is a major tendency for the calculated solutions of $f(r)$ and $g(r)$ to blow up and become numerically unstable. Thus, the matching concept of the shooting method is presented. We only solve from each direction up to a certain "middle" point which we denote as r_m (note, this point does not have to be halfway in distance) where the two solutions are then matched.

To define the difference $\text{Δ}\left(B_0\mathrm{,}{a}_0\right)$ between the two solutions at the matching point, there are a few options [15]. These include simple direct differences $\text{Δ}f\left(B_0\mathrm{,}{a}_0\right)=f_{\text{out}}\left(r_{\text{m}}\right)-f_{\text{in}}\left(r_{\text{m}}\right)\mathrm{,}$ (21) $\text{Δ}g\left(B_0\mathrm{,}{a}_0\right)=g_{\text{out}}\left(r_{\text{m}}\right)-g_{\text{in}}\left(r_{\text{m}}\right)\mathrm{,}$ (22) as well as relative average value based differences $\text{Δ}f\left(B_0\mathrm{,}{a}_0\right)=\frac{2\left(f_{\text{out}}\left(r_{\text{m}}\right)-f_{\text{in}}\left(r_{\text{m}}\right)\right)}{f_{\text{out}}\left(r_{\text{m}}\right)+f_{\text{in}}\left(r_{\text{m}}\right)}\mathrm{,}$ (23) $\text{Δ}g\left(B_0\mathrm{,}{a}_0\right)=\frac{2\left(g_{\text{out}}\left(r_{\text{m}}\right)-g_{\text{in}}\left(r_{\text{m}}\right)\right)}{g_{\text{out}}\left(r_{\text{m}}\right)+g_{\text{in}}\left(r_{\text{m}}\right)}\mathrm{.}$ (24) Here, $B_0$ and $a_0$ are the shooting parameters that must be adjusted such that the inside and outside Dirac spinor solutions match at the chosen matching radius $r_{\text{m}}$, that is, the direct differences in the equations above yield zero. After the trial, simple direct differences were found to converge most consistently and were used in the solver.

To drive the differences to zero, we must define the change in the difference in terms of the parameters. We do this by linearizing through a first order Taylor expansion $\text{Δ}f+\frac{\partial \text{Δ}f}{\partial B}\left(B_{\text{new}}-B_{\text{old}}\right)+\frac{\partial \text{Δ}f}{\partial a_0}\left(a_{\text{0,new}}-a_{\text{0,old}}\right)=\mathrm{0,}$ (25) $\text{Δ}g+\frac{\partial \text{Δ}g}{\partial B}\left(B_{\text{new}}-B_{\text{old}}\right)+\frac{\partial \text{Δ}g}{\partial a_0}\left(a_{\text{0,new}}-a_{\text{0,old}}\right)=0.$ (26) This defines a linear system, which we solve iteratively using Newton’s root-finding method. The resultant solution is considered to be the correct solution after satisfactory convergence, that is, the change in the solution becomes insignificant with further iterations. The overall solution method is shown in Fig. 1.

Fig. 1Full-screen View

Flowchart of the DiracSVT solver scheme. IVP stands for Initial Value Problem

Regarding the implementation of The numerical solver is written in three different languages: Python, MATLAB, and C++. It was developed and tested mainly for Linux, and the experiments were run on both Intel i7-2600 @3.40 GHz and AMD Ryzen 5 3600 @3.60 GHz. It is worth noting that the sensitivity of the solution to the precision of the system solution is very high in some cases. Changing the float precision can induce convergence to very different solutions. Similarly, merely varying the underlying linear algebra routine version (e.g., BLAS/LAPACK) also resulted in similar problems. As such, while the solver has been tested on different architectures, discrepancies remain regarding convergence that the reader should be aware of, resulting from different back-end libraries being utilized.

Numerical results

For example, we calculated the single-particle energies of all known single-particle states and single-hole states around the doubly magic nuclei for all three scenarios and all three language frameworks in which the solver is written. The binding energies against which the parameters were fitted were extracted from the experimental nuclear data. The parameters were fitted to minimize the RMS error between the numerical solutions and the experimental single-particle energies, as specified in Tables 2 and 3, so that they replicate the energies of the spin-orbit partners. The tensor potential is not included in the fitting of the parameters given below because of the scarcity of experimental data that can fully constrain its strength. Its effect is expected to be more significant as it approaches the neutron-rich nuclei around the neutron dripline.

In Table 1, the fitted parameters are given for Scenario 1, which were used to obtain the values in Tables 2 and 3. The difference between the scenarios relies on the value of δso: Scenario 2 is identical to Scenario 1 but takes δso=0 and Scenario 3 takes ${\delta }_{so}=-\delta$. The position was scaled by $1/\hslash c$, acquiring units of MeV^-1, yielding a Coulomb barrier in MeV as per Equation 15.

Some of the results are listed in Tables 2 and 3 for comparison purposes. The solutions from the three different code frameworks were identical within the numerical errors and were in reasonable agreement with the experimentally determined energies from Ref. [21]. The results for the three codes with the same inputs are listed. In some rare cases, one of the solvers can display a solution, while another diverges (marked as–in the tables). This seems likely to be the result of different languages using various routines to invert matrices. Moreover, the solver can sometimes converge to the energy levels of different states (e.g., the proton 2d_3/2 and 2f_7/2 states in ²⁰⁸Pb). This artifact can be easily spotted and is sensitive to the initial conditions and area of integration. We did not correct for comparison purposes. The convergence can be improved by slightly modifying the initial value, matching radius, and maximal radius. The values in Tables 2 and 3 are given for initial guesses equal to the experimental value for standardization. However, it may not always be the most suitable choice. In the case of unsuccessful convergence, this can be tweaked until the converged value improves. For some outliers, the code converges to a state with different nodes, which is easy to identify because the calculated energy differs significantly from the desired value.

The solver also captured the waveform solution of the Dirac equation for the considered particle and state, as shown in Fig. 2. In the present work, we focus on bound states and very low-lying unbound states, as all quantities are assumed to be real. The bound and unbound states appear quite similar, except for the boundary conditions. This also implies that the code is not yet suitable for handling wide resonances.

Fig. 2Full-screen View

Radial wave function for ¹⁶O, with neutron hole in state 1p_1/2. The results from the three codes were indistinguishable. The position was scaled by $1/\hslash c$, acquiring units of MeV^-1

Summary

We presented the Dirac equation solver DiracSVT, which solves the Dirac equation in coordinate space with Scalar, Vector, and tensor nuclear potentials. The code is written in three different languages and uses the shooting method with the Runge-Kutta 4 integration scheme. The above potentials were parameterized in the Saxon form. We show that the potential can reproduce the known single-particle states very well around all the doubly magic nuclei. This code is freely available in Ref. [28].