Uniform error bounds of time-splitting methods for the nonlinear Dirac equation in the nonrelativistic regime without magnetic potential

TL;DR

This paper rigorously analyzes super-resolution properties of Lie-Trotter and Strang splitting methods for the nonlinear Dirac equation in the nonrelativistic regime, revealing their uniform error bounds and super-resolution capabilities despite high oscillations.

Contribution

It provides the first rigorous proof of super-resolution and uniform error bounds for splitting methods applied to the nonlinear Dirac equation without magnetic potential in the nonrelativistic limit.

Findings

Both $S_1$ and $S_2$ exhibit 1/2 order convergence uniformly in $ ext{varepsilon}$.

Non-resonant time steps improve convergence to first and 3/2 order for $S_1$ and $S_2$ respectively.

Super-resolution holds for higher order splitting methods.

Abstract

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 regime with a small parameter $0<\varepsilon\leq 1$ inversely proportional to the speed of light. In this regime, the solution highly oscillates in time with wavelength at $O(\varepsilon^2)$. The splitting methods surprisingly show super-resolution, i.e. the methods can capture the solution accurately even if the time step size $\tau$ is much larger than the sampled wavelength at $O(\varepsilon^2)$. Similar to the linear case, $S_1$ and $S_2$ both exhibit $1/2$ order convergence uniformly with respect to $\varepsilon$. Moreover, if $\tau$ is non-resonant, i.e. $\tau$ is away from certain region determined by $\varepsilon$, $S_1$ would yield an improved uniform first order $O(\tau)$ 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.

Tables

Table 1. Table 3.2: Discrete H 1 superscript 𝐻 1 H^{1} temporal errors e ε , τ ​ ( t = 2 ​ π ) superscript 𝑒 𝜀 𝜏 𝑡 2 𝜋 e^{\varepsilon,\tau}(t=2\pi) for the wave function of the NLDE ( 2.2 ) with resonant time step size, S 2 subscript 𝑆 2 S_{2} method.

$e^{\epsilon ,\tau }\left(t=2\pi \right)$	${\tau }_0=\pi /4$	${\tau }_0/4$	${\tau }_0/4^2$	${\tau }_0/4^3$	${\tau }_0/4^4$	${\tau }_0/4^5$	${\tau }_0/4^6$
${\epsilon }_0=1$	1.17E+1	2.55E-1	1.37E-2	8.49E-4	5.30E-5	3.31E-6	2.07E-7
order	–	2.76	2.11	2.01	2.00	2.00	2.00
${\epsilon }_0/2$	4.63	4.32E-1	7.83E-3	4.84E-4	3.02E-5	1.89E-6	1.18E-7
order	–	1.71	2.89	2.01	2.00	2.00	2.00
${\epsilon }_0/2^2$	4.36	1.50	1.04E-2	6.00E-4	3.73E-5	2.33E-6	1.45E-7
order	–	0.77	3.59	2.05	2.00	2.00	2.00
${\epsilon }_0/2^3$	3.61	8.39E-1	7.79E-1	1.02E-3	5.98E-5	3.72E-6	2.32E-7
order	–	1.05	0.05	4.79	2.05	2.00	2.00
${\epsilon }_0/2^4$	3.51	4.38E-1	4.14E-1	4.02E-1	1.19E-4	6.95E-6	4.32E-7
order	–	1.50	0.04	0.02	5.86	2.05	2.00
${\epsilon }_0/2^5$	3.50	2.44E-1	2.09E-1	2.08E-1	2.05E-1	1.47E-5	8.55E-7
order	–	1.92	0.11	0.00	0.01	6.89	2.05
${\epsilon }_0/2^9$	3.46	1.10E-1	1.45E-2	1.31E-2	1.31E-2	1.31E-2	1.31E-2
order	–	2.49	1.46	0.07	0.00	0.00	0.00
${\epsilon }_0/2^{13}$	3.45	1.08E-1	4.76E-3	9.11E-4	8.21E-4	8.18E-4	8.18E-4
order	–	2.50	2.25	1.19	0.08	0.00	0.00
${\epsilon }_0/2^{17}$	3.45	1.08E-1	4.57E-3	3.18E-4	7.94E-5	7.57E-5	7.57E-5
order	–	2.50	2.28	1.92	1.00	0.03	0.00
	1.17E+1	1.50	7.79E-1	4.02E-1	2.05E-1	1.04E-1	5.21E-2
order	–	1.48	0.47	0.48	0.49	0.49	0.50

Table 2. Table 3.3: Discrete L 1 superscript 𝐿 1 L^{1} temporal errors e ρ ε , τ ​ ( t = 2 ​ π ) superscript subscript 𝑒 𝜌 𝜀 𝜏 𝑡 2 𝜋 e_{\rho}^{\varepsilon,\tau}(t=2\pi) for the probability density of the NLDE ( 2.2 ) with resonant time step size, S 2 subscript 𝑆 2 S_{2} method.

$e_{\rho }^{\epsilon ,\tau }\left(t=2\pi \right)$	${\tau }_0=\pi /4$	${\tau }_0/4$	${\tau }_0/4^2$	${\tau }_0/4^3$	${\tau }_0/4^4$	${\tau }_0/4^5$	${\tau }_0/4^6$
${\epsilon }_0=1$	1.79	3.63E-2	2.04E-3	1.27E-4	7.94E-6	4.96E-7	3.11E-8
order	–	2.81	2.08	2.00	2.00	2.00	2.00
${\epsilon }_0/2$	1.09	4.94E-2	1.56E-3	9.66E-5	6.03E-6	3.77E-7	2.37E-8
order	–	2.23	2.49	2.01	2.00	2.00	2.00
${\epsilon }_0/2^2$	1.37	4.68E-1	2.86E-3	1.61E-4	9.97E-6	6.23E-7	3.87E-8
order	–	0.77	3.68	2.08	2.00	2.00	2.01
${\epsilon }_0/2^3$	1.06	3.87E-1	2.97E-1	3.05E-4	1.74E-5	1.08E-6	6.72E-8
order	–	0.73	0.19	4.96	2.07	2.00	2.00
${\epsilon }_0/2^4$	9.00E-1	2.05E-1	1.89E-1	1.70E-1	3.50E-5	2.01E-6	1.25E-7
order	–	1.07	0.06	0.08	6.12	2.06	2.00
${\epsilon }_0/2^5$	8.28E-1	1.13E-1	9.58E-2	9.49E-2	9.02E-2	4.20E-6	2.43E-7
order	–	1.44	0.12	0.01	0.04	7.20	2.06
${\epsilon }_0/2^9$	7.66E-1	3.01E-2	7.02E-3	6.02E-3	5.97E-3	5.96E-3	5.96E-3
order	–	2.33	1.05	0.11	0.01	0.00	0.00
${\epsilon }_0/2^{13}$	7.63E-1	2.67E-2	1.84E-3	4.39E-4	3.76E-4	3.73E-4	3.73E-4
order	–	2.42	1.93	1.03	0.11	0.01	0.00
${\epsilon }_0/2^{17}$	7.62E-1	2.65E-2	1.62E-3	1.14E-4	2.60E-5	2.27E-5	2.27E-5
order	–	2.42	2.02	1.92	1.06	0.10	0.00
	1.79	4.68E-1	2.97E-1	1.70E-1	9.02E-2	4.64E-2	2.35E-2
order	–	0.97	0.33	0.40	0.46	0.48	0.49

Table 3. Table 3.4: Discrete relative L 1 superscript 𝐿 1 L^{1} temporal errors e 𝐉 ε , τ ​ ( t = 2 ​ π ) superscript subscript 𝑒 𝐉 𝜀 𝜏 𝑡 2 𝜋 e_{\mathbf{J}}^{\varepsilon,\tau}(t=2\pi) for the current density of the NLDE ( 2.2 ) with resonant time step size, S 2 subscript 𝑆 2 S_{2} method.

$e_{𝐉}^{\epsilon ,\tau }\left(t=2\pi \right)$	${\tau }_0=\pi /4$	${\tau }_0/4$	${\tau }_0/4^2$	${\tau }_0/4^3$	${\tau }_0/4^4$	${\tau }_0/4^5$	${\tau }_0/4^6$
${\epsilon }_0=1$	7.11E-1	1.47E-2	8.30E-4	5.16E-5	3.22E-6	2.02E-7	1.26E-8
order	–	2.80	2.07	2.00	2.00	2.00	2.00
${\epsilon }_0/2$	5.93E-1	2.55E-2	8.37E-4	5.18E-5	3.23E-6	2.02E-7	1.27E-8
order	–	2.27	2.46	2.01	2.00	2.00	2.00
${\epsilon }_0/2^2$	5.71E-1	3.34E-1	1.74E-3	9.99E-5	6.22E-6	3.88E-7	2.41E-8
order	–	0.39	3.79	2.06	2.00	2.00	2.00
${\epsilon }_0/2^3$	4.14E-1	2.19E-1	2.06E-1	1.98E-4	1.15E-5	7.18E-7	4.47E-8
order	–	0.46	0.05	5.01	2.05	2.00	2.00
${\epsilon }_0/2^4$	3.58E-1	1.17E-1	1.16E-1	1.13E-1	2.36E-5	1.38E-6	8.56E-8
order	–	0.81	0.01	0.02	6.11	2.05	2.00
${\epsilon }_0/2^5$	3.46E-1	6.07E-2	5.95E-2	5.95E-2	5.88E-2	2.90E-6	1.69E-7
order	–	1.26	0.01	0.00	0.01	7.16	2.05
${\epsilon }_0/2^9$	3.42E-1	1.28E-2	3.85E-3	3.81E-3	3.81E-3	3.81E-3	3.81E-3
order	–	2.37	0.86	0.01	0.00	0.00	0.00
${\epsilon }_0/2^{13}$	3.42E-1	1.24E-2	7.76E-4	2.41E-4	2.38E-4	2.38E-4	2.38E-4
order	–	2.39	2.00	0.84	0.01	0.00	0.00
${\epsilon }_0/2^{17}$	3.42E-1	1.24E-2	7.51E-4	4.68E-5	1.35E-5	1.37E-5	1.37E-5
order	–	2.39	2.02	2.00	0.90	-0.01	0.00
	7.11E-1	3.34E-1	2.06E-1	1.13E-1	5.88E-2	3.00E-2	1.51E-2
order	–	0.55	0.35	0.43	0.47	0.49	0.49

Table 4. Table 3.5: Relative temporal errors e E ε , τ ​ ( t = 2 ​ π ) superscript subscript 𝑒 𝐸 𝜀 𝜏 𝑡 2 𝜋 e_{E}^{\varepsilon,\tau}(t=2\pi) for the energy of the NLDE ( 2.2 ) with resonant time step size, S 2 subscript 𝑆 2 S_{2} method.

$e_E^{\epsilon ,\tau }\left(t=2\pi \right)$	${\tau }_0=\pi /4$	${\tau }_0/4$	${\tau }_0/4^2$	${\tau }_0/4^3$	${\tau }_0/4^4$	${\tau }_0/4^5$	${\tau }_0/4^6$
${\epsilon }_0=1$	1.30E-1	1.94E-3	1.10E-4	6.86E-6	4.29E-7	2.69E-8	1.78E-9
order	–	3.03	2.07	2.00	2.00	2.00	1.96
${\epsilon }_0/2$	3.29E-2	2.16E-3	7.02E-5	4.27E-6	2.67E-7	1.67E-8	1.12E-9
order	–	1.96	2.47	2.02	2.00	2.00	1.95
${\epsilon }_0/2^2$	2.01E-2	2.53E-2	2.54E-4	1.36E-5	8.44E-7	5.25E-8	3.10E-9
order	–	-0.17	3.32	2.11	2.01	2.00	2.04
${\epsilon }_0/2^3$	3.20E-2	1.50E-3	9.21E-3	3.37E-5	1.90E-6	1.18E-7	7.07E-9
order	–	2.21	-1.31	4.05	2.08	2.01	2.03
${\epsilon }_0/2^4$	4.05E-2	4.25E-4	1.65E-3	3.32E-3	4.29E-6	2.45E-7	1.51E-8
order	–	3.29	-0.98	-0.50	4.80	2.06	2.01
${\epsilon }_0/2^5$	4.45E-2	1.50E-3	7.52E-4	8.92E-4	1.29E-3	5.35E-7	3.09E-8
order	–	2.45	0.50	-0.12	-0.27	5.62	2.06
${\epsilon }_0/2^9$	4.65E-2	2.51E-3	1.05E-4	4.42E-5	5.35E-5	5.41E-5	5.42E-5
order	–	2.11	2.29	0.62	-0.14	-0.01	0.00
${\epsilon }_0/2^{13}$	4.66E-2	2.57E-3	1.55E-4	6.54E-6	2.75E-6	3.33E-6	3.36E-6
order	–	2.09	2.02	2.28	0.63	-0.14	-0.01
${\epsilon }_0/2^{17}$	4.66E-2	2.57E-3	1.59E-4	1.03E-5	1.03E-6	4.49E-7	4.49E-7
order	–	2.09	2.01	1.97	1.66	0.60	0.00
	1.30E-1	2.53E-2	9.21E-3	3.32E-3	1.29E-3	5.44E-4	2.45E-4
order	–	1.18	0.73	0.74	0.68	0.62	0.58

Equations

$$i{\partial }_t\Phi =\left(-\frac{i}{\epsilon }\sum _{j=1}^d{\sigma }_j{\partial }_j+\frac{1}{{\epsilon }^2}{\sigma }_3\right)\Phi +V(\mathbf{x})\Phi +\mathbf{F}(\Phi )\Phi ,\mathbf{x}\in {\mathbb{R}}^d,d=1,2,t>0,$$

$$\sigma_{1}=\left(\begin{array}[]{cc}0&1\\ 1&0\\ \end{array}\right),\qquad\sigma_{2}=\left(\begin{array}[]{cc}0&-i\\ i&0\\ \end{array}\right),\qquad\sigma_{3}=\left(\begin{array}[]{cc}1&0\\ 0&-1\\ \end{array}\right).$$

$$\mathbf{F}(\Phi )={\lambda }_1({\Phi }^{\ast }{\sigma }_3\Phi ){\sigma }_3+{\lambda }_2|\Phi {|}^2{I}_2,$$

$$\Phi (t=0,\mathbf{x})={\Phi }_0(\mathbf{x}),\mathbf{x}\in {\mathbb{R}}^d,d=1,2.$$

$$Q^{\epsilon }=-i\epsilon {\sigma }_1{\partial }_x+{\sigma }_3,x\in \mathbb{R},$$

$$i{\partial }_t\Phi (t,x)=\frac{1}{{\epsilon }^2}{Q}^{\epsilon }\Phi (t,x)+V(x)\Phi (t,x)+\mathbf{F}(\Phi (t,x))\Phi (t,x),x\in \mathbb{R},$$

$${\Phi }^{n+1}(x)=e^{-\frac{i\tau }{{\epsilon }^2}{Q}^{\epsilon }}{e}^{-i\tau \left[V(x)+\mathbf{F}({\Phi }^n(x))\right]}{\Phi }^n(x),\text{with }{\Phi }^0(x)={\Phi }_0(x),x\in \mathbb{R}.$$

$$(A)V(x)\in W^{2m+1,\infty }(\mathbb{R}),m\in {\mathbb{N}}^{\ast }.$$

$$(B)\Phi (t,x)\in L^{\infty }([0,T];(H^{2m+1}(\mathbb{R}){)}^2),m\in {\mathbb{N}}^{\ast }.$$

$${\mathbf{e}}^n(x)=\Phi (t_n,x)-{\Phi }^n(x),0\le n\le \frac{T}{\tau },$$

$$\Vert \Phi {\Vert }_{H^1}^2=\Vert \Phi {\Vert }_{L^2}^2+\Vert {\partial }_x\Phi {\Vert }_{L^2}^2,$$

$$\Vert {\mathbf{e}}^n(x){\Vert }_{H^1}\lesssim \tau +\epsilon ,\Vert {\mathbf{e}}^n(x){\Vert }_{H^1}\lesssim \tau +\tau /\epsilon ,0\le n\le \frac{T}{\tau }.$$

$$\Vert {\mathbf{e}}^n(x){\Vert }_{H^1}\lesssim \tau +\underset{0<\epsilon \le 1}{\max }\min \{\epsilon ,\tau /\epsilon \}\lesssim \sqrt{\tau },0\le n\le \frac{T}{\tau }.$$

$$\Vert S_{n,\tau }^{\mathrm{Lie}}(\Phi )-S_{n,\tau }^{\mathrm{Lie}}(\Psi ){\Vert }_{H^1}\le e^{c_1\tau }\Vert \Phi -\Psi {\Vert }_{H^1},$$

$$M_1=\underset{\epsilon \in (0,1]}{\sup }\Vert \Phi (t,x){\Vert }_{L^{\infty }([0,T];(H^1(\mathbb{R}){)}^2)}.$$

$$Q^{\epsilon }=\sqrt{Id-{\epsilon }^2\Delta }{\Pi }_{+}^{\epsilon }-\sqrt{Id-{\epsilon }^2\Delta }{\Pi }_{-}^{\epsilon },$$

$${\Pi }_{+}^{\epsilon }=\frac{1}{2}\left[Id+(Id-{\epsilon }^2\Delta {)}^{-1/2}{Q}^{\epsilon }\right],{\Pi }_{-}^{\epsilon }=\frac{1}{2}\left[Id-(Id-{\epsilon }^2\Delta {)}^{-1/2}{Q}^{\epsilon }\right].$$

$${\Pi }_{\pm }^{\epsilon }={\Pi }_{\pm }^0\pm \epsilon {\mathcal{R}}_1={\Pi }_{\pm }^0\mp i\frac{\epsilon }{2}{\sigma }_1{\partial }_x\pm {\epsilon }^2{\mathcal{R}}_2,{\Pi }_{+}^0=\text{diag}(1,0),{\Pi }_{-}^0=\text{diag}(0,1),$$

$${\Phi }_{\pm }^{\epsilon }(t,x):={\Pi }_{\pm }^{\epsilon }\Phi (t,x).$$

$${\mathcal{D}}^{\epsilon }=\frac{1}{{\epsilon }^2}(\sqrt{Id-{\epsilon }^2\Delta }-Id)=-(\sqrt{Id-{\epsilon }^2\Delta }+Id{)}^{-1}\Delta ,$$

$$e^{\frac{it}{{\epsilon }^2}{Q}^{\epsilon }}=e^{\frac{it}{{\epsilon }^2}(\sqrt{Id-{\epsilon }^2\Delta }{\Pi }_{+}^{\epsilon }-\sqrt{Id-{\epsilon }^2\Delta }{\Pi }_{-}^{\epsilon })}=e^{\frac{it}{{\epsilon }^2}}{e}^{it{\mathcal{D}}^{\epsilon }}{\Pi }_{+}^{\epsilon }+e^{-\frac{it}{{\epsilon }^2}}{e}^{-it{\mathcal{D}}^{\epsilon }}{\Pi }_{-}^{\epsilon }.$$

$${\mathbf{e}}^{n+1}(x)=e^{-\frac{i\tau }{{\epsilon }^2}{Q}^{\epsilon }}{e}^{-i\tau \mathbf{F}({\Phi }^n)}{\mathbf{e}}^n(x)+{\eta }_1^n(x)+e^{-\frac{i\tau }{{\epsilon }^2}{Q}^{\epsilon }}{\eta }_2^n(x),$$

$$f_2^n(s)=$$

$$-i{e}^{\frac{-i2s}{{\epsilon }^2}}\left[{\Pi }_{+}^{\epsilon }\left({\mathbf{g}}_1^n(x){\Phi }_{+}^{\epsilon }(t_n)\right)+{\Pi }_{-}^{\epsilon }\left({\mathbf{g}}_2^n(x){\Phi }_{+}^{\epsilon }(t_n)+{\mathbf{g}}_1^n(x){\Phi }_{-}^{\epsilon }(t_n)\right)\right]$$

$$-i{e}^{\frac{2is}{{\epsilon }^2}}\left[{\Pi }_{-}^{\epsilon }\left(\overline{{\mathbf{g}}_1^n}(x){\Phi }_{-}^{\epsilon }(t_n)\right)+{\Pi }_{+}^{\epsilon }\left({\mathbf{g}}_2^n(x){\Phi }_{-}^{\epsilon }(t_n)+\overline{{\mathbf{g}}_1^n}(x){\Phi }_{+}^{\epsilon }(t_n)\right)\right],$$

$${\mathbf{g}}_1({\Phi }_{+}^{\epsilon }(t_n),{\Phi }_{-}^{\epsilon }(t_n))={\lambda }_1\left(({\Phi }_{-}^{\epsilon }(t_n){)}^{\ast }{\sigma }_3{\Phi }_{+}(t_n)\right){\sigma }_3+{\lambda }_2\left(({\Phi }_{-}^{\epsilon }(t_n){)}^{\ast }{\Phi }_{+}^{\epsilon }(t_n)\right)I_2,$$

$${\mathbf{g}}_2({\Phi }_{+}^{\epsilon }(t_n),{\Phi }_{-}^{\epsilon }(t_n))=\sum _{\sigma =\pm }\left[{\lambda }_1(({\Phi }_{\sigma }^{\epsilon }(t_n){)}^{\ast }{\sigma }_3{\Phi }_{\sigma }^{\epsilon }(t_n)){\sigma }_3+{\lambda }_2|{\Phi }_{\sigma }^{\epsilon }(t_n){|}^2{I}_2\right].$$

$${\mathbf{e}}^{n+1}(x)=e^{-\frac{i\tau }{{\epsilon }^2}{Q}^{\epsilon }}{e}^{-i\tau \mathbf{F}({\Phi }^n)}{\mathbf{e}}^n(x)+{\eta }^n(x),0\le n\le \frac{T}{\tau }-1,x\in \mathbb{R},$$

$${\eta }^n(x)=\Phi (t_{n+1},x)-e^{-\frac{i\tau }{{\epsilon }^2}{Q}^{\epsilon }}{e}^{-i\tau \mathbf{F}({\Phi }^n)}\Phi (t_n,x),x\in \mathbb{R}.$$

$$\Phi (t_n+s,x)=e^{-\frac{is}{{\epsilon }^2}{Q}^{\epsilon }}\Phi (t_n,x)-i{\int }_0^s{e}^{-\frac{i(s-w)}{{\epsilon }^2}{Q}^{\epsilon }}\mathbf{F}(\Phi (t_n+w,x))\Phi (t_n+w,x)dw,0\le s\le \tau ,$$

$${\eta }^n(x)=e^{-\frac{i\tau }{{\epsilon }^2}{Q}^{\epsilon }}\left({\int }_0^{\tau }{f}^n(s)ds-\tau f^n(0)\right)+R_1^n(x)+R_2^n(x),$$

$$f^n(s)=-i{e}^{\frac{is}{{\epsilon }^2}{Q}^{\epsilon }}\left(\mathbf{F}(e^{-\frac{is}{{\epsilon }^2}{Q}^{\epsilon }}\Phi (t_n))e^{-\frac{is}{{\epsilon }^2}{Q}^{\epsilon }}\Phi (t_n,x)\right),R_1^n(x)=e^{-\frac{i\tau }{{\epsilon }^2}{Q}^{\epsilon }}\left({\Lambda }_1^n(x)+{\Lambda }_2^n(x)\right),$$

Full text

1 Introduction

The splitting methods form an important group of methods which are quite accurate and efficient [60]. Actually, they have been widely applied for dealing with highly oscillatory systems such as the Schrödinger/nonlinear Schrödinger equations [1, 8, 9, 24, 25, 58, 70], the Dirac/nonlinear Dirac equations [5, 6, 14, 57], the Maxwell-Dirac system [10, 52], the Zakharov system [12, 13, 44, 53], the Gross-Pitaevskii equation for Bose-Einstein condensation (BEC) [11], the Stokes equation [23], and the Enrenfest dynamics [35], etc.

In this paper, we consider the splitting methods applied to the nonlinear Dirac equation (NLDE) [30, 31, 40, 43, 50, 27, 46, 47, 36, 37, 38, 39, 64, 66, 73] in the nonrelativistic 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 $\Phi:=\Phi(t,{\bf x})=(\phi_{1}(t,{\bf x}),\phi_{2}(t,{\bf x}))^{T}\in\mathbb{C}^{2}$ [6]:

[TABLE]

where $i=\sqrt{-1}$ is the imaginary unit, $t$ is time, ${\bf x}=(x_{1},...,x_{d})^{T}$, $\partial_{j}=\frac{\partial}{\partial x_{j}}$ ($j=1,...,d$), $\varepsilon\in(0,1]$ is a dimensionless parameter inversely proportional to the speed of light, and $V:=V({\bf x})$ is a real-valued function denoting the external electric potential. $\sigma_{1}$, $\sigma_{2}$, $\sigma_{3}$ are the Pauli matrices defined as

[TABLE]

The nonlinearity $\mathbf{F}(\Phi)$ in (1.1) is usually taken as

[TABLE]

with $|\Phi|^{2}=\Phi^{*}\Phi$, where $\lambda_{1}$, $\lambda_{2}\in\mathbb{R}$ are two given real constants, $\Phi^{*}=\overline{\Phi}^{T}$ is the complex conjugate transpose of $\Phi$ and $I_{2}$ is the $2\times 2$ identity matrix. The above choice of nonlinearity is motivated from the so-called Soler model in quantum field theory, e.g. $\lambda_{2}=0$ and $\lambda_{1}\neq 0$ [40, 43, 71], and BEC with a chiral confinement and/or spin-orbit coupling, e.g. $\lambda_{1}=0$ and $\lambda_{2}\neq 0$ [27, 46, 47]. In order to study the dynamics, the initial data is chosen as

[TABLE]

When $\varepsilon=1$ 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, 26, 32, 33, 34, 54] and references therein. Particularly, for the case where $d=1$, $V(x)\equiv 0$, $\lambda_{1}=-1$ and $\lambda_{2}=0$ in the choice of $\mathbf{F}(\Phi)$, the NLDE (1.1) admits explicit soliton solutions [28, 43, 48, 55, 59, 63, 68, 69]. 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, 49, 62], the time-splitting Fourier spectral (TSFP) methods [10, 18, 42, 52] and the Runge-Kutta discontinuous Galerkin methods [51].

On the other hand, when $0<\varepsilon\ll 1$ (the nonrelativistic regime where the wave speed is much smaller than the speed of light), as indicated by previous analysis in [6, 41, 61, 20], the wavelength of the solution in time is at $O(\varepsilon^{2})$. The oscillation of the solution as well as the unbounded and indefinite energy functional w.r.t. $\varepsilon$ [16, 34] cause much burden in the analysis and computation. Indeed, it would require that the time step size $\tau$ to be strictly reliant on $\varepsilon$ to capture the exact solution. Numerical studies in [6] have confirmed this dependence. The error bounds show that $\tau=O(\varepsilon^{3})$ is required for the conservative Crank-Nicolson finite difference (CNFD) method [6], and $\tau=O(\varepsilon^{2})$ is required for the exponential wave integrator Fourier pseudospectral (EWI-FP) method as well as the time-splitting Fourier pseudospectral (TSFP) method [6]. To overcome the restriction, recently, uniform accurate (UA) schemes with two-scale formulation approach [56] or multiscale time integrator pseudospectral method [4, 22] or nested Picard iterative integrators [21] have been designed for the NLDE in the nonrelativistic regime, where the time step size $\tau$ could be independent of $\varepsilon$.

Though the error of the TSFP method (also called $S_{2}$ later in this paper) has a $\tau^{2}/\varepsilon^{4}$ dependence on the small parameter $\varepsilon$ [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 error of $S_{2}$ is independent of $\varepsilon$ and uniform w.r.t. $\varepsilon$. In other words, $S_{2}$ for the NLDE (1.1) in the absence of magnetic potentials displays super-resolution w.r.t. $\varepsilon$.

The super-resolution here suggests independence of the oscillation wavelength. It is even stronger than the ‘super-resolution’ in [29] for the Schrödinger equation in the semiclassical regime, where the restriction on the time steps is still related to the wavelength, but not so strict as the resolution of the oscillation by fixed number of points per wavelength. This 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 ($S_{1}$) and second-order ($S_{2}$) 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 [58], 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 adopted 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 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, $A\lesssim B$ is used in the sense that there exists a generic constant $C>0$ independent of $\varepsilon$ and $\tau$, such that $|A|\leq C\,B$. $A\lesssim_{\delta}B$ has a similar meaning that there exists a generic constant $C_{\delta}>0$ dependent on $\delta$ but independent of $\varepsilon$ and $\tau$, such that $|A|\leq C_{\delta}\,B$.

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 ($d=1$). Extensions to (1.1) in 2D and/or the four component form of the NLDE with $d=1,2,3$ [6] are straightforward.

Denote the free Dirac Hermitian operator

[TABLE]

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

[TABLE]

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

Choose $\tau>0$ as the time step size and $t_{n}=n\tau$ for $n=0,1,...$ as the time steps. Denote $\Phi^{n}(x)$ to be the numerical approximation of $\Phi(t_{n},x)$, where $\Phi(t,x)$ 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) [72], $S_{1}$ can be represented as [6]:

[TABLE]

For simplicity, we also write $\Phi^{n+1}(x):=S_{n,\tau}^{\text{Lie}}(\Phi^{n})$, where $S_{n,\tau}^{\text{Lie}}$ denotes the numerical propagator of the Lie-Trotter splitting.

2.1 A uniform error bound

For any $0<T<T^{*}$, where $T^{*}$ denotes the common maximal existence time of the solution for (1.1) with (1.3) and (1.4) for all $0<\varepsilon\leq 1$, we are going to consider smooth solutions, i.e. we assume the electric potential satisfies

[TABLE]

In addition, we assume the exact solution $\Phi(t,x)$ satisfies

[TABLE]

For the numerical approximation $\Phi^{n}(x)$ obtained from $S_{1}$ (2.3), we introduce the error function

[TABLE]

then the following uniform error bound in $H^{1}$ norm can be established, where the $H^{1}$ norm for function $\Phi(x)=(\phi_{1},\phi_{2})^{T}\in\mathbb{C}^{2}$ is given by

[TABLE]

with $L^{2}$ norm defined as $\|\Phi\|_{L^{2}}=\sqrt{\int_{\mathbb{R}}|\Phi(x)|^{2}\,dx}=\sqrt{\int_{\mathbb{R}}\left(|\phi_{1}(x)|^{2}+|\phi_{2}(x)|^{2}\right)\,dx}$.

Theorem 1**.**

Let $\Phi^{n}(x)$ be the numerical approximation obtained from $S_{1}$ (2.3), then under assumptions $(A)$ and $(B)$ with $m=1$, there exists $0<\tau_{0}\leq 1$ independent of $\varepsilon$ such that the following two error estimates hold for $0<\tau<\tau_{0}$

[TABLE]

Consequently, there is a uniform error bound for $S_{1}$ when $0<\tau<\tau_{0}$

[TABLE]

Remark 2.1**.**

Instead of proving the $L^{2}$ error bounds as in the linear case, in Theorem 1 and the other results in this paper for the 1D problem, we prove the $H^{1}$ error bounds for ${\bf e}^{n}(x)$ due to the fact that $H^{1}(\mathbb{R})$ is an algebra, and the corresponding estimates should be in $H^{2}$ norm for 2D and 3D cases (2D case in the sense of (1.1), and 3D case in the sense of the four-component nonlinear Dirac equation given in [6]) with of course higher regularity assumptions (higher order Sobolev norm estimates need higher regularity of the exact solution).

Remark 2.2**.**

In Theorem 1, the $H^{3}$ regularity ($m=1$ in assumptions (A) and (B)) is assumed for the first order Lie splitting scheme, and this regularity assumption is sharp for the results stated in the theorem. Heuristically, the estimates of the type $\|{\bf e}^{n}(x)\|_{H^{1}}\lesssim\tau+\varepsilon$ hold for $\varepsilon\in(0,1]$, while in the limit $\varepsilon\to 0^{+}$, the NLDE (1.1) converges to the coupled nonlinear Schrödinger equations (CNLSE) after filtering out the nonrelativistic temporal oscillations [6, 41, 61, 20]. Thus, letting $\varepsilon\to 0^{+}$, the estimates $\|{\bf e}^{n}(x)\|_{H^{1}}\lesssim\tau+\varepsilon$ will become the error bounds for the Lie splitting method applied to CNLSE. For the $H^{1}$ error estimates of the Lie splitting in the case of Schrödinger type equations, the regularity requirement of the exact solution should be 2 orders higher [58] (one oder temporal derivative corresponds to two order spatial derivative in Schrödinger type equations), i.e. $H^{3}$ regularity of the exact solution is needed.

For simplicity of the presentation, in the proof for this theorem and other theorems later for NLDE in this paper, we take $V(x)\equiv 0$. Extension to the case where $V(x)\neq 0$ is straightforward [7]. Compared to the linear case [7], the nonlinear term is much more complicated to analyze. As discussed earlier in the introduction, for the linear Dirac equation, the linearity and $L^{2}$ unitary property of the numerical propagator enable the explicit expression (exact) of the error ${\bf e}^{n}(x)$ by the local error (see Lemma 3) without any extra condition on time step $\tau$. Therefore, the error estimates (in $L^{2}$ norm) in [7] are obtained by carefully studying the accumulation of the local errors. However, for the nonlinear Dirac equation case, the approach (highly depend on the linear property) in [7] fails. Different from the linear case, the novelty of the strategy we adopt for the nonlinear case lies in the following aspects: (i) carefully carry out expansions of the nonlinear terms to analyze the local errors and identify the leading temporal oscillations; and (ii) estimate the errors in $H^{1}$ norm where the conditional stabilities of the numerical propagators (see Lemma 2) and the NLDE (1.1) hold, and then control the nonlinear terms by mathematical induction, the uniform error estimates (2.7) and Sobolev inequalities (see (2.35) and the proof after). We emphasis here that our analysis and convergence rate results are valid for the linear Dirac equation, while the approach and error estimates in the linear case [7] can not be applied here for the nonlinear case. Of course, the dependence of the constant on time $T$ in front of the convergence rate is sharper in the linear case in [7] than that in Theorem 2.1.

As mentioned above, a key issue of the error analysis for NLDE is to control the nonlinear term of numerical solution $\Phi^{n}$, and for which we require the following stability lemma [58].

Lemma 2**.**

Suppose $V(x)\in W^{1,\infty}(\mathbb{R})$, and $\Phi(x),\Psi(x)\in(H^{1}(\mathbb{R}))^{2}$ satisfy $\|\Phi\|_{H^{1}},\|\Psi\|_{H^{1}}\leq M$, we have

[TABLE]

where $c_{1}$ depends on $M$ and $\|V(x)\|_{W^{1,\infty}}$.

Proof.

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

Under the assumption (B) ($m\geq 1$), for $\varepsilon\in(0,1]$, we denote $M_{1}>0$ as

[TABLE]

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

The following properties of $Q^{\varepsilon}$ will be frequently used in the analysis. $Q^{\varepsilon}$ is diagonalizable in the phase space (Fourier domain) and can be decomposed as

[TABLE]

where $\Delta=\partial_{xx}$ is the Laplace operator in 1D, $Id$ is the identity operator, and $\Pi_{+}^{\varepsilon}$, $\Pi_{-}^{\varepsilon}$ are projectors defined as

[TABLE]

It is straightforward to verify that $\Pi_{+}^{\varepsilon}+\Pi_{-}^{\varepsilon}=Id$, $\quad\Pi_{+}^{\varepsilon}\Pi_{-}^{\varepsilon}=\Pi_{-}^{\varepsilon}\Pi_{+}^{\varepsilon}=0$, $\quad(\Pi_{\pm}^{\varepsilon})^{2}=\Pi_{\pm}^{\varepsilon}$, and through Taylor expansion, we have [16]

[TABLE]

with $\mathcal{R}_{1}:(H^{m}(\mathbb{R}))^{2}\rightarrow(H^{m-1}(\mathbb{R}))^{2}$ for $m\geq 1$, $\mathcal{R}_{2}:(H^{m}(\mathbb{R}))^{2}\rightarrow(H^{m-2}(R))^{2}$ for $m\geq 2$ being uniformly bounded operators w.r.t. $\varepsilon$. For simplicity of expression, we denote

[TABLE]

In order to characterize the oscillatory features of the solution, noticing $(\sqrt{Id-\varepsilon^{2}\Delta}-Id)(\sqrt{Id-\varepsilon^{2}\Delta}+Id)=-\varepsilon^{2}\Delta$, we denote

[TABLE]

which is a uniformly bounded operator w.r.t $\varepsilon$ from $(H^{m}(\mathbb{R}))^{2}\rightarrow(H^{m-2}(\mathbb{R}))^{2}$ for $m\geq 2$, then the evolution operator $e^{\frac{it}{\varepsilon^{2}}Q^{\varepsilon}}$ can be expressed as

[TABLE]

For simplicity, here we use $\Phi(t):=\Phi(t,x)$, $\Phi^{n}:=\Phi^{n}(x)$ in short.

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

Lemma 3**.**

Let $\Phi^{n}(x)$ ( $0\leq n\leq\frac{T}{\tau}-1$) be obtained from $S_{1}$ (2.3) satisfying $\|\Phi^{n}(x)\|_{H^{1}}\leq M_{1}+1$, under the assumptions of Theorem 1, we have

[TABLE]

with $\|\eta_{1}^{n}(x)\|_{H^{1}}\leq c_{1}\tau^{2}+c_{2}\tau\|{\bf e}^{n}(x)\|_{H^{1}}$, $\eta_{2}^{n}(x)=\int_{0}^{\tau}f_{2}^{n}(s)ds-\tau f_{2}^{n}(0)$, where $c_{1}$ depends on $M_{1}$, $\lambda_{1}$, $\lambda_{2}$ and $\|\Phi(t,x)\|_{L^{\infty}([0,T];(H^{3})^{2})}$; $c_{2}$ depends on $M_{1}$, $\lambda_{1}$, and $\lambda_{2}$. Here

[TABLE]

where $\mathbf{g}_{j}^{n}(x)=\mathbf{g}_{j}(\Phi_{+}^{\varepsilon}(t_{n}),\Phi_{-}^{\varepsilon}(t_{n}))$ and

[TABLE]

Proof.

Through the definition of ${\bf e}^{n}(x)$ (2.4), noticing the formula (2.3), we have

[TABLE]

where $\eta^{n}(x)$ is the “local truncation error” (notice that this is not the usual local truncation error, compared with $\Phi(t_{n+1},x)-S_{n,\tau}^{\text{Lie}}\Phi(t_{n},x)$),

[TABLE]

By Duhamel’s principle, the solution $\Phi(t,x)$ to (2.2) satisfies

[TABLE]

which implies that $\|\Phi(t_{n}+s,x)-e^{-\frac{is}{\varepsilon^{2}}Q^{\varepsilon}}\Phi(t_{n},x)\|_{H^{1}}\lesssim\tau$ ($s\in[0,\tau]$). Setting $s=\tau$ in (2.22), we have from (2.21),

[TABLE]

where

[TABLE]

with

[TABLE]

Noticing (2.9), (2.22), and the fact that $e^{-isQ^{\varepsilon}/\varepsilon^{2}}$ preserves $H^{k}$ norm, it is not difficult to find

[TABLE]

On the other hand, from the definition of F and the fact that $H^{1}(\mathbb{R})$ is an algebra, we have for any $\Phi_{j}(x)=(\phi_{j1}(x),\phi_{j2}(x))^{T}\in\mathbb{C}^{2}$, $j=1,2,3$

[TABLE]

Having the above inequality, using the assumption that $\|\Phi^{n}\|_{H^{1}}\leq M_{1}+1$, and the Taylor expansion in $\Lambda_{1}^{n}(x)$, we get

[TABLE]

It remains to estimate the $f^{n}(s)$ part. Using the decomposition (2.15) and the Taylor exapnsion $e^{i\tau\mathcal{D}^{\varepsilon}}=Id+O(\tau\mathcal{D}^{\varepsilon})$ (in the sense of phase space), we have $e^{\frac{-isQ^{\varepsilon}}{\varepsilon^{2}}}\Phi(t_{n})=e^{\frac{-is}{\varepsilon^{2}}}\Phi_{+}^{\varepsilon}(t_{n})+e^{\frac{is}{\varepsilon^{2}}}\Phi_{-}^{\varepsilon}(t_{n})+O(s)$,

[TABLE]

where for $s\in[0,\tau]$,

[TABLE]

Since $\mathbf{F}$ is of polynomial type, by direct computation, we can further simplify (2.30) to get

[TABLE]

where $f_{2}^{n}(s)$ is given in (2.17) and $\tilde{f}^{n}(s)$ is independent of $s$ as

[TABLE]

with $\mathbf{g}_{1,2}^{n}$ defined in (2.18)-(2.19).

Now, it is easy to verify that $\eta^{n}(x)=\eta_{1}^{n}(x)+\eta_{2}^{n}(x)$ with $\eta_{2}^{n}(x)$ given in Lemma 3 by choosing

[TABLE]

Noticing that $\tilde{f}^{n}(s)$ is independent of $s$ and $\|f_{1}^{n}(s)\|_{H^{1}}\lesssim\tau$, combining (2.27) and (2.29), we can get

[TABLE]

which completes the proof of Lemma 3. □

Now, we proceed to prove Theorem 1.

Proof.

We will prove by induction that the estimates (2.6)-(2.7) hold for all time steps $n\leq\frac{T}{\tau}$ together with

[TABLE]

Since initially $\Phi^{0}=\Phi_{0}(x)$, $n=0$ case is obvious. Assume (2.6)-(2.7) and (2.35) hold true for all $0\leq n\leq p\leq\frac{T}{\tau}-1$, then we are going to prove the case $n=p+1$.

From Lemma 3, we have

[TABLE]

with $\|\eta_{1}^{n}(x)\|_{H^{1}}\lesssim\tau^{2}+\tau\|{\bf e}^{n}(x)\|_{H^{1}}$, ${\bf e}^{0}=0$ and $\eta_{2}^{n}(x)$ given in Lemma 3.

Denote $\mathcal{L}_{n}=e^{-\frac{i\tau}{\varepsilon^{2}}Q^{\varepsilon}}\left(e^{-i\tau\mathbf{F}(\Phi^{n})}-I_{2}\right)$ ($0\leq n\leq p\leq\frac{T}{\tau}-1$), and it is straightforward to calculate

[TABLE]

with $C_{M_{1}}$ only depending on $M_{1}$. Thus we can obtain from (2.36) that for $0\leq n\leq p$,

[TABLE]

Since $\|\eta_{1}^{k}(x)\|_{H^{1}}\lesssim\tau^{2}+\tau\|{\bf e}^{n}(x)\|_{H^{1}}$, $k=0,1,...,n$, and $e^{-is/\varepsilon^{2}Q^{\varepsilon}}$ ($s\in\mathbb{R}$) preserves $H^{1}$ norm, we have from (2.37)

[TABLE]

which leads to

[TABLE]

To analyze $\eta_{2}^{n}(x)=\int_{0}^{\tau}f_{2}^{n}(s)ds-\tau f_{2}^{n}(0)$, using (2.12), we can find $f_{2}^{n}(s)=O(\varepsilon)$, e.g.

[TABLE]

and the other terms in $f_{2}^{n}(s)$ can be estimated similarly. As $\mathcal{R}_{1}:(H^{m})^{2}\to(H^{m-1})^{2}$ is uniformly bounded with respect to $\varepsilon\in(0,1]$, we have (with detailed computations omitted)

[TABLE]

Noticing the assumptions of Theorem 1, we obtain from (2.17)

[TABLE]

which leads to

[TABLE]

On the other hand, using Taylor expansion and the second inequality in (2.42), we have

[TABLE]

Combining (2.43) and (2.44), we arrive at

[TABLE]

Then from (2.40), we get for $n\leq p$

[TABLE]

Using discrete Gronwall’s inequality, we have

[TABLE]

which shows that (2.6)-(2.7) hold for $n=p+1$. It can be checked that all the constants appearing in the estimates depend only on $M_{1},\lambda_{1},\lambda_{2},T$ and $\|\Phi(t,x)\|_{L^{\infty}([0,T];(H^{3})^{2})}$, and

[TABLE]

for some $C=C(M_{1},\lambda_{1},\lambda_{2},T,\|\Phi(t,x)\|_{L^{\infty}([0,T];(H^{3})^{2})})$. Choosing $\tau\leq\frac{1}{C^{2}}$ will justify (2.35) at $n=p+1$, which finishes the induction process, and the proof for Theorem 1 is completed. □

2.2 An improved error bound for non-resonant time steps

The leading term in the NLDE (2.2) is $\frac{1}{\varepsilon^{2}}\sigma_{3}\Phi$, suggesting that the solution behaves almost periodically in time with periods $2k\pi\varepsilon^{2}$ ($k\in\mathbb{N}^{*}$, the periods of $e^{-i\sigma_{3}/\varepsilon^{2}}$). From numerical results, we observe that $S_{1}$ behave much better than the results in Theorem 1 when $4\tau$ (which is derived from the proof) is not close to the leading temporal oscillation periods $2k\pi\varepsilon^{2}$. In fact, for given $0<\delta\leq 1$, define

[TABLE]

then when $\tau\in\mathcal{A}_{\delta}(\varepsilon)$, i.e., when non-resonant time step sizes are chosen, the errors of $S_{1}$ can be improved. To illustrate $\mathcal{A}_{\delta}(\varepsilon)$ (compared to the linear case [7], the region of the resonant steps $\mathcal{A}_{\delta}^{c}(\varepsilon):=\mathbb{R}^{+}\backslash\mathcal{A}_{\delta}(\varepsilon)$ for fixed $\varepsilon$ are doubled due to the cubic nonlinearity), we show in Figure 2.1 for $\varepsilon=1$ and $\varepsilon=0.5$ with fixed $\delta=0.15$.

For $\tau\in\mathcal{A}_{\delta}(\varepsilon)$, we can derive improved uniform error bounds for $S_{1}$ as follows.

Theorem 4**.**

Let $\Phi^{n}(x)$ be the numerical approximation obtained from $S_{1}$ (2.3). If the time step size $\tau$ is non-resonant, i.e. there exists $0<\delta\leq 1$, such that $\tau\in\mathcal{A}_{\delta}(\varepsilon)$, then under the assumptions $(A)$ and $(B)$ with $m=1$, we have an improved uniform error bound for small enough $\tau>0$

[TABLE]

Proof.

First of all, the assumptions of Theorem 1 are satisfied in Theorem 4, so we can directly use the results of Theorem 1. In particular, the numerical solution $\Phi^{n}$ are bounded in $H^{1}$ as $\|\Phi^{n}\|_{H^{1}}\leq M_{1}+1$ (2.35) and Lemma 3 for local truncation error holds.

We start from (2.40). The improved estimates rely on the cancellation phenomenon for the $\eta_{2}^{k}$ term in (2.40). From Lemma 3, (2.17), (2.18) and (2.19), we can write $\eta_{2}^{k}(x)$ as

[TABLE]

where $\mathcal{R}_{j,\pm}(\Phi_{+}^{\varepsilon},\Phi_{-}^{\varepsilon})$ ($j=2,4,\Phi_{+}^{\varepsilon},\Phi_{-}^{\varepsilon}:\mathbb{R}\to\mathbb{C}^{2}$) are as follows

[TABLE]

with $\mathbf{g}_{1},\mathbf{g}_{2}$ given in (2.18)-(2.19) (Lemma 3), and

[TABLE]

It is obvious that $|p_{1}(\tau)|,|p_{2}(\tau)|\leq 2\tau$ and (2.40) implies that

[TABLE]

To proceed, we introduce $\widetilde{\Phi}_{\pm}^{\varepsilon}(t)$ as

[TABLE]

Since $\Phi(t,x)$ solves the NLDE (1.1) (or (2.2)), noticing the properties of $Q^{\varepsilon}$ as in (2.10) and (2.14) and the $L^{2}$ orthogonal projections $\Pi_{\pm}^{\varepsilon}$, it is straightforward to compute that

[TABLE]

and the assumptions of Theorem 1 would yield

[TABLE]

Now, we can deal with the terms involving $\mathcal{R}_{j,\pm}$ ($j=2,4$) in (2.52).

For $\mathcal{R}_{4,-}$: By direct computation, we get $\mathcal{R}_{4,-}(\Phi_{+}^{\varepsilon}(t_{k}),\Phi_{-}^{\varepsilon}(t_{k}))=e^{-\frac{3it_{k}}{\varepsilon^{2}}}\mathcal{R}_{4,-}(\widetilde{\Phi}_{+}^{\varepsilon}(t_{k}),\widetilde{\Phi}_{-}^{\varepsilon}(t_{k}))$. In view of (2.15) and (2.52), we have for $0\leq k\leq n\leq\frac{T}{\tau}-1$,

[TABLE]

Denoting

[TABLE]

and noticing that $\partial_{t}A(t)=-ie^{-it\mathcal{D}^{\varepsilon}}\mathcal{D}^{\varepsilon}\mathcal{R}_{4,-}(\widetilde{\Phi}_{+}^{\varepsilon}(t),\widetilde{\Phi}_{-}^{\varepsilon}(t))+e^{-it\mathcal{D}^{\varepsilon}}\partial_{t}\mathcal{R}_{4,-}(\widetilde{\Phi}_{+}^{\varepsilon}(t),\widetilde{\Phi}_{-}^{\varepsilon}(t))$, we can derive from (2.57) and the fact that $\mathcal{D}^{\varepsilon}:(H^{m})^{2}\to(H^{m-2})^{2}$ is uniformly bounded w.r.t $\varepsilon$,

[TABLE]

Using (2.60), (2.58), $\|A(t)\|_{L^{\infty}([0,T];(H^{1})^{2})}\lesssim 1$, the property that $e^{it\mathcal{D}^{\varepsilon}}$ preserves $H^{1}$ norm, summation by parts formula and triangle inequality, we have

[TABLE]

with

[TABLE]

For $\tau\in\mathcal{A}_{\delta}(\varepsilon)$ (2.49), we have $|1-e^{-\frac{i4\tau}{\varepsilon^{2}}}|=|2\sin(2\tau/\varepsilon^{2})|\geq 2\delta$ and $|\theta_{k}|\leq\frac{2}{2\delta}=1/\delta$, and (2.61) leads to

[TABLE]

For $\mathcal{R}_{2,-}$: Similar to the case $\mathcal{R}_{4,-}$ (slightly different), it is straightforward to show that

[TABLE]

where

[TABLE]

$B(t)$ and $C(t)$ satisfy the same estimates as $A(t)$ (2.60). Therefore, similar procedure will give

[TABLE]

with $\widetilde{\theta}_{k}=\sum_{j=0}^{k}e^{-\frac{i2j\tau}{\varepsilon^{2}}}=\frac{1-e^{-\frac{i2(k+1)\tau}{\varepsilon^{2}}}}{1-e^{\frac{-i2\tau}{\varepsilon^{2}}}},\quad k\geq 0,\quad\widetilde{\theta}_{-1}=0$. For $\tau\in\mathcal{A}_{\delta}(\varepsilon)$ (2.49), we know $|1-e^{\frac{-i2\tau}{\varepsilon^{2}}}|=|2\sin(\tau/\varepsilon^{2})|\geq|4\sin(2\tau/\varepsilon^{2})|\geq 4\delta$ and $|\widetilde{\theta}_{k}|\leq\frac{2}{4\delta}=2/\delta$, which shows

[TABLE]

For $\mathcal{R}_{4,+}$ and $\mathcal{R}_{2,+}$: It is easy to see that the $\mathcal{R}_{4,+}$ and $\mathcal{R}_{2,+}$ terms in (2.54) can be bounded exactly the same as the $\mathcal{R}_{4,-}$ and $\mathcal{R}_{2,-}$ terms, respectively.

Finally, combining (2.54), (2.63), (2.68) and above observations, we have for $\tau\in\mathcal{A}_{\delta}(\varepsilon)$,

[TABLE]

and discrete Gronwall inequality yields $\|{\bf e}^{n+1}(x)\|_{H^{1}}\lesssim\frac{\tau}{\delta}$ ($0\leq n\leq\frac{T}{\tau}-1$) for small enough $\tau\in\mathcal{A}_{\delta}(\varepsilon)$. The proof is completed. □

2.3 Numerical results

To verify our error bounds in Theorems 1 and 4, we show a numerical example here. In this example and all the numerical examples later, we always use Fourier pseudospectral method for spatial discretization.

As a common practice when applying the Fourier pseudospectral method, in our numerical simulations, we truncate the whole space onto a sufficiently large bounded domain $\Omega=(a,b)$, and assume periodic boundary conditions. The mesh size is chosen as $h:=\triangle x=\frac{b-a}{M}$ with $M$ being an even positive integer. Then the grid points can be denoted as $x_{j}:=a+jh$, for $j=0,1,...,M$.

In this example, we choose the electric potential $V(x)\equiv 0$. For the nonlinearity (1.3), we take $\lambda_{1}=1$, $\lambda_{2}=0$, i.e.

[TABLE]

and the initial data $\Phi_{0}=(\phi_{1},\phi_{2})$ in (1.4) is given as

[TABLE]

As only the temporal errors are concerned in this paper, during the computation, the spatial mesh size is always set to be $h=\frac{1}{16}$ so that the spatial errors are negligible.

We first take resonant time steps, that is, for small enough chosen $\varepsilon$, there is a positive $k_{0}$, such that $\tau=\frac{1}{2}k_{0}\varepsilon^{2}\pi$, to check the error bounds in Theorem 1. The bounded computational domain is taken as $\Omega=(-32,32)$, i.e., $a=-32$ and $b=32$. Because the exact solution is unknown, for comparison, we use a numerical ‘exact’ solution generated by the second-order time-splitting method ($S_{2}$), which will be introduced later, with a very fine time step size $\tau_{e}=2\pi\times 10^{-6}$.

To display the numerical results, we introduce the discrete $H^{1}$ errors of the numerical solution. Let $\Phi^{n}=(\Phi_{0}^{n},\Phi_{1}^{n},...,\Phi_{M-1}^{n},\Phi_{M}^{n})^{T}$ be the numerical solution obtained by a numerical method with given $\varepsilon$, time step size $\tau$ as well as the fine mesh size $h$ at time $t=t_{n}$, and $\Phi(t,x)$ be the exact solution, then the discrete $H^{1}$ error is defined as

[TABLE]

where

[TABLE]

with $\mu_{l}$, $\widehat{\Phi}_{l}^{n}\in\mathbb{C}^{2}$ defined as

[TABLE]

and $\Phi^{\prime}(t_{n},x_{j})$ is defined similarly. Then $e^{\varepsilon,\tau}(t_{n})$ should be close to the $H^{1}$ errors in Theorem 1 for fine spatial mesh sizes $h$.

Table 2.1 shows the temporal errors $e^{\varepsilon,\tau}(t=2\pi)$ with different $\varepsilon$ and time step size $\tau$ for $S_{1}$.

The last two rows of Table 2.1 show the largest error of each column for fixed $\tau$. The errors exhibit $1/2$ order convergence, which coincides well with Theorems 1. More specifically, we can observe when $\tau\gtrsim\varepsilon$ (below the lower bolded diagonal line), there is first order convergence, which agrees with the error bound $\|\Phi(t_{n},x)-\Phi^{n}(x)\|_{H^{1}}\lesssim\tau+\varepsilon$. When $\tau\lesssim\varepsilon^{2}$ (above the upper bolded diagonal line), there is also first order convergence, which matches the other error bound $\|\Phi(t_{n},x)-\Phi^{n}(x)\|_{H^{1}}\lesssim\tau+\tau/\varepsilon$.

To support the improved uniform error bound in Theorem 4, we further test the discrete errors using non-resonant time steps, i.e., we choose $\tau\in\mathcal{A}_{\delta}(\varepsilon)$ for some given $\varepsilon$ and fixed $0<\delta\leq 1$. In this case, the bounded computational domain is set as $\Omega=(-16,16)$.

For comparison, the numerical ‘exact’ solution is computed by the second-order time-splitting method ($S_{2}$) with a very small time step size $\tau_{e}=8\times 10^{-6}$.

Figure 2.2 shows the errors $e^{\varepsilon,\tau}(t=4)$ with different $\varepsilon$ and time step size $\tau$ for $S_{1}$.

From the left part of Fig. 2.2, we could see that for each $\varepsilon\in(0,1]$, there is always first order convergence in $\tau$ for non-resonant time steps. From the right part, we find that for fixed time step size $\tau$, i.e., for each line in the figure, the error $e^{\varepsilon,\tau}(t=4)$ does not change much with different $\varepsilon$. This verifies the temporal uniform first order convergence for $S_{1}$ with non-resonant time step size, as stated in Theorem 4.

Through the results of this example, we successfully validate the uniform error bounds for $S_{1}$ in Theorems 1 & 4.

3 Extension to the second-order splitting method

In this section, we extend the results in the previous section to the second-order Strang splitting method.

Applying the discrete-in-time second-order splitting (Strang splitting, $S_{2}$) to (2.2), we have the numerical method as [6, 67]

[TABLE]

with $\Phi^{0}(x)=\Phi_{0}(x)$. We write the numerical propagator for $S_{2}$ as $\Phi^{n+1}(x):=S_{n,\tau}^{\text{Str}}(\Phi^{n})$.

3.1 Uniform error bounds

For the numerical approximation $\Phi^{n}(x)$ obtained from $S_{2}$ (3.1), we introduce the error function as in $S_{1}$

[TABLE]

and the following uniform error bounds hold.

Theorem 5**.**

Let $\Phi^{n}(x)$ be the numerical approximation obtained from $S_{2}$ (3.1), then under the assumptions $(A)$ and $(B)$ with $m=2$, there exists $0<\tau_{0}\leq 1$ independent of $\varepsilon$ such that the following error estimates hold for $0<\tau<\tau_{0}$,

[TABLE]

As a result, there is a uniform error bound for $S_{2}$ for $\tau>0$ small enough

[TABLE]

Proof.

As the proof of the theorem is not difficult to establish by combining the techniques used in proving Theorem 1 and the ideas in the proof of the uniform error bounds for $S_{2}$ in the linear case [7], we only give the outline of the proof here. For simplicity, we assume $V(x)\equiv 0$ and denote $\Phi(t):=\Phi(t,x)$, $\Phi^{n}:=\Phi^{n}(x)$ in short. Similar to the $S_{1}$ case, the $H^{1}$ bound of the numerical solution $\Phi^{n}$ is needed and can be done by using mathematical induction. For simplicity, we will assume the $H^{1}$ bound of $\Phi^{n}$ as in (2.35).

Step 1. Use Taylor expansion and Duhamel’s principle repeatedly to represent the ‘local truncation error’ $\eta^{n}(x)=\Phi(t_{n+1})-e^{-\frac{i\tau}{2\varepsilon^{2}}Q^{\varepsilon}}e^{-i\tau\mathbf{F}(e^{-\frac{i\tau}{2\varepsilon^{2}}Q^{\varepsilon}}\Phi^{n})}e^{-\frac{i\tau}{2\varepsilon^{2}}Q^{\varepsilon}}\Phi(t_{n})$ [6, 58] as

[TABLE]

where $\|R^{n}(x)\|_{H^{1}}\lesssim\tau^{3}+\tau\|\mathbf{e}^{n}(x)\|_{H^{1}}$, $f^{n}(s)$ is the same as that in Lie splitting $S_{1}$ case (2.24) and

[TABLE]

Step 2. For $h^{n}(s)$, using Duhamel’s principle to get

[TABLE]

where $\phi^{n}(s)=e^{-\frac{is}{\varepsilon^{2}}Q^{\varepsilon}}\Phi(t_{n})$, $\hat{f}^{n}(s)=ie^{-\frac{is}{\varepsilon^{2}}Q^{\varepsilon}}\int_{0}^{s}(f^{n}(w)-f^{n}(s))\,dw$, and we could find

[TABLE]

Recalling $\hat{f}^{n}(s)=O(s)$ and (2.30), we get $f^{n}(s)-f^{n}(w)=f^{n}_{2}(s)-f_{2}^{n}(w)+O(s)$ with $f_{2}^{n}(s)$ given in (2.17). Finally, under the assumption of Theorem 5, expanding $e^{-\frac{isQ^{\varepsilon}}{\varepsilon^{2}}}\Phi(t_{n})=e^{-\frac{is}{\varepsilon^{2}}}\Phi_{+}^{\varepsilon}(t_{n})+e^{\frac{is}{\varepsilon^{2}}}\Phi_{-}^{\varepsilon}(t_{n})+O(s)$, we can write the $h^{n}(s)$ term as

[TABLE]

with $\zeta_{1}^{n}(x)$ given as

[TABLE]

By taking $e^{-\frac{isQ^{\varepsilon}}{\varepsilon^{2}}}\approx e^{-\frac{is}{\varepsilon^{2}}}\Pi_{+}+e^{\frac{is}{\varepsilon^{2}}}\Pi_{-}$, it can be proved that $\|\zeta_{1}^{n}(x)\|_{H^{1}}\lesssim\min\{\tau^{2}\varepsilon,\frac{\tau^{3}}{\varepsilon}\}$.

Similarly, $g^{n}(s,w)$ can be written as

[TABLE]

where $\|\mathcal{G}_{3}^{n}(s,w)\|_{H^{1}}\lesssim\tau$, the oscillatory term (in time) $\mathcal{G}_{1}^{n}(s,w)$ simplifies $g^{n}(s,w)$ by using $e^{-\frac{isQ^{\varepsilon}}{\varepsilon^{2}}}\approx e^{-\frac{is}{\varepsilon^{2}}}\Pi_{+}+e^{\frac{is}{\varepsilon^{2}}}\Pi_{-}$ and removing the non-oscillatory terms as in (2.33), $\mathcal{G}_{2}^{n}(s,w)=\mathcal{G}_{2}^{n}(0,0)$ is the non-oscillatory term ($s,w$ independent) similar to (2.33), $\|\mathcal{G}_{1}^{n}(s,w)\|_{H^{1}}\lesssim\varepsilon$. We can prove $\|\partial_{s}\mathcal{G}_{1}^{n}(s,w)\|_{H^{1}}\lesssim 1/\varepsilon$, $\|\partial_{w}\mathcal{G}_{1}^{n}(s,w)\|_{H^{1}}\lesssim 1/\varepsilon$.

Lastly, $f^{n}(s)$ can be decomposed as

[TABLE]

where $\|\mathcal{F}_{3}^{n}(s)\|_{H^{1}}\lesssim\tau^{2}$, the oscillatory term (in time) $\mathcal{F}_{1}^{n}(s)$ simplifies $f^{n}(s)$ by using $e^{-\frac{isQ^{\varepsilon}}{\varepsilon^{2}}}=e^{-\frac{is}{\varepsilon^{2}}}(I_{2}-is\mathcal{D}^{\varepsilon})\Pi_{+}+e^{\frac{is}{\varepsilon^{2}}}(I_{2}+is\mathcal{D}^{\varepsilon})\Pi_{-}+O(s^{2})$ and removing the non-oscillatory terms as in (2.33), $\mathcal{F}_{2}^{n}(s)=\mathcal{F}_{2}^{n}(0)$ is the non-oscillatory term ($s$ independent) similar to (2.33). We can prove $\|\mathcal{F}_{1}^{n}(s)\|_{H^{1}}\lesssim\varepsilon$, $\|\partial_{s}\mathcal{F}_{1}^{n}(s)\|_{H^{1}}\lesssim 1/\varepsilon$, $\|\partial_{ss}\mathcal{F}_{1}^{n}(s)\|_{H^{1}}\lesssim 1/\varepsilon^{3}$.

Denote

[TABLE]

and we have

[TABLE]

where $\kappa^{n}(x)=R^{n}(x)+e^{-\frac{i\tau}{\varepsilon^{2}}Q^{\varepsilon}}\left(\kappa_{1}^{n}(x)+\int_{0}^{\tau}\mathcal{F}_{3}^{n}(s)\,ds-\tau\mathcal{F}_{3}^{n}\left(\frac{\tau}{2}\right)-\int_{0}^{\tau}\int_{0}^{s}\mathcal{G}_{3}^{n}(s,w)\,dwds+\frac{\tau^{2}}{2}\mathcal{G}_{3}^{n}\left(\frac{\tau}{2},\frac{\tau}{2}\right)\right)$ and $\|\kappa^{n}(x)\|_{H^{1}}\lesssim\tau^{3}+\tau\|\mathbf{e}^{n}(x)\|_{H^{1}}$.

Following the idea in $S_{1}$ case (2.38), we have the error equation for $S_{2}$

[TABLE]

where $\widetilde{L}_{n}{\mathbf{e}^{n}}(x)=e^{-\frac{i\tau}{2\varepsilon^{2}}Q^{\varepsilon}}\left(e^{-i\tau\mathbf{F}\left(e^{-\frac{i\tau}{2\varepsilon^{2}}Q^{\varepsilon}}\Phi^{n}\right)}-I_{2}\right)e^{-\frac{i\tau}{2\varepsilon^{2}}Q^{\varepsilon}}$, and $\|\widetilde{L}_{n}{\mathbf{e}^{n}(x)}\|_{H^{1}}\leq e^{c_{M_{1}}\tau}\|\mathbf{e}^{n}(x)\|_{H^{1}}$ ($c_{M_{1}}$ depends on $M_{1}$). For $0\leq n\leq\frac{T}{\tau}-1$, we would have (following (2.54)),

[TABLE]

Under the hypothesis of Theorem 5, we have

[TABLE]

which together with (3.11) gives $\|\zeta_{2}^{n}(x)\|_{H^{1}}\lesssim\min\{\varepsilon\tau,\tau^{3}/\varepsilon^{3}\}$ and $\|\zeta_{3}^{n}(x)\|_{H^{1}}\lesssim\min\{\varepsilon\tau^{2},\tau^{3}/\varepsilon\}$. Since $\|\zeta_{1}^{n}(x)\|_{H^{1}}\lesssim\min\{\tau^{2}\varepsilon,\frac{\tau^{3}}{\varepsilon}\}$, we derive from (3.14) that

[TABLE]

The discrete Gronwall’s inequality gives the desired results in Theorem 5 with the help of mathematical induction. □

For non-resonant time steps, i.e., for $\tau\in\mathcal{A}_{\delta}(\varepsilon)$, similar to $S_{1}$, we can derive improved uniform error bounds for $S_{2}$ as shown in the following theorem.

Theorem 6**.**

Let $\Phi^{n}(x)$ be the numerical approximation obtained from $S_{2}$ (3.1). If the time step size $\tau$ is non-resonant, i.e. there exists $0<\delta\leq 1$, such that $\tau\in\mathcal{A}_{\delta}(\varepsilon)$, then under the assumptions $(A)$ and $(B)$ with $m=2$, the following two error estimates hold for small enough $\tau>0$

[TABLE]

As a result, there is an improved uniform error bound for $S_{2}$ when $\tau>0$ is small enough

[TABLE]

Proof.

As the proof is extended from the techniques used for $S_{1}$ and the proof for improved uniform error bounds for $S_{2}$ in the linear case [7], here we just show the outline of the proof for brevity.

We start from (3.15). Following the strategy in the $S_{1}$ case, the key idea is to extract the leading terms from $\Phi(t,x)$ as (2.55) for estimating $\zeta_{2}^{n}(x)$, and the computations are more or less the same. Recalling (3.11) , noticing $\mathcal{F}_{1}^{n}(s)$ is similar to $f_{2}^{n}(s)$ (2.17) and $\|\zeta_{2}^{n}(x)\|_{H^{1}}\lesssim\min\{\varepsilon\tau,\tau^{2}/\varepsilon\}$, following the computations in the proof of Theorem 4, we would get for $0\leq n\leq\frac{T}{\tau}-1$ and $\tau\in A_{\delta}(\varepsilon)$,

[TABLE]

and the conclusions of Theorem 6 hold by applying the discrete Gronwall inequality to (3.15). □

3.2 Numerical results

In this subsection, we use a numerical example to validate our uniform error bounds in Theorems 5 and 6.

In the example, we choose the nonlinearity and the initial values as (2.70) and (2.71). In order to show that the error estimates still hold for $V\neq 0$, here we take the electric potential

[TABLE]

We first test the errors for resonant time steps, that is, for small enough chosen $\varepsilon$, there is a positive $k_{0}$, such that $\tau=\frac{1}{2}k_{0}\varepsilon^{2}\pi$, to check the error bounds in Theorem 5. In this case, the bounded computational domain is taken as $\Omega=(-32,32)$. The numerical ‘exact’ solution is generated by $S_{2}$ with a very fine time step size $\tau_{e}=2\pi\times 10^{-6}$.

The discrete $H^{1}$ error $e^{\varepsilon,\tau}(t_{n})$ used to show the results is defined in (2.72). It should be close to the $H^{1}$ errors in Theorems 5 here. In addition, we test the performance of $S_{2}$ in approximating the physical observables including probability density, current density, and energy. The discrete $l^{1}$ error for probability density is defined as

[TABLE]

the discrete relative $l^{1}$ error for current density is given by

[TABLE]

where $\mathbf{J}(\Phi^{n})=(J_{1}(\Phi^{n}),J_{2}(\Phi^{n}))^{T}$, with

[TABLE]

and the relative error for energy is defined as

[TABLE]

where

[TABLE]

Tables 3.2 to 3.5 exhibit the corresponding numerical temporal errors $e^{\varepsilon,\tau}(t=2\pi)$, $e_{\rho}^{\varepsilon,\tau}(t=2\pi)$, $e_{\mathbf{J}}^{\varepsilon,\tau}(t=2\pi)$, and $e_{E}^{\varepsilon,\tau}(t=2\pi)$ for $S_{2}$ with different $\varepsilon$ and resonant time step size $\tau$.

In these tables, the last two rows show the largest error of each column for fixed $\tau$. We could observe similar patterns for the errors of the wave function, and the physical observables. Clearly, overall there is $1/2$ order convergence, which agrees well with Theorem 5 for the wave function, and also suggests the same convergence rate for the observables. More specifically, from Tables 3.2 to 3.5, we can see when $\tau\gtrsim\sqrt{\varepsilon}$ (below the lower bolded diagonal line), there is second order convergence, which coincides with the error bound $\tau^{2}+\varepsilon$; when $\tau\lesssim\varepsilon^{2}$ (above the upper bolded diagonal line), we also observe second order convergence, which matches the other error bound $\tau^{2}+\tau^{2}/\varepsilon^{3}$.

Furthermore, to support the improved uniform error bound in Theorems 6, we test the error bounds using non-resonant time step sizes, i.e., we choose $\tau\in\mathcal{A}_{\delta}(\varepsilon)$ for some given $\varepsilon$ and fixed $0<\delta\leq 1$. Similar to the resonant time step case, we also test the errors for physical observables. The bounded computational domain is set as $\Omega=(-16,16)$.

For comparison, the numerical ‘exact’ solution is computed by $S_{2}$ with a very small time step size $\tau_{e}=8\times 10^{-6}$. Spatial mesh size is fixed as $h=1/16$ for all the numerical simulations.

Tables 3.6 to 3.9 show the numerical temporal errors $e^{\varepsilon,\tau}(t=4)$, $e_{\rho}^{\varepsilon,\tau}(t=4)$, $e_{\mathbf{J}}^{\varepsilon,\tau}(t=4)$, and $e_{E}^{\varepsilon,\tau}(t=4)$ with different $\varepsilon$ and non-resonant time step size $\tau$ for $S_{2}$.

The last two rows in Table 3.6 to 3.9 show the largest error of each column for fixed $\tau$, which gives $3/2$ order of uniform convergence, and it is consistent with Theorem 6 for the wave function. We could conclude that for physical observables, the convergence rate is the same. More specifically, in these tables, we can roughly observe the second order convergence when $\tau\gtrsim\varepsilon$ (below the lower bolded diagonal line) or when $\tau\lesssim\varepsilon^{2}$ (above the upper bolded diagonal line), agreeing with the error bound $\tau^{2}+\tau\varepsilon$ and the other error bound $\tau^{2}+\tau^{2}/\varepsilon$, respectively. When $\tau$ is large, the performance of the algorithm for probability density and current density is better than the performance for wave function and energy.

Through the results of this example, we successfully validate the uniform error bounds of $S_{2}$ in Theorems 5 and 6.

Remark 3.1**.**

Through extensive numerical results not shown here for brevity, we found out that the super-resolution property also holds true for higher order time-splitting methods in solving the NLDE. Specifically, the fourth-order compact splitting method for the Dirac equation [14] and the fourth-order partitioned Runge-Kutta splitting method for the NLDE [17, 6] exhibits 1/2 order uniform convergence under resonant time steps, and the uniform order could be improved to 3/2 under non-resonant time steps. The details are omitted here for brevity.

4 Conclusion

We studied the super-resolution property of time-splitting methods for the nonlinear Dirac equation in the nonrelativistic regime without magnetic potential in this paper. The uniform and improved uniform error bounds under non-resonant time step sizes for Lie-Trotter splitting ($S_{1}$) and Strang splitting ($S_{2}$) were rigorously established. For $S_{1}$, there are two independent error bounds $\tau+\varepsilon$ and $\tau+\tau/\varepsilon$, which give a uniform $1/2$ order convergence. Surprisingly, there is an improved uniform first order convergence if the time step sizes are non-resonant. For $S_{2}$, the two different error bounds are $\tau^{2}+\varepsilon$ and $\tau^{2}+\tau^{2}/\varepsilon^{3}$, also resulting in a uniform $1/2$ order convergence. For non-resonant time step sizes, the convergence rates can be improved to $3/2$ for $S_{2}$, with the two independent error bounds as $\tau^{2}+\tau\varepsilon$ and $\tau^{2}+\tau^{2}/\varepsilon$. Numerical results agreed with our theorems and suggested that our estimates are sharp. We remark that super-resolution also holds true for higher order splitting methods. Moreover, although only 1D cases are presented in this paper, these results are valid in higher dimensions, and the proofs can be easily generalized.