On numerical methods for the semi-nonrelativistic limit system of the nonlinear Dirac equation

Solving the nonlinear Dirac equation in the nonrelativistic limit regime numerically is difficult, because the solution oscillates in time with frequency of Oε-2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathscr {O}} \! \left( \varepsilon ^{-2}\right) $$\end{document}, where 0<ε≪1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0<\varepsilon \ll 1$$\end{document} is inversely proportional to the speed of light. Yongyong Cai and Yan Wang have shown, however, that such solutions can be approximated up to an error of Oε2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathscr {O}} \! \left( \varepsilon ^2\right) $$\end{document} by solving the semi-nonrelativistic limit system, which is a non-oscillatory problem. For this system, we construct a two-step method, called the explicit exponential midpoint rule, and prove second-order convergence of the semi-discretization in time. Furthermore, we construct a benchmark method based on standard techniques and compare the efficiency of both methods. Numerical experiments show that the new integrator reduces the computational costs per time step to 40% and within a given runtime improves the accuracy significantly.


Introduction
The Dirac equation describes the relativistic motion of spin-1/2 particles such as, e.g., electrons, positrons, protons, neutrons, and quarks, under the influence of external electromagnetic fields. Since its derivation by Dirac in [10], it has become one of the cornerstones of relativistic quantum mechanics; cf. [26]. Nonlinear versions of the Dirac equation have been proposed to model self-interaction of particles and other phenomena; see, e.g., [11,23,25,27]. In the nonrelativistic limit regime, the Dirac equation involves a small parameter 0 < ε 1 which is inversely proportional to the speed of light, and non-trivial solutions oscillate in time with frequency of O ε −2 . Using traditional numerical methods to approximate such solutions is inefficient, because then the oscillations have to be resolved with a tiny step size, which causes prohibitively large computational costs; cf. [2]. Hence, constructing and analyzing numerical methods for the nonlinear Dirac equation in the nonrelativistic limit regime is a considerable challenge.
For the special case that no magnetic potential is present, it was shown in [3] that the classical Lie-Trotter and Strang splitting with step size τ have an error of O τ 1/2 , but uniformly in ε. For non-resonant step sizes, the accuracy improves to O(τ ) for Lie-Trotter and O τ 3/2 for the Strang splitting. For the general case (with magnetic potential), uniformly accurate methods have been proposed and analyzed in [5,8,18]. Under certain assumptions, the error of the time discretization with the multiscale time integrator pseudospectral method from [5] is bounded by a constant times min{τ 2 +ε 2 , τ 2 ε 2 }, where τ is the step size. This yields second-order convergence if τ ≥ ε or τ ε and, what is more important, first-order convergence uniformly in ε. The nested Picard iterative integrator constructed in [8] converges even with order two in time and uniformly in ε. However, the correct implementation of both methods is not easy, because they are based on complicated expansions and involve a plethora of terms. A different approach was proposed in [18] in one space dimension. The idea is to consider an augmented problem where the slow and fast time scales are distinguished. A formal Chapman-Enskog expansion is used to construct initial data for the augmented problem such that the corresponding solution has three uniformly bounded time derivatives, which paves the way for the construction of a uniformly accurate second-order scheme. However, the price to pay is that the augmented problem involves one additional dimension representing the fast time scale, which increases the numerical work significantly.
The nonlinear Dirac equation in the nonrelativistic limit regime has also been intensively studied in analysis; cf. [7,[20][21][22]24]. It was shown that the solution ψ ε (t, x) ∈ C 4 can be approximated by ψ ε (t, x) = e −itβ/ε 2 ϕ(t, x) + O(ε) (1.1) where β = diag(1, 1, −1, −1) ∈ R 4×4 is a diagonal matrix, and where ϕ is the solution of a nonlinear Schrödinger equation which does not depend on ε and is thus easier to approximate numerically. A precise formulation of this result and its proof are given in [7,Theorem 2.3]. The main result of [7], however, is that a better approximation ψ ε (t, x) = e −it/ε 2 ϕ + ε (t, x) + e it/ε 2 ϕ − ε (t, x) + O ε 2 (1.2) can be obtained, where ϕ ± ε are the solutions of two coupled semilinear PDEs called the semi-nonrelativistic limit system; see [7,Theorem 2.2] or Theorem 2.3 below for details. In contrast to the above-mentioned nonlinear Schrödinger equation, the seminonrelativistic limit system does still involve the parameter ε, but in contrast to the original problem, the solution does not oscillate in time; cf. [7,Theorem 2.2]. Hence, (1.2) offers a way to approximate the highly oscillatory solution of the nonlinear Dirac equation without having to solve a highly oscillatory problem. Of course, one cannot expect the error of this approximation to be smaller than O ε 2 , but in this work we assume that this accuracy is sufficient. If a higher accuracy is required, one has to use the uniformly accurate integrators from [8,18] with a step size τ < ε, which is computationally intense.
Solving the semi-nonrelativistic limit system numerically is much easier than solving the nonlinear Dirac equation in the nonrelativistic limit regime, but it is not straightforward. For example, explicit Runge-Kutta methods suffer from severe CFL conditions, whereas fully implicit methods come at the price of solving a large nonlinear system in every time step. Constructing splitting methods 1 in a straightforward way is not an option, either, due to the particular structure of the semi-nonrelativistic limit system. These problems can be avoided with exponential integrators. Such integrators are typically constructed by applying variation of constants on the interval [t n , t n+1 ] (where t n = nτ are the times where numerical approximations are supposed to be computed) and approximating the convolution integral, e.g. by expanding the nonlinearity in such a way that the integral can be solved analytically. The corresponding techniques are nowadays well-known in the context of dispersive equations and in particular highly oscillatory problems. Because of the special structure of the semi-nonrelativistic limit system, however, such a method requires many forward and backward Fourier transforms per time step, which is the dominating factor in the computational costs. In this paper, we propose a non-standard second-order exponential integrator. The idea is to apply variation of constants over the interval [t n−1 , t n+1 ], which simplifies the treatment of the nonlinearity a lot. This approach leads to a twostep method which we call the explicit exponential midpoint rule. The new method is time-symmetric, simpler to implement, and considerably more efficient than the standard second-order exponential integrator.
In Sect. 2 we introduce the nonlinear Dirac equation in the nonrelativistic limit regime and sketch the derivation of the semi-nonrelativistic limit system as presented in [7]. Moreover, we specify our assumptions and quote a number of important results from [7]. Time-integrators for the semi-nonrelativistic limit system are constructed in Sect. 3. The first method is an exponential integrator which is based on well-known techniques, and which is therefore considered as a benchmark method. The second method is the explicit exponential midpoint rule. For this integrator we carry out a detailed error analysis; cf. Theorem 3.7. In Sect. 4 we test the efficiency of both methods in a numerical experiment. It turns out that our new method reduces the computational costs per time step to about 40% and, within the same runtime as the benchmark method, improves the accuracy by a factor of about 4.6 in L 2 and 6 in H 1 . We explain the reason for these improvements.

Nonlinear Dirac equation in the nonrelativistic limit regime
We consider the nonlinear Dirac equation (NLDE) for x ∈ R 3 and t > 0. In (2.1), ψ ε := ψ ε (t, x) ∈ C 4 is the complex-valued vector wave function with initial data ψ 0 = ψ 0 (x) ∈ C 4 . The parameter ε > 0 is inversely proportional to the speed of light and thus is very small in the nonrelativistic limit regime. Furthermore, T ε and W denote the free Dirac operator and the electromagnetic potential, respectively, given by where V (t, x) ∈ R is the electric scalar potential and A(t, x) = (A 1 (t, x), ..., A 3 (t, x)) T is the magnetic vector potential. The Dirac matrices are determined by the Pauli matrices Finally, F is the nonlinearity given by where v * = v T denotes the conjugate transpose and |v| = √ v * v the Euclidean norm of a vector v, respectively. This type of nonlinearity is motivated by numerous applications in physics and describes self-interaction of Dirac fermions; see, e.g., [11,23,25,27] and the references in [2,3,5,8,18]. For simplicity, we assume that λ = 0 henceforth, but all results and proofs can be adapted to the the case λ = 0.
Throughout this paper, we will use the following notation: v H m denotes the standard Sobolev norm of a scalar-valued function v ∈ H m (R 3 ), whilst for a C 4 -valued The following assumptions regarding the initial data and the potential W will be made. Recall that W is determined by V and A j via (2.2). Assumption 2.1 Let 0 < T 0 < ∞ be an arbitrary fixed time. For some m ≥ 2, we assume that The following theorem quoted from [ with uniform estimates where C is independent of ε.
The original formulation of this result in [7] is slightly more general and applies also to the case where the initial data in (2.1) depend on ε to some extent. Solving (2.1) numerically is a challenging task, because typical solutions oscillate in time with frequency of O ε −2 due to the term − i ε 2 T ε ψ ε (t, x) on the right-hand side. Applying traditional time-integrators such as, e.g., Runge-Kutta or standard multistep methods is inefficient, because such methods only achieve an acceptable accuracy if the ratio of the step size and the highest frequency is small; see, e.g., [2]. One possibility to solve this problem is to construct special integrators which do not suffer from such a severe step size restriction. This has been done in [5,8,18], but the implementation of such methods is quite involved. In this paper we pursue a different goal. In [7] it was rigorously shown that in the limit ε → 0 the solution of the NLDE (2.1) can be approximated up to O ε 2 by solving a non-oscillatory system of PDEs known as the semi-nonrelativistic limit system; cf. (2.10) below. Since this accuracy is good enough in many applications, our main goal is to construct a particularly efficient method for the semi-nonrelativistic limit system. This is done in Sect. 3. Before that, we briefly outline the derivation of the semi-nonrelativistic limit system given in [7].

Transformed Dirac equation
In this and the next subsection we summarize the main results from [7].
By performing an eigenspace decomposition in Fourier space, the operator T ε can be decomposed as [4,Eq. (1.22)] with the scalar operator Λ ε = √ I d − ε 2 Δ and the two projection operators Π ± ε given by can easily be checked; cf. [7]. As a mapping from (H m (R 3 )) 4 to (H m (R 3 )) 4 the projectors Π ± ε are uniformly bounded w.r.t. ε; cf. [4, Lemma 2.1]. The decomposition (2.3) allows us to filter out the main part of the temporal oscillations in a solution ψ ε of the NLDE (2.1). This is achieved by considering the functions instead of ψ ε . Substituting (2.4) into the NLDE (2.1) shows that φ + ε and φ − ε are the solution of the two coupled PDEs with the differential operator cf. [7,Sect. 2.1]. From the solution φ ± ε of (2.5) we can reconstruct the solution ψ ε of the NLDE (2.1) by In Fourier space, application of D ε to a function v corresponds to multiplication of the Fourier transform of v at ξ ∈ R 3 with This yields the bound which means that the operator is uniformly bounded w.r.t. ε for all n ∈ N 0 . Hence, the first time derivative of a solution φ ± ε is uniformly bounded w.r.t. ε, which is not true for a solution ψ ε of the NLDE due to the factor 1/ε 2 on the right-hand side. In this sense (2.5) is better suited for numerical purposes than the original form (2.1) of the NLDE. However, solving (2.5) with standard methods still suffers from severe step size restrictions, because the solution of (2.5) still oscillates with the same frequency as the original problem, albeit with smaller amplitude. In the next subsection, these oscillations are completely removed at the cost of an approximation error.

Semi-nonrelativistic limit system
Omitting the terms containing highly oscillatory exponential functions in (2.5) (including those in the nonlinearity g) yields the semi-nonrelativistic limit system [7, Eq. (2.14)] (2.10) Well-posedness of (2.10) and regularity of solutions of (2.10) has been shown in [7]. Furthermore, the authors proved that solutions of (2.10) provide approximations to a solution of the original problem (2.1):

Theorem 2.3 [7, Theorem 2.2] Under the assumptions (A) and (B)
, there is a time T 2 ∈ (0, T 0 ] such that for any ε ∈ (0, 1], the semi-nonrelativistic limit system (2.10) admits a unique solution Moreover, ϕ ± ε remain in the eigenspaces associated with Π ± ε , respectively. If in addition the assumption V , In [7] this theorem is formulated in a more general way which, however, exceeds our demands.
The inequality (2.11) implies that a solution of the NLDE (2.1) can be approximated up to O ε 2 using a solution of the semi-nonrelativistic limit system. In this paper, we consider the case where ε is small enough such that this approximation is satisfactory. Thus, instead of developing time-integrators for the NLDE (2.1) or its transformed version (2.5), we can focus on the simpler semi-nonrelativistic limit system (2.10). Solutions of (2.10) are not affected by oscillations, because there is neither a factor ε −2 on the right-hand side (in contrast to (2.1)) nor oscillating exponentials (in contrast to (2.5)). The solution only depends on ε because the projectors Π ± ε and the differential operator D ε do, but in a non-critical way. In spite of these advantages, solving (2.10) with standard methods is still not a good option. If explicit Runge-Kutta or multistep methods are used, then the spatial discretization of the differential operator D ε causes severe CFL conditions, whereas a time step with an implicit method is somewhat costly due to the nonlinearity and the projectors Π ± ε . Applying a splitting method to (2.10) in a straightforward way is not feasible, because the sub-problems involving the projectors Π ± ε cannot be propagated exactly or particularly efficiently. These disadvantages can be avoided by exponential integrators. Two such methods are presented and compared in the next section.

Time integration methods for the semi-nonrelativistic limit system
Our goal now is to compute approximations ϕ ± n ≈ ϕ ± ε (t n ) of the solution of the semi-nonrelativistic limit system at discrete times t n = nτ , where τ > 0 is the step size. We propose two exponential integrators which converge with order two in τ . The first one is constructed by applying variation of constants over the interval [t n , t n+1 ], approximating the integrand in a suitable way and computing the resulting integrals exactly. This strategy is, of course, not new, and the related techniques have been used for various types of PDEs, in particular in the context of highly oscillatory problems. We consider this first method only as a benchmark method, and for this reason we refrain from an extensive error analysis. The main contribution of this paper is the second time-integrator. The crucial idea is to use variation of constants over the time interval [t n−1 , t n+1 ] instead of [t n , t n+1 ], which makes the approximation of the resulting integrals much easier. This leads to a novel exponential two-step method called the explicit exponential midpoint rule (EEMR). This time-integrator is timesymmetric, simple to implement, and considerably more efficient than the benchmark method. We present a detailed error analysis for the EEMR and explain the speed-up observed in numerical examples.
Assumptions and notation. In the end, we want to obtain error bounds in L 2 . Our methods will rely on the uniform boundedness of the second time derivative of ϕ ± in L 2 . According to Theorem 2.3, this is given under the following assumptions that we will assume for the rest of the paper: In order to prove convergence of the methods, we will also require the assumption Assumptions (I) and (II) coincide with Assumption 2.1 for m = 4.
To increase readability, we define the function space for T = min{T 1 , T 2 } from Theorem 2.3, which then states that if assumptions (I)-(III) hold, then ϕ ± ∈ S T with uniform bounds in ε.
From now on we assume that ε is small but fixed. We can thus omit the index ε in our notation such that the semi-nonrelativistic limit system (2.10) reads All bounds presented below are uniform in ε in the sense that the constants do not depend on ε.
Let τ ∈ (0, T ) be the time step size and let t n = nτ , n = 0, 1, ..., T /τ . To improve readability, we omit the spatial variable x on the solution and the potential in the following. For a function f ε = f ε (s), we write f ε = O(s p ) for some p ∈ N 0 to express that f ε (s) L 2 ≤ Cs p for s → 0 with some constant C which does not dependent on s and ε.
We will repeatedly use that there is a constant C such that 4 . These inequalities follow from the Sobolev embedding

The benchmark method
Using variation of constants, we can express the solution ϕ ± of the semi-nonrelativistic limit system (3.1) at time t n + τ as given by The operators e ∓i(τ −s)D and Π ± are both bounded in (L 2 (R 3 )) 4 . Thus, in order to obtain a third-order approximation (in τ ) to the integrals, we need a second-order approximation (in s) to the integrands. Under assumption (IV), W (t n + s) can be replaced by the Taylor expansion 4 under assumptions (I)-(III), we can also expand being the first time derivative of ϕ ± at time t n . It is obtained by evaluating the right-hand side of the PDE (3.1) at time t n : Before we continue by inserting (3.5) and (3.6) into (3.4), let us quickly comment on an alternative approach to construct a second-order approximation to ϕ ± (t n + s).
Using variation of constants once again, but now over a time interval of length s, and fixing ϕ ± as well as W at time t n inside the integrals yields This approach does only rely on boundedness of the first time derivative of ϕ ± and thus is, at first glance, feasible under lower regularity assumptions on the potential W and the initial data. Unfortunately, inserting (3.8) into (3.4) leads to integrals which cannot be computed analytically. In order to avoid this problem, we could use the formal approximations Using these approximations in (3.8), however, yields exactly the same second-order approximation to ϕ ± (t n + s) as (3.6) together with (3.7). Now we continue the construction of the benchmark method. Inserting (3.5) and (3.6) into I ± 1 yields is given by The functions p 1 and p 2 are defined by When inserting (3.6) into I ± 2 , we can additionally drop all O s 2 -terms in the integrand that arise due to the nonlinearity. Overall, we obtain where I ± 2 = I ± 2 (ϕ + (t n ), ϕ − (t n ), t n ) is given by and Θ ± from (3.7). A third-order approximation to ϕ ± (t n + τ ) is obtained by simply replacing the integrals I ± 1 and I ± 2 in (3.4) by their approximations I ± 1 and I ± 2 : This approximation suggests a numerical method with local error of order O τ 3 which, however, would not be stable. The reason for this instability is the term ∓iDϕ ± (t n ) which appears in Θ ± , cf. (3.7), and thus also in ζ ± . A bound for the norm of D that is independent of ε can only be established when interpreting D as mapping from H 2 to L 2 , cf. (2.9). Hence, the L 2 -norm of Θ ± and thus ζ ± , I ± 1 and I ± 2 can only be bounded using the H 2 -norm of ϕ ± (t n ), which would not be sufficient for stability. This is why we replace D in Θ ± by a filtered version as, e.g., in [6,8]. It is not difficult to show that for every τ > 0, D(τ ) is a bounded operator from L 2 to L 2 with D(τ ) ≤ 1 τ , and that 4 under assumptions (I) and (II), it follows that replacing D by D(τ ) in Θ ± and hence also in ζ ± causes an error of O(τ ). But in I ± 1 and I ± 2 , the terms including Θ ± or ζ ± are multiplied by a factor τ 2 . Thus, substituting D(τ ) for D in the right-hand side of (3.10) causes only an additional error of O τ 3 and hence does not affect the overall approximation error. All in all, this yields the numerical method with the numerical flow (3.13) I ± 1 and I ± 2 correspond to I ± 1 and I ± 2 , respectively, but with D replaced by D(τ ) in Θ ± and ζ ± , i.e.
For an efficient implementation, the two integrals I ± 1 and γ I ± 2 can be combined to Under assumptions (I)-(IV) the local error in L 2 is bounded by Cτ 3 by construction.
With well-known techniques, it can be shown that under assumptions (I)-(IV) there are constants τ 0 > 0 and C such that for all step sizes τ ∈ (0, τ 0 ] the bound for the global error holds. We omit the proof, because our focus is not on the benchmark method. The step size restriction τ ≤ τ 0 is required to obtain uniform boundedness of the numerical approximations in H 2 (R 3 ), which is required for stability; for the EEMR this issue is discussed in the proof of Theorem 3.7.

Remark 3.1
The method (3.12)-(3.13) is certainly not new. We have described the construction only for the convenience of the reader and in order to keep the paper self-contained. In fact, (3.12)-(3.13) coincides with a "part" of the multiscale method for the NLDE (2.1) which has been proposed in [5]. The idea is, roughly speaking, to make the ansatz i.e. to decompose the solution of (2.1) into the part provided by the semi-nonrelativistic limit system plus a rest r (t, x) which, according to (2.11), is only O ε 2 . Substituting this ansatz into the NLDE and replacing ∂ t ϕ ± by (3.1) yields a PDE for r (t, x) with a rather complicated right-hand side. Then, a numerical method for ϕ ± and r is constructed in [5]. Within this method the part which approximates ϕ ± is almost identical to what we call the benchmark method. The only differences are that instead of (3.11) a different filter is used in [5], and that the authors consider the full discretization in time and space.

Remark 3.2
In this work we only consider time discretizations. For a full discretization in time and space on the torus, the benchmark method (3.12)-(3.13) can be combined with a Fourier pseudospectral method, such that ϕ ± n is approximated by a trigonometric polynomial. All operators involving spatial derivatives (which includes the projectors Π ± ) are applied in Fourier space, whereas pointwise multiplications of functions such as, e.g., W (t n )ϕ ± n or W (t n ) Θ ± correspond to entry-wise multiplications of vectors. In order to compute all terms required for one time step, the fast Fourier transform (FFT) or its inverse has to be applied quite a number of times, and in spite of the efficiency of the FFT, this causes the dominating part of the numerical work.

Explicit exponential midpoint rule
We will now propose and analyze a new exponential integrator which converges with order two under the same regularity assumptions as the benchmark method, but which is conceptually simpler, easier to implement and significantly faster. The new integrator is time-symmetric, in contrast to the benchmark method. For time-dependent potentials, the new method does not require evaluations of ∂ t W , which is convenient in situations where no explicit formula for W (t, x) is available.

Construction
We again use variation of constants to express the solution ϕ ± of the semi-nonrelativistic limit system at time t n + τ , but now over a time interval of length 2τ . This yields Under the assumptions (I)-(III), we know that ϕ ± ∈ S T with uniformly bounded (w.r.t. ε) derivatives according to Theorem 2.3. If additionally assumption (IV) is fulfilled, then according to the estimates (3.2) and (3.3) the same holds for the functions Π ± W ϕ ± and γ Π ± ϕ + 2 + ϕ − 2 ϕ ± which appear in the integrands of I ± 1 and I ± 2 . For a function v ∈ S T , a third-order approximation to integrals of the form τ −τ e ∓i(τ −s)D v(t n + s) ds is obtained by fixing v at the midpoint t n , as the following lemma confirms.

Lemma 3.3 Let v ∈ S T and τ ∈ (0, T ). Then
for some constant C that only depends on ∂ t v H 2 and ∂ tt v L 2 .
2 Note that I ± 1 and I ± 2 are different from the integral terms which were denoted with I ± 1 and I ± 2 in the previous section. Many other objects which appeared in the previous section such as, e.g., I ± 1 , I ± 2 , Φ ± τ , Φ τ etc., will be re-defined in a different way in this section.

Proof
Since v is twice continuously differentiable, Taylor's theorem yields with the remainders For R 2 , we can derive the bound We analyze the norm of R 1 in Fourier space. Recall that application of D corresponds to multiplication with δ ε (ξ ) in Fourier space, cf. (2.7). This yields where we used the bound (2.7) on δ ε (ξ ) in the last step. Overall, we have Since v ∈ S T , the assertion follows from (3.15) and (3.16).
After fixing v at the midpoint as in the previous lemma, the remaining integral can be computed as with p 1 from (3.9). Applying the lemma to the integrals I ± 1 and I ± 2 in (3.14) thus yields given by Omitting the O τ 3 -terms and replacing exact solutions with approximations ϕ ± n ≈ ϕ ± (t n ) leads to the integrator with the numerical flow For an efficient implementation, we can again combine the two integrals I ± 1 and γ I ± When using a Fourier pseudospectral method for space discretization on the torus, only one FFT per time step is required in the computation of I ± before being able to apply the operators p 1 (∓2iτ D) Π ± . If the approximations obtained in the two previous steps are saved in physical as well as in Fourier space, only one inverse FFT is necessary for retransforming the result of (3.19)-(3.20) into physical space.
We call this method the explicit exponential midpoint rule (EEMR), because it can be regarded as the exponential counterpart of the classical explicit midpoint rule. Since (3.19)-(3.20) is a two-step method, the first approximation ϕ ± 1 ≈ ϕ ± (t 1 ) has to be computed with a starting step. Only an accuracy of O τ 2 is required for ϕ ± 1 , and this can be achieved easily by using variation of constants over a time span of length τ as in the benchmark method and then approximating the integrals via the rectangle rule. After omitting the O τ 2 -terms and replacing ϕ ± (t 1 ) with ϕ ± 1 we obtain To improve the accuracy of ϕ ± 1 , one can replace (3.21) by η ∈ N such steps with step size τ/η.

Remark 3.4
The idea to construct exponential multistep methods by applying variation of constants over the time-interval [(n − )τ, (n + 1)τ ] for some ≥ 1 has already been used in [9,14,15], but the PDEs and the methods considered in these references are completely different. The exponential multistep methods reviewed in [12,Section 2.5] are of Adams type, which is different from what we propose here.

Error analysis
Our goal is to prove that the EEMR (3.19)-(3.21) is indeed second-order convergent under the regularity assumptions (I)-(IV), which are also required for the benchmark method. For this purpose we reformulate the EEMR as a one-step method by introducing the vectors We then have u n+1 = Φ τ (u n , t n ) for n ∈ N with the numerical flow with Φ ± τ defined by (3.20). For vectors of the form v = (v 1 , ..., v 4 ) T with four functions v 1 , ..., v 4 ∈ (L 2 (R 3 )) 4 , we define the norm and analogously for the H m -norm for some m > 0. The following bounds for the local error and the starting step are an immediate consequence of the construction of the EEMR.

Under assumptions (I)-(IV) there is a constant C E 2 such that the inequality
holds for some constant C E 2 .
Proof By definition of u, Φ τ and |||·||| L 2 we have Since all approximations made during the construction of the method are O τ 3 , the bound (3.23) follows. In a similar way, the bound for can be shown with standard arguments.
Next, we discuss stability. In order to simplify presentation, we will henceforth assume that the electromagnetic potential W does not depend on time. In this case, the semi-nonrelativistic limit system (2.10) is autonomous, and as a consequence the numerical flows Φ ± τ defined in (3.20) and Φ τ defined in (3.22) do not depend on t, either. This allows us to omit the last variable in Φ ± τ and Φ τ , which makes the following equations easier to read. We stress, however, that under Assumption (I)-(IV) the following proofs could be extended to a time-dependent W at the cost of a more involved notation.

Under assumptions (I)-(II) there is a constant C S such that the stability estimate
holds for all τ > 0. The constant C S depends on V H 2 , A j H 2 , j = 1, 2, 3, v ± 1 H 2 and w ± 1 H 2 , but not on τ .
Proof It follows from (3.22) and (3.20) that (3.27) Since e ∓2iτ D is an isometry in L 2 , the first term on the right-hand side equals |||v − w||| L 2 . Now we insert the definition (3.17) of I ± 1 and use that p 1 (∓2iτ D) and the projectors Π ± are bounded operators in L 2 . Applying the estimate (3.2) for the product with the potential W yields the inequality for a constant C which depends on V H 2 and A j H 2 , j = 1, 2, 3, but not on τ . Now we will prove such a bound for the term I ± 2 containing the nonlinearity. The estimates (3.2) and (3.3) imply the inequality v + with a constant C which depends on v ± 1 H 2 and w ± 1 H 2 . Together with the arguments mentioned above, we obtain for some constant C which depends on v ± 1 H 2 and w ± 1 H 2 . Combining (3.27), (3.28) and (3.30) proves the assertion.
We are now in a position to prove second-order convergence for the EEMR.
Theorem 3.7 (Global error of the EEMR) Assume that assumptions (I) and (II) hold and that W does not depend on t. Let τ > 0 be the step size and let ϕ ± n be the approximations obtained by (3.19) and (3.21) with step size τ and initial data ϕ ± 0 = ϕ ± (0) = Π ± [ψ 0 ]. Then, there are constants C and τ 0 > 0 such that the global error bound holds for all τ ∈ (0, τ 0 ].

Remark 3.8
Under the assumptions (I)-(IV) the theorem remains true for a timedependent potential W = W (t).

Remark 3.9
Under stronger regularity assumptions on the initial data and the potentials, one could of course obtain error bounds in higher-order Sobolev spaces. More precisely, to obtain an identical error bound in H m , m ≥ 0, one would require the initial data ψ 0 and the potentials V , A j to be in H m+4 with ) recursively for n ∈ N, such that Φ n τ (v) denotes the result of n steps of the EEMR in the one-step formulation with initial data v.
In order to prove the global error bound, we combine the local error bounds (3.23) and (3.24) with the stability estimate (3.26) in the classical construction known as Lady Windermere's fan. Using a telescopic sum, we have for n = 1, 2, ..., T /τ (3.31) At this point, we would like to control the term by applying the stability estimate (3.26) k times. A minor technical difficulty is the fact that the constant C S in (3.26) depends on the H 2 -norms of the two functions involved, which in our situation are Φ j τ (u(t n−k )) and Φ j+1 τ (u(t n−k−1 )), respectively, with j = 0, . . . , k − 1. In order to obtain a corresponding bound for (3.32) with a constant which does not depend on n, k or τ , we need that there are constants τ 0 and C such that max j, =0,..., T /τ holds for some constant c s which depends on |||v||| H 2 and |||w||| H 2 . Note that in contrast to (3.26), the H 2 -norm is used on both sides of (3.34). The estimate can be shown by adjusting the proof of Lemma 3.6. Secondly, one has to prove that under the assumptions of Lemma 3.5, the local error bound holds for some constant C E 3 and for all τ > 0 and n = 1, 2, ..., T /τ . In contrast to (3.23), the local error is measured in the H 2 -norm instead of L 2 , but the power of τ is only 2 instead of 3. The proof of (3.35) is straightforward. Having established (3.34) and (3.35), one can prove (3.33) by induction. Since the procedure is essentially the same as, e.g., in [19] or [16,Section 8], we omit this part. Now we return to (3.31). Combining (3.33) with the stability estimate (3.26) yields under the condition τ ≤ τ 0 that with a constant C which only depends on V H 2 , A j H 2 , j = 1, 2, 3, and on C from (3.33). Applying the local error bound (3.23) as well as the bound (3.24) for the starting step and using that (n − 1)τ = t n−1 ≤ T shows that The EEMR for the semi-nonrelativistic limit system is thus indeed second-order accurate. Both the solution and the numerical approximation still depend on ε (see the remark at the beginning of this section), but all bounds are uniform in ε.

Numerical experiments
In this section, we present numerical results to illustrate our error analysis, and to compare the efficiency of the EEMR and the benchmark method. For simplicity, we conduct our experiments in one space dimension, where the NLDE can be reduced to the system x) ∈ C 2 with initial data ψ ε (0, x) = ψ 0 (x) ∈ C 2 (see e.g. [2]). Here, the differential operator T , the potential W and the nonlinearity F are given by For simplicity, we omit the˜in the following, and we chose γ = 1. The properties of the semi-nonrelativistic limit system, the construction of the numerical methods presented above as well as the obtained error results can be formulated for this reduced system in exactly the same manner.
As is common practice [1][2][3]13], we truncate the whole space problem to a bounded interval Ω = [a, b] which is large enough such that the truncation error is negligible. We impose periodic boundary conditions and discretize Ω through the grid points For all following error plots, we compute approximations to solutions at time t end = 1 using the presented methods, and compare them to reference solutions computed via Matlab's ode45 routine using the same spatial grid and very small tolerances. The error is always computed in the L 2 -norm, approximated by v 2 Performing various numerical experiments, we have observed that an even number of time steps gives slightly better results for the two-step method. Therefore, in all following plots, the step sizes are chosen such that the number of time steps t end /τ is even. The quality of approximations after many time steps can be further improved In the experiments below, we always use η = 3 substeps in the starting step (3.21) to obtain an approximation at time t 1 . Increasing η further did not improve the results significantly.
Accuracy. Figure 1shows the global error of both methods in dependency of the step size τ in logarithmic axes. On the left-hand side, the approximations obtained by both methods are compared to a reference solution ϕ ± of the semi-nonrelativistic limit system (3.1). The solid lines represent the benchmark method, the coloured dashed lines the EEMR. Different values of ε are depicted through different colours, but the six lines coincide almost. The results confirm that both methods are second-order accurate in the step size τ , and that the error constants do not depend on ε. For a fixed step size τ , both methods yield approximations of nearly the same accuracy in this example. On the right-hand side of Fig. 1, the approximations are compared to a reference solution of the NLDE (2.1). Here the numerical approximations ϕ ± n are interpreted as approximations to a solution φ ± of the transformed Dirac equation (2.5). Then, according to (2.6) the function approximates a solution of the original problem (2.1). Since the overall error is composed of two parts: -the approximation error of the numerical methods in comparison to the exact solution of the semi-nonrelativistic limit system, which is of order O(τ 2 ), -the difference between solutions of the semi-nonrelativistic limit system and the transformed NLDE, which, for a fixed time t end , is of order O(ε 2 ).
The overall approximation error is thus of order O(τ 2 ) + O(ε 2 ). Consequently, the O(τ 2 )-terms are dominating for large step sizes τ > ε, and we observe second-order convergence w.r.t. τ until τ = ε (this value is indicated by the vertical dashed-dotted lines). For τ < ε however, the O(ε 2 )-terms are dominating and no further convergence is achieved when the step size τ is reduced. Thus, applying the two methods to the semi-nonrelativistic limit system allows us to compute very accurate approximations to the highly oscillatory solution of the original NLDE in the non-relativistic regime, where ε is very small.
Efficiency. Whilst the experiments above suggest that both methods perform equally well, they do not take into account the computational effort required for one time step in each method. In the EEMR, we have symmetry of the integration interval of the integrals I ± 1 and I ± 2 in the variation of constants formula. This is why the required accuracy was achieved by essentially only keeping the constant term of the Taylor expansions of the integrands. In contrast to that, in the benchmark method, the linear terms of such Taylor expansions had to be taken into account as well. Those terms include several pointwise multiplications of space-dependent functions (with the potential W or with the functions ϕ ± itself) as well as applications of the projectors Π ± . Whilst the former has to be done in physical space, the latter can only be done in Fourier space. Consequently, computing those linear terms requires additional (inverse) Fourier transforms, which are the dominating operations in computational costs; cf. Remark 3.2. One time step of the benchmark method is thus significantly more expensive than of the EEMR. In an efficient implementation, one time step of the EEMR can be done using one Fourier and one inverse Fourier transform (where we count one transformation of a function v = (v + , v − ) T , v ± ∈ (L 2 (Ω)) 2 , into or out of Fourier space as one transform). One time step of the benchmark method, however, requires three Fourier and two inverse Fourier transforms. Hence, the computational costs of the benchmark method is about 5 2 times larger than one time step of the EEMR. Figure 2 shows the results of numerical experiments comparing the efficiency of both methods. For the plot on the left hand side, both methods have been applied using different step sizes and their computation time has been measured. To lower the impact of background processes, the average of multiple runs is taken and the different step sizes are used in random order. For each method, a reference line has been added that fits best to the measurement points. Comparing the constants of those lines, one can see that for a given computation time, the error of the benchmark method is about 4.6 times larger than the error of the EEMR, despite the fact that for this specific data, the EEMR performs a little worse for a fixed time step size than the benchmark method. On the right hand side, the number of Fourier and inverse Fourier transforms is counted. Again, for a fixed number of Fourier transforms, the error of the benchmark method is about 4.7 times larger than the error of the EEMR.
If the error is measured in the discrete H 1 -norm instead (cf. Remark 3.9), then the accuracy of both methods for a fixed time step is almost indistinguishable for this model problem (data not shown). For a fixed computation time, the error of the benchmark method in the discrete H 1 -norm is six times larger than the error of EEMR (data not shown).
The reason for the observed factors is the following. Let w j (N j ) be the numerical work (measured in runtime or in number of FFTs) required for N j steps with the benchmark method ( j = 1) or the EEMR ( j = 2), respectively. For a given N we have w 1 (N ) ≈ 5 2 cN and w 2 (N ) ≈ cN with some constant c > 0. If we fix the numerical work w > 0 we can thus perform N 1 ≈ 2w 5c steps with the benchmark method and N 2 ≈ w c with the EEMR. According to Fig. 1 the error err j (N j ) of both methods is err j (N j ) ≈ C j N −2 j with some error constant C j . Hence, the errors for a fixed numerical work w are This implies that err 1 (N 1 ) ≈ C 1 C 2 25 4 err 2 (N 2 ). For the error in the L 2 -norm, the ratio C 1 C 2 ≈ 0.8 can be calculated from Fig. 1, which yields err 1 (N 1 ) ≈ 5err 2 (N 2 ). For the error in H 1 , we have that C 1 C 2 ≈ 1, which yields err 1 (N 1 ) ≈ 6.25err 2 (N 2 ). Although the values 5 and 6.25 are a bit larger than the observed factors 4.7 and 6, they predict approximately the improvement obtained with EEMR.
Funding Open Access funding enabled and organized by Projekt DEAL. This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -Project-ID 258734477 -SFB 1173.

Conflict of interest
The authors have no competing interests to declare that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.