Super-resolution of time-splitting methods for the Dirac equation in the nonrelativistic limit regime

Super-resolution of time-splitting methods for the Dirac equation in the nonrelativistic limit regime

Weizhu Bao Department of Mathematics, National University of Singapore, Singapore 119076 (matbaowz@nus.edu.sg    Yongyong Cai Beijing Computational Science Research Center, Beijing 100193, P. R. China (yongyong.cai@csrc.ac.cn)    Jia Yin NUS Graduate School for Integrative Sciences and Engineering (NGS), National University of Singapore, Singapore 117456 (e0005518@u.nus.edu)
Abstract

We establish error bounds of the Lie-Trotter splitting () and Strang splitting () for the Dirac equation in the nonrelativistic limit regime in the absence of external magnetic potentials, with a small parameter inversely proportional to the speed of light. In this limit regime, the solution propagates waves with wavelength in time. Surprisingly, we find out that the splitting methods exhibit super-resolution, in the sense of breaking the resolution constraint under the Shannon’s sampling theorem, i.e. the methods can capture the solutions accurately even if the time step size is much larger than the sampled wavelength at . shows order convergence uniformly with respect to , by establishing that there are two independent error bounds and . Moreover, if is non-resonant, i.e. is away from certain region determined by , would yield an improved uniform first order error bound. In addition, we show is uniformly convergent with 1/2 order rate for general time step size and uniformly convergent with order rate for non-resonant time step size. Finally, numerical examples are reported to validate our findings.

Key words. Dirac equation, super-resolution, nonrelativistic limit regime, time-splitting, uniform error bound

1 Introduction

The splitting technique introduced by Trotter in 1959 [46] has been widely applied in analysis and numerical simulation [2, 9, 10, 19, 20], especially in computational quantum physics. In the Hamiltonian system and general ordinary differential equations (ODEs), the splitting approach has been shown to preserve the structural/geometric properties [31, 47] and are superior in many applications. Developments of splitting type methods in solving partial differential equations (PDEs) include utilization in Schrödinger/nonlinear Schrödinger equations [2, 9, 10, 19, 20, 37, 45], Dirac/nonlinear Dirac equations [7, 8, 14, 36], Maxwell-Dirac system [11, 32], Zakharov system [12, 13, 28, 34, 35], Stokes equation [18], and Enrenfest dynamics [25], etc.

When dealing with oscillatory problems, the splitting method usually performs much better than traditional numerical methods [31, 9]. For instance, in order to obtain “correct” observables of the Schrödinger equation in the semiclassical limit regime, the time-splitting spectral method requires much weaker constraints on time step size and mesh size than the finite difference methods [9]. Similar properties have been observed for the nonlinear Schrödinger equation (NLSE)/Gross-Pitaevskii equation (GPE) in the semiclassical limit regime [2] and the Enrenfest dynamics [25]. However, in general, splitting methods still suffer from the mesh size/time step constraints related to the high frequencies in the aforementioned problems, i.e. they need to obey the resolution constraint determined by the Shannon’s sampling theorem [43] – in order to resolve a wave one needs to use a few grid points per wavelength. In this paper, we report a surprising finding that the splitting methods are uniformly accurate (w.r.t. the rapid oscillations), when applied to the Dirac equation in the nonrelativistic limit regime without external magnetic field. This fact reveals that there is no mesh size/time step restriction for splitting methods in this situation, e.g. the splitting methods have super-resolution, which is highly nontrivial. In the rest of the paper, we will discuss the oscillatory Dirac equation in the nonrelativistic limit regime, with conventional time splitting numerical approach and its super-resolution properties.

Proposed by British physicist Paul Dirac in 1928 [23], the Dirac equation has now been extensively applied in the study of the structures and/or dynamical properties of graphene, graphite, and other two-dimensional (2D) materials [1, 26, 39, 40], as well as the relativistic effects of molecules in super intense lasers, e.g., attosecond lasers [16, 27]. Mathematically, the -dimensional () Dirac equation with external electro-magnetic potentials [7, 14] for the complex spinor vector field can be written as

 i∂tΨ=(−iεd∑j=1αj∂j+1ε2β)Ψ+(V(t,x)I4−d∑j=1Aj(t,x)αj)Ψ,x∈Rd,t>0, \hb@xt@.01(1.1)

with initial value

 Ψ(t=0,x)=Ψ0(x),x∈Rd, \hb@xt@.01(1.2)

where , is time, is the spatial coordinate vector, (), and () are the given real-valued electric and magnetic potentials, respectively, is a dimensionless parameter inversely proportional to the speed of light. There are two important regimes for the Dirac equation (LABEL:eq:dirac4): the relativistic case (wave speed is comparable to the speed of light) and the nonrelativistic limit case (wave speed is much less than the speed of light). is the identity matrix for , and the matrices , , and are

 α1=(0σ1σ10),α2=(0σ2σ20),α3=(0σ3σ30),β=(I200−I2), \hb@xt@.01(1.3)

where , , are the Pauli matrices

 σ1=(0110),σ2=(0−ii0),σ3=(100−1). \hb@xt@.01(1.4)

In the relativistic regime , extensive analytical and numerical studies have been carried out for the Dirac equation (LABEL:eq:dirac4) in the literature. In the analytical aspect, for the existence and multiplicity of bound states and/or standing wave solutions, we refer to [21, 22, 24, 29, 30, 42] and references therein. In the numerical aspect, many accurate and efficient numerical methods have been proposed and analyzed [3, 38], such as the finite difference time domain (FDTD) methods [4, 41], time-splitting Fourier pseudospectral (TSFP) method [7, 32], exponential wave integrator Fourier pseudospectral (EWI-FP) method [7], and the Gaussian beam method [48], etc.

In the nonrelativistic limit regime, as , the Dirac equation (LABEL:eq:dirac4) converges to Pauli equation [15, 33] or Schrödinger equation [5, 15], and the solution propagates waves with wavelength O() in time and O(1) in space, respectively. The highly oscillatory nature of the solution in time brings severe difficulties in numerical computation in the nonrelativistic limit regime, i.e. when . In fact, it would cause the time step size strictly dependent on in order to capture the solution accurately. Rigorous error estimates were established for the finite difference time domain method (FDTD), exponential wave integrator Fourier pseudospectral method (EWI-FP) and time-splitting Fourier pseudospectral method (TSFP) in this parameter regime [7]. The error bounds suggested for FDTD and for EWI-FP and TSFP. A new fourth-order compact time-splitting method () was recently put forward to improve the efficiency and accuracy [14]. Moreover, a uniformly accurate multiscale time integrator pseudospectral method was proposed and analyzed for the Dirac equation in the nonrelativistic limit regime, where the errors are uniform with respect to [6], allowing for -independent time step .

From the analysis in [7], the error bounds for second order Strang splitting TSFP (also called as later in this paper) depends on the small parameter as . Surprisingly, through our extensive numerical experiments, we find out that if the magnetic potentials for in (LABEL:eq:dirac4), the errors of TSFP are then independent of and uniform w.r.t. , i.e., for Dirac equation (LABEL:eq:dirac4) without magnetic potentials has super-resolution w.r.t. . In such case, (LABEL:eq:dirac4) reduces to

 i∂tΨ(t,x)=(−iεd∑j=1αj∂j+1ε2β+V(t,x)I4)Ψ(t,x),x∈Rd,d=1,2,3,t>0, \hb@xt@.01(1.5)

with the initial value given in (LABEL:eq:initial). In lower dimensions (), the four component Dirac equation (LABEL:eq:dirac4_nomag) can be reduced to the following two-component form for [7]:

 i∂tΦ(t,x)=(−iεd∑j=1σj∂j+1ε2σ3+V(t,x)I2)Φ(t,x),x∈Rd,d=1,2,t>0, \hb@xt@.01(1.6)

with initial value

 Φ(t=0,x)=Φ0(x),x∈Rd,d=1,2. \hb@xt@.01(1.7)

The two component form (LABEL:eq:dirac2_nomag) is widely used in lower dimensions due to its simplicity compared to the four component form (LABEL:eq:dirac4_nomag).

Our extensive numerical studies and theoretical analysis show that for first-order, second-order, and even higher order time-splitting Fourier pseudospectral methods, there are always uniform error bounds w.r.t. . In other words, the splitting methods can capture the solutions accurately even if the time step size is much larger than the sampled wavelength at , i.e. they exhibit super-resolution in the sense of breaking the resolution constraint under the Shannon’s sampling theorem [43]. This super-resolution property of the splitting methods makes them more efficient and reliable for solving the Dirac equation without magnetic potentials in the nonrelativisitc limit regime, compared to other numerical approaches in the literature. In the sequel, we will study rigorously the super-resolution phenomenon for first-order () and second-order () time-splitting methods, and present numerical results to validate the conclusions.

The rest of the paper is organized as follows. In section 2, we review the first and second order time-splitting methods for the Dirac equation in the nonrelativistic limit regime without magnetic potential, and state the main results. In section 3 and section 4 respectively, detailed proofs for the uniform error bounds and improved uniform error bounds are presented. Section 5 is devoted to numerical tests, and finally, some concluding remarks are drawn in section 6. Throughout the paper, we adopt the standard Sobolev spaces and the corresponding norms. Meanwhile, is used with the meaning that there exists a generic constant independent of and , such that . has a similar meaning that there exists a constant dependent on but independent of and , such that .

2 Time-splitting methods and main results

In this section, we recall the first and second order time-splitting methods applied to the Dirac equation and state the main results of this paper. For simplicity of presentation, we only carry out the splitting methods and corresponding analysis for (LABEL:eq:dirac2_nomag) in 1D (). Generalization to (LABEL:eq:dirac4_nomag) and/or higher dimensions is straightforward and results remain valid without modifications.

2.1 Time-splitting methods

Denote the Hermitian operator

 Tε=−iεσ1∂x+σ3,x∈R, \hb@xt@.01(2.1)

and the Dirac equation (LABEL:eq:dirac2_nomag) in 1D can be written as

 i∂tΦ(t,x)=1ε2TεΦ(t,x)+V(t,x)Φ(t,x),x∈R, \hb@xt@.01(2.2)

with initial value

 Φ(0,x)=Φ0(x),x∈R. \hb@xt@.01(2.3)

Choose to be the time step size and for as the time steps. Denote as the numerical approximation of , where is the exact solution to (LABEL:eq:dirac_nomag) with (LABEL:eq:dirac_initial), then the first-order and second-order time-splitting methods can be expressed as follows.

First-order splitting (Lie-Trotter splitting). The discrete-in-time first-order splitting () is written as [46]

 Φn+1(x)=e−iτε2Tεe−i∫tn+1tnV(s,x)dsΦn(x),with Φ0(x)=Φ0(x),x∈R. \hb@xt@.01(2.4)

Second-order splitting (Strang splitting). The discrete-in-time second-order splitting () is written as [44]

 Φn+1(x)=e−iτ2ε2Tεe−i∫tn+1tnV(s,x)dse−iτ2ε2TεΦn(x),with Φ0(x)=Φ0(x),x∈R. \hb@xt@.01(2.5)

Then the main results of this paper can be summarized below.

2.2 Uniform error bounds

For any , we are going to consider smooth enough solutions, i.e. we assume the electric potential satisfies

 (A)V(t,x)∈Wm,∞([0,T];L∞(R))∩L∞([0,T];W2m+m∗,∞(R)),m∈N∗,m∗∈{0,1}.

In addition, we assume the exact solution satisfies

 (B)Φ(t,x)∈L∞([0,T],(H2m+m∗(R))2),m∈N∗,m∗∈{0,1}.

We remark here that if the initial value , then condition is implied by condition .

For the numerical approximation obtained from (LABEL:eq:Lie-trotter) or (LABEL:eq:Strang), we introduce the error function

 en(x)=Φ(tn,x)−Φn(x),0≤n≤Tτ, \hb@xt@.01(2.6)

then the following error estimates hold.

Theorem 2.1

Let be the numerical approximation obtained from (LABEL:eq:Lie-trotter), then under the assumptions and with and , we have the following error estimates

 ∥en(x)∥L2≲τ+ε,∥en(x)∥L2≲τ+τ/ε,0≤n≤Tτ. \hb@xt@.01(2.7)

As a result, there is a uniform error bound for

 ∥en(x)∥L2≲τ+max0<ε≤1min{ε,τ/ε}≲√τ,0≤n≤Tτ. \hb@xt@.01(2.8)
Theorem 2.2

Let be the numerical approximation obtained from (LABEL:eq:Strang), then under the assumptions and with and , we have the following error estimates

 ∥en(x)∥L2≲τ2+ε,∥en(x)∥L2≲τ2+τ2/ε3,0≤n≤Tτ. \hb@xt@.01(2.9)

As a result, there is a uniform error bound for

 ∥en(x)∥L2≲τ2+max0<ε≤1min{ε,τ2/ε3}≲√τ,0≤n≤Tτ. \hb@xt@.01(2.10)
Remark 2.1

The error bounds in Theorem LABEL:thm:lie can be expressed as

 ∥en(x)∥L2≤C1∥Φ0(x)∥H2(τ+max0<ε≤1min{ε,τ/ε}),0≤n≤Tτ, \hb@xt@.01(2.11)

and the error estimates in Theorem LABEL:thm:strang can be restated as

 ∥en(x)∥L2≤C2∥Φ0(x)∥H4(τ2+max0<ε≤1min{ε,τ2/ε3}),0≤n≤Tτ, \hb@xt@.01(2.12)

where and are constants depending only on and .

We note that higher order time-splitting methods also have similar super-resolution property, but for simplicity, we only focus on and here. Remark LABEL:remark:err could be easily derived by examining the proofs of Theorems LABEL:thm:lie & LABEL:thm:strang, and the details will be skipped.

2.3 Improved uniform error bounds for non-resonant time steps

In the Dirac equation (LABEL:eq:dirac2_nomag) or (LABEL:eq:dirac4_nomag), the leading term is or , which suggests the solution exhibits almost periodicity in time with periods (, the periods of and ). From numerical results, we observe the errors behave much better compared to the results in Theorems LABEL:thm:lie& LABEL:thm:strang, when is away from the leading temporal oscillation periods . In fact, for given , define

 Aδ(ε):=∞⋃k=0[ε2kπ+ε2arcsinδ,ε2(k+1)π−ε2arcsinδ],0<ε≤1, \hb@xt@.01(2.13)

and the errors of and can be improved compared to the previous subsection when . To illustrate , we show in Figure LABEL:fig:axis for and with fixed .

For , we can derive improved uniform error bounds for the two splitting methods as shown in the following two theorems.

Theorem 2.3

Let be the numerical approximation obtained from (LABEL:eq:Lie-trotter). If the time step size is non-resonant, i.e. there exists , such that , under the assumptions and with and , we have an improved uniform error bound

 ∥en(x)∥L2≲δτ,0≤n≤Tτ. \hb@xt@.01(2.14)
Theorem 2.4

Let be the numerical approximation obtained from (LABEL:eq:Strang). If the time step size is non-resonant, i.e. there exists , such that , under the assumptions and with and , we assume an extra regulariry and then the following two error estimates hold

 ∥en(x)∥L2≲δτ2+τε,∥en(x)∥L2≲δτ2+τ2/ε,0≤n≤Tτ. \hb@xt@.01(2.15)

As a result, there is an improved uniform error bound for

 ∥en(x)∥L2≲δτ2+max0<ε≤1min{τε,τ2/ε}≲δτ3/2,0≤n≤Tτ. \hb@xt@.01(2.16)
Remark 2.2

In Theorems LABEL:thm:lie2 and LABEL:thm:strang2, the constants in the error estimates depend on and the proof in the paper suggests that the constants are bounded from above by and with some common factor independent of and . The optimality of the uniform error bounds in Theorems LABEL:thm:lie2 and LABEL:thm:strang2 will be verified by numerical examples presented in section 5.

3 Proof of Theorems LABEL:thm:lie and LABEL:thm:strang

In this section, we prove the uniform error bounds for the splitting methods and . As is diagonalizable in the phase space (Fourier domain), it can be decomposed as [6, 7, 15]

 Tε=√Id−ε2ΔΠε+−√Id−ε2ΔΠε−, \hb@xt@.01(3.1)

where is the Laplace operator in 1D and is the identity operator. and are projectors defined as

 Πε+=12[I2+(Id−ε2Δ)−1/2Tε],Πε−=12[I2−(Id−ε2Δ)−1/2Tε]. \hb@xt@.01(3.2)

It is straightforward to see that , and , . Furthermore, through Taylor expansion, we have [15]

 Πε+=Π0++εR1=Π0+−iε2σ1∂x+ε2R2,Π0+=diag(1,0), \hb@xt@.01(3.3) Πε−=Π0−−εR1=Π0−+iε2σ1∂x−ε2R2,Π0−=diag(0,1), \hb@xt@.01(3.4)

where for , , and for , are uniformly bounded operators with respect to .

To help capture the features of solutions, denote

 Dε=1ε2(√Id−ε2Δ−Id)=−(√Id−ε2Δ+Id)−1Δ, \hb@xt@.01(3.5)

where is a uniformly bounded operator with respect to from to for , then we have the decomposition for the unitary evolution operator as

 eitε2Tε=eitε2(√Id−ε2ΔΠε+−√Id−ε2ΔΠε−)=eit/ε2eitDεΠε++e−it/ε2e−itDεΠε−. \hb@xt@.01(3.6)

For the ease of the proof, we first introduce the following two lemmas for the Lie-Trotter splitting (LABEL:eq:Lie-trotter) and the Strang splitting (LABEL:eq:Strang), respectively. For simplicity, we denote , and in short.

Lemma 3.1

Let be the numerical approximation obtained from the Lie-Trotter splitting (LABEL:eq:Lie-trotter), then under the assumptions and with and , we have

 en+1(x)=e−iτε2Tεe−i∫tn+1tnV(s,x)dsen(x)+ηn1(x)+ηn2(x),0≤n≤Tτ−1, \hb@xt@.01(3.7)

with , , where

 fn2(s)=ei2s/ε2eisDεΠε+(V(tn)Πε−eisDεΦ(tn))+e−i2s/ε2e−isDεΠε−(V(tn)Πε+e−isDεΦ(tn)). \hb@xt@.01(3.8)

Proof. From the definition of , noticing the Lie-Trotter splitting formula (LABEL:eq:Lie-trotter), we have

 en+1(x)=e−iτε2Tεe−i∫tn+1tnV(s,x)dsen(x)+ηn(x),0≤n≤Tτ−1,x∈R, \hb@xt@.01(3.9)

where is the local truncation error defined as

 ηn(x)=Φ(tn+1,x)−e−iτε2Tεe−i∫tn+1tnV(s,x)dsΦ(tn,x),x∈R. \hb@xt@.01(3.10)

Noticing (LABEL:eq:dirac_nomag), applying Duhamel’s principle, we derive

 Φ(tn+1,x)=e−iτε2TεΦ(tn,x)−i∫τ0e−i(τ−s)ε2TεV(tn+s,x)Φ(tn+s,x)ds, \hb@xt@.01(3.11)

while Taylor expansion gives

 e−iτε2Tεe−i∫tn+1tnV(s,x)dsΦ(tn,x)= \hb@xt@.01(3.12)

Combining (LABEL:eq:Duhamel), (LABEL:eq:tl:s1) and (LABEL:eq:loc1), we get

 ηn(x)= τie−iτε2TεV(tn,x)Φ(tn,x)−i∫τ0e−i(τ−s)ε2Tε(V(tn,x)e−isε2TεΦ(tn,x))ds+2∑j=1Rnj(x), \hb@xt@.01(3.13)

where

 Rn1(x) Rn2(x) =−i∫τ0e−i(τ−s)ε2Tε(V(tn)λn4(s,x))ds−i∫τ0e−i(τ−s)ε2Tε(λn3(s,x)Φ(tn+s,x))ds,

with

 λn1(x)=e−i∫tn+1tnV(s,x)ds−(1−i∫tn+1tnV(s,x)ds), \hb@xt@.01(3.14) λn2(x)=−i∫tn+1tnV(s,x)ds+iτV(tn,x),λn3(s,x)=V(tn+s,x)−V(tn,x),0≤s≤τ, \hb@xt@.01(3.15) λn4(s,x)=−i∫s0e−i(s−w)ε2Tε(V(tn+w,x)Φ(tn+w,x))dw,0≤s≤τ. \hb@xt@.01(3.16)

It is easy to see that for ,

 ∥λn1(x)∥L∞≲τ2∥V(t,x)∥2L∞(L∞),∥λn4(s,x)∥L∞([0,τ];(L2)2)≲τ∥V(t,x)∥L∞(L∞)∥Φ(t,x)∥L∞((L2)2), ∥λn2(x)∥L∞≲τ2∥∂tV(t,x)∥L∞(L∞),∥λn3(s,x)∥L∞([0,τ];L∞)≲τ∥∂tV(t,x)∥L∞(L∞).

As a consequence, we obtain the following bounds for ,

 ∥Rn1(x)∥L2 ≲(∥λn1(x)∥L∞+∥λn2(x)∥L∞)∥Φ(tn)∥L2≲τ2, \hb@xt@.01(3.17) ∥Rn2(x)∥L2 ≲τ(∥V(tn)∥L∞∥λn4(s,x)∥L∞([0,τ];(L2)2)+∥λn3(s,x)∥L∞([0,τ];L∞)∥Φ∥L∞((L2)2)) \hb@xt@.01(3.18) ≲τ2.

Recalling given in Lemma LABEL:lemma:lie, we introduce

 fn(s):=fn(s,x)=eisε2Tε(V(tn,x)e−isε2TεΦ(tn,x))=fn1(s)+fn2(s),0≤s≤τ, \hb@xt@.01(3.19)

with given in (LABEL:eq:f2nlie) and from the decomposition (LABEL:eq:dec:ev) as

 fn1(s) =eisDεΠε+(V(tn)e−isDεΠε+Φ(tn))+e−isDεΠε−(V(tn)eisDεΠε−Φ(tn)),

and then (LABEL:eq:S1_eta) can be written as

 \hb@xt@.01(3.20)

Now, it is easy to verify that with given in Lemma LABEL:lemma:lie if we let

 ηn1(x)=−ie−iτε2Tε(∫τ0fn1(s)ds−τfn1(0))+Rn1(x)+Rn2(x). \hb@xt@.01(3.21)

Noticing that

 ∥∥∥e−iτε2Tε(∫τ0fn1(s)ds−τfn1(0))∥∥∥L2≲ τ2∥∂sfn1(⋅)∥L∞([0,τ];(L2)2) ≲ τ2(∥V(tn)∥L∞∥Φ(tn)∥H2+∥V(tn)∥W2,∞∥Φ(tn)∥H2),

recalling the regularity assumptions and , combining (LABEL:eq:S1_R1) and (LABEL:eq:S1_R3) , we can get

 ∥ηn1(x)∥L2≤∥Rn1(x)∥L2+∥Rn2(x)∥L2+∥∥∥e−iτε2Tε(∫τ0fn1(s)ds−τfn1(0))∥∥∥L2≲τ2,

which completes the proof of Lemma LABEL:lemma:lie.

Lemma 3.2

Let be the numerical approximation obtained from the Strang splitting (LABEL:eq:Strang), then under the assumptions and with and , we have

 en+1(x)=e−iτ2ε2Tεe−i∫tn+1tnV(s,x)dse−iτ2ε2Tεen(x)+ηn1(x)+ηn2(x)+ηn3(x),0≤n≤Tτ−1, \hb@xt@.01(3.22)

with

 ∥ηn1(x)∥L2≲τ3,ηn2(x)=−ie−iτε2Tε(∫τ0fn2(s)ds−τfn2(τ/2)), \hb@xt@.01(3.23) ηn3(x)=−e−iτε2Tε(∫τ0∫s04∑j=2gnj(s,w)dwds−τ224∑j=2gnj(τ/2,τ/2)), \hb@xt@.01(3.24)

where

 fn2(s) =ei2sε2eisDεΠε+(V(tn+s)eisDεΠε−Φ(tn))+e−i2sε2e−isDεΠε−(V(tn+s)e−isDεΠε+Φ(tn)), \hb@xt@.01(3.25) gn2(s,w) =ei2w/ε2eisDεΠε+(V(tn)e−i(s−w)DεΠε+(V(tn)eiwDεΠε−Φ(tn))) +e−i2w/ε2e−isDεΠε−(