Solution to the Dirac equation using the finite difference method

Solution to the Dirac equation using the finite difference method 1 first-author https://orcid.org/0000-0002-5635-5594 1 1 corresp Corresponding author, hength@ahu.edu.cn Corresponding author, hength@ahu.edu.cn hength@ahu.edu.cn https://orcid.org/0000-0003-1642-920X School of Physics and Materials Science, Anhui University, Hefei 230601, China School of Physics and Materials Science, Anhui University, Hefei 230601, China In this study, single-particle energy was examined using the finite difference method by taking 208 Pb as an example. If the first derivative term in the spherical Dirac equation is discretized using a three-point formula, a one-to-one correspondence occurs between the physical and spurious states. Although these energies are exactly the same, the wave functions of the spurious states exhibit a much faster staggering than those of the physical states. Such spurious states can be eliminated when applying the finite difference method by introducing an extra Wilson term into the Hamiltonian. Furthermore, it was also found that the number of spurious states can be reduced if we improve the accuracy of the numerical differential formula. The Dirac equation is then solved in a momentum space in which there is no differential operator, and we found that the spurious states can be completely avoided in the momentum space, even without an extra Wilson term. Finite difference method Spurious states Momentum space 1 Introduction 1 In recent decades, with a new generation of radioactive beam technologies and an upgrade in the detection equipment, the structure of exotic nuclei far from the β -stability line has become a hot research topic. The nuclei near the neutron or proton drip-lines exhibit numerous exotic physical phenomena, such as a neutron halo [ 1 1 ] or changes in the shell structure [ 2 2 ]. The relativistic density functional theory is a microscopic theoretical model widely applied in nuclear physics, and has provided satisfactory results [ 3 3 - 5 5 3-5 ], including nuclear binding energies [ 6 6 - 9 9 6-9 ], β -decay rates [ 10 10 - 12 12 10-12 ], low-lying spectra [ 13 13 , 14 14 ], and charge-exchange excitations [ 15 15 , 16 16 ]. In addition, numerous nuclear properties have been predicted, and have been widely used in r -process calculations [ 17 17 - 20 20 17-20 ] as well in as the dating of the cosmic age [ 21 21 - 23 23 21-23 ]. Moreover, the relativistic Hartree Bogoliubov theory has an important role in describing exotic nuclei [ 24 24 - 27 27 24-27 ]. Moreover, the density-dependent relativistic Hartree-Fock theory provides a satisfactory description of finite nuclei and nuclear matter [ 28 28 - 30 30 28-30 ]. A new and effective interaction, PKA1, has also been proposed that can eliminate the spurious shell structure and recover the N =64 subshell closure. In addition, the descriptions of single-particle spectra are improved by including tensor correlations [ 29 29 ]. Better systematics have been found in regions from the stable to neutron-rich sides [ 30 30 ]. Within this relativ istic framework, numerous methods for investigating the bound [ 27 27 , 31 31 - 33 33 31-33 ] and resonant [ 34 34 - 37 37 34-37 ] states have been developed. Therefore, the finite difference method was studied within the relativistic framework and the spurious states applied in the corresponding numerical calculations are discussed in detail herein. For the analysis of numerous physical phenomena, it is necessary to deduce the wave function from the Dirac equation. In relativistic quantum mechanics, the Dirac equation can usually be solved either analytically or numerically. However, the Hamiltonians for numerous quantum systems are generally complicated, and thus the analytic solutions of the Dirac equation are possibly not retained, and thus the application of advanced algorithms is required [ 38 38 , 39 39 ]. Based on the development of computer facilities, the base expansion method has been used to solve the Dirac equation in a multidimensional shape space [ 40 40 ] or a three-dimensional coordinate lattice space [ 41 41 ]. The numerical solutions are helpful in promoting the development of computational physics, and our understanding of an increasing number of physical phenomena. The common expansion methods include the momentum, harmonic oscillator, transformed harmonic oscillator, Woods-Saxon, and Laguerre bases, based upon which the Dirac Hamiltonian matrix is obtained by expanding the wave function. In addition to a base expansion method, a shooting method combined with the Runge-Kutta formula is commonly used in solving the Dirac equation within a coordinate space [ 3 3 , 28 28 ]. The finite difference method is another simple and effective way to solve the differential equation numerically because of its easy coding and debugging. Because the calculation program used with this method is widely suitable, it can deal with almost all forms of potential energy functions. It can then be used to apply more accurate values. However, the finite difference method applied for solving the Dirac equation will lead to spurious states [ 42 42 - 45 45 42-45 ]. Within a coordinate space, this is similar to that applied in a lattice QCD [ 46 46 ] caused by a numerical derivative formula. In a lattice QCD, to solve the problem of a spurious state, Wilson introduced the quadratic term of momentum [ 47 47 ], which contributes additionally to the dispersion relationship. Thus, the energy of the spurious states will be increased, and the degeneracy of the spurious and physical states will be eliminated. Within the framework of the relativistic mean field, to realize a self-consistent solution in a three-dimensional case, the Wilson term has been improved [ 41 41 ] through the introduction of an even higher order term of the momentum, allowing the spurious states to be excluded in the case of little influence on the physical states. The remainder of this paper is organized as follows. In Sec. 2, the theoretical framework for solving the Dirac equation is briefly introduced. Two methods are given to deal with the spurious states. In Sec. 3, we numerically apply the finite difference method to a spherical nucleus of 208 Pb and analyze the numerical details and results. In Sec. 4, a summary and some perspectives regarding the proposed approach are described. 2 Theoretical framework 1 The Dirac equation successfully merges quantum mechanics with special relativity. It provides a natural description of the electron spin, predicts the existence of antimatter, and is able to accurately reduce the spectrum of a hydrogen atom. In this study, the Dirac equation plays a significant role in describing the motion of the nucleons. In addition, from the numerical side it is regarded as a differential equation. Within the framework of relativistic mean-field theory, a nucleon is described as a Dirac particle moving within the field of vector potential V ( r ) and scalar potential S ( r ) [ 3 3 , 4 4 , 48 48 ]. The Dirac equation for the nucleon is as follows: <math display="block" id="M1"><mrow><mo stretchy="false">[</mo><mi>α</mi><mo>⋅</mo><mstyle mathvariant="bold" mathsize="normal"><mi>p</mi></mstyle><mo>+</mo><mi>V</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo stretchy="false">)</mo><mo>+</mo><mi>β</mi><mo stretchy="false">(</mo><mi>M</mi><mo>+</mo><mi>S</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mo stretchy="false">]</mo><msub><mi>ψ</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo stretchy="false">)</mo><mo>=</mo><msub><mi>ε</mi><mi>i</mi></msub><msub><mi>ψ</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo stretchy="false">)</mo><mn>,</mn></mrow></math> where α and β are Dirac matrices, ε i is the energy eigenvalue, and M is the nucleon mass. For spherical nuclei, the radial and angular parts of the eigen wave function can be separated into the following: <math display="block" id="M2"><mrow><msub><mi>ψ</mi><mrow><mi>α</mi><mi>κ</mi><mi>m</mi></mrow></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mn>,</mn><mi>s</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mn>1</mn><mi>r</mi></mfrac><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><mi>i</mi><msubsup><mi>G</mi><mi>α</mi><mi>κ</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><msubsup><mi>Y</mi><mi>l</mi><mrow><mi>j</mi><mi>m</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>θ</mi><mn>,</mn><mi>φ</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mrow><mo>−</mo><msubsup><mi>F</mi><mi>α</mi><mi>κ</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><msubsup><mi>Y</mi><mrow><mover accent="true"><mi>l</mi><mo stretchy="true">˜</mo></mover></mrow><mrow><mi>j</mi><mi>m</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>θ</mi><mn>,</mn><mi>φ</mi><mo stretchy="false">)</mo></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mn>,</mn><mi/><mi>j</mi><mo>=</mo><mi>l</mi><mo>±</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mn>.</mn></mrow></math> Based on Eq. 2, the one-dimensional radial part of the Dirac equation is written in matrix form as follows: <math display="block" id="M3"><mrow><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><mi>Σ</mi><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mrow></mtd><mtd><mrow><mfrac><mi>κ</mi><mi>r</mi></mfrac><mo>−</mo><mfrac><mi>d</mi><mrow><mi>d</mi><mi>r</mi></mrow></mfrac></mrow></mtd></mtr><mtr><mtd><mrow><mfrac><mi>d</mi><mrow><mi>d</mi><mi>r</mi></mrow></mfrac><mo>+</mo><mfrac><mi>κ</mi><mi>r</mi></mfrac></mrow></mtd><mtd><mrow><mi>Λ</mi><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><mi>G</mi><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mrow><mi>F</mi><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mo>=</mo><mi>E</mi><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><mi>G</mi><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mrow><mi>F</mi><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mn>,</mn></mrow></math> where <math display="block" id="M4"><mrow><mi>E</mi><mo>=</mo><mi>ε</mi><mo>−</mo><mi>M</mi><mn>,</mn></mrow></math> <math display="block" id="M5"><mrow><mi>Λ</mi><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>=</mo><mi>V</mi><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>−</mo><mi>S</mi><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>−</mo><mn>2</mn><mi>M</mi><mn>,</mn></mrow></math> <math display="block" id="M6"><mrow><mi>Σ</mi><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>=</mo><mi>V</mi><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>+</mo><mi>S</mi><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mn>.</mn></mrow></math> Here, the equation can be solved using the finite difference method within the coordinate space. The derivative term is replaced by the numerical derivative formula. In Eq. 3, the Hamiltonian matrix is symmetric. The Dirac equation (1) can then be solved and discussed based on the eigen solution of the Hamiltonian matrix as follows: <math display="block" id="M7"><mrow><msub><mi>H</mi><mi>κ</mi></msub><mo>=</mo><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><mi>Σ</mi><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mrow></mtd><mtd><mrow><mfrac><mi>κ</mi><mi>r</mi></mfrac><mo>−</mo><mfrac><mi>d</mi><mrow><mi>d</mi><mi>r</mi></mrow></mfrac></mrow></mtd></mtr><mtr><mtd><mrow><mfrac><mi>d</mi><mrow><mi>d</mi><mi>r</mi></mrow></mfrac><mo>+</mo><mfrac><mi>κ</mi><mi>r</mi></mfrac></mrow></mtd><mtd><mrow><mi>Λ</mi><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mn>.</mn></mrow></math> Taking the three-point numerical differential formula as an example, the first derivative term is approximately as follows: <math display="block" id="M8"><mrow><mfrac><mi>d</mi><mrow><mi>d</mi><mi>r</mi></mrow></mfrac><mo>=</mo><mfrac><mn>1</mn><mrow><mn>2</mn><mi>h</mi></mrow></mfrac><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mn>0</mn></mtd><mtd><mn>1</mn></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd><mtd><mn>0</mn></mtd><mtd><mn>1</mn></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow><mn>...</mn></mrow></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd><mtd><mn>0</mn></mtd><mtd><mn>1</mn></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mn>0</mn></mtd><mtd><mn>1</mn></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mn>,</mn></mrow></math> where h is the step size. The boundary conditions of the wave function can be obtained from Eq. 8. The results obtained through the finite difference method include spurious states. To eliminate these states, the Wilson term is introduced [ 41 41 ]. For the radial Dirac equation, the Wilson term can be expressed as follows: <math display="block" id="M9"><mrow><msub><mi>H</mi><mi>w</mi></msub><mo>=</mo><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><msup><mrow><mo stretchy="false">(</mo><mo>−</mo><mn>1</mn><mo stretchy="false">)</mo></mrow><mi>m</mi></msup><msup><mi>h</mi><mrow><mn>2</mn><mi>m</mi><mo>−</mo><mn>1</mn></mrow></msup><mi>R</mi><msub><mo>∂</mo><mi>r</mi></msub><msup><mrow/><mrow><mn>2</mn><mi>m</mi></mrow></msup></mrow></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mrow><msup><mrow><mo stretchy="false">(</mo><mo>−</mo><mn>1</mn><mo stretchy="false">)</mo></mrow><mrow><mi>m</mi><mo>+</mo><mn>1</mn></mrow></msup><msup><mi>h</mi><mrow><mn>2</mn><mi>m</mi><mo>−</mo><mn>1</mn></mrow></msup><mi>R</mi><msub><mo>∂</mo><mi>r</mi></msub><msup><mrow/><mrow><mn>2</mn><mi>m</mi></mrow></msup></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mn>,</mn></mrow></math> where h is the step size, and R and m are parameters. The Wilson term adopts the second and fourth derivatives of m of 1 and 2, respectively. The Hamiltonian matrix with the Wilson term becomes the following: <math display="block" id="M10"><mrow><mi>H</mi><mo>=</mo><msub><mi>H</mi><mi>κ</mi></msub><mo>+</mo><msub><mi>H</mi><mi>w</mi></msub><mn>.</mn></mrow></math> To clarify fully the influence of the spurious states, we also solve the Dirac equation within the momentum space in which the derivative term does not appear. The Dirac equation (1) can be expressed in a momentum representation as follows: <math display="block" id="M11"><mrow><mstyle displaystyle="true"><mrow><mo>∫</mo></mrow></mstyle><mi>d</mi><mstyle mathvariant="bold" mathsize="normal"><msup><mi>k</mi><mo>′</mo></msup></mstyle><mo>〈</mo><mstyle mathvariant="bold" mathsize="normal"><mi>k</mi></mstyle><mo>|</mo><mi>H</mi><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><msup><mi>k</mi><mo>′</mo></msup></mstyle><mo>〉</mo><mi>ψ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><msup><mi>k</mi><mo>′</mo></msup></mstyle><mo stretchy="false">)</mo><mo>=</mo><mi>ε</mi><mi>ψ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>k</mi></mstyle><mo stretchy="false">)</mo><mn>,</mn></mrow></math> where H is the Hamiltonian operator. In [ 36 36 ], the Hamiltonian matrix is given as follows: <math display="block" id="M12"><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover></mstyle><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><msubsup><mi>A</mi><mrow><mi>i</mi><mi>j</mi></mrow><mo>+</mo></msubsup></mrow></mtd><mtd><mrow><msub><mi>B</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>B</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></mtd><mtd><mrow><msubsup><mi>A</mi><mrow><mi>i</mi><mi>j</mi></mrow><mo>−</mo></msubsup></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><mi>f</mi><mo stretchy="false">(</mo><msub><mi>k</mi><mi>j</mi></msub><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mrow><mi>g</mi><mo stretchy="false">(</mo><msub><mi>k</mi><mi>j</mi></msub><mo stretchy="false">)</mo></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mo>=</mo><mi>ε</mi><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><mi>f</mi><mo stretchy="false">(</mo><msub><mi>k</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mrow><mi>g</mi><mo stretchy="false">(</mo><msub><mi>k</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mn>,</mn></mrow></math> where <math id="M13"><mrow><msubsup><mi>A</mi><mrow><mi>i</mi><mi>j</mi></mrow><mo>±</mo></msubsup><mo>=</mo><mo>±</mo><mi>M</mi><msub><mi>δ</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>+</mo><msqrt><mrow><msub><mi>w</mi><mi>i</mi></msub><msub><mi>w</mi><mi>j</mi></msub></mrow></msqrt><msub><mi>k</mi><mi>i</mi></msub><msub><mi>k</mi><mi>j</mi></msub><msup><mi>V</mi><mo>±</mo></msup><mo stretchy="false">(</mo><msub><mi>k</mi><mi>i</mi></msub><mn>,</mn><msub><mi>k</mi><mi>j</mi></msub><mo stretchy="false">)</mo></mrow></math> and B ij =- k i δ ij . Similar to the description provided in the previous subsection, we apply a diagonalization method to the Dirac Hamiltonian in the momentum representation, achieving satisfactory results. 3 Results and discussion 1 Based on the theoretical formalism mentioned in the above section, we investigated the single-particle energy levels of 208 Pb and dealt with the appearance of spurious states. The Dirac equation is solved using the finite difference method in a coordinate representation. A radial coordinate r that is 20 fm in size is discretized using N = 200 grid points. The single-particle energy levels and wave functions of 208 Pb are shown in Fig. 1 Fig. 1 along with the quantum numbers κ = 1 and κ = -1. Here, the first derivative term is expanded using a three-point numerical derivative formula. The two sets of energy levels with κ = ± 1 are as labeled in the figure. Both physical and spurious states are involved in the left and right panels. In Fig. 1 Fig. 1 , the wave functions, which oscillate extremely rapidly, are those of the spurious states. The amplitude of the spurious state is the same as that of the corresponding physical state, but the wave functions of the spurious states have opposite signs on adjacent mesh points. Owing to an improved calculation accuracy, the five-point derivative formula is adapted from this point. We found that the energy levels of the physical states remain almost unchanged, whereas the eigen energies of such spurious states move appreciably. In addition, an extra Wilson term can also be introduced into the Hamiltonian matrix; thus, the energies of the spurious states are pushed up again. Here, it was found that m = 2 can increase the energy levels of the spurious states faster as compared with m = 1. Owing to the increased energy levels, the spurious states are cleared from the low-energy region in Fig. 2 Fig. 2 relative to the left panel of Fig. 1 Fig. 1 . The states with the three lowest energy levels are the physical states. Moreover, the value of parameter R of the Wilson term affects the energy of the spurious state. For example, by adapting the value R = 0.005, the energy levels of the spurious states enter a continuous spectrum, whereas the energies and wave functions of the physical states remain the same. The energy spectrum in 208 Pb achieved by solving the Dirac equation is shown in Fig. 3 Fig. 3 . The eigen states 2s 1/2 and 3s 1/2 are represented by a star and circle, respectively, and the spurious states are represented by a square and triangle. For consistency in the calculation accuracy, the derivative in the Wilson term is computed using the seven-point formula with the parameter ( m = 2), and the first derivative term in the matrix H k is computed using the five-point formula. It can easily be seen from the figure that, for the physical states, the eigen energies are represented by the horizontal lines, and for the spurious states, the energy levels are linearly monotonous depending on the parameter R . Because parameter R reaches 0.001, the energies of the two spurious states are over 0 MeV, and thus are completely separated from the physical states. The appearance of a spurious state is caused by the treatment of the first derivative when using the numerical derivative formula. As can be seen from Fig. 1 Fig. 1 , when using the three-point formula, each pair of physical and spurious states has the same energy level. It is then necessary to remove the spurious states from the energy region we are focused upon. In addition, adding the Wilson term into the Hamiltonian matrix is another way to improve the calculation accuracy. Instead of the three-point formula, we show the results given by the five-point formula in Figs. 2 Figs. 2 and 3 3 . To study the effect of improving the calculation accuracy, we even generated a solution using the nine-point formula without the Wilson term. The wave functions obtained are given in Fig. 4 Fig. 4 by solving the Dirac equation. Comparing the left panel of Fig. 1 Fig. 1 , the energy levels of the spurious state are improved, particularly the energy level of the first spurious state, which increase from -51.559 to -27.576 MeV. However, the energy levels of the physical states remain mostly unchanged within the error considered. Improving the calculation accuracy is one way to increase the energy levels of the spurious states. Another way is introducing the Wilson term into the Hamiltonian matrix. Adjusting the values of parameters R and m , the spurious states can be pushed into the continuum spectrum. However, the spurious states do not completely disappear, which is a result of the numerical derivative formula. The Dirac equation can also be studied for a momentum representation in which the Hamiltonian matrix does not contain the first derivative term. The momentum operator is regarded as an ordinary parameter. The corresponding wave functions obtained in the momentum space are shown in Fig. 5 Fig. 5 . For comparison with the results given by the finite difference method, the wave functions were transformed into the coordinate space. The spurious states are cleared entirely and the wave functions for the physical states are the same as those in the above figures. The single-particle energies of 208 Pb obtained by solving the Dirac equation within the coordinate space using the finite difference method, as well as within the momentum space, are listed in Table 1 Table 1 . In the solution using the finite difference method, the first derivative term is discretized through the three-, five-, seven-, or nine-point differential formulas. Among them, the Wilson term is added into the Hamiltonian matrix when the five-point formula is used. The relative parameters m and R are set to 2 and 0.005, respectively. The energy levels of the physical states are almost the same when applying the different formulas with or without the Wilson term. The energy values in the parentheses are for the spurious states. The data are not equal in every column and the spurious states are mixed with the physical states. The data in the parentheses are the energy levels of the spurious state. All quantities are in MeV 1 1 Three-point formula 1 1 Five-point formula 1 1 Wilson term 1 1 Seven-point formula 1 1 Nine-point formula 1 1 In momentum space 1 1 -58.111(-51.559) 1 1 -58.108(-42.067) 1 1 -58.107 1 1 -58.107 1 1 -58.107 1 1 -58.107 1 1 -40.084(-29.068) 1 1 -40.051 1 1 -40.047 1 1 -40.051 1 1 -40.051 1 1 -40.051 1 1 -18.082(-7.723) 1 1 -17.979(-4.862) 1 1 -17.973 1 1 -17.979 1 1 -17.978 1 1 -17.987 1 1 -0.302 1 1 -0.223 1 1 -0.161 1 1 -0.222 1 1 -0.215 1 1 -0.242 4 Summary and perspectives 1 In this study, we first solved the Dirac equation within the coordinate space by applying the finite difference method. We employed a three-point numerical derivative formula to replace the first derivative term in the Hamiltonian matrix. Using 208 Pb as an example, the results prove that energy levels with quantum numbers of κ = ±1 decrease for the physical and spurious states. Improving the accuracy of the numerical differential formula and introducing the Wilson term into the Hamiltonian matrix are two ways to treat the spurious states. The five-point formula using the Wilson formula at m = 2 and R = 0.005 can be applied to push the spurious states into the continuum spectrum. The more accurate nine-point formula was adapted, and thus only one energy level of the spurious state reached a value of less than 0 MeV. Finally, we solved the Dirac equation within the momentum space in which there is no differential operator and the spurious states can be completely avoided. References I. Tanihata , Nuclear structure studies from reaction induced by radioactive nuclear beams . Prog. Part. Nucl. Phys . 35 , 505 ( 1995 ). https://doi.org/10.1016/0146-6410(95)00046-L https://doi.org/10.1016/0146-6410(95)00046-L Nuclear structure studies from reaction induced by radioactive nuclear beams A. Ozawa , T. Kobayashi , T. Suzuki et al ., New Magic Number, N=16, near the Neutron Drip Line . Phys. Rev. Lett 84 , 5493 ( 2000 ). https://doi.org/10.1103/PhysRevLett.84.5493 https://doi.org/10.1103/PhysRevLett.84.5493 New Magic Number, N=16, near the Neutron Drip Line J. Meng , H. Toki , S.G. Zhou et al ., Relativistic continuum Hartree Bogoliubov theory for ground-state properties of exotic nuclei . Prog. Part. Nucl. Phys . 57 , 470 ( 2006 ). https://doi.org/10.1016/j.ppnp.2005.06.001 https://doi.org/10.1016/j.ppnp.2005.06.001 Relativistic continuum Hartree Bogoliubov theory for ground-state properties of exotic nuclei D. Vretenar , A.V. Afanasjev , G.A. Lalazissis et al ., Relativistic Hartree-Bogoliubov theory: static and dynamic aspects of exotic nuclear structure . Phys. Rep . 409 , 101 ( 2005 ). https://doi.org/10.1016/j.physrep.2004.10.001 https://doi.org/10.1016/j.physrep.2004.10.001 Relativistic Hartree-Bogoliubov theory: static and dynamic aspects of exotic nuclear structure H.Z. Liang , J. Meng , S.G. Zhou , Hidden pseudospin and spin symmetries and their origins in atomic nuclei . Phys. Rep . 570 , 1 ( 2015 ). https://doi.org/10.1016/j.physrep.2014.12.005 https://doi.org/10.1016/j.physrep.2014.12.005 Hidden pseudospin and spin symmetries and their origins in atomic nuclei L.S. Geng , H. Toki , J. Meng , Masses, Deformations and Charge Radii-Nuclear Ground-State Properties in the Relativistic Mean Field Model . Prog. Theor. Phys . 113 , 785 ( 2005 ). https://doi.org/10.1143/PTP.113.785 https://doi.org/10.1143/PTP.113.785 Masses, Deformations and Charge Radii-Nuclear Ground-State Properties in the Relativistic Mean Field Model X.M. Hua , T.H. Heng , Z.M. Niu et al ., Comparative study of nuclear masses in the relativistic mean-field model . Sci. China. Phys. Mech . 55 , 2414 ( 2012 ). https://doi.org/10.1007/s11433-012-4943-y https://doi.org/10.1007/s11433-012-4943-y Comparative study of nuclear masses in the relativistic mean-field model Q.S. Zhang , Z.M. Niu , Z.P. Li et al ., Global dynamical correlation energies in covariant density functional theory: Cranking approximation . Front. Phys . 9 , 529 ( 2014 ). https://doi.org/10.1007/s11467-014-0413-5 https://doi.org/10.1007/s11467-014-0413-5 Global dynamical correlation energies in covariant density functional theory: Cranking approximation X.W. Xia , Y. Lim , P.W. Zhao et al ., The limits of the nuclear landscape explored by the relativistic continuum Hartree-Bogoliubov theory . At. Data Nucl. Data Tables 121-122 , 1 ( 2018 ). https://doi.org/10.1016/j.adt.2017.09.001 https://doi.org/10.1016/j.adt.2017.09.001 The limits of the nuclear landscape explored by the relativistic continuum Hartree-Bogoliubov theory Z.M. Niu , Y.F. Niu , H.Z. Liang et al ., β-decay half-lives of neutron-rich nuclei and matter flow in the r-process . Phys. Lett. B 723 , 172 ( 2013 ). https://doi.org/10.1016/j.physletb.2013.04.048 https://doi.org/10.1016/j.physletb.2013.04.048 β-decay half-lives of neutron-rich nuclei and matter flow in the r-process Z.M. Niu , Y.F. Niu , Q. Liu et al ., Nuclear β+/EC decays in covariant density functional theory and the impact of isoscalar proton-neutron pairing . Phys. Rev. C 87 , 051303(R) ( 2013 ). https://doi.org/10.1103/PhysRevC.87.051303 https://doi.org/10.1103/PhysRevC.87.051303 Nuclear β+/EC decays in covariant density functional theory and the impact of isoscalar proton-neutron pairing T. Marketin , L. Huther , and G. Martínez-Pinedo , Large-scale evaluation of β-decay rates of r-process nuclei with the inclusion of first-forbidden transitions . Phys. Rev. C 93 , 025805 ( 2016 ). https://doi.org/10.1103/PhysRevC.93.025805 https://doi.org/10.1103/PhysRevC.93.025805 Large-scale evaluation of β-decay rates of r-process nuclei with the inclusion of first-forbidden transitions Z.P. Li , T. Nikšić , D. Vretenar et al ., Microscopic analysis of nuclear quantum phase transitions in the N≈90 region . Phys. Rev. C 79 , 054301 ( 2009 ). https://doi.org/10.1103/PhysRevC.79.054301 https://doi.org/10.1103/PhysRevC.79.054301 Microscopic analysis of nuclear quantum phase transitions in the N≈90 region Z.P. Li , C.Y. Li , J. Xiang et al ., Enhanced collectivity in neutron-deficient Sn isotopes in energy functional based collective Hamiltonian . Phys. Lett. B 717 , 470 ( 2012 ). https://doi.org/10.1016/j.physletb.2012.09.061 https://doi.org/10.1016/j.physletb.2012.09.061 Enhanced collectivity in neutron-deficient Sn isotopes in energy functional based collective Hamiltonian N. Paar , P. Ring , T. Nikšić et al ., Quasiparticle random phase approximation based on the relativistic Hartree-Bogoliubov model . Phys. Rev. C 67 , 034312 ( 2003 ). https://doi.org/10.1103/PhysRevC.67.034312 https://doi.org/10.1103/PhysRevC.67.034312 Quasiparticle random phase approximation based on the relativistic Hartree-Bogoliubov model Z.M. Niu , Y.F. Niu , H.Z. Liang et al ., Self-consistent relativistic quasiparticle random-phase approximation and its applications to charge-exchange excitations . Phys. Rev. C 95 , 044301 ( 2017 ). https://doi.org/10.1103/PhysRevC.95.044301 https://doi.org/10.1103/PhysRevC.95.044301 Self-consistent relativistic quasiparticle random-phase approximation and its applications to charge-exchange excitations B. Sun , F. Montes , L. S. Geng et al ., Application of the relativistic mean-field mass model to the r-process and the influence of mass uncertainties . Phys. Rev. C 78 , 025806 ( 2008 ). https://doi.org/10.1103/PhysRevC.78.025806 https://doi.org/10.1103/PhysRevC.78.025806 Application of the relativistic mean-field mass model to the r-process and the influence of mass uncertainties X.D. Xu , B. Sun , Z.M. Niu et al ., Reexamining the temperature and neutron density conditions for r-process nucleosynthesis with augmented nuclear mass models . Phys. Rev. C 87 , 015805 ( 2013 ). https://doi.org/10.1103/PhysRevC.87.015805 https://doi.org/10.1103/PhysRevC.87.015805 Reexamining the temperature and neutron density conditions for r-process nucleosynthesis with augmented nuclear mass models J.S. Zheng , N.Y. Wang , Z.Y. Wang et al ., Mass predictions of the relativistic mean-field model with the radial basis function approach . Phys. Rev. C 90 , 014303 ( 2014 ). https://doi.org/10.1103/PhysRevC.90.014303 https://doi.org/10.1103/PhysRevC.90.014303 Mass predictions of the relativistic mean-field model with the radial basis function approach Z. Li , Z.M. Niu , B.H. Sun , Influence of nuclear physics inputs and astrophysical conditions on r-process . Sci. China. Phys. Mech 62 , 982011 ( 2019 ). https://doi.org/10.1007/s11433-018-9355-y https://doi.org/10.1007/s11433-018-9355-y Influence of nuclear physics inputs and astrophysical conditions on r-process Z.M. Niu , B. Sun , and J. Meng , Influence of nuclear physics inputs and astrophysical conditions on the Th/U chronometer . Phys. Rev. C 80 , 065806 ( 2009 ). https://doi.org/10.1103/PhysRevC.80.065806 https://doi.org/10.1103/PhysRevC.80.065806 Influence of nuclear physics inputs and astrophysical conditions on the Th/U chronometer J. Meng , Z.M. Niu , H.Z. Liang et al ., Selected issues at the interface between nuclear physics and astrophysics as well as the standard model . Sci. China: Phys. Mech. Astron . 54 , S119 ( 2011 ). https://doi.org/10.1007/s11433-011-4439-1 https://doi.org/10.1007/s11433-011-4439-1 Selected issues at the interface between nuclear physics and astrophysics as well as the standard model W.H. Zhang , Z.M. Niu , F. Wang et al ., Uncertainties of nucleo-chronometers from nuclear physics inputs . Acta Phys. Sin . 61 , 112601 ( 2012 ). https://doi.org/10.7498/aps.61.112601 https://doi.org/10.7498/aps.61.112601 Uncertainties of nucleo-chronometers from nuclear physics inputs J. Meng , P. Ring , Relativistic Hartree-Bogoliubov Description of the Neutron Halo in 11Li . Phys. Rev. Lett . 77 , 3963 ( 1996 ). https://doi.org/10.1103/PhysRevLett.77.3963 https://doi.org/10.1103/PhysRevLett.77.3963 Relativistic Hartree-Bogoliubov Description of the Neutron Halo in 11Li J. Meng , P. Ring , Giant Halo at the Neutron Drip Line . Phys. Rev. Lett . 80 , 460 ( 1998 ). https://doi.org/10.1103/PhysRevLett.80.460 https://doi.org/10.1103/PhysRevLett.80.460 Giant Halo at the Neutron Drip Line J. Meng , H. Toki , J.Y. Zeng et al ., Giant halo at the neutron drip line in Ca isotopes in relativistic continuum Hartree-Bogoliubov theory . Phys. Rev. C 65 , 041302(R) ( 2002 ). https://doi.org/10.1103/PhysRevC.65.041302 https://doi.org/10.1103/PhysRevC.65.041302 Giant halo at the neutron drip line in Ca isotopes in relativistic continuum Hartree-Bogoliubov theory S.G. Zhou , J. Meng , P. Ring et al ., Neutron halo in deformed nuclei . Phys. Rev. C 82 , 011301(R) ( 2010 ). https://doi.org/10.1103/PhysRevC.82.011301 https://doi.org/10.1103/PhysRevC.82.011301 Neutron halo in deformed nuclei W.H. Long , N.V. Giai , J. Meng , Density-dependent relativistic Hartree - Fock approach . Phys. Lett. B 640 , 150 ( 2006 ). https://doi.org/10.1016/j.physletb.2006.07.064 https://doi.org/10.1016/j.physletb.2006.07.064 Density-dependent relativistic Hartree - Fock approach W.H. Long , H. Sagawa , N.V. Giai et al ., Shell structure and ρ-tensor correlations in density dependent relativistic Hartree-Fock theory . Phys. Rev. C 76 , 034314 ( 2007 ). https://doi.org/10.1103/PhysRevC.76.034314 https://doi.org/10.1103/PhysRevC.76.034314 Shell structure and ρ-tensor correlations in density dependent relativistic Hartree-Fock theory W.H. Long , P. Ring , N.V. Giai et al ., Relativistic Hartree-Fock-Bogoliubov theory with density dependent meson-nucleon couplings . Phys. Rev. C 81 , 024308 ( 2010 ). https://doi.org/10.1103/PhysRevC.81.024308 https://doi.org/10.1103/PhysRevC.81.024308 Relativistic Hartree-Fock-Bogoliubov theory with density dependent meson-nucleon couplings J. Meng , Relativistic continuum Hartree-Bogoliubov theory with both zero range and finite range Gogny force and their application . Nucl. Phys. A 635 , 3 ( 1998 ). https://doi.org/10.1016/S0375-9474(98)00178-X https://doi.org/10.1016/S0375-9474(98)00178-X Relativistic continuum Hartree-Bogoliubov theory with both zero range and finite range Gogny force and their application Z.X. Ren , S.Q. Zhang , and J. Meng , Solving Dirac equations on a 3D lattice with inverse Hamiltonian and spectral methods . Phys. Rev. C 95 , 024313 ( 2017 ). https://doi.org/10.1103/PhysRevC.95.024313 https://doi.org/10.1103/PhysRevC.95.024313 Solving Dirac equations on a 3D lattice with inverse Hamiltonian and spectral methods S.G. Zhou , J. Meng , and P. Ring , Spherical relativistic Hartree theory in a Woods-Saxon basis . Phys. Rev. C 68 , 034323 ( 2003 ). https://doi.org/10.1103/PhysRevC.68.034323 https://doi.org/10.1103/PhysRevC.68.034323 Spherical relativistic Hartree theory in a Woods-Saxon basis J.Y. Guo , X.Z. Fang , P. Jiao et al ., Application of the complex scaling method in relativistic mean-field theory . Phys. Rev. C 82 , 034318 ( 2010 ). https://doi.org/10.1103/PhysRevC.82.034318 https://doi.org/10.1103/PhysRevC.82.034318 Application of the complex scaling method in relativistic mean-field theory Z.L. Zhu , Z.M. Niu , D.P. Li et al ., Probing single-proton resonances in nuclei by the complex-scaling method . Phys. Rev. C 89 , 034307 ( 2014 ). https://doi.org/10.1103/PhysRevC.89.034307 https://doi.org/10.1103/PhysRevC.89.034307 Probing single-proton resonances in nuclei by the complex-scaling method N. Li , M. Shi , J.Y. Guo et al ., Probing Resonances of the Dirac Equation with Complex Momentum Representation . Phys. Rev. Lett . 117 , 062502 ( 2016 ). https://doi.org/10.1103/PhysRevLett.117.062502 https://doi.org/10.1103/PhysRevLett.117.062502 Probing Resonances of the Dirac Equation with Complex Momentum Representation M. Shi , Z.M. Niu , and H.Z. Liang , Combination of complex momentum representation and Green’s function methods in relativistic mean-field theory . Phys. Rev. C 97 , 064301 ( 2018 ). https://doi.org/10.1103/PhysRevC.97.064301 https://doi.org/10.1103/PhysRevC.97.064301 Combination of complex momentum representation and Green’s function methods in relativistic mean-field theory T.S. Kosmas , I.E. Lagaris , On the muon - nucleus integrals entering the neutrinoless μ- e- conversion rates . J. Phys. G: Nucl. Part. Phys . 28 , 2907 ( 2002 ). https://doi.org/10.1088/0954-3899/28/12/302 https://doi.org/10.1088/0954-3899/28/12/302 On the muon - nucleus integrals entering the neutrinoless μ- e- conversion rates K. Ryuichiro , K. Masafumi , O. Yasuhiro , Detailed calculation of lepton flavor violating muon-electron conversion rate for various nuclei . Phys. Rev. D . 66 , 096002 ( 2002 ). https://doi.org/10.1103/PhysRevD.66.096002 https://doi.org/10.1103/PhysRevD.66.096002 Detailed calculation of lepton flavor violating muon-electron conversion rate for various nuclei B.N. Lu , J. Zhao , E.G. Zhao , Multidimensionally-constrained relativistic mean-field models and potential-energy surfaces of actinide nuclei . Phys. Rev. C 89 , 014323 ( 2014 ). https://doi.org/10.1103/PhysRevC.89.014323 https://doi.org/10.1103/PhysRevC.89.014323 Multidimensionally-constrained relativistic mean-field models and potential-energy surfaces of actinide nuclei Y. Tanimura , K. Hagino , H.Z. Liang , 3D mesh calculations for covariant density functional theory . Prog. Theor. Exp. Phys . 7 , 073D01 ( 2015 ). https://doi.org/10.1093/ptep/ptv083 https://doi.org/10.1093/ptep/ptv083 3D mesh calculations for covariant density functional theory V.M. Shabaev , I.I. Tupitsyn , V.A. Yerokhin et al ., Dual Kinetic Balance Approach to Basis-Set Expansions for the Dirac Equation . Phys. Rev. Lett 93 , 130405 ( 2004 ). https://doi.org/10.1103/PhysRevLett.93.130405 https://doi.org/10.1103/PhysRevLett.93.130405 Dual Kinetic Balance Approach to Basis-Set Expansions for the Dirac Equation G. Pestka , Spurious roots in the algebraic Dirac equation . Phys. Scripta . 68 , 254 2003 . https://doi.org/10.1016/S0009-2614(03)01061-3 https://doi.org/10.1016/S0009-2614(03)01061-3 Spurious roots in the algebraic Dirac equation C. Müller , N. Grün , W. Scheid , Finite element formulation of the Dirac equation and the problem of fermion doubling . Phys. Lett. A 242 , 245 ( 1998 ). https://doi.org/10.1016/S0375-9601(98)00218-7 https://doi.org/10.1016/S0375-9601(98)00218-7 Finite element formulation of the Dirac equation and the problem of fermion doubling S. Zhao , On the spurious solutions in the high-order finite difference methods for eigenvalue problems . Comput. Method. Appl. M . 196 , 5031 ( 2007 ). https://doi.org/10.1016/j.cma.2007.06.018 https://doi.org/10.1016/j.cma.2007.06.018 On the spurious solutions in the high-order finite difference methods for eigenvalue problems S. Salomonson , P. Öster , Relativistic all-order pair functions from a discretized single-particle Dirac Hamiltonian . Phys. Rev. A 40 , 5548 ( 1989 ). https://doi.org/10.1103/PhysRevA.40.5548 https://doi.org/10.1103/PhysRevA.40.5548 Relativistic all-order pair functions from a discretized single-particle Dirac Hamiltonian K.G. Wilson , Proceedings of the First Half of the 1975 International School of Subnuclear Physics Erice , Sicily , July 11 - August 1, 1975 p 69 B.D. Serot , J.D. Walecka , The Relativistic Nuclear Many Body Problem . Adv. Nucl. Phys . 16 , 1 ( 1986 ). https://doi.org/10.1007/978-1-4684-5179-5_8 https://doi.org/10.1007/978-1-4684-5179-5_8 The Relativistic Nuclear Many Body Problem 10.1007/s41365-020-0728-6 This work was partly supported by the National Natural Science Foundation of China (No.11875070), and the Natural Science Foundation of Anhui Province (No. 1908085MA16). 2019-09-12 2019-12-26 2020-01-05 2020-01-31 2020-02-01 2021-04-10 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2020 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-020-0728-6 http://dx.doi.org/10.1007/s41365-020-0728-6 . 2020 Nuclear Science and Techniques 2 31