Time-Compact Scheme for the One-Dimensional Dirac Equation

Abstract

Based on the Lie-algebra, a new time-compact scheme is proposed to solve the one-dimensional Dirac equation. This time-compact scheme is proved to satisfy the conservation of discrete charge and is unconditionally stable. The time-compact scheme is of fourth-order accuracy in time and spectral order accuracy in space. Numerical examples are given to test our results.

[ProQuest: [...] denotes non US-ASCII text; see PDF]

Academic Editor:Cengiz Çinar

School of Applied Science, Beijing Information Science and Technology University, Beijing 100192, China

Received 26 July 2016; Accepted 4 September 2016

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The Dirac equation [1, 2], which is a spinor field equation and is applied to conquer the difficulty of the negative probability of Klein-Gordon equation, was proposed by the famous British physicist Dirac in 1928. In order to consider spin degree of freedom of electron, Dirac introduced multicomponent wave function and defined positive definite density of probability. The property of electron with high speed was discussed by using the Dirac equation. The fine structure of hydrogen atom was given when the Dirac equation was used to study the energy level distribution of hydrogen atom. In Dirac theory, spin-1/2 and intrinsic magnetic moment of electron can be obtained. The Dirac equation predicts the existence of the antiparticle partner to the electron. Dirac equation was widely studied after producing of the graphene in the lab in 2003 [3, 4].

In this paper, we consider the one-dimensional Dirac equation [5] [figure omitted; refer to PDF] where t is time, i=-1, and the magnetic potential _A1 (x) and electronic potential V(x) are real. The complex-valued wave function of spinor field is Φ(t,x)=^{(_[varphi]1 (t,x),_[varphi]2 (t,x))T} ∈R. _σ1 and _σ3 are Pauli matrices [figure omitted; refer to PDF] [straight epsilon] ([straight epsilon][: =]v/c) is a dimensionless parameter, v is the velocity of electron, and c is the light speed. When 0<v<<c, that is, the velocity of electron is far less than light velocity, there is 0<[straight epsilon]<<1, and (1) changes into the nonrelativistic model. On the contrary, it is relativistic model.

The Dirac equation (1) is dispersive and time symmetric. Here we introduce its position density _ρj (j=1,2) and the total density ρ which are defined as [figure omitted; refer to PDF] where ^{Φ[low *]} is conjugate transpose of Φ. Then, the charge of Dirac equation (1) is given as [figure omitted; refer to PDF]

For Dirac equation (1), it is hard to find the exact solution for general condition of electromagnetic potential. But for the special condition, such as coulomb problem, there was only nontrivial exact solution of Dirac equation [6]. The particle dynamics in relativistic quantum mechanics are described by using Dirac equation. Solution of Dirac equation is important in describing the nuclear shell structure [7]. In order to overcome the difficulty of solving analytically Dirac equation, many numerical methods, such as Crank-Nicolson finite difference method, time-splitting method [5], and pseudospectral method [8], have been used to solve numerically Dirac equation. To our best knowledge, the existing numerical methods for solving the Dirac equation have at most second-order accuracy in time. In this paper, we give the time-compact scheme which uses the fewer time steps to reach the fourth-order accuracy to solve the one-dimensional Dirac equation. This scheme may be extended to three-dimensional Dirac equation.

The arrangement of the rest for this paper is organized as follows. The time-compact scheme with fourth-order accuracy is presented in Section 2. Numerical experiment is given to test the accuracy order and conservation of discrete charge in Section 3. Some conclusions are drawn in Section 4.

2. Time-Compact Scheme and Analysis

Recently, there has been growing interest in high-order compact method for solving partial differential equation, especially the time-compact methods [9-14], which can enhance the accuracy order in time. In addition, time-compact scheme is efficient for dealing with high frequency oscillation problem.

In this section, we will use the time-compact scheme with fourth-order accuracy to solve the Dirac equation (1) numerically; that is, we will discretize the equation by adopting time-splitting method in time and pseudospectral method in space. The aim is to enhance the accuracy order to the fourth order in time.

In practical computation, the computational domain is Ω=(-a,a) (a>0), which is large enough. We divide the interval into M equal parts; then there are h=2a/M and _xj =-a+jh; here j=0,1,2,...,M. Choose time step τ; then _tn =nτ; here n=0,1,2,.... That is, [figure omitted; refer to PDF] with the homogeneous boundary condition [figure omitted; refer to PDF] and the initial condition [figure omitted; refer to PDF]

Setting [figure omitted; refer to PDF] (5) can be rewritten in form of [figure omitted; refer to PDF]

Obviously, (9) is a functional differential equation. So, we solve (1) in (_tn ,_tn +τ) and obtain the formal solution [figure omitted; refer to PDF] The key of solving the solution of (10) is to give the approximate value of the operator ^e-iτ(T+Z) .

For operators T and Z, if they are commutative operators, the operator ^e-iτ(T+Z) can be written as follows: [figure omitted; refer to PDF] In the general case, the operators T and Z are noncommutative, and the equation above is not tenable. Assuming that X and Y are noncommutative operators, according to the Baker-Campbell-Hausdorff formula (BCH) [11], ^eXeYeX can be expressed in the form of a single exponential function as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] Here we use the notation of the commutator [figure omitted; refer to PDF]

In order to solve the operator ^e-iτ(T+Z) , based on Lie-algebra and [12], one can obtain an approximate factorization in the form of [figure omitted; refer to PDF] Denote κ=-iτ; then [figure omitted; refer to PDF] According to (12), we have [figure omitted; refer to PDF] Then [figure omitted; refer to PDF] with [figure omitted; refer to PDF] In order to let κ_W3 match κ(T+Z), we must choose _a1 +_a2 =1/2 and 2_b1 +_b2 =1. Consider that the third-order term (T,T,Z) has the partial derivative of the unknown function (i.e., (T,T,Z)=(1/^{[straight epsilon]2} )_VxxI2 +(1/^{[straight epsilon]2} )_A1xxσ1 -(2i/^{[straight epsilon]3} )_A1σ3∂x +(2i/^{[straight epsilon]3} )_Vxσ1σ3 -(4/^{[straight epsilon]4} )_A1σ1 ), so we eliminate it by requiring 6^a22_b1 =1/4. When _a1 =1/6, _a2 =1/3, _b1 =3/8, and _b2 =1/4, the other third-order term (Z,Z,T) can be removed to yield the time-compact scheme [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

Since [figure omitted; refer to PDF] we have [figure omitted; refer to PDF] Therefore, [figure omitted; refer to PDF] So we can obtain [figure omitted; refer to PDF] From time t=_tn to time t=_tn+1 , for the operators T, Z, and Z~, we have three independent differential equations as follows: [figure omitted; refer to PDF]

According to the form of (26)-(28), (26) is discretized by the Fourier spectral method in space, (27) and (28) are functional differential system, and we can solve them analytically.

Next we begin to solve (10). In the first step, in (_tn ,_tn +τ/6), applying the pseudospectral method to solve (26) with initial value ^Φjn [approximate]Φ(_tn ,_xj ), we can obtain the solution ^Φj(1) via the inverse discrete Fourier transform as follows: [figure omitted; refer to PDF] where _μl =πl/a, j=0,1,2,...,M, and the coefficient _ψl (t) is obtained by using the discrete Fourier transform [figure omitted; refer to PDF] Substituting (29) into (26) and combining with the initial condition _ψl (_tn )=^ψln , taking t=_tn +τ/6, we have [figure omitted; refer to PDF] where [figure omitted; refer to PDF] Substituting ^ψl(1) (t) into (29), we can obtain [figure omitted; refer to PDF]

In the second step, we solve (27) in (_tn ,_tn +3τ/8) by using variable separation method, the initial value is ^Φj(1) , and we get the solution ^Φj(2) [figure omitted; refer to PDF] where G(x)=V(x)_I2 -_A1 (x)_σ1 =B_Λ1^{B[low *]} , [figure omitted; refer to PDF]

The third step is to solve (26) in [_tn +τ/6,_tn +τ/2] according to the first step with initial value ^Φj(2) . Thus, we can get the solution ^Φj(3) .

In the fourth step, we solve (28) by using the variable separation method in [_tn +3τ/8,_tn +5τ/8] with the initial value ^Φj(3) solved from the third step and write the solution as ^Φj(4) [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

The methods of the fifth, sixth, and seventh steps are obtained as the third, second, and first steps, respectively.

The time-compact scheme for solving the Dirac equation at (_tn ,_tn+1 ) can be obtained as follows: [figure omitted; refer to PDF] where j=0,1,2,...,M-1, l=-M/2,-M/2+1,...,M/2-1, and [figure omitted; refer to PDF] The above scheme is explicit, and it is easy to iterate.

Theorem 1.

The time-compact scheme (38a)-(38g) conserves the charge in the discretized level; that is, [figure omitted; refer to PDF]

Proof.

Introduce the definition of the discrete inner product; that is, [figure omitted; refer to PDF] where u=(_u1 ,_u2 ) and v=(_v1 ,_v2 ). Then we define the _l2 norm as follows: [figure omitted; refer to PDF] Equation (38a) is a Fourier transformation. From Parseval's equality combining with (30), one obtains [figure omitted; refer to PDF] For (38b), we have [figure omitted; refer to PDF] Similarly, the proof of (38c), (38e), and (38g) and (38d) and (38f) is obtained as the proof of (38a) and (38b), respectively. That is, [figure omitted; refer to PDF] So, one can obtain the conservation of discrete charge [figure omitted; refer to PDF]

3. Numerical Example

In this section, we test the order of accuracy and stability of the time-compact scheme. In order to test the accuracy, we choose the electromagnetic potentials in (5) as [figure omitted; refer to PDF] and the initial condition as [figure omitted; refer to PDF]

We solve problem (5)-(7) numerically under the condition Ω=(-8,8) from t=0 to t=1 by using the time-compact scheme. The numerical results calculated by the time-compact scheme (TCS) with fourth-order accuracy and Strang splitting method (SSM) [5] with second-order accuracy are listed in Table 1.

Table 1: (Example 1) comparison of errors, orders, and CPU times for different schemes with h=0.0156 at t=0.5.

τ	e ∞	Order	CPU time (s)
TCS	SSM	TCS	SSM	TCS	SSM
0.100000	6.2862 e - 05	2.8575 e - 02	--	--	7.310	3.300
0.050000	3.7013 e - 06	7.0465 e - 03	3.96	2.01	14.565	6.427
0.025000	2.2801 e - 07	1.7576 e - 03	4.00	2.00	29.199	12.680
0.012500	1.4218 e - 08	4.3919 e - 04	4.02	2.02	58.864	25.789
0.006250	9.0579 e - 10	1.0977 e - 04	4.10	2.00	116.640	51.960
0.003125	--	2.7441 e - 05	--	2.00	--	95.155

From Table 1, we can know that they have different accuracy order under the same condition, and the error of the time-compact scheme is much smaller than the error of Strang splitting method on the same row.

In every time step, the Strang splitting method with second-order accuracy needs three steps, and the time-compact scheme with fourth-order accuracy needs seven steps; that is, the total number of steps of the second-order accuracy scheme is 3/7 of the fourth-order accuracy scheme. But, from the column of error and time, we can find that the fourth-order accuracy scheme has much bigger step size and costs less time than the second-order accuracy when their errors have the same magnitude (such as ^10-5 ).

In order to test convergence of the algorithm with fourth-order O(^τ4 ), in the calculation, we verify the fourth order by taking enough small h to examine the temporal fourth order. Figure 1 takes log[...](_{(^en )∞} )-log[...](τ) of enough small h. We can know that the scheme in Section 2 is converged in maximum modulus [14], the convergence order for O(^τ4 ).

Figure 1: Fourth-order accuracy analysis in time for h=0.0156 and τ=1/160,1/80,1/40,1/20, and 1/10.

[figure omitted; refer to PDF]

The discrete charge calculated by the time-compact scheme is given at different time. As the calculated results have shown, one can see that the time-compact scheme conserves the discrete charge.

From Theorem 1, we make the conclusion that the time-compact scheme is unconditionally stable.

When [straight epsilon] decreases, high frequency oscillation increases. The time evolution of Re(_[varphi]1 (t,0)) for different [straight epsilon] is shown in Figure 2.

Figure 2: Oscillation of Dirac equation (1) with different [straight epsilon]=0.4, 0.2, and 0.1.

[figure omitted; refer to PDF]

4. Conclusion

Based on Lie-algebra, the time-compact scheme is presented for solving the one-dimensional linear Dirac equation. Then we test whether the time-compact scheme has fourth-order accuracy in time and is proved to keep the conservation of discrete charge. From the numerical results, the time-compact scheme performs much better than the Strang splitting method in the error analysis, in terms of accuracy and efficiency. In addition, the time-compact scheme is unconditionally stable, and numerical experiment is presented to discuss the changes of the frequency oscillation with different [straight epsilon]. In the past years, much work has been done to investigate the dynamical properties of nonlinear Schrödinger equations and other physically important nonlinear wave equations [15-18]. Thus, in the future, we will develop the time-compact scheme to study the nonlinear Dirac equation [19, 20], which is a model of self-interacting Dirac fermions in quantum field theory.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China (no. 11671044) and in part by the Beijing Municipal Education Commission under Grant no. PXM2016_014224_000028.

References

[1] P. A. M. Dirac, "The quantum theory of the electron," Proceedings of the Royal Society of London A: Mathematical and Physical Sciences , vol. 117, no. 778, pp. 610-624, 1928.

[2] P. A. M. Dirac, "A theory of electrons and protons," Proceedings of the Royal Society of London A , vol. 126, no. 801, pp. 360-365, 1930.

[3] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, A. A. Firsov, "Two-dimensional gas of massless Dirac fermions in graphene," Nature , vol. 438, no. 7065, pp. 197-200, 2005.

[4] D. A. Abanin, S. V. Morozov, L. A. Ponomarenko, R. V. Gorbachev, A. S. Mayorov, M. I. Katsnelson, K. Watanabe, T. Taniguchi, K. S. Novoselov, L. S. Levitov, A. K. Geim, "Giant nonlocality near the dirac point in graphene," Science , vol. 332, no. 6027, pp. 328-330, 2011.

[5] W. Z. Bao, Y. Y. Cai, X. W. Jia, Q. L. Tang, "Numerical methods and comparison for the Dirac equation in the nonrelativistic limit regime," https://arxiv.org/abs/1504.02881

[6] H. Bahlouli, E. B. Choubabi, A. Jellal, "Solution of one-dimensional Dirac equation via poincare map," Europhysics Letters , vol. 95, pp. 1305-1323, 2011.

[7] M. Hamzavi, A. A. Rajabi, "Solution of Dirac equation with Killingbeck potential by using wave function ansatz method under spin symmetry limit," Communications in Theoretical Physics , vol. 55, no. 1, pp. 35-37, 2011.

[8] H. Wang, "A time-splitting spectral method for computing dynamics of spinor F=1 Bose-Einstein condensates," International Journal of Computer Mathematics , vol. 84, pp. 925-944, 2007.

[9] J. L. Qiu, "Positive solutions for a nonlinear periodic boundary-value problem with a parameter," Electronic Journal of Differential Equations , vol. 2012, no. 133, pp. 1-10, 2012.

[10] R. D. Ruth, "A canonical integration technique," IEEE Transactions on Nuclear Science , vol. 30, no. 4, pp. 2669-2671, 1983.

[11] M. Suzuki, "General theory of higher-order decomposition of exponential operators and symplectic integrators," Physics Letters A , vol. 165, no. 5-6, pp. 387-395, 1992.

[12] S. A. Chin, "Symplectic integrators from composite operator factorizations," Physics Letters A , vol. 226, no. 6, pp. 344-348, 1997.

[13] M. Suzuki, "Hybrid exponential product formulas for unbounded operators with possible applications to Monte Carlo simulations," Physics Letters A , vol. 201, no. 5-6, pp. 425-428, 1995.

[14] H. Yoshida, "Construction of higher order symplectic integrators," Physics Letters A , vol. 150, no. 5-7, pp. 262-268, 1990.

[15] D.-S. Wang, X.-H. Hu, J. Hu, W. M. Liu, "Quantized quasi-two-dimensional Bose-Einstein condensates with spatially modulated nonlinearity," Physical Review A , vol. 81, no. 2, 2010.

[16] D.-S. Wang, S. J. Yin, Y. Tian, Y. Liu, "Integrability and bright soliton solutions to the coupled nonlinear Schrödinger equation with higher-order effects," Applied Mathematics and Computation , vol. 229, pp. 296-309, 2014.

[17] W. Z. Bao, Y. Y. Cai, H. Q. Wang, "Efficient numerical methods for computing ground states and dynamics of dipolar Bose-Einstein condensates," Journal of Computational Physics , vol. 229, no. 20, pp. 7874-7892, 2010.

[18] D.-S. Wang, X. Wei, "Integrability and exact solutions of a two-component Korteweg-de Vries system," Applied Mathematics Letters , vol. 51, pp. 60-67, 2016.

[19] F. W. Hehl, B. K. Datta, "Nonlinear spinor equation and asymmetric connection in general relativity," Journal of Mathematical Physics , vol. 12, no. 7, pp. 1334-1339, 1971.

[20] N. J. Poplawski, "Nonsingular Dirac particles in spacetime with torsion," Physics Letters B , vol. 690, no. 1, pp. 73-77, 2010.