# Uniform error bounds of time-splitting methods for the nonlinear Dirac equation in the nonrelativistic limit regime

Super-resolution of the Lie-Trotter splitting (S_1) and Strang splitting (S_2) is rigorously analyzed for the nonlinear Dirac equation without external magnetic potentials in the nonrelativistic limit regime with a small parameter 0<ε≤ 1 inversely proportional to the speed of light. In this limit regime, the solution highly oscillates in time with wavelength at O(ε^2) in time. The splitting methods surprisingly show super-resolution, in the sense of breaking the resolution constraint under the Shannon's sampling theorem, i.e. the methods can capture the solution accurately even if the time step size τ is much larger than the sampled wavelength at O(ε^2). Similar to the linear case, S_1 and S_2 both exhibit 1/2 order convergence uniformly with respect to ε. Moreover, if τ is non-resonant, i.e. τ is away from certain region determined by ε, S_1 would yield an improved uniform first order O(τ) error bound, while S_2 would give improved uniform 3/2 order convergence. Numerical results are reported to confirm these rigorous results. Furthermore, we note that super-resolution is still valid for higher order splitting methods.

## Authors

• 17 publications
• 6 publications
• 7 publications
05/21/2021

### Spatial resolution of different discretizations over long-time for the Dirac equation with small potentials

We compare the long-time error bounds and spatial resolution of finite d...
12/07/2021

### Improved uniform error bounds on time-splitting methods for the long-time dynamics of the Dirac equation with small potentials

We establish improved uniform error bounds on time-splitting methods for...
09/18/2021

### Improved uniform error bounds for the time-splitting methods for the long-time dynamics of the Schrödinger/nonlinear Schrödinger equation

We establish improved uniform error bounds for the time-splitting method...
01/29/2020

### Uniform error bounds of a time-splitting spectral method for the long-time dynamics of the nonlinear Klein-Gordon equation with weak nonlinearity

We establish uniform error bounds of a time-splitting Fourier pseudospec...
09/22/2021

### Error bounds of fourth-order compact finite difference methods for the Dirac equation in the massless and nonrelativistic regime

We establish the error bounds of fourth-order compact finite difference ...
10/27/2021

### Scattering and uniform in time error estimates for splitting method in NLS

We consider the nonlinear Schrödinger equation with a defocusing nonline...
12/17/2021

### An adaptive splitting method for the Cox-Ingersoll-Ross process

We propose a new splitting method for strong numerical solution of the C...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The splitting methods form an important group of methods which are quite accurate and efficient [57]. Actually, they have been widely applied for dealing with highly oscillatory systems such as the Schrödinger/nonlinear Schrödinger equations [1, 8, 9, 22, 23, 55, 67], the Dirac/nonlinear Dirac equations [5, 6, 14, 54], the Maxwell-Dirac system [10, 49], the Zakharov system [12, 13, 41, 50], the Gross-Pitaevskii equation for Bose-Einstein condensation (BEC) [11], the Stokes equation [21], and the Enrenfest dynamics [32], etc.

In this paper, we consider the splitting methods applied to the nonlinear Dirac equation (NLDE) [27, 28, 37, 40, 47, 25, 43, 44, 33, 34, 35, 36, 61, 63, 70] in the nonrelativistic limit regime without magnetic potential. In one or two dimensions (1D or 2D), the equation can be represented in the two-component form with wave function [6]:

 (1.1) i∂tΦ=(−iεd∑j=1σj∂j+1ε2σ3)Φ+V(x)Φ+F(Φ)Φ,x∈Rd,d=1,2,t>0,

where is the imaginary unit, is time, , (), is a dimensionless parameter inversely proportional to the speed of light, and is a real-valued function denoting the external electric potential. , , are the Pauli matrices defined as

 (1.2) σ1=(0110),σ2=(0−ii0),σ3=(100−1).

The nonlinearity in (1.1) is usually taken as

 (1.3) F(Φ)=λ1(Φ∗σ3Φ)σ3+λ2|Φ|2I2,

where , are two given real constants, is the complex conjugate transpose of and is the identity matrix. The above choice of nonlinearity is motivated from the so-called Soler model in quantum field theory, e.g. and [37, 40, 68], and BEC with a chiral confinement and/or spin-orbit coupling, e.g. and [25, 43, 44]. In order to study the dynamics, the initial data is chosen as

 (1.4) Φ(t=0,x)=Φ0(x),x∈Rd,d=1,2.

When in (1.1), which corresponds to the classical regime of the nonlinear Dirac equation, there have been comprehensive analytical and numerical results in the literatures. In the analytical aspect, for the existence and multiplicity of bound states and/or standing wave solutions, we refer to [2, 3, 15, 24, 29, 30, 31, 51] and references therein. Particularly, for the case where , , , and in the choice of , the NLDE (1.1) admits explicit soliton solutions [26, 40, 45, 52, 56, 60, 65, 66]. In the numerical aspect, many accurate and efficient numerical methods have been proposed and analyzed, such as the finite difference time domain (FDTD) methods [19, 46, 59], the time-splitting Fourier spectral (TSFP) methods [10, 18, 39, 49] and the Runge-Kutta discontinuous Galerkin methods [48].

On the other hand, when (the nonrelativistic limit regime where the wave speed is much smaller than the speed of light), as indicated by previous analysis in [6, 38, 58, 20], the wavelength of the solution in time is at . The oscillation of the solution as well as the unbounded and indefinite energy functional w.r.t. [16, 31] cause much burden in the analysis and computation. Indeed, it would require that the time step size to be strictly reliant on to capture the exact solution, as suggested by the Shannon’s sampling theorem [62]. Numerical studies in [6] have confirmed this dependence. The error bounds show that is required for the conservative Crank-Nicolson finite difference (CNFD) method, and is required for the exponential wave integrator Fourier pseudospectral (EWI-FP) method as well as the time-splitting Fourier pseudospectral (TSFP) method. To overcome the restriction, recently, uniform accurate (UA) schemes with two-scale formulation approach [53] or multiscale time integrator pseudospectral method [20] have been designed for the NLDE in the nonrelativistic limit regime, where the time step size could be independent of .

Though the TSFP method (also called later in this paper) has a dependence on the small parameter [6], under the specific case where there is a lack of magnetic potential, as in (1.1), we find out through our recent extensive numerical experiments that the errors of is independent of and uniform w.r.t. . In other words, for the NLDE (1.1) in the absence of magnetic potentials displays super-resolution w.r.t. .

The super-resolution property for the time-splitting methods makes them superior in solving the NLDE in the absence of magnetic potentials in the nonrelativistic regime as they are more efficient and reliable as well as simple compared to other numerical methods in the literature. In this paper, the super-resolution for the first-order () and second-order () time-splitting methods will be rigorously analyzed, and numerical results will be presented to validate the conclusions. We remark that similar results have been analyzed for the Dirac equation [7]

, where the linearity enables us to explicitly track the error exactly and make estimation at the target time step without using Gronwall type arguments. However, in the nonlinear case, it is impossible to follow the error propagation exactly and estimations have to be done at each time step. As a result, Gronwall arguments will be involved together with the mathematical induction to control the nonlinearity and to bound the numerical solution. In particular, instead of the previously adopted Lie calculus approach

[55], Taylor expansion and Duhamel principle are employed to study the local error of the splitting methods, which can identify how temporal oscillations propagate numerically. In other words, the techniques adapted to establish uniform error bounds of the time-splitting methods for the NLDE are completely different with those used for the Dirac equation [7].

The rest of the paper is organized as follows. In section 2, we establish uniform error estimates of the first-order time-splitting method for the NLDE without magnetic potentials in the nonrelativistic limit regime and report numerical results to confirm our uniform error bounds. Similar results are presented for the second-order time-splitting method in section 3 with a remark on extension to higher order splitting methods. Some conclusions are drawn in section 4. Throughout the paper, we adopt the standard Sobolev spaces and the corresponding norms. Meanwhile, is used in the sense that there exists a generic constant independent of and , such that . has a similar meaning that there exists a generic constant dependent on but independent of and , such that .

## 2 Uniform error bounds of the first-order Lie-Trotter splitting method

For simplicity of notations and without loss of generality, here we only consider (1.1) in 1D (). Extensions to (1.1) in 2D and/or the four component form of the NLDE with [6] are straightforward.

Denote the free Dirac Hermitian operator

 (2.1) Tε=−iεσ1∂x+σ3,x∈R,

then the NLDE (1.1) in 1D can be written as

 (2.2) i∂tΦ(t,x)=1ε2TεΦ(t,x)+V(x)Φ(t,x)+F(Φ(t,x))Φ(t,x),x∈R,

with nonlinearity (1.3) and the initial condition (1.4).

Choose as the time step size and for as the time steps. Denote to be the numerical approximation of , where is the exact solution of (2.2) with (1.3) and (1.4), then through applying the discrete-in-time first-order splitting (Lie-Trotter splitting) [69], can be represented as:

 (2.3) Φn+1(x)=e−iτε2Tεe−iτ[V(x)+F(Φn(x))]Φn(x),with Φ0(x)=Φ0(x),x∈R.

For simplicity, we also write , where denotes the numerical propagator of the Lie-Trotter splitting.

### 2.1 A uniform error bound

For any , where denotes the common maximal existence time of the solution for (1.1) with (1.4) for all , we are going to consider smooth solutions, i.e. we assume the electric potential satisfies

 (A)V(x)∈W2m+1,∞(R),m∈N∗.

In addition, we assume the exact solution satisfies

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

For the numerical approximation obtained from (2.3), we introduce the error function

 (2.4) en(x)=Φ(tn,x)−Φn(x),0≤n≤Tτ,

then the following uniform error bound can be established. Let be the numerical approximation obtained from (2.3), then under assumptions and with , there exists independent of such that the following two error estimates hold for

 (2.5) ∥en(x)∥H1≲τ+ε,∥en(x)∥H1≲τ+τ/ε,0≤n≤Tτ.

Consequently, there is a uniform error bound for when

 (2.6) ∥en(x)∥H1≲τ+max0<ε≤1min{ε,τ/ε}≲√τ,0≤n≤Tτ.
###### Remark 2.1

In Theorem 2.1 and the other results in this paper for the 1D case, we prove the error bounds for due to the fact that is an algebra, and the corresponding estimates should be in norm for 2D and 3D cases with of course higher regularity assumptions.

For simplicity of the presentation, in the proof for this theorem and other theorems later for NLDE in this paper, we take . Extension to the case where is straightforward [7]. Compared to the linear case [7], the nonlinear term is much more complicated to analyze. A key issue of the error analysis for NLDE is to control the nonlinear term of numerical solution , and for which we require the following stability lemma [55]. Suppose , and satisfy , we have

 (2.7) ∥SLien,τ(Φ)−SLien,τ(Ψ)∥H1≤ec1τ∥Φ−Ψ∥H1,

where depends on and . The proof is quite similar to the nonlinear Schrödinger equation case in [55] and we omit it here for brevity.

Under the assumption (B) (), for , we denote as

 (2.8) M1=supε∈(0,1]∥Φ(t,x)∥L∞([0,T];(H1(R))2).

Based on (2.8) and Lemma 2.1, one can control the nonlinear term once the hypothesis of the lemma is fulfilled. Making use of the fact that is explicit, together with the uniform error estimates in Theorem 2.1, we can use mathematical induction to complete the proof.

The following properties of will be frequently used in the analysis. is diagonalizable in the phase space (Fourier domain) and can be decomposed as

 (2.9) Tε=√Id−ε2ΔΠε+−√Id−ε2ΔΠε−,

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

 (2.10) Πε+=12[Id+(Id−ε2Δ)−1/2Tε],Πε−=12[Id−(Id−ε2Δ)−1/2Tε].

It is straightforward to verify that , , , and through Taylor expansion, we have [16]

 (2.11) Πε±=Π0±±εR1=Π0±∓iε2σ1∂x±ε2R2,Π0+=diag(1,0),Π0−=% diag(0,1),

with for , for being uniformly bounded operators w.r.t. .
In order to characterize the oscillatory features of the solution, denote

 (2.12) Dε=1ε2(√Id−ε2Δ−Id)=−(√Id−ε2Δ+Id)−1Δ,

which is a uniformly bounded operator w.r.t from for , then the evolution operator can be expressed as

 (2.13) eitε2Tε=eitε2(√Id−ε2ΔΠε+−√Id−ε2ΔΠε−)=eitε2eitDεΠε++e−itε2e−itDεΠε−.

For simplicity, here we use , in short.

Now we are ready to introduce the following lemma for proving Theorem 2.1.

Let ( ) be obtained from (2.3) satisfying , under the assumptions of Theorem 2.1, we have

 (2.14) en+1(x)=e−iτε2Tεe−iτF(Φn)en(x)+ηn1(x)+e−iτε2Tεηn2(x),

with , , where depends on , , and ; depends on , , and . Here

 fn2(s)= −ie−4isε2Πε−(gn1(x)Πε+Φ(tn))−ie4isε2Πε+(¯¯¯¯¯¯gn1(x)Πε−Φ(tn)) −ie−i2sε2[Πε+(gn1(x)Πε+Φ(tn))+Πε−(gn2(x)Πε+Φ(tn)+gn1(x)Πε−Φ(tn))] (2.15)

where with and

 (2.16) g1(Φ+(tn),Φ−(tn))=λ1((Φ−(tn))∗σ3Φ+(tn))σ3+λ2((Φ−(tn))∗Φ+(tn))I2, (2.17) g2(Φ+(tn),Φ−(tn))=∑σ=±[λ1((Φσ(tn))∗σ3Φσ(tn))σ3+λ2|Φσ(tn)|2I2].

Through the definition of (2.4), noticing the formula (2.3), we have

 (2.18) en+1(x)=e−iτε2Tεe−iτF(Φn)en(x)+ηn(x),0≤n≤Tτ−1,x∈R,

where is the “local truncation error” (notice that this is not the usual local truncation error, compared with ),

 (2.19) ηn(x)=Φ(tn+1,x)−e−iτε2Tεe−iτF(Φn)Φ(tn,x),x∈R.

By Duhamel’s principle, the solution to (2.2) satisfies

 (2.20) Φ(tn+s,x)=e−isε2TεΦ(tn,x)−i∫s0e−i(s−w)ε2TεF(Φ(tn+w,x))Φ(tn+w,x)dw,0≤s≤τ,

which implies that (). Setting in (2.20), we have from (2.19),

 (2.21) ηn(x)=e−iτε2Tε(∫τ0fn(s)ds−τfn(0))+Rn1(x)+Rn2(x),

where

 (2.22) fn(s)=−ieisε2Tε(F(e−isε2TεΦ(tn))e−isε2TεΦ(tn,x)),Rn1(x)=e−iτε2Tε(Λn1(x)+Λn2(x)), (2.23) Rn2(x)=−i∫τ0e−i(τ−s)ε2Tε[F(Φ(tn+s))Φ(tn+s)−F(Φ(e−isε2TεΦ(tn)))e−isε2TεΦ(tn)]ds,

with

 (2.24) Λn1(x)=−(e−iτF(Φn)−(I2−iτF(Φn)))Φ(tn),Λn2(x)=(iτ(F(Φn)−F(Φ(tn))))Φ(tn).

Noticing (2.20), the assumption that , and the fact that preserves norm, it is not difficult to find

 (2.25) ∥Rn2(x)∥H1≲(M1+1)2∫τ0∥Φ(tn+s,x)−e−isε2TεΦ(tn,x)∥H1ds≲τ2.

On the other hand, using Taylor expansion in and the local Lipschitz property of , we get

 (2.26) ∥Rn1(x)∥H1 ≲τ2∥Φn∥2H1∥Φ(tn)∥H1+τ(M1+1)2∥Φn−Φ(tn)∥H1≲τ2+τ∥en(x)∥H1.

It remains to estimate the part. Using the decomposition (2.13) and the Taylor exapnsion , we have (),

 (2.27) fn(s)=−i∑σ=±eσisε2Πεσ{F(e−isε2Φ+(tn)+eisε2Φ−(tn))(e−isε2Φ+(tn)+eisε2Φ−(tn))}+fn1(s),

where for ,

 (2.28) ∥fn1(s)∥H1≲τ∥Φ(tn)∥3H3≲τ.

Since is of polynomial type, by direct computation, we can further simplify (2.27) to get

 (2.29) fn(s)=fn1(s)+fn2(s)+~fn(s),0≤s≤τ,

where is given in (2.15) and is independent of as

 (2.30) ~fn(s)≡−i[Πε+(gn2(x)Πε+Φ(tn)+gn1(x)Πε−Φ(tn))+Πε−(gn2(x)Πε−Φ(tn)+¯¯¯¯¯¯gn1(x)Πε+Φ(tn))],

with defined in (2.16)-(2.17).

Now, it is easy to verify that with given in Lemma 2.1 by choosing

 (2.31) ηn1(x)=e−iτε2Tε(∫τ0(fn1(s)+~fn(s))ds−τ(fn1(0)+~fn(0)))+Rn1(x)+Rn2(x).

Noticing that is independent of and , combining (2.25) and (2.26), we can get

 ∥ηn1(x)∥H1≤2∑j=1∥Rnj(x)∥H1+∥∥∥∫τ0fn1(s)ds−τfn1(0)∥∥∥H1≲τ∥en(x)∥H1+τ2,

which completes the proof of Lemma 2.1.

Now, we proceed to prove Theorem 2.1.  We will prove by induction that the estimates (2.5)-(2.6) hold for all time steps together with

 (2.32) ∥Φn∥H1≤M1+1.

Since initially , case is obvious. Assume (2.5)-(2.6) and (2.32) hold true for all , then we are going to prove the case .

From Lemma 2.1, we have

 (2.33) en+1(x)=e−iτε2Tεe−iτF(Φn)en(x)+ηn1(x)+e−iτε2Tεηn2(x),0≤n≤m,

with , and given in Lemma 2.1.

Denote (), and it is straightforward to calculate

 (2.34) ∥LnΨ(x)∥H1≤CM1τ∥Ψ∥H1,∀Ψ∈(H1(R))2,

with only depending on . Thus we can obtain from (2.33) that for ,

 en+1(x) =e−iτε2Tεen(x)+ηn1(x)+e−iτε2Tεηn2(x)+Lnen(x) =e−2iτε2Tεen−1(x)+e−iτε2Tε(ηn−11(x)+e−iτε2Tεηn−12(x)+Ln−1en−1) +(ηn1(x)+e−iτε2Tεηn2(x)+Lnen) =... (2.35) =e−i(n+1)τTε/ε2e0(x)+n∑k=0e−i(n−k)τε2Tε(ηk1(x)+e−iτε2Tεηk2(x)+Lkek(x)).

Since , , and () preserves norm, we have from (2.34)

 (2.36) ∥∥ ∥∥n∑k=0e−i(n−k)τε2Tε(ηk1(x)+Lkek)∥∥ ∥∥H1≲n∑k=0τ2+n∑k=0τ∥ek(x)∥H1≲τ+τn∑k=0∥ek(x)∥H1,

 (2.37) ∥en+1(x)∥H1≲τ+τn∑k=0∥ek(x)∥H1+∥∥ ∥∥n∑k=0e−i(n−k+1)τε2Tεηk2(x)∥∥ ∥∥H1,n≤m.

To analyze , using (2.11), we can find , e.g.

 (Πε+Φ(tn))∗σ3(Πε−Φ(tn))= −ε(Πε+Φ(tn))∗σ3(R1Φ(tn))+ε(R1Φ(tn))∗σ3(Πε−Φ(tn)),

and the other terms in can be estimated similarly. As is uniformly bounded with respect to , we have (with detailed computations omitted)

 (2.38) ∥fn2(⋅)∥L∞([0,τ];(H1)2)≲ε∥Φ(tn)∥3H2≲ε.

Noticing the assumptions of Theorem 2.1, we obtain from (2.15)

 (2.39) ∥fn2(⋅)∥L∞([0,τ];(H1)2)≲ε,∥∂s(fn2)(⋅)∥L∞([0,τ];(H1)2)≲ε/ε2=1/ε,