The splitting methods form an important group of methods which are quite accurate and efficient . 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) , the Stokes equation , and the Enrenfest dynamics , 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 :
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
The nonlinearity in (1.1) is usually taken as
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
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 .
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 . Numerical studies in  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  or multiscale time integrator pseudospectral method  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 , 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 
, 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, 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 .
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
Denote the free Dirac Hermitian operator
then the NLDE (1.1) in 1D can be written as
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) , can be represented as:
For simplicity, we also write , where denotes the numerical propagator of the Lie-Trotter splitting.
2.1 A uniform error bound
In addition, we assume the exact solution satisfies
For the numerical approximation obtained from (2.3), we introduce the error function
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
Consequently, there is a uniform error bound for when
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 . Compared to the linear case , 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 . Suppose , and satisfy , we have
where depends on and . The proof is quite similar to the nonlinear Schrödinger equation case in  and we omit it here for brevity.
Under the assumption (B) (), for , we denote as
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
where is the Laplace operator in 1D, is the identity operator, and , are projectors defined as
It is straightforward to verify that , , , and through Taylor expansion, we have 
with for , for being uniformly bounded operators w.r.t. .
In order to characterize the oscillatory features of the solution, denote
which is a uniformly bounded operator w.r.t from for , then the evolution operator can be expressed as
For simplicity, here we use , in short.
Now we are ready to introduce the following lemma for proving Theorem 2.1.
with , , where depends on , , and ; depends on , , and . Here
where with and
where is the “local truncation error” (notice that this is not the usual local truncation error, compared with ),
By Duhamel’s principle, the solution to (2.2) satisfies
Noticing (2.20), the assumption that , and the fact that preserves norm, it is not difficult to find
On the other hand, using Taylor expansion in and the local Lipschitz property of , we get
It remains to estimate the part. Using the decomposition (2.13) and the Taylor exapnsion , we have (),
where for ,
Since is of polynomial type, by direct computation, we can further simplify (2.27) to get
where is given in (2.15) and is independent of as
Now, it is easy to verify that with given in Lemma 2.1 by choosing
which completes the proof of Lemma 2.1.
Denote (), and it is straightforward to calculate
with only depending on . Thus we can obtain from (2.33) that for ,
Since , , and () preserves norm, we have from (2.34)
which leads to
To analyze , using (2.11), we can find , e.g.
and the other terms in can be estimated similarly. As is uniformly bounded with respect to , we have (with detailed computations omitted)
which leads to
On the other hand, using Taylor expansion and the second inequality in (2.39), we have