# Time-Compact Scheme for the One-Dimensional Dirac Equation

Authors: Jun-Jie Cao; Xiang-Gui Li; Jing-Liang Qiu; Jing-Jing Zhang
Journal: Discrete Dynamics in Nature and Society
 (2016)

---

## 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.


---

## Body

## 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](1)i∂tΦt,x=-iεσ1∂x+1ε2σ3Φt,x+VxI2-A1xσ1Φt,x,x∈R,wheretis time, i=-1, and the magnetic potential A1x and electronic potential Vx are real. The complex-valued wave function of spinor field is Φt,x=ϕ1t,x,ϕ2t,xT∈R.σ1 and σ3are Pauli matrices(2)σ1=0110,σ3=100-1,I2=1001.εε≔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<ε≪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(3)ρjt,x=ϕjt,x2,j=1,2,ρt,x=∑j=12ρjt,x=Φ∗Φ,where Φ∗ is conjugate transpose of Φ. Then, the charge of Dirac equation (1) is given as(4)Φt,·2≔∫RΦt,x2dx=∫R∑j=1ϕjt,x2dx.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 Section2. 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,aa>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,(5)i∂tΦt,x=-iεσ1∂x+1ε2σ3Φt,x+VxI2Φt,x-A1xσ1Φt,x,x∈-a,a,t>0,with the homogeneous boundary condition(6)Φt,-a=Φt,a,t≥0and the initial condition(7)Φ0,x=Φ0x,x∈Ω¯.Setting(8)T=-iεσ1∂x+1ε2σ3,Z=VxI2-A1xσ1, (5) can be rewritten in form of(9)i∂Φt,x∂t=T+ZΦt,x.Obviously, (9) is a functional differential equation. So, we solve (1) in tn,tn+τ and obtain the formal solution(10)Φtn+1,x=e-iτT+ZΦtn,x.The key of solving the solution of (10) is to give the approximate value of the operator e-iτT+Z.For operatorsT and Z, if they are commutative operators, the operator e-iτT+Z can be written as follows: (11)e-iτT+Z=e-iτTe-iτZ=e-iτZe-iτT.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 theBaker-Campbell-Hausdorff formula (BCH) [11], eXeYeX can be expressed in the form of a single exponential function as(12)eXeYeX=eW,where(13)W=2X+Y+16Y,Y,X-16X,X,Y+7360X,X,X,X,Y+1360Y,Y,X+190X,Y,Y,Y,X+145Y,X,X,X,Y-160X,X,Y,Y,X+130Y,Y,X,X,Y+⋯.Here we use the notation of the commutator(14)X,Y≔XY-YX,X,X,Y≔X,X,Y.In order to solve the operatore-iτT+Z, based on Lie-algebra and [12], one can obtain an approximate factorization in the form of(15)e-iτT+Z=e-iτa1Te-iτb1Ze-iτa2Te-iτb2Ze-iτa2Te-iτb1Ze-iτa1T.Denoteκ=-iτ; then(16)eκa2Teκb2Zeκa2T=eκW1,eκb1ZeκW1eκb1Z=eκW2,eκa1TeκW2eκa1T=eκW3.According to (12), we have(17)κW1=2κa2T+κb2Z+16κ3b2Z,b2Z,a2T-16κ3a2T,a2T,b2Z+Oκ5,κW2=2κb1Z+κW1+16κW1,κW1,κb1Z-16κb1Z,κb1Z,κW1+Oκ5,κW3=2κa1T+κW2+16κW2,κW2,κa1T-16κa1T,κa1T,κW2+Oκ5.Then(18)eκa1Teκb1Zeκa2Teκb2Zeκa2Teκb1Zeκa1T=eκW3,with(19)κW3=2κa1+a2T+κ2b1+b2Z+16κ3a2b22-2b1a2b2-2b12a2+a12b1+b22Z,Z,T-16κ3a22b2-4a22b1+2a12b1+b2+a122b1+b2T,T,Z+κ5.In order to letκW3match κT+Z, we must choose a1+a2=1/2 and 2b1+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/ε2VxxI2+1/ε2A1xxσ1-(2i/ε3)A1σ3∂x+2i/ε3Vxσ1σ3-4/ε4A1σ1), so we eliminate it by requiring 6a22b1=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(20)e-iτT+Z=e-iτ1/6Te-iτ3/8Ze-iτ1/3Te-iτ1/4Z~e-iτ1/3Te-iτ3/8Ze-iτ1/6T+Oτ5,where (21)Z~=Z-τ2Z,T,Z48.Since(22)T,ZΦ=TZΦ-ZTΦ=-iεσ1∂x+1ε2σ3VI2-A1σ1ΦVI2-A1σ1-iεσ1∂x+1ε2σ3Φ=-iεσ1∂xVI2Φ-A1σ1Φ+1ε2σ3VI2Φ-A1σ1Φ-VI2-iεσ1Φx+1ε2σ3Φ+A1σ1-iεσ1Φx+1ε2σ3Φ=-iεσ1VxΦ+iεI2A1xΦ+21ε2σ1σ3A1Φ,we have(23)T,Z=-iεσ1Vx+iεI2A1x+21ε2σ1σ3A1.Therefore,(24)Z,T,ZΦ=ZT,Z-T,ZZΦ=ZT,ZΦ-T,ZZΦ=VI2-A1σ1-iεσ1Vx+iεI2A1x+21ε2σ1σ3A1Φ--iεσ1Vx+iεI2A1x+21ε2σ1σ3A1VI2-A1σ1Φ=VI2-A1σ1-iεσ1VxΦ+iεI2A1xΦ+21ε2σ1σ3A1Φ--iεσ1Vx+iεI2A1x+21ε2σ1σ3A1VI2Φ-A1σ1Φ=-41ε2σ3A1t,x2Φ.So we can obtain(25)Z,T,Z=-41ε2σ3A1x2,Z~=VxI2-A1xσ1+τ212ε2σ3A1x2.From time t=tn to time t=tn+1, for the operators T, Z, and Z~, we have three independent differential equations as follows:(26)i∂Φt,x∂t=-iεσ1∂x+1ε2σ3Φt,x=TΦt,x,(27)i∂Φt,x∂t=VxI2-A1xσ1Φt,x=ZΦt,x,(28)i∂Φt,x∂t=VxI2-A1xσ1+τ212ε2σ3A1x2Φt,x=Z~Φt,x.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≈Φtn,xj, we can obtain the solution Φj(1) via the inverse discrete Fourier transform as follows:(29)Φj1=1M∑l=-M/2M/2-1ψl1texp⁡iμlxj+a,where μl=πl/a, j=0,1,2,…,M, and the coefficient ψlt is obtained by using the discrete Fourier transform(30)ψlt=∑j=0M-1Φt,xjexp⁡-iμlxj+a,l=-M2,-M2+1,…,M2-1.Substituting (29) into (26) and combining with the initial condition ψltn=ψln, taking t=tn+τ/6, we have(31)ψl1t=exp⁡-iτ61εμlσ1+1ε2σ3ψln=exp⁡-iτ6ε2Γψln=Qlexp⁡-iτ6ε2DlQl∗ψln,where(32)Γ=1εμlεμl-1,D1=λl00λl,Ql=1+λl2λl1+λl-εμl2λl1+λlεμl2λl1+λl1+λl2λl1+λl,λl=1+εμl2.Substituting ψl(1)t into (29), we can obtain(33)Φj1=1M∑l=-M/2M/2-1Qlexp⁡-iτ16DlQl∗ψlnexp⁡iμlxj+a.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)(34)Φj2=exp⁡-i3τ8VxI2-A1xσ1Φj1=exp⁡-i3τ8GxΦj1=Bexp⁡-i3τ8Λ1B∗Φj1,where Gx=VxI2-A1xσ1=BΛ1B∗,(35)
G
x
=
V
x
-
A
1
x
-
A
1
x
V
x
,
B
=
1
2
1
2
-
1
2
1
2
,
Λ
1
=
V
x
-
A
1
x
0
0
V
x
+
A
1
x
.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)(36)Φj4=exp⁡-iτ4VxI2-A1xσ1+τ212ε2σ3A1x2Φj3=exp⁡-iτ48ε2HxΦj3=Pexp⁡-iτ48ε2Λ2P∗Φj3,where (37)Hx=PΛ2P∗,Λ2=12ε2Vx+A1xδ0012ε2Vx-A1xδ,P=12ε22δδ-τ2A1x12ε22δδ+τ2A1xτ2A1x-δ2δδ-τ2A1xτ2A1x+δ2δδ+τ2A1x,δ=τ4A1x2+144ε4.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 attn,tn+1 can be obtained as follows:(38a)Φj1=1M∑l=-M/2M/2-1Qlexp⁡-iτ16DlQl∗ψlnexp⁡iμlxj+a,(38b)Φj2=Bexp⁡-i3τ8Λ1B∗Φj1,(38c)Φj3=1M∑l=-M/2M/2-1Qlexp⁡-iτ13DlQl∗ψl2exp⁡iμlxj+a,(38d)Φj4=Pexp⁡-iτ48ε2Λ2P∗Φj3,(38e)Φj5=1M∑l=-M/2M/2-1Qlexp⁡-iτ13DlQl∗ψl4exp⁡iμlxj+a,(38f)Φj6=Bexp⁡-i3τ8Λ1B∗Φj5,(38g)Φjn+1=1M∑l=-M/2M/2-1Qlexp⁡-iτ16DlQl∗ψl6exp⁡iμlxj+a,where j=0,1,2,…,M-1, l=-M/2,-M/2+1,…,M/2-1, and(39)ψlk=∑j=0M-1Φjkexp⁡-iμlxj+a,k=0,1,…,7.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,(40)Φnl22≔h∑j=0M-1Φjn2≡h∑j=0M-1Φj02=Φ0l22.Proof.
Introduce the definition of the discrete inner product; that is,(41)u,v=u∗v=u1v1+u2v2,where u=u1,u2 and v=v1,v2. Then we define the l2 norm as follows:(42)u2=u,u.Equation (38a) is a Fourier transformation. From Parseval’s equality combining with (30), one obtains (43)Φj1l2=Φjnl2.For (38b), we have(44)Φj2l2=Φj2,Φj2=Φj2∗Φj2=Bexp⁡-i3τ8Λ1B∗Φj1∗Bexp⁡-i3τ8Λ1B∗Φj1=Φj1∗Bexp⁡-i3τ8Λ1∗B∗Bexp⁡-i3τ8Λ1B∗Φj1=Φj1∗Bexp⁡i3τ8Λ1exp⁡-i3τ8Λ1B∗Φj1=Φj1∗Φj1=Φj1l2.Similarly, the proof of (38c), (38e), and (38g) and (38d) and (38f) is obtained as the proof of (38a) and (38b), respectively. That is,(45)Φj2l2=Φj1l2,Φj3l2=Φj2l2,Φj4l2=Φj3l2,Φj5l2=Φj4l2,Φj6l2=Φj5l2,Φjn+1l2=Φj6l2.So, one can obtain the conservation of discrete charge(46)Φjn+1l2=Φjnl2.


## 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(47)Vx=1+x21+x2,A1x=1-x1+x2,x∈Ωand the initial condition as(48)Φ0x=ϕ1,ϕ2T=exp⁡-x22,exp⁡-x-122T,x∈Ω.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 withh=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.155From Table1, 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 is3/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-orderOτ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 forh=0.0156 and τ=1/160,1/80,1/40,1/20, and 1/10.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 Theorem1, we make the conclusion that the time-compact scheme is unconditionally stable.Whenε decreases, high frequency oscillation increases. The time evolution of Reϕ1t,0 for different ε is shown in Figure 2.Figure 2
Oscillation of Dirac equation (1) with different ε=0.4, 0.2, and 0.1.


## 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ε. 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.



---

*Source: 3670139-2016-10-16.xml*