1 Introduction↩︎

The Dirac equation, formulated by Paul Dirac in 1928, describes the relativistic quantum dynamics of spin-1/2 particles. The equation not only accounts for the fine structure of the hydrogen spectrum, but also successfully predicted the existence of antimatter, thus having a profound impact on particle physics [1]. In condensed matter systems, the massless Dirac equation governs the low‑energy electronic behavior of graphene, explaining its unusual linear dispersion and the associated quantum Hall effect [2], [3]. Furthermore, the Dirac equation provides the fundamental framework for understanding the topological surface states and associated dissipationless edge currents in topological insulators and topological semimetals [4], [5]. In the field of quantum simulation, it has been employed to model physical phenomena in curved spacetime [6]. In strong‑field physics, the time‑dependent Dirac equation is essential for simulating highly relativistic electron dynamics under ultra‑intense laser fields (e.g., attosecond lasers), including highly non-linear quantum electrodynamical processes such as electron-positron pair production [7]. Among its various regimes, the Dirac equation in the semiclassical limit is of particular interest due to its distinctive physical and mathematical features, as it describes the transition from quantum to classical dynamics.

According to the framework established in [8], [9], the \(d\)-dimensional (\(d = 1,2,3\)) semiclassical Dirac equation admits the following formulation:

\[\label{def-Dirac-complete} i \epsilon \partial_t \boldsymbol{\boldsymbol{\Psi}} = \left[ \left( -i \epsilon \sum_{j=1}^d \alpha_j \partial_j + \beta \right) + \left( V(t, \boldsymbol{x}) I_4 - \sum_{j=1}^dA_j(t, \boldsymbol{x}) \alpha_j \right) \right] \boldsymbol{\Psi},\tag{1}\] where \(i = \sqrt{-1}\) is the imaginary unit, \(\boldsymbol{x} = (x_1, \cdots,x_d)\) represents the spatial coordinate, \(t\) is time, \(I_4\) denotes the \(4\times 4\) identity matrix and the 4 \(\times\) 4 matrices \(\alpha_j\) (\(j = 1\), \(\cdots\), \(d\)), \(\beta\) are defined as follows:

\[\alpha_j = \begin{bmatrix} \boldsymbol{0} & \sigma_j \\ \sigma_j & \boldsymbol{0} \end{bmatrix}, \qquad j = 1, \cdots, d,\qquad \beta = \begin{bmatrix} I_2 & \boldsymbol{0} \\ \boldsymbol{0} & - I_2 \end{bmatrix}.\] In the above expression, \(I_2\) denotes the 2 \(\times\) 2 identity matrix and \(\sigma_j\) (\(j = 1\), \(\cdots\), \(d\)) represent the Pauli matrices which are given as

\[\sigma_1 = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}, \qquad \sigma_2 = \begin{bmatrix} 0 & -i \\ i & 0 \end{bmatrix}, \qquad \sigma_3 = \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix}.\]

Moreover, in the equation 1 , \(\boldsymbol{\Psi}\in\mathbb{C}^4\) represents the four-component spinor wave function, and \(V, A_j \in \mathbb{R}\) (\(j = 1\), \(\cdots\), \(d\)) represent the electric and magnetic potentials, respectively. The parameter \(\epsilon \in (0,1)\) is the scaled Planck constant, which influences the scale of oscillations in the solution.

In the cases \(d=1,2\), as given in [10], we can reduce the Dirac equation 1 to the equation of the two-component wave function \(\Phi\) :

\[\label{def-Dirac-simple} i \epsilon \partial_t \boldsymbol{\Phi} = \left[ \left( -i \epsilon \sum_{j=1}^d \sigma_j \partial_j + \sigma_3 \right) + \left( V(t, \boldsymbol{x}) I_2 - \sum_{j=1}^dA_j(t, \boldsymbol{x}) \sigma_j \right) \right] \boldsymbol{\Phi},\tag{2}\]

where \(\Phi = \left( \Psi_1, \Psi_4 \right)^T\) or \(\left( \Psi_2, \Psi_3 \right)^T \in \mathbb{C}^2\). For simplicity, we only deal with the one- or two-dimensional (1D or 2D) semiclassical Dirac equation with the reduced form 2 in the following discussion.

In the classical regime where \(\epsilon = 1\), there have been a variety of efficient numerical methods for the Dirac equation, such as the finite difference time domain (FDTD) methods [11], [12] and exponential wave integrator Fourier pseudospectral (EWI-FP) method [8]. However, in the semiclassical regime, these methods may suffer from less efficiency due to the high frequency oscillation of the solution. Specificially, their accuracy typically depends on the discretization parameters resolving the small wavelength associated with \(\epsilon\). In the semiclassical regime, i.e. \(0 < \epsilon \ll 1\), the solution oscillates with the wavelength \(O( \epsilon)\) in time and space. This high-frequency oscillation requires fine mesh size \(h = O(\epsilon)\) and \(\tau = O(\epsilon)\) , which makes direct numerical resolution prohibitively expensive. This motivates the study of numerical methods with improved \(\epsilon\)-scalability. Several efficient numerical methods for the semiclassical Dirac equation including LFFD, SIFD1 and SIFD2 [13]. Most FDTD schemes are conditionally stable, requiring a coupling between the spatial and temporal mesh sizes. In contrast, the CNFD method is unconditionally stable but entailing a significant computational cost as it requires solving a linear system at each time step. To address the above difficulties, we apply the time-splitting methods to the semiclassical Dirac equation, which is used to solve the Schrödinger equation in [14]–[16], the Stokes equation [17] and the Dirac equation in different regimes [9], [18], [19]. Specifically, we consider the first-order (\(S_1\)), the second-order (\(S_2\)) and the compact fourth-order (\(S_{4c}\)) splitting methods. Among these methods, \(S_1\) and \(S_2\) are mass-preserving, while \(S_\text{4c}\) also conserves mass in 1D. Compared to FDTD schemes, the time-splitting methods are highly efficient and exhibit higher-order accuracy.

The rest of this paper is structured as follows: In Section 2, we introduce the time-splitting algorithms, including \(S_1\), \(S_2\) and \(S_\text{4c}\), for solving the semiclassical Dirac equation and show the mass conservation property. In Section 3, we give the rigorous error estimates of the time-splitting methods under some regularity conditions and prove the result for \(S_1\). We conduct numerical experiments in Section 4 to verify the error estimates established. Finally, we conclude in Section 5.

2 Time-splitting Spectral Methods↩︎

In this section, we apply the time-splitting methods to the semiclassical Dirac equation and prove the property of mass conservation. To employ spectral methods with Fourier basis which enables the efficient use of the Fast Fourier Transform (FFT), we impose periodic boundary conditions on the Dirac equation 3 . This assumption is reasonable since the wave function usually vanishes on the boundary of a sufficiently large domain. For the sake of notational simplicity and clarity of presentation, we restrict our discussion to the one-dimensional case of 2 , for which the numerical schemes and the corresponding analysis are presented in detail. The extension to 1 and higher-dimensional settings is standard and straightforward, and all the results established in this work remain valid without essential modifications.

Moreover, in order to avoid unnecessary technicalities, we confine ourselves to the case where the potential functions are independent of time. The more general situation involving explicitly time-dependent potentials can be handled by means of the time-ordering technique; we refer the reader to [19] for a detailed treatment.

2.1 Operator splitting for the Dirac equation↩︎

For the Dirac equation in 1D \(\left( d = 1 \right)\), we consider the following initial-boundary value problem:

\[\begin{align} i \epsilon \partial_t \Phi(t,\mathbf{x}) &= \left[ \left( -i \epsilon \sigma_1 \partial_x + \sigma_3 \right) + \left( V(x) I_2 - A_1(x) \sigma_1 \right) \right] \Phi(t,x) \notag \\ &\qquad (t,x) \in \mathbb{R}^+ \times \Omega, \label{def95Dirac95simplified}\\ \Phi(t, a) &= \Phi(t, b), \quad \partial_{x}\Phi(t, a) = \partial_{x}\Phi(t, b), \quad t > 0, \\ \Phi(0,x) &= \Phi_0(x), \quad x \in \bar{\Omega}. \end{align}\tag{3}\]

where \(\Phi, \Phi_0 \in \mathbb{C}^2, x \in \Omega := (a,b), t \in \left[0,T_\text{end}\right]\). To simplify the notation, we define:

\[\begin{align} & T = -\sigma_1 \partial_x - \frac{i}{\epsilon} \sigma_3, \tag{4}\\ & W = -\frac{i}{\epsilon} V(x) I_2 + \frac{i}{\epsilon} A_1 (x) \sigma_1. \tag{5} \end{align}\]

Then the Dirac equation can be written in the compact form \[\label{def95dirac95operato} \partial_t \Phi = (T + W)\Phi.\tag{6}\] Formally, its solution is given by \[\Phi(t,x)=e^{t(T+W)}\Phi(0,x).\] Applying the operator splitting technique, the evolution operator can be approximated by \[\Phi(t,x)\approx \prod_{j=1}^{n} e^{a_j t T} e^{b_j t W}\Phi(0,x),\] where \(n\geq 1\), and the coefficients \[a_j,b_j\in\mathbb{R}, \qquad j=1,2,\ldots,n,\] are parameters to be determined. The coefficients are chosen such that the resulting splitting scheme attains a prescribed order of accuracy with respect to the time step.

2.2 The full-discretization of the splitting methods↩︎

Let \(h = \frac{b - a}{M}\) denote the mesh size with \(M\) being an even integer and \(\tau>0\) denote the time step. We denote the grid points and the time steps as:

\[\begin{align} x_j = a + j h,\quad j = 0,1,\dots,M-1; \quad t_n = n \tau,\quad n = 0,1,2,\dots \end{align}\]

Let \(\Phi_{j}^{ n} \approx \Phi(t_n, x_j)\) be the numerical solution at the grid point \(x_j\) and the time step \(t_n\), \(V_j := V( x_j)\) and \(A_{1,j} := A_1(x_j)\) denote the electric and magnetic potentials, respectively.

The Fourier interpolation function of \(f\) based on the given grid points can be formulated as follows:

\[\begin{align} & f_I(x) = \sum_{l = -\frac{M}{2}}^{\frac{M}{2}-1}\hat{f}_l e^{i \mu_{l}(x - a)}, \quad\hat{f}_l = \frac{1}{M} \sum_{j = 0}^{M - 1} f(x_j) e^{-i \mu_l(x_j - a)}, \tag{7}\\ & \mu_l = \frac{2 \pi l}{b - a},\quad l = - \frac{M}{2}, \cdots, \frac{M}{2} - 1. \tag{8} \end{align}\]

In this paper, we define the \(l^2\)-norm for a function as:

\[\label{l942-norm} \left\| \Phi(t, \cdot) \right\|_{L^2} := \sqrt{\int_{\Omega} \left| \Phi(t,x)^2 \right| dx}.\tag{9}\]

2.2.1 The first-order time-splitting method (\(S_1\))↩︎

In the first-order splitting, from time \(t = t_n\) to time \(t = t_{n+1}\), we solve the Dirac equation in two sub-steps. The approximate solution of the Dirac equation 6 can be computed by the following iterative formula:

\[\begin{align} & \Phi(t_{n+1}, \cdot) = e^{\tau (W +T) } \Phi(t_{n}, \cdot) \approx e^{\tau W } e^{\tau T} \Phi(t_{n}, \cdot), \quad n \geq 0. \end{align}\]

Specifically, in the first sub-step, we compute \(e^{\tau T}\Phi(t_n, \cdot)\) in the phase space using the Fourier spectral method. After that, in the second sub-step, we integrate the resulting system exactly in time for the \(W\) operator. More precisely, we denote by \(\Phi^n=\left(\Phi_0^n,\Phi_1^n,\cdots,\Phi_{M-1}^n\right)\) the numerical solution at the temporal grid point \(t=t_n\). The proposed fully discrete scheme is formulated as follows: for each \(n\geq 0\),

\[\begin{align} &\Phi^{n, *}_j = \sum_{l = -\frac{M}{2}}^{\frac{M}{2}-1} e^{\tau \left( -i \sigma_1 \mu_l - \frac{i}{\epsilon} \sigma_3 \right)} \hat{\Phi}^n_l e^{i \mu_l (x_{j} - a)}, \tag{10}\\ &\Phi^{n+1}_j = e^{\tau \left[ -\frac{i}{\epsilon} V( x_j) I_2 + \frac{i}{\epsilon} A_1 ( x_j) \sigma_1 \right] } \Phi^{n, *}_j, \tag{11} \end{align}\]

where \(\mu_l\) is defined as the equality 8 , the initial condition is given by \(\Phi^0=\left( \Phi^0_0,\Phi^0_1,\cdots,\Phi^0_{M-1} \right)=\left( \Phi^0(x_0),\Phi^0(x_1),\cdots,\Phi^0(x_{M-1}) \right)\) and \(\hat{\Phi}^n_j\) gives the Fourier coefficients of \(\Phi_j^n\)

\[\label{sp1953} \hat{\Phi}^n_j = \frac{1}{M} \sum_{j = 0}^{M-1} \Phi^n_j e^{-i \mu_l(x_j - a)}.\tag{12}\] Obviously, the main computational cost of this algorithm comes from FFT, which scales as \(O(M\log M)\).

2.2.2 The second-order time-splitting method (\(S_2\))↩︎

To improve the accuracy of the approximate solution, the \(S_2\) algorithm is designed to narrow the gap between the exact solution operator and the splitting operator [20]. Specifically, from time \(t = t_n\) to time \(t = t_{n + 1}\), we split the solution operator as follows:

\[\begin{align} & \Phi(t_{n+1}, \cdot) = e^{\tau (W +T) } \Phi(t_{n}, \cdot) \approx e^{\frac{\tau}{2} W } e^{\tau T} e^{\frac{\tau}{2} W } \Phi(t_{n}, \cdot). \end{align}\] Similar to \(S_1\), the approximate solution can be obtained by the following steps for \(n\geq 0\):

\[\begin{align} & \Phi^{n, *}_j = e^{\frac{1}{2} \tau \left[ -\frac{i}{\epsilon} V( x) I_2 + \frac{i}{\epsilon} A_1 (x) \sigma_1 \right] } \Phi^n_j, \tag{13} \\ & \Phi ^{n,**}_j = \sum_{l = -\frac{M}{2}}^{\frac{M}{2}-1} e^{\tau \left( -i \sigma_1 \mu_l - \frac{i}{\epsilon} \sigma_3 \right)} \hat{\Phi}^{n, *}_l e^{i \mu_l(x_j - a)}, \tag{14} \\ & \Phi_j ^{n + 1} = e^{\frac{1}{2} \tau \left[ -\frac{i}{\epsilon} V( x) I_2 + \frac{i}{\epsilon} A_1 (x) \sigma_1 \right] } \Phi^{n, **}_j, \tag{15} \end{align}\] where the Fourier coefficients of \(\Phi^{n, *}_j\) are defined as

\[\label{sp2954} \hat{\Phi}^{n, *}_l = \frac{1}{M} \sum_{j = 0}^{M-1} \Phi^{n, *}_j e^{-i \mu_l (x_j - a)}.\tag{16}\]

We remark here that in practical computation, \(S_2\) can be more effectively implementated through the following observation:

\[\label{consolidation95S2} \begin{align} e^{T_\text{end}(T + W)} &\approx e^{\frac{\tau}{2} W} e^{\tau T} e^{\frac{\tau}{2} W} \cdots e^{\frac{\tau}{2} W} e^{\tau T} e^{\frac{\tau}{2} W} \\ & = e^{\frac{\tau}{2} W} e^{\tau T} e^{\tau W} e^{\tau T} e^{\tau W} \cdots e^{\tau T} e^{\tau W} e^{\frac{\tau}{2} W}. \end{align}\tag{17}\]

It is apparent that under this framework, the computational cost for \(S_2\) is nearly identical to that of \(S_1\), yet it cab yield a substantial reduction in the approximation error.

2.2.3 The compact fourth-order time-splitting method (\(S_\text{4c}\))↩︎

We can further reduce the error by enhancing the accuracy of the splitting scheme. In addition, to avoid the negative time steps [9], [21], we use \(S_\text{4c}\) to obtain the approximate solution [9] by using the following formula for \(n\geq 0\):

\[\label{def95s4c} \Phi^{n + 1} = e^{\frac{1}{6} \tau W } e^{\frac{\tau}{2} T} e^{\frac{2}{3} \tau W + \frac{\tau^3}{72} \left[ W, \left[ T, W \right] \right]} e^{\frac{\tau}{2} T} e^{\frac{1}{6} \tau W} \Phi^n.\tag{18}\]

In this scheme, we need to compute an additional operator \(e^{\frac{2}{3} \tau W + \frac{\tau^3}{72} \left[ W, \left[T, W \right]\right]}.\) Similar to the computation in [9], we can find out that this operator contains partial derivatives \(\partial_{x_1}, \cdots, \partial_{x_d}\) in 2D and 3D. As a result, in higher dimensions, to evaluate this operator and maintain fourth-order accuracy, we use the fourth-order Runge–Kutta method (RK4) in time and spectral discretization in space. We remark here that a simplification similar to 17 can be applied to \(S_\text{4c}\). For the sake of brevity, the detailed algorithm for \(S_\text{4c}\) is omitted here.

2.3 Mass Conservation↩︎

Let \(X_M= \Bigl\{ U=(U_0,U_1,\cdots,U_M) \;\big|\; U_j\in\mathbb{C}^2,\;j=0,1,\cdots,M,\;U_0=U_M \Bigr\}\), and adopt the periodic extension \(U_{-1}=U_{M-1}, U_{M+1}=U_1\), whenever such indices are involved. On the space \(X_M\), the discrete \(l_2\)-norm is defined by \[\|U\|_{l_2} = \sqrt{h\sum_{j=0}^{M-1}|U_j|^2}, \qquad U\in X_M.\]

We establish the discrete mass conservation property for the time-splitting scheme \(S_1\).

Similar conservation properties hold for \(S_2\), as well as \(S_\text{4c}\) in 1D. For brevity, we only give the mass conservation result and omit the proof here.

3 Error Estimates↩︎

For simplicity, in this paper we only show the error estimates of 2 for the case \(d = 1\). The generalization to the Dirac equation 1 and higher dimensions will follow directly with few modifications.

To obtain a proper error estimate, our analysis assumes that the exact solution and the electromagnetic potentials satisfy certain regularity conditions:

\[\begin{align} (A) \qquad &\Phi \in C^2 \left( \left[0,T_\text{end}\right]; W^{m,\infty}_p \left(\Omega \right) \right), \quad V, A_1 \in C^m\left(\Omega \right), \quad m \geq 4, \end{align}\] where \(W^{m,\infty}_p \left(\Omega \right) = \left\{ f | f \in W^{m,\infty} \left( \Omega \right), \partial_x^k f(a) = \partial_x^k f(b), k = 0,\cdots,m - 1 \right\}\), and \(m\) can be any integer with \(m \geq 4\), which will be related to the numerical accuracy in space. In practice, the electromagnetic potentials \(V\) and \(A_1\) may not be inherently periodic, but they nearly attain constant values at the boundaries of \(\Omega\). Under this circumstance, they can be reasonably treated as periodic over \(\Omega\). Consequently, the wave function \(\Phi\) will also be periodic.

In addition, the derivatives of the wave function and the electromagnetic potentials should satisfy the following conditions:

\[\begin{align} (B) \qquad & \left \| \frac{\partial^{r + k}}{\partial t^{r} \partial x^k} \Phi \right \|_{L^\infty([0,T_{end}]; (L^\infty(\Omega))^2)} \leq \frac{C_{r + k}}{\epsilon^{r + k}}\\ &\left \| \frac{\partial^{ k}}{\partial x^k} V \right \|_{L^\infty((L^\infty(\Omega))^2)} \leq D_{k}, \quad\left \| \frac{\partial^{k}}{\partial x^{k}} A_1 \right \|_{L^\infty((L^\infty(\Omega))^2)} \leq D_{k}, \\ & \qquad \qquad \qquad0 \leq r \leq 2,\, 0 \leq k \leq m. \end{align}\] where the constants \(C_{r + k}\) and \(D_{k}\) are independent of \(\epsilon\).

3.1 The Main Results↩︎

Based on the assumptions above, we can get the main results of the error estimates. To simplify the notation, we denote the numerical interpolation function of \(\Phi^n\) by \(\Phi_I^n\) at time \(t = t_n\). Specifically, we have: \[\label{def95numerical95interpolation} \Phi^n_I(x) = \sum_{l = -\frac{M}{2}}^{\frac{M}{2}-1} \hat{\Phi}^n_l e^{i \mu_l(x - a)}, \quad \hat{\Phi}^n_l = \frac{1}{M} \sum_{j = 0}^{M-1} \Phi^n_j e^{-i \mu_l (x - a)},\tag{19}\] with \(\mu_l = \frac{2 \pi l}{b - a}\), \(l=-\frac{M}{2}, -\frac{M}{2}+1, \cdots, \frac{M}{2}-1\). Then we can get the following theorems for the error estimates.

Furthermore, similar error estimates for the observables, such as the total probability density and the current density can be established by using the control relation between error estimates for the observables and the error estimate for the wave function \(\Phi\) [14], where the total probability density \(\rho\) and the current density \(J\) are defined as follows:

\[\begin{align} & \rho(t,x) = \Phi(t,x)^* \Phi(t,x), \tag{20}\\ & J(t,x) = \Phi(t,x)^* \sigma_1 \Phi(t,x). \tag{21} \end{align}\] Similar to \(\Phi_I^n(x)\), we can define the numerical interpolation functions \(\rho_I^n(\cdot)\) and \(J_I^n(\cdot)\): \[\begin{align} &\rho^n_I(x) = \sum_{l = -\frac{M}{2}}^{\frac{M}{2}-1} \hat{\rho}^n_l e^{i \mu_l(x - a)}, \quad \hat{\rho}^n_l = \frac{1}{M} \sum_{j = 0}^{M-1} \rho^n_j e^{-i \mu_l (x - a)},\\ & J^n_I(x) = \sum_{l = -\frac{M}{2}}^{\frac{M}{2}-1} \hat{J}^n_l e^{i \mu_l(x - a)}, \quad \hat{J}^n_l = \frac{1}{M} \sum_{j = 0}^{M-1} J^n_j e^{-i \mu_l (x - a)}. \end{align}\]

From the above theorems and corollaries, we can obtain that if \(\delta\) is the desired error bound, the associated discretization parameters shall satisfy the following conditions:

\[\begin{align} & S_1: \frac{\tau}{\epsilon^2} \lesssim \delta , \qquad \frac{h}{\epsilon} \lesssim \left(\delta \tau\right) ^{\frac{1}{m}}, \tag{22} \\ & S_2: \frac{\tau^2}{\epsilon^3} \lesssim \delta , \qquad \frac{h}{\epsilon} \lesssim \left(\delta \tau\right) ^{\frac{1}{m}}, \tag{23} \\ & S_\text{4c}: \frac{\tau^4}{\epsilon^5} \lesssim \delta , \qquad \frac{h}{\epsilon} \lesssim \left(\delta \tau\right) ^{\frac{1}{m}}. \tag{24} \end{align}\]

3.2 Proof of the Main Results↩︎

Before proving the Theorems 1 - 3, we first prove the following lemmas.

Next, we introduce an error estimate for the Fourier interpolation.

We refer the proof of the above Lemma 5 to [23] and the details are omitted here for brevity.
We are now ready to prove Theorem 1. The proof of Theorems 2 and 3 are similar, and thus omitted here.
Proof for Theorem 1

For convenience, we introduce the auxiliary functions \[\Gamma_n(\cdot)=e^{\tau W}e^{\tau T}\Phi(t_n,\cdot), \qquad \gamma_n(\cdot)=e^{\tau T}\Phi(t_n,\cdot).\] Here, \(\Gamma_n\) denotes one exact splitting step starting from the exact solution at time \(t_n\). Consequently, the difference between \(\Phi(t_{n+1},\cdot)\) and \(\Gamma_n(\cdot)\) measures the local error generated by the time-splitting procedure. By Lemma 4, this local truncation error satisfies \[\label{fisrt95term95thm95proof} \|\Phi(t_{n+1},\cdot)-\Gamma_n(\cdot)\|_{L^2} \le C_1\frac{\tau^2}{\epsilon^2}.\tag{25}\]

To derive the full discretization error estimate, we decompose the error into three different parts corresponding respectively to the temporal splitting error, the Fourier interpolation error, and the propagated numerical error. More precisely, by the triangle inequality, \[\begin{align} \label{error95decomposition} \|\Phi(t_{n+1},\cdot)-\Phi_I^{n+1}(\cdot)\|_{L^2} \le & \|\Phi(t_{n+1},\cdot)-\Gamma_n(\cdot)\|_{L^2} \nonumber\\ & +\|\Gamma_n(\cdot)-\Gamma_{n,I}(\cdot)\|_{L^2} +\|\Gamma_{n,I}(\cdot)-\Phi_I^{n+1}(\cdot)\|_{L^2}. \end{align}\tag{26}\]

We first estimate the interpolation error term. Applying Lemma 5, we have \[\|\Gamma_n(\cdot)-\Gamma_{n,I}(\cdot)\|_{L^2} \lesssim M^{-m}\|\partial_x^m\Gamma_n(\cdot)\|_{L^2}.\] Therefore, it remains to control the \(m\)-th order spatial derivative of \(\Gamma_n\). Recalling the definition of \(\Gamma_n\) and using Leibniz’s rule, we obtain \[\begin{align} \left\|\partial_x^m\Gamma_n(\cdot) \right\|_{L^2} &= \left\| \partial_x^m \left( e^{\tau W}e^{\tau T}\Phi(t_n,\cdot) \right) \right\|_{L^2} \\ &\le \sum_{q=0}^m \binom{m}{q} \left\|\partial_x^{m-q}e^{\tau W}\right\|_{L^\infty} \, \left\|\partial_x^q(e^{\tau T}\Phi(t_n,\cdot))\right\|_{L^2}. \end{align}\] where \(\left\|\partial_x^{m-q}e^{\tau W}\right\|_{L^\infty}\) represents the standard matrix \(L^{\infty}\)-norm of \(\partial_x^{m-q}e^{\tau W}\). Under the regularity assumptions imposed on \(W\) and the exact solution \(\Phi\), each term on the right-hand side is bounded, which leads to \[\|\partial_x^m\Gamma_n(\cdot)\|_{L^2} \lesssim \epsilon^{-m}.\] Substituting the above estimate into the Fourier interpolation inequality yields \[\label{second95term95thm95proof} \|\Gamma_n(\cdot)-\Gamma_{n,I}(\cdot)\|_{L^2} \lesssim \left( \frac{h}{\epsilon(b-a)} \right)^m.\tag{27}\]

Next, we turn to the last term in 26 , which corresponds to the propagation of the numerical error from the previous time level. To simplify the notation, define \[\Gamma^n= (\Gamma_n(x_0),\ldots,\Gamma_n(x_{M-1})), \qquad \gamma^n= (\gamma_n(x_0),\ldots,\gamma_n(x_{M-1})),\] and \[\Phi_{\mathrm{exact}}^n= (\Phi(t_n,x_0),\ldots,\Phi(t_n,x_{M-1})).\] Using the definition of the numerical scheme together with Parseval’s identity, we derive \[\begin{align} \label{third95term95thm95proof} \|\Gamma_{n,I}(\cdot)-\Phi_I^{n+1}(\cdot)\|_{L^2} &= \|\Gamma^n-\Phi^{n+1}\|_{l^2} \nonumber\\ &= \|\gamma^n-\Phi^{n,*}\|_{l^2} = \|\Phi_{\mathrm{exact}}^n-\Phi^n\|_{l^2} \nonumber\\ &= \|\Phi_I(t_n,\cdot)-\Phi_I^n(\cdot)\|_{L^2}. \end{align}\tag{28}\]

The above identity shows that the error at time \(t = t_{n+1}\) can be controlled by the numerical error at the previous step. Applying the triangle inequality once again gives \[\|\Phi_I(t_n,\cdot)-\Phi_I^n(\cdot)\|_{L^2} \le \|\Phi_I(t_n,\cdot)-\Phi(t_n,\cdot)\|_{L^2} + \|\Phi(t_n,\cdot)-\Phi_I^n(\cdot)\|_{L^2}.\] Moreover, Lemma 5 implies that \[\|\Phi_I(t_n,\cdot)-\Phi(t_n,\cdot)\|_{L^2} \lesssim \left( \frac{h}{\epsilon(b-a)} \right)^m.\]

Combining the estimates 25 , 27 , and 28 , we arrive at the recursive inequality \[\|\Phi(t_{n+1},\cdot)-\Phi_I^{n+1}(\cdot)\|_{L^2} \le C_1\frac{\tau^2}{\epsilon^2} + E_{m,1} \left( \frac{h}{\epsilon(b-a)} \right)^m + \|\Phi(t_n,\cdot)-\Phi_I^n(\cdot)\|_{L^2}.\] Finally, applying induction with respect to \(n\) immediately yields the desired estimate ?? . \(\square\)

Since the control coefficients of the splitting error and the Fourier interpolation error are independent of \(\epsilon\), the control coefficients of ?? ,?? ,?? are also independent of \(\epsilon\).

4 Numerical Results↩︎

In this section, we conduct the numerical experiments of the time-splitting spectral methods for the semiclassical Dirac equation to verify the error estimates given in Section 3. To evaluate the numerical errors, we use the following relative errors of the wave function \(\Phi\), the total probability density function \(\rho\) and the current density function \(\boldsymbol{J}\):

\[\begin{align} & e_{\Phi}^r (t_n) = \frac{\left\| \Phi(t_n, \cdot) - \Phi_I^n(\cdot) \right\|_{l^2}}{\left\| \Phi(t_n, \cdot) \right\|_{l^2}}, \tag{29} \\ & e_{\rho}^r (t_n) = \frac{\left\| \rho(t_n, \cdot) - \rho_I^n(\cdot) \right\|_{l^2}}{\left\| \rho(t_n, \cdot) \right\|_{l^2}}, \tag{30} \\ & e_{\boldsymbol{J}}^r (t_n) = \frac{\left\| \boldsymbol{J} (t_n, \cdot) - \boldsymbol{J}_I^n(\cdot) \right\|_{l^2}}{\left\| \boldsymbol{J}(t_n, \cdot) \right\|_{l^2}}, \tag{31} \end{align}\] where the total probability density function \(\rho\) and the current density function \(\boldsymbol{J}\) are defined in 20 and 21 .

4.1 An example in 1D↩︎

First, we consider a one-dimensional case with \(d = 1\). The electromagnetic potentials are given by \(V(x) = \frac{1 - x}{1 + x^2}, A_1(x) = \frac{(x + 1)^2}{1 + x^2}\) and the initial condition is chosen as follows:

\[\begin{align} & \phi_1(0,x) = \frac{1}{2} e^{-4x^2 + i \frac{S_0(x)}{\epsilon}} \left( 1 + \sqrt{1 + S_0'(x)^2} \right), \\ & \phi_2(0,x) = \frac{1}{2} e^{-4x^2 + i \frac{S_0(x)}{\epsilon}} S_0'(x), \end{align}\] where \(S_0(x) = \frac{1}{40} \left(1 + \cos(2 \pi x) \right)\) and \(\epsilon \in (0,1]\).

In our numerical experiments, we take the spatial domain \(\Omega = (-16,16)\) and we compute the solutions until \(t = 2\). Moreover, as there is no explicit expression for the exact solution, we use the numerical solution obtained via the \(S_\text{4c}\) method with time step size \(\tau = \frac{0.1}{4^{8}}\) and spatial grid size \(h = \frac{1}{2^{10}}\) as the “exact” solution for comparison.

Firstly, we test the dependence of the error on time step sizes and present the relative errors of \(S_1\) and \(S_2\) in Tables ¿tbl:tab:error95fun95tau9536S95136? and ¿tbl:tab:error95fun95tau9536S95236?. Tables ¿tbl:tab:error95fun95tau9536S95136? and ¿tbl:tab:error95fun95tau9536S95236? clearly show the first order convergence in time for \(S_1\) and the second order convergence in time for \(S_2\), when the spatial error is negligible. Moreover, the elements of the diagonal lines of the two tables exhibit that the \(\epsilon\)-scalability of \(S_1\) and \(S_2\) are \(\tau = O(\epsilon^2)\) and \(\tau = O(\epsilon^{\frac{3}{2}})\), respectively.

In addition, the relative error of the wave function obtained via the \(S_\text{4c}\) method in Table ¿tbl:tab:error95fun95tau95S4c? demonstrates the fourth-order accuracy in time. Compared to \(S_1\) and \(S_2\), \(S_\text{4c}\) achieves higher accuracy with less stringent time step size requirements owing to the higher \(\epsilon\)-scalability.

Furthermore, we compute the relative temporal errors of the total probability density \(e^r_{\rho}\) and the current density \(e^r_{\boldsymbol{J}}\) for \(S_1\) and \(S_2\). For simplicity, we omit the results of \(S_2\) here. The data of Table ¿tbl:tab:error95position95tau9536S95136? and ¿tbl:tab:error95curr95tau9536S95136? confirms the error bounds in Corollary [corollary-sp1].

And to make the dependence of various errors on time step sizes more explicit, we present the following convergence analysis plot 1 (a) and 1 (b) with the dimensionless parameter \(\epsilon = 1\), \(\frac{1}{16}\) respectively.

Figures 1 (a) and 1 (b) display the relationships between the relative errors \(e_{\Phi}^r,\;e_{\rho}^r,\;\text{and}\;e_{J}^r\) and the time step size \(\tau\) for \(\epsilon = 1\) and \(\epsilon = \frac{1}{16}\), respectively. It can be observed from these figures that the relative errors \(e_{\Phi}^r,\;e_{\rho}^r,\;\text{and}\;e_{J}^r\) exhibit the same convergence order for different values of \(\epsilon\), indicating a consistent numerical behavior of the scheme with respect to different physical quantities.

Furthermore, the \(S_2\) method achieves a higher convergence rate than the \(S_1\) method. To attain an accuracy of order \(10^{-6}\), smaller values of \(\epsilon\) require finer time step sizes \(\tau\). In addition, the \(S_1\) method demands significantly smaller time step sizes than the \(S_2\) method to achieve the same level of accuracy, demonstrating the superior efficiency and accuracy balance of the \(S_2\) scheme.

Finally, Table ¿tbl:tab:error95fun95h9536S95136? reports the spatial convergence behavior of the \(S_1\) method with the time step size fixed at \(\tau = \frac{0.1}{4^{8}}\). It can be observed that, as the spatial mesh size \(h\) decreases, the relative error exhibits a typical spectral convergence behavior. This result confirms the high spatial accuracy of the proposed scheme.

4.2 An example in 2D↩︎

We now consider an example in 2D with the electromagnetic potentials: \[\begin{align} & V(\boldsymbol{x}) = \cos \left( \frac{4 \pi}{\sqrt{3}} x_1 \right) + \cos \left[ \frac{4 \pi}{\sqrt{3}} \left( \frac{1}{2} x_1 + \frac{\sqrt{3}}{2} x_2 \right) \right] + \cos \left[ \frac{4 \pi}{\sqrt{3}} \left( \frac{1}{2} x_1 - \frac{\sqrt{3}}{2} x_2 \right) \right], \\ & A_1(\boldsymbol{x}) = 0 \quad, \qquad A_2(\boldsymbol{x}) = x_1, \end{align}\] where \(\boldsymbol{x} = (x_1, x_2) \in \Omega = (-10,10) \times (-10,10)\). And the initial condition is given as follows: \[\phi_1(0,\boldsymbol{x}) = e^{-\frac{x_1^2 + x_2^2}{2}}, \quad \phi_2(0,\boldsymbol{x}) = e^{-\frac{\left(x_1 - 1 \right)^2 + x_2^2}{2}}.\]

Similar to the example in 1D, we use the numerical solution obtained via the \(S_\text{4c}\) method with \(\tau = \frac{0.1}{4^{5}}\) and \(h = \frac{1}{32}\) as the “exact” solution. And we employ the relative errors introduced in 29 to evaluate the results.

To exhibit the temporal error, we always take the mesh size \(h = \frac{1}{32}\) in the numerical experiments such that the spatial error could be neglected. The Table ¿tbl:tab:error95fun95tau9536S95136952D? and ¿tbl:tab:error95fun95tau9536S95236952D? present the temporal errors of \(S_1\) and \(S_2\). We can clearly identify that \(S_1\) is first-order accurate in time, while \(S_2\) achieves second-order accuracy. Moreover, the numerical results in the Table ¿tbl:tab:error95fun95tau95S4c952D? demonstrate that the fourth-order accuracy is preserved in 2D, despite the introduction of additional errors arising from the implementation of the RK4 method.

5 Conclusion↩︎

This paper investigates time-splitting methods for the semiclassical Dirac equation. Through rigorous theoretical analysis, we prove that both the first-order and second-order time-splitting schemes preserve the mass conservation property of the Dirac equation. Furthermore, rigorous error estimates and convergence analysis are established for the proposed schemes. Theoretical results show that the proposed time-splitting methods achieve spectral accuracy in space, with temporal error bounds of \(O\left(\tau/\varepsilon^2\right), O\left(\tau^2/\varepsilon^3\right), O\left(\tau^4/\varepsilon^5\right)\) for \(S_1\), \(S_2\), and \(S_{\mathrm{4c}}\), respectively. Numerical experiments are presented to validate the theoretical results and demonstrate the effectiveness of the proposed methods.

Data Availability Statement↩︎

The implementation code and numerical datasets supporting this study have been deposited in a Zenodo repository. Access to the repository has been provided in the submission system.

References↩︎