Single-cone real-space finite difference scheme for the time-dependent Dirac equation

René Hammer rene.hammer@uni-graz.at Walter Pötz walter.poetz@uni-graz.at and Anton Arnold anton.arnold@tuwien.ac.at Institut für Physik, Karl-Franzens-Universität Graz, Universitätsplatz 5, 8010 Graz, Austria Institut für Analysis und Scientific Computing, TU-Wien, Wiedner Hauptstr. 8, 1040 Wien, Austria

Abstract

A finite difference scheme for the numerical treatment of the (3+1)D Dirac equation is presented. Its staggered-grid intertwined discretization treats space and time coordinates on equal footing, thereby avoiding the notorious fermion doubling problem. This explicit scheme operates entirely in real space and leads to optimal linear scaling behavior for the computational effort per space-time grid-point. It allows for an easy and efficient parallelization. A functional for a norm on the grid is identified. It can be interpreted as probability density and is proved to be conserved by the scheme. The single-cone dispersion relation is shown and exact stability conditions are derived. Finally, a single-cone scheme for the two-component (2+1)D Dirac equation, its properties, and a simulation of scattering at a Klein step are presented.

keywords:

Dirac equation , leap-frog , staggered grid , fermion doubling , FDTD, Klein step

1 Introduction

1.1 The Dirac Equation and Numerical Schemes

Even more than $80$ years after its presentation by P. A. M. Dirac in 1928, the Dirac equation has not lost its fascination and significance in physics [1]. The Dirac equation (in quantized form) has been of fundamental importance to the development of modern field theories and many-particle physics [2, 3, 4, 5]. Interestingly, even today analytic solutions are very rare¹¹1E.g., spin- $1/2$ particle in a homogeneous magnetic field [6], the Dirac oscillator [7], an electromagnetic plane wave [8], or the Coulomb potential [9]. [6, 7, 8, 9]. In general, solutions have to be obtained numerically. They have gained in relevance by rapidly increasing computational resources, as well as the development of efficient numerical schemes. In a single particle picture, the (3+1)D Dirac equation is applicable whenever the external electromagnetic fields are strong enough to accelerate a spin- $1/2$ particle to relativistic speeds, but many-particle effects and electron-positron pair production can be neglected [10, 11]. This regime is reached in the study of light-matter interaction with the availability of short-pulse laser light in the (sub-)femtosecond range and intensities in excess of $10^{18}$ W/cm², which corresponds to the relativistic threshold²²2The electron is accelerated to relativistic speeds during one laser cycle. [10, 11]. Much of the physics in strong laser fields has been understood within a classical treatment of the relativistic electron. More recently, a numerical treatment of the quantum wave packet dynamics has become feasible [12, 13]. For an electron in a plane wave field, a wave-packet description of an electron reveals a drift of the wave packet in the direction of light propagation along with its spreading and shearing [12]. For this investigation a (2+1)D FFT-split-operator code was used. In such an approach, the propagation induced by the momentum part of the Hamiltonian is computed in momentum space, and the remainder in real-space, using fast Fourier transformation between the two representations [12]. The computational effort scales like $\mathcal{O}(N\ln N)$ where $N$ is the number of grid-points. An efficient code using operator splitting in real space, in combination with the exact characteristic solutions in each space direction, was introduced for the (3+1)D case of the Dirac equation recently [13]. It leads to the highly efficient operations count of $\mathcal{O}(N)$ .

In condensed matter physics, relativistic effects frequently are well accounted for by corrections to the Pauli equation derived from the Dirac equation [2]. Recently, however, (topological) condensed matter systems supporting effective Dirac and Majorana fermions have become a vivid playground for this community [14, 15, 17]. In particular, metallic surface states on topological insulator surfaces display 2D Dirac cone dispersion [16, 17, 18, 19]. A dynamic analysis of such Dirac fermions, e.g., in presence of effective electromagnetic fields calls for numerical schemes which faithfully represent the low-energy excitation spectrum. We have recently developed and applied such schemes for the (1+1)D and (2+1)D effective 2-component Dirac equation. In the (1+1)D case, we have presented a single-cone lattice scheme for which exactly absorbing open boundaries were derived [25]. In (2+1)D, we have first developed a staggered grid scheme with one additional artificial cone which, however, is able to preserve the linear dispersion of the free $m=0$ Dirac spin- $1/2$ particle along $x$ - and $y$ -axis [26]. More recently, a staggered grid single-cone scheme was developed for the two-component (2+1)D Dirac equation and used in studies of Dirac fermion dynamics on textured TI surfaces [27, 28]. This scheme, its properties, and its generalization to the four-spinor-component (3+1)D Dirac equation are topic of this paper. It operates entirely in real space and, due to symmetric staggering, treats space and time on equal terms. This scheme avoids the notorious fermion doubling characteristic of direct discretization of the Dirac derivative operators on a real space lattice [20]. To our knowledge, this is the first multi-D finite difference scheme with this property. This is achieved by redistributing the spinor components on a grid, staggered in space and time, such that individual spinor components sit on different (adjacent) time sheets. The proposed scheme, as will be shown below, shows an $\mathcal{O}(N)$ scaling behavior.

1.2 The Fermion Doubling Problem

Real-space finite-difference schemes for the Dirac equation have been plagued by the fermion doubling problem³³3It manifests itself in a non-monotonic dispersion relation leading to additional Dirac cones, in addition to the one at $k=0$ of the underlying continuum Hamiltonian. For $d$ discretized spatial dimensions one ends up with up to $2^{d}$ cones (“fermion flavors”).. It was shown rigorously by Nielsen and Ninomiya that the discretization of the Dirac equation on a regular space grid forbids a single chirally invariant fermion flavor without breaking one or more of the following assumptions: translational invariance, locality, and Hermiticity [20]. Obtaining additional spurious solutions due to discretization is not a problem specific to the Dirac equation, but can occur in all the cases where one discretizes a first derivative operator on a grid. Standard symmetric finite-difference approximations for the first derivatives are applied to preserve Hermiticity. They leave out the central grid point which, in turn, can take on arbitrary values without changing the specific value which the finite difference expression yields. For example, both a constant function and a function oscillating with the maximum frequency which can be resolved on the grid lead to zero for the central finite difference expression. Already the simplest model, the advection equation in 1D shows a non-monotonic behavior with a second minimum in the dispersion when a central approximation for the first spatial derivative is used (e.g., the forward-time central-space (FTCS) method). In other words, a symmetric first derivative in space utilizes twice the lattice spacing of the underlying grid thereby, in the language of solid state physics, shrinking the effective Brillouin zone in this direction by a factor of two.

Fermion doubling on a grid can be avoided by following basically two main strategies: (i) the incorporation or (ii) the complete avoidance of the central grid point in the scheme. As to the former, the use of a one-sided finite difference operator leads to so-called upwind schemes. For the Dirac equation this seems to work only in 1D, when a unitary time-evolution (conserving the norm) is to be maintained [29]. In this particular case it is equivalent to distributing the spinor components over a staggered grid [29, 25]. Unfortunately in higher dimensions this strategy either leads to non-Hermiticity and a non-unitary scheme or again to the introduction of fermion doubling. We note that the claim for a recently presented coordinate space operator splitting scheme was not troubled by fermion doubling could not be verified by us [13]. In contrast, we find that monotonic dispersion behavior is limited to $K=c\Delta_{t}/\Delta_{x}\leq 1/2$ for the (2+1)D and (3+1)D case within this scheme. However, in its present form it can only be executed for $K\in\mathbb{N}$ , because characteristic solutions for advection equations, e.g. $f(x,t+\Delta_{t})=f(x+K\Delta_{x},t)$ , are used. It might be possible to repair the scheme by using the reservoir technique, where the main idea is to wait $n$ time steps before updating the solution, which then allows one to use a smaller $\tilde{\Delta}_{t}=\Delta_{t}/n$ [30].

The presence of the central grid point in the scheme can also be enforced by adding artificial terms to the scheme, such as the momentum-dependent mass term suggested by Wilson, which lifts (splits) the spurious cones at high momentum, but retains a $k=0$ cone in good approximation of the continuum dispersion [32]. This strategy comes with its price. In case of the Dirac equation, it spoils chiral invariance for the physical mass $m=0$ case⁴⁴4For $m\neq 0$ chiral symmetry already is broken in the continuum problem. [32]. Alternatively, one can violate one of the other premises of the Nielsen-Ninomiya no-go theorem, for example, by breaking translational invariance in an extra dimension. This leads to domain-wall fermions [33, 20]. Interestingly, nature uses this trick in topological insulators [15].

The other option is to retain the centered approximation and to ensure that the disturbing central grid point cannot cause a problem: get rid of the point - get rid of the problem. This is achieved by the use of a staggered (”checkerboard”) grid whereby the grid is subdivided into subsets, with a particular spinor component defined on one of them, but not on the other(s). For the Dirac equation the idea of a staggered grid may be tracked back to Kogut and Susskind in the context of lattice QCD [31]. Here we utilize a staggered grid in space and time which treats time and space on equal footing, in accordance with the covariance of the underlying continuum model. Given initial conditions on a time slice, time propagation is executed in a leap-frog manner where the upper spinor components are computed from the lower ones and vice versa in an explicit recursive scheme.

The paper is organized as follows. In Section 2 we present the (3+1)D scheme and proceed with a discussion of its properties, such as norm conservation, stability, and the scaling behavior. In Section 3 we present this analysis for the corresponding scheme applied to (2+1)D. The gauge-invariant inclusion of an external electromagnetic potential and its consequences are outlined in Section 4. In Section 5 we use the (2+1)D scheme to show the scattering of an initial Gaussian wave packet at a Klein step. Summary and conclusions are given in Section 6.

2 Numerical Scheme for the (3+1)D Dirac Equation

Refer to caption — Figure 1: (color online). Space-time stencil for the first sub-step with explicit leap-frog time stepping for the space-time staggered finite difference scheme of the (3+1)D Dirac equation: The spinor components $A$ , $B$ , $C$ , and $D$ are defined on special positions of the space-time grid in order to allow a centered approximation of all first derivatives without inducing fermion doubling. The time propagation is sketched in horizontal direction. The z-axis is given the vertical direction. In the first step the two spinor components $A$ and $B$ , initially placed on time sheet $t^{j_{o}-1/2}$ , are propagated to time sheet $t^{j_{o}+1/2}$ , using components $C$ and $D$ located on time sheet $t^{j_{o}}$ . For component $A$ , the spinor components entering the centered symmetric difference quotients representing the partial derivatives are indicated in dark (blue) color. In the second half-step components $C$ and $D$ are propagated from time sheet $t^{j_{o}}$ to time sheet $t^{j_{o}+1}$ (not shown in the figure) using centered differences of components $A$ and $B$ on time sheet $t^{j_{o}+1/2}$ .

The single particle Dirac equation offers a relativistic description of spin $1/2$ particles capturing their particle-antiparticle and spin degree of freedom. A common representation is the (3+1)D Schrödinger form [34, 2]

i\hbar c\frac{\partial}{\partial x_{o}}{\vec{\psi}}(x_{\mu})=\left[-i\hbar c\sum_{j=1,2,3}\alpha_{j}\frac{\partial}{\partial x_{j}}+mc^{2}\beta+V(x_{\mu})\mbox{$1\hskip-3.69885pt1$}_{4}\right]{\vec{\psi}}(x_{\mu})~.

(1)

Here we use the $4\times 4$ Dirac matrices in Pauli-Dirac (“standard”) form [34]

\alpha_{j}=\sigma_{x}\otimes\sigma_{j}=\left(\begin{array}[]{cc}0&\sigma_{j}\\ \sigma_{j}&0\end{array}\right)

(2)

and

\beta=\sigma_{z}\otimes\mbox{$1\hskip-3.69885pt1$}_{2}=\left(\begin{array}[]{cc}\mbox{$1\hskip-3.69885pt1$}_{2}&0\\ 0&-\mbox{$1\hskip-3.69885pt1$}_{2}\end{array}\right)~,

(3)

expressed by the Pauli matrices $\sigma_{j}$ and the $2\times 2$ unit matrices $\mbox{$1\hskip-3.69885pt1$}_{2}$ . The first element in the direct product form accounts for the particle-antiparticle subspace and the second one represents the spin degree of freedom. As usual, $(x_{\mu})=(x_{o}=ct,x_{1}=x,x_{2}=y,x_{3}=z)$ is the space-time four vector in Minkowski space.
Writing the Dirac spinor as

{\vec{\psi}}(x_{\mu})=\left(\begin{array}[]{c}\psi_{1}(x_{\mu})\\ \psi_{2}(x_{\mu})\\ \psi_{3}(x_{\mu})\\ \psi_{4}(x_{\mu})\end{array}\right)=\left(\begin{array}[]{c}A(x_{\mu})\\ B(x_{\mu})\\ C(x_{\mu})\\ D(x_{\mu})\end{array}\right)

(4)

one observes that, in the standard representation, components $A$ and $B$ , as well as $C$ and $D$ , decouple. The remaining couplings are of two types: the mass and potential terms couple the time derivative of a particular component with itself, secondly, the time derivative of each of the upper two components couples to spatial derivatives of the lower two components. As an illustration we now write out the first line of (1):

i\hbar c\frac{\partial A}{\partial x_{o}}=-i\hbar c\left[\frac{\partial D}{\partial x_{1}}-i\frac{\partial D}{\partial x_{2}}+\frac{\partial C}{\partial x_{3}}\right]+(m^{2}c+V)\,A~.

(5)

With $j_{o}$ , $j_{x}$ , $j_{y}$ , and $j_{z}$ $\in\mathbb{Z}$ , we define the time step progression in units $\Delta_{t}$ , with $\Delta_{o}=c\Delta_{t}$ , for the spinor in two steps as shown in Fig. 1. First $A$ and $B$ are propagated from time-sheet $t-\Delta_{t}/2$ (i.e. index $j_{o}-1/2$ ) to $t+\Delta_{t}/2$ (i.e. index $j_{o}+1/2$ ), followed by the propagation of $C$ and $D$ from time-sheet $t$ ( $j_{o}$ ) to $t+\Delta_{t}$ ( $j_{o}+1$ ). One needs $x$ -, $y$ -, and $z$ -derivatives of $C$ and $D$ on time sheet $j_{0}$ , for the former, and of $A$ and $B$ on time sheet $j_{o}+1/2$ , for the latter update. Each spinor component is placed onto every other time sheet: $A$ and $B$ are put on half-integer time sheets $j_{o}\pm 1/2$ and $C$ and $D$ are put on integer time sheets, such as $j_{o}$ and $j_{o}+1$ . This allows one to compute symmetric time derivatives with step width $\Delta_{t}$ . If one then puts $A$ on a rectangular spatial grid, $j_{x},j_{y},j_{z}$ the implementation of symmetric $x$ - and $y$ -derivatives of $D$ requires the latter to be put on a face-centered rectangular (fcr) sublattice with grid points $(j_{x},j_{y}+1/2)$ , $(j_{x}+1/2,j_{y})$ in the $x$ - $y$ planes. In order to be able, in turn, to provide symmetric $x$ -and $y$ -derivatives of $A$ when $D$ is updated, the sublattice of $A$ has to be extended from simple rectangular to fcr (adding lattice points $j_{x}+1/2,j_{y}+1/2$ ). Since it is necessary to compute the $z$ -derivative of $C$ to obtain $A$ at the new time sheet, $C$ has to be given at the $j_{z}+1/2$ sites, sharing its $x$ - and $y$ -positions with $A$ . $A$ and $D$ may share the same $j_{z}$ grid points. Similarly, $B$ and $D$ share grid points on the $x$ - and $y$ - axis, while $B$ and $C$ share their grid points along the $z$ -axis. This determines all lattice sites on which individual spinor components need to be defined as follows (again with $j_{o}$ , $j_{x}$ , $j_{y}$ and $j_{z}$ $\in\mathbb{Z}$ )

$\displaystyle A(x_{\mu})$	$\displaystyle~\rightarrow~A^{j_{o}-1/2,j_{z}}_{j_{x},j_{y}},\quad A^{j_{o}-1/2,j_{z}}_{j_{x}+1/2,j_{y}+1/2}~,$
$\displaystyle B(x_{\mu})$	$\displaystyle~\rightarrow~B^{j_{o}-1/2,j_{z}+1/2}_{j_{x}+1/2,j_{y}},\quad B^{j_{o}-1/2,j_{z}+1/2}_{j_{x},j_{y}+1/2}~,$
$\displaystyle C(x_{\mu})$	$\displaystyle~\rightarrow~C^{j_{o},j_{z}+1/2}_{j_{x},j_{y}},\quad C^{j_{o},j_{z}+1/2}_{j_{x}+1/2,j_{y}+1/2}~,$
$\displaystyle D(x_{\mu})$	$\displaystyle~\rightarrow~D^{j_{o},j_{z}}_{j_{x}+1/2,j_{y}},\quad D^{j_{o},j_{z}}_{j_{x},j_{y}+1/2}~~.$	(6)

As the attentive reader probably has already recognized, the construction of the grid above was influenced by the use of the standard representation and the requirement of symmetric derivatives using single-lattice constant discretization with a fully symmetric structure in time and space. One can summarize: $A$ and $B$ ( $C$ and $D$ ) live on the same time sheets, $A$ and $D$ ( $B$ and $C$ ) on the same $z$ -sheets, and $A$ and $C$ ( $B$ and $D$ ) are defined on the same $x-y$ positions respectively. Let us discuss the number of grid-points available (occupied and unoccupied) vs. the number of grid-points where a spinor component is actually defined. The number of grid-points altogether (occupied and unoccupied) is $2^{4}=16$ times the number of grid-points available on a space-time-grid without the half-integer positions. Each component is defined on two times the number of integer positions. Thus we can conclude that $1/8$ of all available grid-points are occupied by one concrete spinor component (e.g. $A$ ). Remember that for an non-staggered grid, in three spatial dimensions, one obtains $2^{3}=8$ Dirac cones. Here, as we will see below 2.1, one has only one Dirac cone. Due to the staggering of the lattice (here also in time!), all partial derivatives can be performed in centered fashion about the grid points without modification of the ”lattice constants”.

With $j_{o}$ , $j_{x}$ , $j_{y}$ , and $j_{z}$ $\in\mathbb{Z}$ the scheme looks as follows. For the update of the first spinor component $A$ on the $(j_{x},j_{y})$ integer grid-points one has (cp. (5))

	$\displaystyle\frac{A^{j_{o}+1/2,j_{z}}_{j_{x},j_{y}}-A^{j_{o}-1/2,j_{z}}_{j_{x},j_{y}}}{\Delta_{o}}+\frac{C^{j_{o},j_{z}+1/2}_{j_{x},j_{y}}-C^{j_{o},j_{z}-1/2}_{j_{x},j_{y}}}{\Delta_{z}}+\frac{D^{j_{o},j_{z}}_{j_{x}+1/2,j_{y}}-D^{j_{o},j_{z}}_{j_{x}-1/2,j_{y}}}{\Delta_{x}}$
	$\displaystyle-i\frac{D^{j_{o},j_{z}}_{j_{x},j_{y}+1/2}-D^{j_{o},j_{z}}_{j_{x},j_{y}-1/2}}{\Delta_{y}}=\frac{1}{i\hbar c}(mc^{2}+V)^{j_{o},j_{z}}_{j_{x},j_{y}}\frac{A^{j_{o}+1/2,j_{z}}_{j_{x},j_{y}}+A^{j_{o}-1/2,j_{z}}_{j_{x},j_{y}}}{2}~.$		(7)

For the half-integer $(j_{x}+1/2,j_{y}+1/2)$ grid-points the update is

	$\displaystyle\frac{A^{j_{o}+1/2,j_{z}}_{j_{x}+1/2,j_{y}+1/2}-A^{j_{o}-1/2,j_{z}}_{j_{x}+1/2,j_{y}+1/2}}{\Delta_{o}}+\frac{C^{j_{o},j_{z}+1/2}_{j_{x}+1/2,j_{y}+1/2}-C^{j_{o},j_{z}-1/2}_{j_{x}+1/2,j_{y}+1/2}}{\Delta_{z}}+\frac{D^{j_{o},j_{z}}_{j_{x}+1,j_{y}+1/2}-D^{j_{o},j_{z}}_{j_{x},j_{y}+1/2}}{\Delta_{x}}$
	$\displaystyle-i\frac{D^{j_{o},j_{z}}_{j_{x}+1/2,j_{y}+1}-D^{j_{o},j_{z}}_{j_{x}+1/2,j_{y}}}{\Delta_{y}}=\frac{1}{i\hbar c}(mc^{2}+V)^{j_{o},j_{z}}_{j_{x}+1/2,j_{y}+1/2}\frac{A^{j_{o}+1/2,j_{z}}_{j_{x}+1/2,j_{y}+1/2}+A^{j_{o}-1/2,j_{z}}_{j_{x}+1/2,j_{y}+1/2}}{2}~.$		(8)

For $B$ one has on the $(j_{x}+1/2,j_{y})$ sub-grid

	$\displaystyle\frac{B^{j_{o}+1/2,j_{z}}_{j_{x}+1/2,j_{y}}-B^{j_{o}-1/2,j_{z}}_{j_{x}+1/2,j_{y}}}{\Delta_{o}}-\frac{D^{j_{o},j_{z}+1/2}_{j_{x}+1/2,j_{y}}-D^{j_{o},j_{z}-1/2}_{j_{x}+1/2,j_{y}}}{\Delta_{z}}+\frac{C^{j_{o},j_{z}}_{j_{x}+1,j_{y}}-C^{j_{o},j_{z}}_{j_{x},j_{y}}}{\Delta_{x}}$
	$\displaystyle+i\frac{C^{j_{o},j_{z}}_{j_{x}+1/2,j_{y}+1/2}-C^{j_{o},j_{z}}_{j_{x}+1/2,j_{y}-1/2}}{\Delta_{y}}=\frac{1}{i\hbar c}(mc^{2}+V)^{j_{o},j_{z}}_{j_{x}+1/2,j_{y}}\frac{B^{j_{o}+1/2,j_{z}}_{j_{x}+1/2,j_{y}}+B^{j_{o}-1/2,j_{z}}_{j_{x}+1/2,j_{y}}}{2}~,$		(9)

and on the $(j_{x},j_{y}+1/2)$ sub-grid one uses

	$\displaystyle\frac{B^{j_{o}+1/2,j_{z}}_{j_{x},j_{y}+1/2}-B^{j_{o}-1/2,j_{z}}_{j_{x},j_{y}+1/2}}{\Delta_{o}}-\frac{D^{j_{o},j_{z}+1/2}_{j_{x},j_{y}+1/2}-D^{j_{o},j_{z}-1/2}_{j_{x},j_{y}+1/2}}{\Delta_{z}}+\frac{C^{j_{o},j_{z}}_{j_{x}+1/2,j_{y}+1/2}-C^{j_{o},j_{z}}_{j_{x}-1/2,j_{y}+1/2}}{\Delta_{x}}$
	$\displaystyle+i\frac{C^{j_{o},j_{z}}_{j_{x},j_{y}+1}-C^{j_{o},j_{z}}_{j_{x},j_{y}}}{\Delta_{y}}=\frac{1}{i\hbar c}(mc^{2}+V)^{j_{o},j_{z}}_{j_{x},j_{y}+1/2}\frac{B^{j_{o}+1/2,j_{z}}_{j_{x},j_{y}+1/2}+B^{j_{o}-1/2,j_{z}}_{j_{x},j_{y}+1/2}}{2}~.$		(10)

The update for the $C$ -components living on $(j_{x},j_{y})$ writes

	$\displaystyle\frac{C^{j_{o}+1,j_{z}}_{j_{x},j_{y}}-C^{j_{o},j_{z}}_{j_{x},j_{y}}}{\Delta_{o}}+\frac{A^{j_{o}+1/2,j_{z}+1/2}_{j_{x},j_{y}}-A^{j_{o}+1/2,j_{z}-1/2}_{j_{x},j_{y}}}{\Delta_{z}}+\frac{B^{j_{o}+1/2,j_{z}}_{j_{x}+1/2,j_{y}}-B^{j_{o}+1/2,j_{z}}_{j_{x}-1/2,j_{y}}}{\Delta_{x}}$
	$\displaystyle-i\frac{B^{j_{o}+1/2,j_{z}}_{j_{x},j_{y}+1/2}-B^{j_{o}+1/2,j_{z}}_{j_{x},j_{y}-1/2}}{\Delta_{y}}=\frac{1}{i\hbar c}(-mc^{2}+V)^{j_{o}+1/2,j_{z}}_{j_{x},j_{y}}\frac{C^{j_{o}+1,j_{z}}_{j_{x},j_{y}}+C^{j_{o},j_{z}}_{j_{x},j_{y}}}{2}~,$		(11)

and for the $C$ -component on $(j_{x}+1/2,j_{y}+1/2)$

	$\displaystyle\frac{C^{j_{o}+1,j_{z}}_{j_{x}+1/2,j_{y}+1/2}-C^{j_{o},j_{z}}_{j_{x}+1/2,j_{y}+1/2}}{\Delta_{o}}+\frac{A^{j_{o}+1/2,j_{z}+1/2}_{j_{x}+1/2,j_{y}+1/2}-A^{j_{o}+1/2,j_{z}-1/2}_{j_{x}+1/2,j_{y}+1/2}}{\Delta_{z}}+\frac{B^{j_{o}+1/2,j_{z}}_{j_{x}+1,j_{y}+1/2}-B^{j_{o}+1/2,j_{z}}_{j_{x},j_{y}+1/2}}{\Delta_{x}}$
	$\displaystyle-i\frac{B^{j_{o}+1/2,j_{z}}_{j_{x}+1/2,j_{y}+1}-B^{j_{o}+1/2,j_{z}}_{j_{x}+1/2,j_{y}}}{\Delta_{y}}=\frac{1}{i\hbar c}(-mc^{2}+V)^{j_{o}+1/2,j_{z}}_{j_{x}+1/2,j_{y}+1/2}\frac{C^{j_{o}+1,j_{z}}_{j_{x}+1/2,j_{y}+1/2}+C^{j_{o},j_{z}}_{j_{x}+1/2,j_{y}+1/2}}{2}~.$		(12)

Finally, for the $D$ component on the $(j_{x}+1/2,j_{y})$ sub-grid the update is

	$\displaystyle\frac{D^{j_{o}+1,j_{z}}_{j_{x}+1/2,j_{y}}-D^{j_{o},j_{z}}_{j_{x}+1/2,j_{y}}}{\Delta_{o}}-\frac{B^{j_{o}+1/2,j_{z}+1/2}_{j_{x}+1/2,j_{y}}-B^{j_{o}+1/2,j_{z}-1/2}_{j_{x}+1/2,j_{y}}}{\Delta_{z}}+\frac{A^{j_{o}+1/2,j_{z}}_{j_{x}+1,j_{y}}-A^{j_{o}+1/2,j_{z}}_{j_{x},j_{y}}}{\Delta_{x}}$
	$\displaystyle+i\frac{A^{j_{o}+1/2,j_{z}}_{j_{x}+1/2,j_{y}+1/2}-A^{j_{o}+1/2,j_{z}}_{j_{x}+1/2,j_{y}-1/2}}{\Delta_{y}}=\frac{1}{i\hbar c}(-mc^{2}+V)^{j_{o}+1/2,j_{z}}_{j_{x}+1/2,j_{y}}\frac{D^{j_{o}+1,j_{z}}_{j_{x}+1/2,j_{y}}+D^{j_{o},j_{z}}_{j_{x}+1/2,j_{y}}}{2}~,$		(13)

and, for $(j_{x},j_{y}+1/2)$ ,

	$\displaystyle\frac{D^{j_{o}+1,j_{z}}_{j_{x},j_{y}+1/2}-D^{j_{o},j_{z}}_{j_{x},j_{y}+1/2}}{\Delta_{o}}-\frac{B^{j_{o}+1/2,j_{z}+1/2}_{j_{x},j_{y}+1/2}-B^{j_{o}+1/2,j_{z}-1/2}_{j_{x},j_{y}+1/2}}{\Delta_{z}}+\frac{A^{j_{o}+1/2,j_{z}}_{j_{x}+1/2,j_{y}+1/2}-A^{j_{o}+1/2,j_{z}}_{j_{x}-1/2,j_{y}+1/2}}{\Delta_{x}}$
	$\displaystyle+i\frac{A^{j_{o}+1/2,j_{z}}_{j_{x},j_{y}+1}-A^{j_{o}+1/2,j_{z}}_{j_{x},j_{y}}}{\Delta_{y}}=\frac{1}{i\hbar c}(-mc^{2}+V)^{j_{o}+1/2,j_{z}}_{j_{x},j_{y}+1/2}\frac{D^{j_{o}+1,j_{z}}_{j_{x},j_{y}+1/2}+D^{j_{o},j_{z}}_{j_{x},j_{y}+1/2}}{2}~.$		(14)

A more compact representation of the scheme, in analogy to the continuum form of the Dirac equation (1) in standard representation, can be given for the spinor components ${\psi}_{k}(x_{\mu}^{(k)})$ as follows

i\hbar cD_{o}{\psi}_{k}(x_{\mu}^{(k)})=\sum_{l=1}^{4}\left[\frac{\hbar c}{i}\sum_{j=1,2,3}\alpha_{j}D_{j}+\left(mc^{2}\beta+V(x_{\mu}^{(k)})\mbox{$1\hskip-3.69885pt1$}_{4}\right)T\right]_{kl}\psi_{l}(x_{\mu}^{(k)})~;~x_{\mu}^{(k)}\in G_{k},k=1,2,3,4.

(15)

Here we use the symmetric difference operators

D_{\nu}\psi(x_{\mu})=\frac{1}{\Delta_{\nu}}\left[\psi(x_{\mu}+\frac{\Delta_{\nu}}{2}{\hat{e}}_{\nu})-\psi(x_{\mu}-\frac{\Delta_{\nu}}{2}{\hat{e}}_{\nu})\right];~\nu=o,x,y,z

and the time-average operator

T\psi(x_{\mu})=\frac{1}{2}\left[\psi(x_{\mu}+\frac{\Delta_{o}}{2}{\hat{e}}_{o})+\psi(x_{\mu}-\frac{\Delta_{o}}{2}{\hat{e}}_{o})\right].

${\hat{e}}_{\nu}$ denotes the unit vector along the $\nu$ -axis in Minkowski space. The grids $G_{k}$ are defined as follows

$\displaystyle G_{1}:$	$\displaystyle~(j_{o},j_{x},j_{y},j_{z})~\cup~(j_{o},j_{x}+1/2,j_{y}+1/2,j_{z})~,$
$\displaystyle G_{2}:$	$\displaystyle~(j_{o},j_{x}+1/2,j_{y},j_{z}+1/2)~\cup~(j_{o},j_{x},j_{y}+1/2,j_{z}+1/2)~,$
$\displaystyle G_{3}:$	$\displaystyle~(j_{o}+1/2,j_{x},j_{y},j_{z}+1/2)~\cup~(j_{o}+1/2,j_{x}+1/2,j_{y}+1/2,j_{z}+1/2)~,$
$\displaystyle G_{4}:$	$\displaystyle~(j_{o}+1/2,j_{x}+1/2,j_{y},j_{z})~\cup~(j_{o}+1/2,j_{x},j_{y}+1/2,j_{z})~,$	(16)

for integer $j_{\nu}$ .

2.1 Dispersion relation

There is no fermion doubling for this scheme. For constant coefficients (potential and mass), this can be shown by considering the equivalent to continuum plane-wave solutions for frequency $\omega$ and $k$ -vector $(k_{x},k_{y},k_{z})$ ,

{\vec{\psi}}_{k_{\mu}}(x_{\mu})=\left(\begin{array}[]{c}\tilde{A}(k_{j})\\ \tilde{B}(k_{j})\\ \tilde{C}(k_{j})\\ \tilde{D}(k_{j})\end{array}\right)e^{-i\omega t}e^{i\sum_{j}k_{j}x_{j}}~.

(17)

On the lattice, given by the scheme, plane-wave solutions take the form

{\vec{\psi}}_{k_{\mu}}(x_{\mu})=\left(\begin{array}[]{c}\tilde{A}(k_{j})e^{i\omega\Delta_{t}/2}\\ \tilde{B}(k_{j})e^{i\omega\Delta_{t}/2}e^{ik_{x}\Delta_{x}/2}e^{ik_{z}\Delta_{z}/2}\\ \tilde{C}(k_{j})e^{ik_{z}\Delta_{z}/2}\\ \tilde{D}(k_{j})e^{ik_{x}\Delta_{x}/2}\end{array}\right)e^{-i\omega j_{t}\Delta_{t}}e^{i\sum_{j}k_{j}j_{j}\Delta_{j}}~,

(18)

where $\omega=ck_{o}$ , $\Delta_{o}=c\Delta_{t}$ , and $j_{t}=j_{o}$ . ⁵⁵5The relative shift of the sub-grids e.g. in the $x$ - $y$ -plane can be established by a translation by one-half lattice spacing in $x$ - or $y$ -direction. In Eq. (18) the translation in $x$ -direction is chosen. This does not induce an asymmetry in the scheme. In each of the four equations the common $e$ -factor for the center position in space and time can be canceled. This neither changes the zeros of the characteristic determinant nor the spinor-component ratios of the eigensolutions. Insertion into the scheme above translates any discrete $x_{o}=ct$ derivative into a multiplication by $-\frac{2i}{\Delta_{o}}\sin(\frac{\omega\Delta_{t}}{2})$ and every discrete partial position $j=x,y,z$ derivative into a multiplication by a factor $\frac{2i}{\Delta_{j}}\sin(\frac{k_{j}\Delta_{j}}{2})$ . A time average leads to a factor $\cos(\frac{\omega\Delta_{t}}{2})$ .
For fixed $\Delta_{\mu}$ , we therefore obtain the dispersion for the lattice from the continuum solution by the substitutions

E=\hbar\omega\rightarrow E(\omega)=\frac{2\hbar}{\Delta_{o}}\sin(\frac{\omega\Delta_{t}}{2})~,

(19)

p_{j}=\hbar k_{j}\rightarrow P_{j}(k_{j})=\frac{2\hbar}{\Delta_{j}}\sin(\frac{k_{j}\Delta_{j}}{2})~,

(20)

and

mc^{2}\rightarrow m(\omega)c^{2}=mc^{2}\cos(\frac{\omega\Delta_{t}}{2})~.

(21)

With these translation rules, the plane-wave solutions on the lattice can be obtained from the continuum solutions [34, 2]. In particular, the energy dispersion reads

E(\omega)=\pm\sqrt{m(\omega)^{2}c^{4}+|c{\bf P}|^{2}}~,

(22)

were $\mathbf{P}=(P_{1},P_{2},P_{3})$ is the momentum vector with the components defined in Eq. 20. This establishes the dispersion relation between $\omega$ and $k_{j}$ . Note that on the lattice the $k_{j}$ ’s are restricted to the values $\big{(}-\frac{\pi}{\Delta_{j}},+\frac{\pi}{\Delta_{j}}\big{]}$ and $\omega$ to $\big{(}-\frac{\pi}{\Delta_{o}},+\frac{\pi}{\Delta_{o}}\big{]}$ . Equation (22) can be solved for $\omega$ by taking the square and rearranging it to the final result for the dispersion relation

\hbar\omega=\pm\frac{2\hbar}{\Delta_{t}}\arcsin\big{[}X(k_{j})\big{]}\quad\mbox{with}\quad X(k_{j})=\sqrt{\frac{(mc^{2})^{2}+|c{\bf P}|^{2}}{(mc^{2})^{2}+(\frac{2\hbar}{\Delta_{t}})^{2}}}~.

(23)

Moreover,

E(\omega)=\pm\frac{2\hbar}{\Delta_{t}}X(k_{j})

(24)

and

m(\omega)c^{2}=mc^{2}\sqrt{1-X(k_{j})^{2}}~.

(25)

The spinor eigenket structure, as for the continuum model, is [2]

{\vec{\psi}}(k_{j})^{(1)}=N\left(\begin{array}[]{c}1\\ 0\\ \frac{cP_{z}}{E(\omega)+m(\omega)c^{2}}\\ \frac{c(P_{x}+iP_{y})}{E(\omega)+m(\omega)c^{2}}\end{array}\right)~,\qquad{\vec{\psi}}(k_{j})^{(2)}=N\left(\begin{array}[]{c}0\\ 1\\ \frac{c(P_{x}-iP_{y})}{E(\omega)+m(\omega)c^{2}}\\ -\frac{cP_{z}}{E(\omega)+m(\omega)c^{2}}\end{array}\right),

{\vec{\psi}}(k_{j})^{(3)}=N\left(\begin{array}[]{c}\frac{-cP_{z}}{\mid E(\omega)\mid+m(\omega)c^{2}}\\ -\frac{c(P_{x}+iP_{y})}{\mid E(\omega)\mid+m(\omega)c^{2}}\\ 1\\ 0\\ \end{array}\right)~,\qquad{\vec{\psi}}(k_{j})^{(4)}=N\left(\begin{array}[]{c}-\frac{c(P_{x}-iP_{y})}{\mid E(\omega)\mid+m(\omega)c^{2}}\\ \frac{cP_{z}}{\mid E(\omega)\mid+m(\omega)c^{2}}\\ 0\\ 1\end{array}\right)~.

(26)

$N$ depends on the choice of normalization.

The dispersion relation is shown in Fig. 2 for mass $m=0$ and two values of $k_{z}$ . The scheme displays a single cone centered around $k=0$ and preserves Kramers ( $k_{j}\rightarrow-k_{j}$ ) and time-reversal symmetry. Further symmetries depend on the ratios between the grid spacings $\Delta_{\nu}$ .

2.2 Stability

A preliminary discussion of stability can be based on the dispersion relation associated with plane wave solutions (in absence of an external potential), Eq. (23). Instability occurs when $|X(k_{j})|$ becomes larger than one, making the two solutions for $\omega$ to become a pair of complex conjugate numbers. Imposing

\Bigg{|}\frac{(mc^{2})^{2}+(2\hbar c)^{2}\sum_{j=x}^{z}\frac{1}{\Delta_{j}^{2}}}{(mc^{2})^{2}+(2\hbar c)^{2}\frac{1}{\Delta_{o}^{2}}}\Bigg{|}\leq 1~,

one arrives at the condition

r_{x}^{2}+r_{y}^{2}+r_{z}^{2}\leq 1~,

for $r_{j}=\Delta_{o}/\Delta_{j}$ . Thus, for equal spatial differences $\Delta_{x}=\Delta_{y}=\Delta_{z}=\Delta$ this gives the condition $r=\Delta_{o}/\Delta\leq 1/\sqrt{3}$ .

In order to discuss stability of this scheme in more generality, we seek a quantity (norm) which is conserved under the discrete time evolution. For this purpose we introduce several definitions and abbreviations to simplify notation. Let us recall on which grid positions the spinor components are defined

$\displaystyle A_{j}~\mbox{with}~j$	$\displaystyle~\in~(j_{x},j_{y},j_{z})~\cup~(j_{x}+1/2,j_{y}+1/2,j_{z})~,$
$\displaystyle B_{j}~\mbox{with}~j$	$\displaystyle~\in~(j_{x}+1/2,j_{y},j_{z}+1/2)~\cup~(j_{x},j_{y}+1/2,j_{z}+1/2)~,$
$\displaystyle C_{j}~\mbox{with}~j$	$\displaystyle~\in~(j_{x},j_{y},j_{z}+1/2)~\cup~(j_{x}+1/2,j_{y}+1/2,j_{z}+1/2)~,$
$\displaystyle D_{j}~\mbox{with}~j$	$\displaystyle~\in~(j_{x}+1/2,j_{y},j_{z})~\cup~(j_{x},j_{y}+1/2,j_{z})~,$	(27)

where $j_{x}$ , $j_{y}$ and $j_{z}$ are integer numbers. In what follows we consider a time step $\Delta_{t}$ where spinor $A$ and $B$ are propagated from $j_{o}-1/2$ to $j_{o}+1/2$ , and $C$ and $D$ are propagated from $j_{o}$ to $j_{o}+1$ , using the short-hand notation $A_{j}^{-},B_{j}^{-},C_{j}^{-},D_{j}^{-}\rightarrow A_{j}^{+},B_{j}^{+},C_{j}^{+},D_{j}^{+}$ . Furthermore we define a scalar product between spinor components $F^{j_{o}}_{j}$ and $G^{j_{o}^{\prime}}_{j+j^{\prime}}$ on the lattice in which a sum over all spatial sites $j$ is performed

(F^{j_{o}},G^{j_{o}^{\prime}})_{j^{\prime}}=\sum_{j}F^{j_{o}}_{j}G^{j_{o}^{\prime}*}_{j+j^{\prime}}=\sum_{j}F^{j_{o}}_{j-j^{\prime}}G^{j_{o}^{\prime}*}_{j}.

(28)

This short-hand notation involves a summation over all the spatial sublattice sites $j$ on the time sheet $j_{o}$ which belong to spinor component $F$ . Also, $j^{\prime}$ is limited to spatial shifts which connect the spatial sublattice of $F$ to the one associated with spinor component $G$ . Note that this scalar product depends on $j_{o}$ , $j_{o}^{\prime}$ , the relative position vector $j^{\prime}$ , and the two spinor components involved. Here, however, we only need to consider products for up to nearest neighbor sites. For simplicity, we write $(F^{j_{o}},G^{j_{o}^{\prime}})_{0}=(F^{j_{o}},G^{j_{o}^{\prime}})$ .
$j_{i}^{\pm}$ is defined as the vector shifting in the spatial direction $i$ by one half grid-spacing. Eg. $j_{x}^{\pm}:=(j_{x}\pm 1/2,j_{y},j_{z})$ . Finally we define for any spinor component

(\delta_{i}F)_{j}:=F_{j_{i}^{+}}-F_{j_{i}^{-}}~.

(29)

Now we introduce the conserved functional:

Lemma 1.

The quadratic functional

$\displaystyle E_{\mathbf{r}}$	$\displaystyle:=(A,A)+(B,B)+(C,C)+(D,D)$
	$\displaystyle~~-r_{z}\Re\Big{[}(\delta_{z}C,A)-(\delta_{z}D,B)\Big{]}-r_{x}\Re\Big{[}(\delta_{x}D,A)+(\delta_{x}C,B)\Big{]}$
	$\displaystyle~~+r_{y}\Re\Big{[}i(\delta_{y}D,A)-i(\delta_{y}C,B)\Big{]}$	(30)

is conserved under time propagation by the scheme (15). Here we used the notation $\mathbf{r}=(r_{x},r_{y},r_{z})$ .

Proof.

The proof follows the strategy from the (1+1)D paper [25]. We begin by writing the scheme in short-hand notation using the definitions given above:

$\displaystyle A^{+}_{j}-A^{-}_{j}+r_{z}(C^{-}_{j_{z}^{+}}-C^{-}_{j_{z}^{-}})+r_{x}(D^{-}_{j_{x}^{+}}-D^{-}_{j_{x}^{-}})-ir_{y}(D^{-}_{j_{y}^{+}}-D^{-}_{j_{y}^{-}})$	$\displaystyle=\frac{\Delta_{o}(mc^{2}+V)_{j}}{2i\hbar}(A^{+}_{j}+A^{-}_{j})~,$	(31)
$\displaystyle B^{+}_{j}-B^{-}_{j}-r_{z}(D^{-}_{j_{z}^{+}}-D^{-}_{j_{z}^{-}})+r_{x}(C^{-}_{j_{x}^{+}}-C^{-}_{j_{x}^{-}})+ir_{y}(C^{-}_{j_{y}^{+}}-C^{-}_{j_{y}^{-}})$	$\displaystyle=\frac{\Delta_{o}(mc^{2}+V)_{j}}{2i\hbar}(B^{+}_{j}+B^{-}_{j})~,$	(32)
$\displaystyle C^{+}_{j}-C^{-}_{j}+r_{z}(A^{+}_{{j}_{z}^{+}}-A^{+}_{{j}_{z}^{-}})+r_{x}(B^{+}_{{j}_{x}^{+}}-B^{+}_{{j}_{x}^{-}})-ir_{y}(B^{+}_{{j}_{y}^{+}}-B^{+}_{{j}_{y}^{-}})$	$\displaystyle=\frac{\Delta_{o}(-mc^{2}+V)^{+}_{j}}{2i\hbar}(C^{+}_{j}+C^{-}_{j})~,$	(33)
$\displaystyle D^{+}_{j}-D^{-}_{j}-r_{z}(B^{+}_{j_{z}^{+}}-B^{+}_{j_{z}^{-}})+r_{x}(A^{+}_{j_{x}^{+}}-A^{+}_{j_{x}^{-}})+ir_{y}(A^{+}_{j_{y}^{+}}-A^{+}_{j_{y}^{-}})$	$\displaystyle=\frac{\Delta_{o}(-mc^{2}+V)^{+}_{j}}{2i\hbar}(D^{+}_{j}+D^{-}_{j})~,$	(34)

where $V=V^{j_{o}}$ , $V^{+}=V^{j_{o}+1/2}$ , and analogously for $m$ .
Each of the Eqs. (31) to (34) is used to obtain an identity for the real part of a scalar product. Eq. (31) is multiplied by $(A^{+}_{j}+A^{-}_{j})^{*}$ and summed over all lattice sites $j$ . Similarly Eq. (32) is multiplied by $(B^{+}_{j}+B^{-}_{j})^{*}$ and so on. Adding up the real part of these four equations one obtains

	$\displaystyle(A^{+},A^{+})-(A^{-},A^{-})+(B^{+},B^{+})-(B^{-},B^{-})+(C^{+},C^{+})-(C^{-},C^{-})+(D^{+},D^{+})-(D^{-},D^{-})$
	$\displaystyle+\Re\left\{r_{z}\left[(A^{+}+A^{-},\delta_{z}C^{-})+(C^{+}+C^{-},\delta_{z}A^{+})-(B^{+}+B^{-},\delta_{z}D^{-})-(D^{+}+D^{-},\delta_{z}B^{+})\right]\right.$
	$\displaystyle\left.+r_{x}\left[(A^{+}+A^{-},\delta_{x}D^{-})+(D^{+}+D^{-},\delta_{x}A^{+})+(B^{+}+B^{-},\delta_{x}C^{-})+(C^{+}+C^{-},\delta_{x}B^{+})\right]\right\}$
	$\displaystyle-r_{y}\Im\left\{(A^{+}+A^{-},\delta_{y}D^{-})-(D^{+}+D^{-},\delta_{y}A^{+})-(B^{+}+B^{-},\delta_{y}C^{-})+(C^{+}+C^{-},\delta_{y}B^{+})\right\}=0~.$		(35)

Whenever possible, terms have been grouped in pairs of the form $\Re\{(F^{+}+F^{-},\delta_{k}G^{-})+(G^{+}+G^{-},\delta_{k}F^{+})\}$ or $\Im\{(F^{+}+F^{-},\delta_{k}G^{-})-(G^{+}+G^{-},\delta_{k}F^{+})\}$ . Each of the second terms in these pairs may be rewritten as follows (i.e. by summation by parts)

(G^{+}+G^{-},\delta_{k}F^{+})=\sum_{j}(G^{+}_{j}+G^{-}_{j})~(F^{+}_{j_{k}^{+}}-F^{+}_{j_{k}^{-}})^{*}=-(\delta_{k}G^{+},F^{+})-(\delta_{k}G^{-},F^{+})~,

We now observe that the contribution arising from mixed terms in time $+$ and $-$ cancel when taking, respectively, the real and imaginary part, and Eq. (35) takes the form $E^{+}-E^{-}=0$ with $E$ given in Eq. (30). ∎

Next we utilize the conserved functional (30) to prove stability of the numerical scheme for arbitrary space- and time-dependent mass and potential.

Proposition 2.

Let $r_{x}+r_{y}+r_{z}<1$ (e.g. using $r=r_{x}=r_{y}=r_{z}<1/3$ ). Then, the spinor elements stay finite and are bounded for all time by the estimate

(A,A)+(B,B)+(C,C)+(D,D)\leq\frac{E^{0}_{\mathbf{r}}}{1-r_{x}-r_{y}-r_{z}}~.

(36)

$E^{0}_{\mathbf{r}}$ is the squared norm (30) of the initial data $A^{-1/2}$ , $B^{-1/2}$ , $C^{0}$ , and $D^{0}$ .

Proof.

We start estimating the functional $E_{\mathbf{r}}$ :

$\displaystyle E_{\mathbf{r}}=E^{0}_{\mathbf{r}}$	$\displaystyle=(A,A)+(B,B)+(C,C)+(D,D)$
	$\displaystyle~~-r_{z}\Re\Big{[}(\delta_{z}C,A)-(\delta_{z}D,B)\Big{]}-r_{x}\Re\Big{[}(\delta_{x}D,A)+(\delta_{x}C,B)\Big{]}$
	$\displaystyle~~+r_{y}\Re\Big{[}i(\delta_{y}D,A)-i(\delta_{y}C,B)\Big{]}$
	$\displaystyle\geq(A,A)+(B,B)+(C,C)+(D,D)$
	$\displaystyle~~-r_{z}\Big{\|}\Re\Big{[}(\delta_{z}C,A)\Big{]}\Big{\|}-r_{z}\Big{\|}\Re\Big{[}(\delta_{z}D,B)\Big{]}\Big{\|}-r_{x}\Big{\|}\Re\Big{[}(\delta_{x}D,A)\Big{]}\Big{\|}-r_{x}\Big{\|}\Re\Big{[}(\delta_{x}C,B)\Big{]}\Big{\|}$
	$\displaystyle~~-r_{y}\Big{\|}\Re\Big{[}(\delta_{y}D,A)\Big{]}\Big{\|}-r_{y}\Big{\|}\Re\Big{[}(\delta_{y}C,B)\Big{]}\Big{\|}$
	$\displaystyle\geq(A,A)+(B,B)+(C,C)+(D,D)$
	$\displaystyle~~-r_{z}\Big{\|}\Re\Big{[}(C_{j_{z}^{+}},A)\Big{]}\Big{\|}-r_{z}\Big{\|}\Re\Big{[}(D_{j_{z}^{+}},B)\Big{]}\Big{\|}-r_{x}\Big{\|}\Re\Big{[}(D_{j_{x}^{+}},A)\Big{]}\Big{\|}-r_{x}\Big{\|}\Re\Big{[}(C_{j_{x}^{+}},B)\Big{]}\Big{\|}$
	$\displaystyle~~-r_{y}\Big{\|}\Re\Big{[}(D_{j_{y}^{+}},A)\Big{]}\Big{\|}-r_{y}\Big{\|}\Re\Big{[}(C_{j_{y}^{+}},B)\Big{]}\Big{\|}$
	$\displaystyle~~-r_{z}\Big{\|}\Re\Big{[}(C_{j_{z}^{-}},A)\Big{]}\Big{\|}-r_{z}\Big{\|}\Re\Big{[}(D_{j_{z}^{-}},B)\Big{]}\Big{\|}-r_{x}\Big{\|}\Re\Big{[}(D_{j_{x}^{-}},A)\Big{]}\Big{\|}-r_{x}\Big{\|}\Re\Big{[}(C_{j_{x}^{-}},B)\Big{]}\Big{\|}$
	$\displaystyle~~-r_{y}\Big{\|}\Re\Big{[}(D_{j_{y}^{-}},A)\Big{]}\Big{\|}-r_{y}\Big{\|}\Re\Big{[}(C_{j_{y}^{-}},B)\Big{]}\Big{\|}$
	$\displaystyle\geq(A,A)+(B,B)+(C,C)+(D,D)$
	$\displaystyle~~-\big{(}r_{x}+r_{y}+r_{z}\big{)}\bigg{(}(A,A)+(B,B)+(C,C)+(D,D)\bigg{)}~.$	(37)

In the last line we used the inequality $2|\Re\left\{(a,b)\right\}|\leq\left\|a\right\|^{2}+\left\|b\right\|^{2}$ . And this yields the estimate (36). ∎

By comparison with the result from the “reality condition” ( $\Im\{\omega\}=0$ ) of the dispersion (obtained for constant mass and potential), one observes that the condition in Proposition 2 for arbitrary space- and time-dependent coefficients gives a narrower bound for stability (i.e. $r<1/3$ vs. $r\leq 1/\sqrt{3}$ ). This is due to the coarse estimates in the short proof above. And, in fact, we shall now improve the stability condition to $r\leq 1/\sqrt{3}$ for the case $r=r_{x}=r_{y}=r_{z}$ :

Proposition 3.

Let $r=1/\sqrt{3}$ . Then the scheme (15) for arbitrary space- and time-dependent mass and potential is stable. It satisfies the following estimate for all time:

\|(\tilde{A},\tilde{B})\|^{2}+\|(\tilde{C},\tilde{D})\|^{2}\leq 2E^{0},

where $E^{0}$ is the “energy” of the initial data, as in (36). And on the l.h.s. we use the following grid-averaged norm of the spinor components:

	$\displaystyle\\|(\tilde{C},\tilde{D})\\|^{2}:=$	$\displaystyle\sum_{j}\left\|\frac{C_{j+(0,0,\frac{1}{2})}+D_{j+(\frac{1}{2},0,0)}-iD_{j+(0,\frac{1}{2},0)}}{2\sqrt{3}}+\frac{C_{j-(0,0,\frac{1}{2})}+D_{j-(\frac{1}{2},0,0)}-iD_{j-(0,\frac{1}{2},0)}}{2\sqrt{3}}\right\|^{2}$		(38)
		$\displaystyle\!\!\!\!\!+\sum_{j}\left\|\frac{D_{j+(\frac{1}{2},0,1)}-C_{j+(1,0,\frac{1}{2})}-iC_{j+(\frac{1}{2},\frac{1}{2},\frac{1}{2})}}{2\sqrt{3}}+\frac{D_{j+(\frac{1}{2},0,0)}-C_{j+(0,0,\frac{1}{2})}-iC_{j+(\frac{1}{2},-\frac{1}{2},\frac{1}{2})}}{2\sqrt{3}}\right\|^{2}~,$

with $j\in\mathbb{Z}^{3}\cup(\mathbb{Z}+1/2)^{2}\times\mathbb{Z}$ in the sum. $\|(\tilde{A},\tilde{B})\|^{2}$ is defined analogously, subject to an index shift detailed in (2.2).

This stability proof is more lengthy and hence deferred to A.

3 Numerical Scheme for the (2+1)D Dirac Equation

In some cases it is not necessary to solve the full (3+1)D Dirac equation. For example, under translational invariance in one spatial direction (when mass and potential terms are independent of this space coordinate) one can perform a partial Fourier transform within the scheme above, essentially leading to the substitution Eq. (20). The spin degree of freedom may be unimportant or it may be locked to the orbital motion, as is the case for the effective model for topological insulator (TI) surface states [15]. Then, a two component spinor in a (2+1)D model is adequate to capture the underlying physics. In some cases, low-dimensional models have been used simply because they are easier to handle.

Once more, we use the Schrödinger form of the (2+1)D Dirac equation

i\hbar c\partial_{t}\mbox{\boldmath$\psi$}(x,y,t)=\bigg{[}-i\hbar c\sigma_{x}\frac{\partial}{\partial x}-i\hbar c\sigma_{y}\frac{\partial}{\partial y}+\sigma_{z}m(x,y,t)+\mathbb{1}_{2}V(x,y,t)\bigg{]}\mbox{\boldmath$\psi$}(x,y,t)~,

(39)

where $\mbox{\boldmath$\psi$}(x,y,t)\equiv(u,v)\in\mathbb{C}^{2}$ is a $2$ -component spinor. For topological insulators, the two components are associated with spin= $1/2$ , whereby the physical spin quantization axis is $\mathbf{S}\propto\mathbf{\hat{z}}\times\mbox{\boldmath$\sigma$}$ . Here, $\mathbf{\hat{z}}$ is the unit-normal vector to the surface and $c$ constitutes the effective velocity of the quasi-particles in the effective model. $V\in\mathbb{R}$ is the scalar potential.

Earlier on and mainly for the treatment of rectangular structures we have presented a scheme which shows perfect dispersion along the coordinate axes for the free mass-less case, but with a second Dirac cone equally shared by the four corners of the first Brillouin zone [26, 28]. It is therefore ideally suited for TI surface states with rectangular structuring. Here a different staggering of the grid in space is applied to the $u$ - and $v$ -component of the spinor $\mbox{\boldmath$\psi$}=(u,v)$ which eliminates the second cone. It is represented graphically in Fig. 3 and follows the procedure for the (3+1)D case: A progression in time by $\Delta_{o}=c\Delta_{t}$ involves two initial, $t=j_{o}-1/2,j_{o}$ , and two final time sheets, $t=j_{o}+1/2,j_{o}+1$ . Spinor component $u$ lives on a face-centered rectangular (fcr) lattice (fcr-u) on the half-integer time sheets. Spinor component $v$ lives on the fcr lattice (fcr-v) formed by the set of midpoint positions of lines connecting nearest neighbors of the simple rectangular sub-lattice of $u$ , however, shifted vertically onto the integer time sheets. In other words, fcr-u is converted into fcr-v by a shift by $(i+1/2)\Delta_{x}{\hat{e}}_{x}$ (or $(i+1/2)\Delta_{y}{\hat{e}}_{y}$ ) followed by $\pm(j+1/2)\Delta_{o}{\hat{e}}_{o},i,j\in\mathbb{Z}$ , and vice versa.

The time progression by $\Delta_{o}$ is executed in two steps. First $u$ is propagated from time sheet $j_{o}-1/2$ to $j_{o}+1/2$ (both supporting fcr-u) by forming symmetric $x$ - and $y$ -derivatives of $v$ using fcr-v on time sheet $j_{o}$ . In the second step, $v$ on fcr-v on time sheet $j_{o}$ is propagated to fcr-v on time sheet $j_{o}+1$ using symmetric $x$ - and $y$ -derivatives of $u$ living on fcr-u on time sheet $j_{o}+1/2$ . According to the discussion above, one has to propagate $u_{j_{x},j_{y}}^{j_{o}-1/2}$ , $u_{j_{x}+1/2,j_{y}+1/2}^{j_{o}-1/2}$ , $v_{j_{x}+1/2,j_{y}}^{j_{o}}$ , and $v_{j_{x},j_{y}+1/2}^{j_{o}}$ for $j_{\nu}\in\mathbb{Z}$ . Using second order accurate symmetric approximations for the derivatives one arrives at the following scheme for this grid

	$\displaystyle\frac{u_{j_{x},j_{y}}^{j_{o}+1/2}-u_{j_{x},j_{y}}^{j_{o}-1/2}}{\Delta_{t}}=$	$\displaystyle\bigg{(}\frac{m+V}{i\hbar c}\bigg{)}_{j_{x},j_{y}}^{j_{o}}\frac{u_{j_{x},j_{y}}^{j_{o}+1/2}+u_{j_{x},j_{y}}^{j_{o}-1/2}}{2}$		(40)
	$\displaystyle-$	$\displaystyle\frac{(v_{j_{x}+1/2,j_{y}}^{j_{o}}-v_{j_{x}-1/2,j_{y}}^{j_{o}})}{\Delta_{x}}+i\frac{(v_{j_{x},j_{y}+1/2}^{j_{o}}-v_{j_{x},j_{y}-1/2}^{j_{o}})}{\Delta_{y}}~,$

	$\displaystyle\frac{v_{j_{x}-1/2,j_{y}}^{j_{o}+1}-v_{j_{x}-1/2,j_{y}}^{j_{o}}}{\Delta_{t}}=-$	$\displaystyle\bigg{(}\frac{m-V}{i\hbar c}\bigg{)}_{j_{x}-1/2,j_{y}}^{j_{o}+1/2}\frac{v_{j_{x}-1/2,j_{y}}^{j_{o}+1}+v_{j_{x}-1/2,j_{y}}^{j_{o}}}{2}$		(41)
	$\displaystyle-$	$\displaystyle\frac{(u_{j_{x},j_{y}}^{j_{o}+1/2}-u_{j_{x}-1,j_{y}}^{j_{o}+1/2})}{\Delta_{x}}-i\frac{(u_{j_{x}-1/2,j_{y}+1/2}^{j_{o}+1/2}-u_{j_{x}-1/2,j_{y}-1/2}^{j_{o}+1/2})}{\Delta_{y}}~,$

and analogous equations for the other sub-grid

	$\displaystyle\frac{u_{j_{x}-1/2,j_{y}-1/2}^{j_{o}+1/2}-u_{j_{x}-1/2,j_{y}-1/2}^{j_{o}-1/2}}{\Delta_{t}}=$	$\displaystyle\bigg{(}\frac{m+V}{i\hbar c}\bigg{)}_{j_{x}-1/2,j_{y}-1/2}^{j_{o}}\frac{u_{j_{x}-1/2,j_{y}-1/2}^{j_{o}+1/2}+u_{j_{x}-1/2,j_{y}-1/2}^{j_{o}-1/2}}{2}$		(42)
	$\displaystyle-$	$\displaystyle\frac{(v_{j_{x},j_{y}-1/2}^{j_{o}}-v_{j_{x}-1,j_{y}-1/2}^{j_{o}})}{\Delta_{x}}+i\frac{(v_{j_{x}-1/2,j_{y}}^{j_{o}}-v_{j_{x}-1/2,j_{y}-1}^{j_{o}})}{\Delta_{y}}~,$

	$\displaystyle\frac{v_{j_{x},j_{y}-1/2}^{j_{o}+1}-v_{j_{x},j_{y}-1/2}^{j_{o}}}{\Delta_{t}}=-$	$\displaystyle\bigg{(}\frac{m-V}{i\hbar c}\bigg{)}_{j_{x},j_{y}-1/2}^{j_{o}+1/2}\frac{v_{j_{x},j_{y}-1/2}^{j_{o}+1}+v_{j_{x},j_{y}-1/2}^{j_{o}}}{2}$		(43)
		$\displaystyle-\frac{(u_{j_{x}+1/2,j_{y}-1/2}^{j_{o}+1/2}-u_{j_{x}-1/2,j_{y}-1/2}^{j_{o}+1/2})}{\Delta_{x}}-i\frac{(u_{j_{x},j_{y}}^{j_{o}+1/2}-u_{j_{x},j_{y}-1}^{j_{o}+1/2})}{\Delta_{y}}~.$

Note that the $u$ -component defined for the discrete time indices $j_{o}-1/2\in\mathbb{Z}$ ‘lives’ on the discrete space gridpoints $(j_{x},j_{y})\in\mathbb{Z}^{2}$ and on $(j_{x}-1/2,j_{y}-1/2)\in\mathbb{Z}^{2}$ , whereas the $v$ -component defined for $j_{o}\in\mathbb{Z}$ is defined for space indices $(j_{x}-1/2,j_{y})\in\mathbb{Z}^{2}$ and $(j_{x},j_{y}-1/2)\in\mathbb{Z}^{2}$ .
For constant coefficients $m$ and $V$ von Neumann stability analysis reveals the stability of the finite difference scheme [36]. Moreover, the growth factor shows that an imaginary potential can be utilized to model an absorbing boundary layer. Fourier transformation in the spatial coordinates gives

	$\displaystyle\underbrace{\left(\begin{array}[]{cc}\frac{1}{\Delta_{t}}-\frac{m+V}{i\hbar c}&0\\ \frac{2i}{\Delta_{x}}\sin{\frac{k_{x}\Delta_{x}}{2}}+i\frac{2i}{\Delta_{y}}\sin{\frac{k_{y}\Delta_{y}}{2}}&\frac{1}{\Delta_{t}}+\frac{m-V}{i\hbar c}\end{array}\right)}_{=:S}\left(\begin{array}[]{c}\tilde{u}^{+}\\ \tilde{v}^{+}\end{array}\right)$		(48)
	$\displaystyle+\underbrace{\left(\begin{array}[]{cc}-\frac{1}{\Delta_{t}}-\frac{m+V}{i\hbar c}&\frac{2i}{\Delta_{x}}\sin{\frac{k_{x}\Delta_{x}}{2}}-i\frac{2i}{\Delta_{y}}\sin{\frac{k_{y}\Delta_{y}}{2}}\\ 0&-\frac{1}{\Delta_{t}}+\frac{m-V}{i\hbar c}\end{array}\right)}_{=:T}\left(\begin{array}[]{c}\tilde{u}^{-}\\ \tilde{v}^{-}\end{array}\right)=0~.$		(53)

One defines the amplification matrix $G=-S^{-1}T=S^{-1}S^{*}~$ . Its eigenvalues are the growth factors (written for $\hbar=1,c=1$ and $\Delta_{x}=\Delta_{y}=\Delta$ )

\lambda_{\pm}=P/2\pm\sqrt{(P/2)^{2}-Q}~,

(54)

with

P=\mbox{tr}[G]=-2\Big{[}(m^{2}-V^{2})\Delta_{t}^{2}-4c^{2}(1-\Delta_{t}^{2})-4c^{2}\Delta_{t}^{2}(\Delta_{y}^{2}\cos k_{x}\Delta_{x}+\Delta_{x}^{2}\cos k_{y}\Delta_{y})/(\Delta_{x}\Delta_{y})^{2}\Big{]}/N~,

(55)

Q=\mbox{det}[G]=(4c^{2}+(m^{2}-V^{2})\Delta_{t}^{2}-4icV\Delta_{t})/N

(56)

and

N=4c^{2}+(m^{2}-V^{2})\Delta_{t}^{2}+4icV\Delta_{t}~.

(57)

For $\Delta_{t}/\Delta\leq 1/\sqrt{2}$ and $m,V\in\mathbb{R}$ the absolute value of the growth factors is $1$ . Whereas for $V\in\mathbb{C}$ it depends on the sign of the imaginary part whether their absolute value is greater or lesser than $1$ . The latter can be used for absorbing boundary layers. The size of the absolute value of the growth (or in this case damping) factors is exemplified in Fig. 4.

The free-particle dispersion relation is revealed directly by using a plane wave ansatz
$u^{j_{o}+1}_{j_{x}+1,j_{y}+1}=e^{i(\omega\Delta_{t}-k_{x}\Delta_{x}-k_{y}\Delta_{y})}u^{j_{o}}_{j_{x},j_{y}}$ (and similarly for $v$ ) leading to the substitutions Eqs. (19)-(21) and the dispersion relation

\hbar\omega=\pm\frac{2\hbar}{\Delta_{t}}\arcsin\Bigg{\{}\sqrt{\frac{1}{(mc^{2})^{2}+2\hbar c/\Delta_{t}}\bigg{[}(mc^{2})^{2}+\Big{(}\frac{2\hbar}{\Delta_{x}}\sin\frac{k_{x}\Delta_{x}}{2}\Big{)}^{2}+\Big{(}\frac{2\hbar}{\Delta_{y}}\sin\frac{k_{y}\Delta_{y}}{2}\Big{)}^{2}\bigg{]}}\Bigg{\}}~.

(58)

The naive discrete expression for the norm $\left\|\psi\right\|:=\sqrt{\sum_{j_{x},j_{y}}|u_{j_{x},j_{y}}|^{2}+|v_{j_{x},j_{y}}|^{2}}$ in general oscillates around its mean value and is conserved in time only on average. We use a multiplication technique to identify an exactly conserved functional and prove stability for arbitrary space and time dependent coefficients, as for the (3+1)D case above. Let us again introduce a short-hand notation using

	$\displaystyle u_{j}~\mbox{with}~j$	$\displaystyle~\in~(j_{x},j_{y})~\cup~(j_{x}+1/2,j_{y}+1/2)~,$
	$\displaystyle v_{j}~\mbox{with}~j$	$\displaystyle~\in~(j_{x}+1/2,j_{y})~\cup~(j_{x},j_{y}+1/2)~.$

We define spatial difference operators as $(\delta_{x}f^{j_{o}})_{j_{x},j_{y}}=f_{j_{x}+1/2,j_{y}}^{j_{o}}-f_{j_{x}-1/2,j_{y}}^{j_{o}}$ , $(\delta_{y}f^{j_{o}})_{j_{x},j_{y}}=f_{j_{x},j_{y}+1/2}^{j_{o}}-f_{j_{x},j_{y}-1/2}^{j_{o}}$ and $(\delta_{\pm}f^{j_{o}})_{j_{x},j_{y}}=(\delta_{x}f^{j_{o}})_{j_{x},j_{y}}\pm i(\delta_{y}f^{j_{o}})_{j_{x},j_{y}}$ . We define the inner product $(u^{j_{o}},v^{j_{o}^{\prime}})_{j^{\prime}}:=\sum_{j}u_{j}^{j_{o}}v^{*j_{o}^{\prime}}_{j+j^{\prime}}=\sum_{j}u_{j-j^{\prime}}^{j_{o}}v^{*j_{o}^{\prime}}_{j}$ on $l^{2}(\mathbb{Z}^{2};\mathbb{C})$ and the notation $\left\|u^{j_{o}}\right\|^{2}:=(u^{j_{o}},u^{j_{o}})$ , with the sum over $j$ running over all spatial lattice points on the time sheet $j_{o}$ . $j^{\prime}$ again denotes a displacement vector connecting the two spatial sublattices of $u$ and $v$ .

In analogy to (30) we now define a conserved functional for the scheme (40)-(43):

Lemma 4.

Let $r=r_{x}=r_{y}$ . Then the functional

E_{r}^{j_{o}}:=\left\|u^{j_{o}+1/2}\right\|^{2}+\left\|v^{j_{o}+1}\right\|^{2}+r\Re\Big{[}(\delta_{-}u^{j_{o}+1/2},v^{j_{o}+1})\Big{]}=\mbox{const}=E_{r}^{0}~

(60)

is conserved under time propagation.

Proof.

The real part of the inner product of (40) with $({u}^{j_{o}+1/2}+{u}^{j_{o}-1/2})$ gives

\left\|u^{j_{o}+1/2}\right\|^{2}-\left\|u^{j_{o}-1/2}\right\|^{2}+r\Re\Big{[}(\delta_{+}v^{j_{o}},u^{j_{o}+1/2}+u^{j_{o}-1/2})\Big{]}=0~.

(61)

Analogously, (41) is multiplied by $(v^{j_{o}+1}+v^{j_{o}})^{*}$ and again the real part is taken to give

\left\|v^{j_{o}+1}\right\|^{2}-\left\|v^{j_{o}}\right\|^{2}+r\Re\Big{[}(\delta_{-}u^{j_{o}+1/2},v^{j_{o}+1}+v^{j_{o}})\Big{]}=0~.

(62)

Performing a summation by parts with vanishing ”boundary terms” at infinity gives

\Re\Big{[}(\delta_{x}v^{j_{o}},u^{j_{o}+1/2}+u^{j_{o}-1/2})\Big{]}=-\Re\Big{[}(\delta_{x}u^{j_{o}+1/2}+\delta_{x}u^{j_{o}-1/2},v^{j_{o}})\Big{]}~,

(63)

\Re\Big{[}(i\delta_{y}v^{j_{o}},u^{j_{o}+1/2}+u^{j_{o}-1/2})\Big{]}=\Re\Big{[}i~(u^{j_{o}+1/2}+u^{j_{o}-1/2},\delta_{y}v^{j_{o}})^{*}\Big{]}=\Re\Big{[}(i\delta_{y}u^{j_{o}+1/2}+i\delta_{y}u^{j_{o}-1/2},v^{j_{o}})\Big{]}~.

(64)

Finally, adding Eq. (61) and Eq. (62) leads to

\left\|u^{j_{o}+1/2}\right\|^{2}+\left\|v^{j_{o}+1}\right\|^{2}+r\Re\Big{[}(\delta_{-}u^{j_{o}+1/2},v^{j_{o}+1})\Big{]}=\left\|u^{j_{o}-1/2}\right\|^{2}+\left\|v^{j_{o}}\right\|^{2}+r\Re\Big{[}(\delta_{-}u^{j_{o}-1/2},v^{j_{o}})\Big{]}~.

(65)

∎

Similar to the (3+1)D case in Section 2.2, we can use this result to prove stability for an arbitrary space- and time-dependent mass and potential.

Proposition 5.

Let $r_{x}+r_{y}<1$ (e.g. using $r_{x}=r_{y}<1/2$ ). Then, the (2+1)D–scheme is stable and satisfies the estimate

\left\|u^{j_{o}+1/2}\right\|^{2}+\left\|v^{j_{o}+1}\right\|^{2}\leq\frac{E^{0}_{\mathbf{r}}}{1-r_{x}-r_{y}}

for all time.

Proof.

We start estimating the conserved functional

	$\displaystyle E^{0}_{\mathbf{r}}=E^{j_{o}}_{\mathbf{r}}:=\left\\|u^{j_{o}+1/2}\right\\|^{2}+\left\\|v^{j_{o}+1}\right\\|^{2}$	$\displaystyle+r_{x}\Re\Big{[}(\delta_{x}u^{j_{o}+1/2},v^{j_{o}+1})\Big{]}-r_{y}\Re\Big{[}(i\delta_{y}u^{j_{o}+1/2},v^{j_{o}+1})\Big{]}$
	$\displaystyle\geq\left\\|u^{j_{o}+1/2}\right\\|^{2}+\left\\|v^{j_{o}+1}\right\\|^{2}$	$\displaystyle-r_{x}~\Big{\|}\Re\Big{[}(\delta_{x}u^{j_{o}+1/2},v^{j_{o}+1})\Big{]}\Big{\|}-r_{y}~\Big{\|}\Re\Big{[}(i\delta_{y}u^{j_{o}+1/2},v^{j_{o}+1})\Big{]}\Big{\|}$
	$\displaystyle\geq\left\\|u^{j_{o}+1/2}\right\\|^{2}+\left\\|v^{j_{o}+1}\right\\|^{2}$	$\displaystyle-\big{(}r_{x}+r_{y}\big{)}\Big{(}\left\\|u^{j_{o}+1/2}\right\\|^{2}+\left\\|v^{j_{o}+1}\right\\|^{2}\Big{)}~,$

where we have used the inequality $2|\Re\left\{(a,b)\right\}|\leq\left\|a\right\|^{2}+\left\|b\right\|^{2}$ . And this yields the claimed result. ∎

Again, the condition in Proposition 5 is too restrictive for constant mass and potential. From the reality condition for the free-particle dispersion, one has $r_{x}^{2}+r_{y}^{2}\leq 1$ . For the special case $r=r_{x}=r_{y}$ this less restrictive stability condition reads $r\leq 1/\sqrt{2}$ . As in Proposition 3 we shall now prove that this stability condition actually also holds for an arbitrary space- and time-dependent mass and potential.

To this end we first define an averaged norm for the spinor components on the grid, pertaining to the special value $r=1/\sqrt{2}$ :

\left\|\tilde{u}\right\|^{2}:=\sum_{j_{x},j_{y}}\bigg{|}\frac{u_{j_{x}+1,j_{y}}-iu_{j_{x}+1/2,j_{y}+1/2}}{2\sqrt{2}}+\frac{u_{j_{x},j_{y}}-iu_{j_{x}+1/2,j_{y}-1/2}}{2\sqrt{2}}\bigg{|}^{2}~.

(66)

Here and in the following equations the summation runs over $j=(j_{x},j_{y})\in\mathbb{Z}^{2}\cup(\mathbb{Z}+1/2)^{2}$ .

Lemma 6.

$\left\|\tilde{u}\right\|^{2}$ is a norm on $l^{2}\big{[}\mathbb{Z}^{2}\cup(\mathbb{Z}+1/2)^{2}\big{]}$ .

Proof.

Assume $\left\|\tilde{u}\right\|^{2}=0$ , which leads to

u_{j_{x}+1,j_{y}}-iu_{j_{x}+1/2,j_{y}+1/2}+u_{j_{x},j_{y}}-iu_{j_{x}+1/2,j_{y}-1/2}=0~.

(67)

This is a leap-frog scheme $u_{j_{x}+1/2,j_{y}-1/2},u_{j_{x},j_{y}}\rightarrow u_{j_{x}+1,j_{y}},u_{j_{x}+1/2,j_{y}+1/2}$ . To show that is has no non-trivial solution in $l^{2}$ we observe that the real part of $u\big{|}_{\mathbb{Z}^{2}}$ is only coupled to the imaginary part of $u\big{|}_{(\mathbb{Z}+1/2)^{2}}$ :

\Re u_{j_{x},j_{y}}+\Im u_{j_{x}+1/2,j_{y}-1/2}=-(\Re u_{j_{x}+1,j_{y}}+\Im u_{j_{x}+1/2,j_{y}+1/2})=\Re u_{j_{x}+1,j_{y}+1}+\Im u_{j_{x}+3/2,j_{y}+1/2}=\cdots

(68)

The $l^{2}$ -summabilty of $u$ implies Eq. (68) $~=0$ giving

\Re u_{j_{x},j_{y}}=-\Im u_{j_{x}+1/2,j_{y}-1/2}~.

(69)

Now we shift (67) by $(\frac{1}{2},-\frac{1}{2})$ and take the imaginary part:

\Im u_{j_{x}+1/2,j_{y}-1/2}-\Re u_{j_{x}+1,j_{y}-1}=-(\Im u_{j_{x}+3/2,j_{y}-1/2}-\Re u_{j_{x}+1,j_{y}})=\Im u_{j_{x}+3/2,j_{y}+1/2}-\Re u_{j_{x}+2,j_{y}}=\cdots=0~,

(70)

again due to the $l^{2}$ -summabilty. Combining (69) and (70) with their integer grid-shifts yields

\Re u_{j_{x},j_{y}}=-\Im u_{j_{x}+1/2,j_{y}-1/2}=\Re u_{j_{x}+1,j_{y}-1}=-\Im u_{j_{x}+3/2,j_{y}-3/2}=\cdots~.

(71)

By the $l^{2}$ -summabilty, Eq. (71) $=0$ .

Analogously, the imaginary part of $u\big{|}_{\mathbb{Z}^{2}}$ is only coupled to the real part of $u\big{|}_{(\mathbb{Z}+1/2)^{2}}$ . The same arguments as before lead to

\Im u_{j_{x},j_{y}}=\Re u_{j_{x}+1/2,j_{y}-1/2}=\Im u_{j_{x}+1,j_{y}-1}=\Re u_{j_{x}+3/2,j_{y}-3/2}=\cdots=0~.

(72)

Hence $\tilde{u}=0$ . ∎

Now we can formulate our stability result:

Proposition 7.

Let $r=r_{x}=r_{y}=1/\sqrt{2}$ hold in (40)-(43). Then this (2+1)D–scheme is stable and satisfies for all time

\left\|\tilde{u}\right\|^{2}+\left\|\tilde{v}\right\|^{2}\leq 2E^{0}~.

Here, $E^{0}$ is the ”energy”

\displaystyle E:=\sum_{j_{x},j_{y}}|u_{j_{x},j_{y}}|^{2}+|v_{j_{x}+1/2,j_{y}}|^{2}+\frac{1}{\sqrt{2}}\Re\bigg{[}\big{(}u_{j_{x}+1,j_{y}}-u_{j_{x},j_{y}}-iu_{j_{x}+1/2,j_{y}+1/2}+iu_{j_{x}+1/2,j_{y}-1/2},v_{j_{x}+1/2,j_{y}}\big{)}\bigg{]}

(73)

of the initial data $u^{-1/2},\,v^{0}$ (use $r=1/\sqrt{2}$ in Eq. (60)).

Proof.

We use the same strategy as in the (1+1)D version of the scheme [25]. We first define the auxiliary quantity

$\displaystyle\tilde{E}:$	$\displaystyle=\frac{1}{2}\sum_{j_{x},j_{y}}\bigg{\|}\frac{u_{j_{x}+1,j_{y}}-iu_{j_{x}+1/2,j_{y}+1/2}}{\sqrt{2}}+v_{j_{x}+1/2,j_{y}}\bigg{\|}^{2}+\frac{1}{2}\sum_{j_{x},j_{y}}\bigg{\|}\frac{u_{j_{x},j_{y}}-iu_{j_{x}+1/2,j_{y}-1/2}}{\sqrt{2}}-v_{j_{x}+1/2,j_{y}}\bigg{\|}^{2}$	(74)
	$\displaystyle=\sum_{j}\|u_{j}\|^{2}+\|v_{j}\|^{2}$
	$\displaystyle~~+\frac{1}{2\sqrt{2}}\sum_{j_{x},j_{y}}\bigg{[}~(u_{j_{x}+1,j_{y}}-iu_{j_{x}+1/2,j_{y}+1/2})~v_{j_{x}+1/2,j_{y}}^{}+(u_{j_{x}+1,j_{y}}-iu_{j_{x}+1/2,j_{y}+1/2})^{}v_{j_{x}+1/2,j_{y}}$
	$\displaystyle\qquad\qquad\qquad-(u_{j_{x},j_{y}}-iu_{j_{x}+1/2,j_{y}-1/2})~v_{j_{x}+1/2,j_{y}}^{}-(u_{j_{x},j_{y}}-iu_{j_{x}+1/2,j_{y}-1/2})^{}v_{j_{x}+1/2,j_{y}}\bigg{]}$
	$\displaystyle~~+\frac{i}{4}\sum_{j_{x},j_{y}}\bigg{[}u_{j_{x}+1,j_{y}}u_{j_{x}+1/2,j_{y}+1/2}^{}-u_{j_{x}+1,j_{y}}^{}u_{j_{x}+1/2,j_{y}+1/2}$
	$\displaystyle\qquad\qquad\quad+u_{j_{x},j_{y}}u_{j_{x}+1/2,j_{y}-1/2}^{}-u_{j_{x},j_{y}}^{}u_{j_{x}+1/2,j_{y}-1/2}\bigg{]}~,$

where the last summation equals zero because the terms cancel, e.g. the first term with $j_{x}=j_{y}=0$ cancels with the last one having $j_{x}=j_{y}=1/2$ . Thus, one gets $\tilde{E}=E$ .

Using $\frac{1}{4}\left\|a_{1}+a_{2}\right\|^{2}\leq\frac{1}{4}\big{(}\left\|a_{1}+b\right\|+\left\|a_{2}-b\right\|\big{)}^{2}\leq\frac{1}{2}\left\|a_{1}+b\right\|^{2}+\frac{1}{2}\left\|a_{2}-b\right\|^{2}~$ gives $\left\|\tilde{u}\right\|^{2}\leq\tilde{E}=E^{0}.$ The symmetry of the scheme and of $E$ in $u$ and $v$ yields also: $\left\|\tilde{v}\right\|^{2}\leq\tilde{E}=E^{0}$ . ∎

4 Gauge-Invariant Introduction of a Vector Potential

In the previous exposition, only the presence of a scalar potential has been considered. The treatment of a Dirac fermion under the influence of a general electromagnetic field, however, calls for the use of the electromagnetic four-vector potential, i.e. both scalar and vector potential. The incorporation of the vector potential ${\bf{\cal A}}(x,y,z,t)$ poses no problem for the two schemes presented above. It can be accomplished in the same fashion as discussed recently in conjunction with our (1+1)D and (2+1)D schemes, following the standard Peierls substitution for the spinor components $\psi$ on the lattice [37, 38, 25, 26],

\displaystyle\psi^{t_{j},z_{j}}_{x_{j},y_{j}}\rightarrow{\hat{\psi}}^{t_{j},z_{j}}_{x_{j},y_{j}}:=\psi^{t_{j},z_{j}}_{x_{j},y_{j}}\exp\{-ia^{t_{j},z_{j}}_{x_{j},y_{j}}\}~.

(75)

Here $\psi$ denotes a spinor component and the real phase $a^{t_{j},z_{j}}_{x_{j},y_{j}}$ is defined as the line integral of the vector potential ${\bf{\cal A}}$ , starting at an arbitrary, but fixed position $(x_{o},y_{o},z_{o})$ and ending on the lattice point $(x_{j},y_{j},z_{j})$ ,

a^{t_{j},z_{j}}_{x_{j},y_{j}}=\frac{q}{\hbar c}\int_{(x_{o},y_{o},z_{o})}^{(x,y,z)}d{\bf s}\cdot{\bf{\cal A}}({\bf s},t)\mid_{x=x_{j},y=y_{j},z=z_{j},t=t_{j}}~.

$q$ is the fermion electric charge. Here the grid notation for the more general (3+1)D case has been adopted. The (2+1)D case can be handled analogously, as discussed in detail for our earlier (2+1)D scheme [26].

The implications of the insertion of the Peierls phase-shifted wave functions into the respective numerical schemes can be summarized as follows [26]:

1.

Since the connectivity of grid points is determined by the structure of the $\alpha-$ or, respectively, the $\sigma-$ matrices, requirements on the grid structure are not affected.
2.

Since every spinor component is multiplied merely by a phase factor, all stability and convergence estimates, as well as the definition of the norm, can be carried over simply by replacing spinor components by their Peierls transformed $\psi\rightarrow{\hat{\psi}}$ .

The effect of the substitution on difference quotients (derivative terms)

\frac{\psi_{1}-\psi_{2}}{\Delta}\rightarrow\frac{e^{-ia_{1}}\psi_{1}-e^{-ia_{2}}\psi_{2}}{\Delta}~.

can be implemented by using the product rule for “differentiation on the lattice”

\frac{e^{-ia_{1}}\psi_{1}-e^{-ia_{2}}\psi_{2}}{\Delta}=f_{+}(a_{1},a_{2})\frac{\psi_{1}-\psi_{2}}{\Delta}+\frac{e^{-ia_{1}}-e^{-ia_{2}}}{\Delta}\frac{\psi_{1}+\psi_{2}}{2}~,

with the definition $f_{\pm}(a_{1},a_{2}):=(e^{-ia_{1}}\pm e^{-ia_{2}})/2$ . For the case of slow variation of the vector potential over the grid, the difference quotient for the exponential may be approximated by using the “chain rule” for the derivative on the grid

\frac{e^{-ia_{1}}-e^{-ia_{2}}}{\Delta}=f_{+}(a_{1},a_{2})\frac{i(a_{2}-a_{1})}{\Delta}+\frac{1}{\Delta}O((a_{1}-a_{2})^{3})~.

(76)

All symmetric averages of the structure

\frac{\psi_{1}+\psi_{2}}{2}\rightarrow\frac{e^{-ia_{1}}\psi_{1}+e^{-ia_{2}}\psi_{2}}{2}~.

may be rewritten into

\frac{e^{-ia_{1}}\psi_{1}+e^{-ia_{2}}\psi_{2}}{2}=f_{+}(a_{1},a_{2})\frac{\psi_{1}+\psi_{2}}{2}+f_{-}(a_{1},a_{2})\frac{\psi_{1}-\psi_{2}}{2}\approx f_{+}(a_{1},a_{2})\frac{\psi_{1}+\psi_{2}}{2}~,

(77)

whereby the $f_{-}$ term may be neglected for vector potentials which are sufficiently smooth on the grid.

For such smooth vector potentials, approximations Eqs. (76) and (77) lead us to the following simple rules:

Every partial time-derivative renormalizes the associated scalar potential term (which occurs in the same equation) according to

V^{t_{j}(,z_{j})}_{x_{j},y_{j}}\rightarrow V^{t_{j}(,z_{j})}_{x_{j},y_{j}}-\frac{q}{c}\int_{x_{o},y_{o}(,z_{o})}^{x_{j},y_{j}(,z_{j})}d{\bf s}\cdot\frac{{\bf{\cal A}}({\bf s},t_{j}+\Delta_{t}/2)-{\bf{\cal A}}({\bf s},t_{j}-\Delta_{t}/2)}{\Delta_{t}}~.

Every spatial $j$ -derivative term transforms according to

\frac{\psi_{j+1/2}-\psi_{j-1/2}}{\Delta_{j}}\rightarrow\frac{\psi_{j+1/2}-\psi_{j-1/2}}{\Delta_{j}}-\frac{iq}{\hbar c}{\cal A}_{j}\frac{\psi_{j+1/2}+\psi_{j-1/2}}{2}~.

The substitutions into the schemes above are straight-forward [25, 26]. Since the resulting expressions at the various stages of approximation are rather lengthy they are not given here in more detail.

5 Simulation of wave packet scattering at a Klein step

Although the (2+1)D scheme has been used in the simulation of Dirac fermions in electromagnetic textures, here we give an example for illustration [28, 27, 40]. In particular we consider the scattering of a Gaussian wave packet at a Klein step. The latter consists of a step of the scalar potential $V$ from zero in the left half-space of Figs. 6 and 7 to a finite value $V>0$ in the right half-space. We consider the fermion mass $m=0$ and mean energies $E>0$ of the wave packet of less than the step hight. This causes a resonant energy overlap of electronic states in the left half-space with hole-states in the right half-space. In this energy window there is a flip in the sign of the phase velocity $v_{ph}=\frac{\omega}{k}$ relative to the group velocity $v_{g}=\frac{\partial\omega}{\partial k}$ ( $\omega=E/\hbar$ ). In Fig. 6 we show an example for normal incidence of the wave packet onto the potential step, while Fig. 7 shows snapshots of the propagation for incidence under $37^{\circ}$ to the $y$ –axis. The simulation shows the nearly complete penetration of the wave packet into and through the barrier under normal incidence (Klein paradox) due to the choice of parameters. For general angle of incidence one has both a reflected and transmitted wave packet whereby the latter displays a negative index of refraction, as evident in Fig. 7 and is the characteristic of a meta-material. We note the faithful representation of the propagation direction of both the wave packet and the phase (indicated in color) due to the absence of spurious additional cones in spite of the presence of an abrupt potential step. Absorbing boundary layers have been used to suppress spurious backscattering from the simulation boundaries [27, 40].

6 Summary and Conclusions

We have presented a real-space finite-difference method for both the four-component (3+1)D and the two-component (2+1)D Dirac equation. In both cases, the scheme is capable of handling general time-dependent electromagnetic potential and mass terms. Stability is proven for this general case and an exactly conserved functional is identified. It can be interpreted as a probability density.

The proposed schemes share and combine several decisive advantages over previous methods. Most importantly, the fermion doubling problem is avoided by staggering the spinor components symmetrically both in space and time. The strictly monotonic dispersion relation features a single Dirac cone. The scheme is second order accurate and, because it is explicit, shows an optimal linear scaling behavior. It lends itself to a convenient introduction of absorbing boundary conditions via boundary regions with imaginary scalar potential contributions, following a strategy discussed elsewhere [26]. Furthermore, this scheme which we have developed explicitly for one, two, and three spatial dimensions can be extended to higher dimensions.

The code readily is parallelized. This can be achieved, for example, by dividing the problem into spatial sub-domains, each computed on one CPU and communicating only the boundary grid-points. Although Maxwell’s equations have a different structure, the method presented here for the Dirac equation has much in common with the finite difference time domain (FDTD) method. Today it is applied with great success in computational electromagnetics [21, 22]. This close relationship also allows for the use of tools, initially developed for the FDTD time domain method, like the perfectly matched layer (PML), in its various variants, acting as an absorbing (open) boundary for finite-domain simulations [23, 24].

Since the presented schemes allow for the incorporation of the full electromagnetic four-vector potential, self-consistent dynamics (e.g., self-energy corrections) mediated by the electromagnetic interaction is readily achieved by solving the electromagnetic potential equations in parallel. This allows the utilization of all the advantages associated with the present approach and eliminates the need for the discretization of non-linear extensions of the Dirac equation. On the other hand, the present method may prove to be useful also for the discretization of such nonlinear Dirac equations, in particular, where nonlinearities arise from other types of interaction [39].

In summary, the combination of these favorable properties makes this approach highly suitable for the numerical treatment of Dirac fermion dynamics in space- and time-dependent external fields. As such, it should be useful to a variety of fields in physics, ranging from elementary particle, atomic, molecular, to condensed matter and astrophysics.

Acknowledgment: This work has been supported by the Austrian Science Foundation under Project No. I395-N16.

Appendix A Stability proof for the (3+1)D leap-frog staggered-grid scheme

In this appendix we provide the proof of Proposition 3. To this end we shall need the conserved functional (30) for the case $r=r_{x}=r_{y}=r_{z}=1/\sqrt{3}$ :

$\displaystyle E=$	$\displaystyle\sum_{j}\|A_{j}\|^{2}+\|B_{j+(\frac{1}{2},0,\frac{1}{2})}\|^{2}+\|C_{j+(0,0,\frac{1}{2})}\|^{2}+\|D_{j+(\frac{1}{2},0,0)}\|^{2}$	(78)
	$\displaystyle-\frac{1}{\sqrt{3}}\Re\bigg{[}\Big{(}C_{j+(0,0,\frac{1}{2})}-C_{j-(0,0,\frac{1}{2})}+D_{j+(\frac{1}{2},0,0)}-D_{j-(\frac{1}{2},0,0)}-iD_{j+(0,\frac{1}{2},0)}+iD_{j-(0,\frac{1}{2},0)},A_{j}\Big{)}\bigg{]}$
	$\displaystyle-\frac{1}{\sqrt{3}}\Re\bigg{[}\Big{(}-D_{j+(\frac{1}{2},0,1)}+D_{j+(\frac{1}{2},0,0)}+C_{j+(1,0,\frac{1}{2})}-C_{j+(0,0,\frac{1}{2})}+iC_{j+(\frac{1}{2},\frac{1}{2},\frac{1}{2})}-iC_{j+(\frac{1}{2},-\frac{1}{2},\frac{1}{2})},B_{j+(\frac{1}{2},0,\frac{1}{2})}\Big{)}\bigg{]}~,$

with $j\in\mathbb{Z}^{3}\cup(\mathbb{Z}+1/2)^{2}\times\mathbb{Z}$ in the sum.

Proof of Proposition 3.

Let us first define the auxiliary quantity

$\displaystyle\tilde{E}:$	$\displaystyle=\frac{1}{2}\sum_{j}\bigg{\|}\frac{C_{j+(0,0,\frac{1}{2})}+D_{j+(\frac{1}{2},0,0)}-iD_{j+(0,\frac{1}{2},0)}}{\sqrt{3}}-A_{j}\bigg{\|}^{2}$	(79)
	$\displaystyle+\frac{1}{2}\sum_{j}\bigg{\|}\frac{C_{j-(0,0,\frac{1}{2})}+D_{j-(\frac{1}{2},0,0)}-iD_{j-(0,\frac{1}{2},0)}}{\sqrt{3}}+A_{j}\bigg{\|}^{2}$
	$\displaystyle+\frac{1}{2}\sum_{j}\bigg{\|}\frac{D_{j+(\frac{1}{2},0,1)}-C_{j+(1,0,\frac{1}{2})}-iC_{j+(\frac{1}{2},\frac{1}{2},\frac{1}{2})}}{\sqrt{3}}+B_{j+(\frac{1}{2},0,\frac{1}{2})}\bigg{\|}^{2}$
	$\displaystyle+\frac{1}{2}\sum_{j}\bigg{\|}\frac{D_{j+(\frac{1}{2},0,0)}-C_{j+(0,0,\frac{1}{2})}-iC_{j+(\frac{1}{2},-\frac{1}{2},\frac{1}{2})}}{\sqrt{3}}-B_{j+(\frac{1}{2},0,\frac{1}{2})}\bigg{\|}^{2}~.$

It satisfies

$\displaystyle\tilde{E}$	$\displaystyle=\sum_{j}\|A_{j}\|^{2}+\|B_{j+(\frac{1}{2},0,\frac{1}{2})}\|^{2}+\|C_{j+(0,0,\frac{1}{2})}\|^{2}+\|D_{j+(\frac{1}{2},0,0)}\|^{2}$	(80)
	$\displaystyle~~+\frac{1}{2\sqrt{3}}\sum_{j}\bigg{[}-\big{(}C_{j+(0,0,\frac{1}{2})}+D_{j+(\frac{1}{2},0,0)}-iD_{j+(0,\frac{1}{2},0)}\big{)}A_{j}^{*}$
	$\displaystyle\qquad\qquad\quad~-\big{(}C_{j+(0,0,\frac{1}{2})}+D_{j+(\frac{1}{2},0,0)}-iD_{j+(0,\frac{1}{2},0)}\big{)}^{*}A_{j}$
	$\displaystyle\qquad\qquad\quad~+\big{(}C_{j-(0,0,\frac{1}{2})}+D_{j-(\frac{1}{2},0,0)}-iD_{j-(0,\frac{1}{2},0)}\big{)}A_{j}^{*}$
	$\displaystyle\qquad\qquad\quad~+\big{(}C_{j-(0,0,\frac{1}{2})}+D_{j-(\frac{1}{2},0,0)}-iD_{j-(0,\frac{1}{2},0)}\big{)}^{*}A_{j}$
	$\displaystyle\qquad\qquad\quad~+\big{(}D_{j+(\frac{1}{2},0,1)}-C_{j+(1,0,\frac{1}{2})}-iC_{j+(\frac{1}{2},\frac{1}{2},\frac{1}{2})})B_{j+(\frac{1}{2},0,\frac{1}{2})}^{*}$
	$\displaystyle\qquad\qquad\quad~+\big{(}D_{j+(\frac{1}{2},0,1)}-C_{j+(1,0,\frac{1}{2})}-iC_{j+(\frac{1}{2},\frac{1}{2},\frac{1}{2})})^{*}B_{j+(\frac{1}{2},0,\frac{1}{2})}$
	$\displaystyle\qquad\qquad\quad~-\big{(}D_{j+(\frac{1}{2},0,0)}-C_{j+(0,0,\frac{1}{2})}-iC_{j+(\frac{1}{2},-\frac{1}{2},\frac{1}{2})})B_{j+(\frac{1}{2},0,\frac{1}{2})}^{*}$
	$\displaystyle\qquad\qquad\quad~-\big{(}D_{j+(\frac{1}{2},0,0)}-C_{j+(0,0,\frac{1}{2})}-iC_{j+(\frac{1}{2},-\frac{1}{2},\frac{1}{2})})^{*}B_{j+(\frac{1}{2},0,\frac{1}{2})}\bigg{]}$
	$\displaystyle~~+\frac{1}{6}\sum_{j}\bigg{[}C_{j+(0,0,\frac{1}{2})}\big{(}D_{j+(\frac{1}{2},0,0)}^{}+iD_{j+(0,\frac{1}{2},0)}^{}\big{)}+C_{j+(0,0,\frac{1}{2})}^{*}\big{(}D_{j+(\frac{1}{2},0,0)}-iD_{j+(0,\frac{1}{2},0)}\big{)}$
	$\displaystyle\qquad\qquad+C_{j-(0,0,\frac{1}{2})}\big{(}D_{j-(\frac{1}{2},0,0)}^{}+iD_{j-(0,\frac{1}{2},0)}^{}\big{)}+C_{j-(0,0,\frac{1}{2})}^{*}\big{(}D_{j-(\frac{1}{2},0,0)}-iD_{j-(0,\frac{1}{2},0)}\big{)}$
	$\displaystyle\qquad\qquad+D_{j+(\frac{1}{2},0,1)}\big{(}-C_{j+(1,0,\frac{1}{2})}^{}+iC_{j+(\frac{1}{2},\frac{1}{2},\frac{1}{2})}^{}\big{)}+D_{j+(\frac{1}{2},0,1)}^{*}\big{(}-C_{j+(1,0,\frac{1}{2})}-iC_{j+(\frac{1}{2},\frac{1}{2},\frac{1}{2})}\big{)}$
	$\displaystyle\qquad\qquad+D_{j+(\frac{1}{2},0,0)}\big{(}-C_{j+(0,0,\frac{1}{2})}^{}+iC_{j+(\frac{1}{2},-\frac{1}{2},\frac{1}{2})}^{}\big{)}+D_{j+(\frac{1}{2},0,0)}^{*}\big{(}-C_{j+(0,0,\frac{1}{2})}-iC_{j+(\frac{1}{2},-\frac{1}{2},\frac{1}{2})}\big{)}\bigg{]}$
	$\displaystyle~~+\frac{i}{6}\sum_{j}\bigg{[}D_{j+(\frac{1}{2},0,0)}D_{j+(0,\frac{1}{2},0)}^{}-D_{j+(\frac{1}{2},0,0)}^{}D_{j+(0,\frac{1}{2},0)}+D_{j-(\frac{1}{2},0,0)}D_{j-(0,\frac{1}{2},0)}^{}-D_{j-(\frac{1}{2},0,0)}^{}D_{j-(0,\frac{1}{2},0)}\bigg{]}$
	$\displaystyle~~+\frac{i}{6}\sum_{j}\bigg{[}-C_{j+(1,0,\frac{1}{2})}C_{j+(\frac{1}{2},\frac{1}{2},\frac{1}{2})}^{}+C_{j+(1,0,\frac{1}{2})}^{}C_{j+(\frac{1}{2},\frac{1}{2},\frac{1}{2})}-C_{j+(0,0,\frac{1}{2})}C_{j+(\frac{1}{2},-\frac{1}{2},\frac{1}{2})}^{}+C_{j+(0,0,\frac{1}{2})}^{}C_{j+(\frac{1}{2},-\frac{1}{2},\frac{1}{2})}\bigg{]}$
	$\displaystyle=E~,$

since the last and the penultimate sum vanish because they are telescopic sums: In both sums the first (resp. second) term with $j=(0,0,0)$ cancels the forth (resp. third) term with $j=(\frac{1}{2},\frac{1}{2},0)$ . In the antepenultimate sum there occurs direct cancelation in two pairs of terms (involving $C_{j+(0,0,\frac{1}{2})}D_{j+(\frac{1}{2},0,0)}^{*}$ and $C_{j+(0,0,\frac{1}{2})}^{*}D_{j+(\frac{1}{2},0,0)}$ ) and the remaining terms are again telescopic sums.

Using $\frac{1}{4}\left\|a_{1}+a_{2}\right\|^{2}\leq\frac{1}{4}\big{(}\left\|a_{1}+b\right\|+\left\|a_{2}-b\right\|\big{)}^{2}\leq\frac{1}{2}\left\|a_{1}+b\right\|^{2}+\frac{1}{2}\left\|a_{2}-b\right\|^{2}~$ gives for the grid-averaged norm (defined in (38))

\|(\tilde{C},\tilde{D})\|^{2}\leq\tilde{E}=E~.

The symmetry of the scheme and of $E$ w.r.t. $(C,D)$ and $(A,B)$ yields also: $\|(\tilde{A},\tilde{B})\|^{2}\leq\tilde{E}=E$ . ∎

References

[1] P. A. M. Dirac, The Quantum Theory of the Electron, Proc. R. Soc. Lond. A 117 (1928) 610-624.
[2] J. J. Sakurai, Advanced Quantum Mechanics, Pearson Education India (2006).
[3] L. H. Ryder, Quantum Field Theory, Cambridge University Press (1996).
[4] C. Itzykson and J. B. Zuber, Quantum Field Theory, Dover Publications Inc. New York (2005).
[5] M. Srednicki, Quantum Field Theory, Cambridge University Press (2007).
[6] I. I. Rabi, Das freie Elektron im homogenen Magnetfeld nach der Diracschen Theorie, Z. Phys. 49 (1928) 507-511.
[7] K. Nikolsky, Das Oszillatorproblem nach der Diracschen Theorie, Z. Phys. 62 (1930) 677 681
[8] D. M. Wolkow, Über eine Klasse von Lösungen der Diracschen Gleichung, Z. Phys. 94 (1935) 250 260.
[9] H. A. Bethe, E. E. Salpeter, Quantum Mechanics of One- and Two-Electron Atoms, Plenum Publishing Corporation New York (1977) 63 71.
[10] Y. I. Salamin, S. X. Hu, K. Z. Hatsagortsyan, and C. H Keitel, Relativistic High-Power Laser-Matter Interactions, Phys. Rep. 427 (2006) 41-155, and references therein.
[11] A. Di Piazza, C. Müller, K. Z. Hatsagortsyan, and C. H. Keitel, Extremly High-Intensity Laser Interactions with Fundamental Quantum Systems, Rev. Mod. Phys. 84 (2012) 1177-1228, and references therein.
[12] G. R. Mocken and C. H. Keitel, FFT-Split-Operator Code for Solving the Dirac Equation in 2+1 Dimensions, Comp. Phys. Comm. 178 (2008) 868-882.
[13] F. Fillion-Gourdeau, E. Lorin, and A. D. Bandrauk, Numerical Solution of the Time-Dependent Dirac Equation in Coordinate Space without Fermion-Doubling, Comp. Phys. Comm. 183 (2012) 1403-1415.
[14] A. C. Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The Electronic Properties of Graphene, Rev. Mod. Phys. 81 (2009) 109-162, and references therein.
[15] X. L. Qi and S. C. Zhang, Topological Insulators and Superconductors, Reviews of Modern Physics 83 (2011) 1057-1110, and references therein.
[16] D. Hsieh, Y. Xia, D. Qian, L. Wray, F. Meier, J. H. Dil, J. Osterwalder, L. Patthey, A. V. Fedorov, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava, and M. Z. Hasan, Observation of Time-Reversal-Protected Single-Dirac-Cone Topological-Insulator States in Bi2Te3 and Sb2Te3, Phys. Rev. Lett. 103 (2009) 146401-1 - 146401-4.
[17] M. Z. Hasan and C. L. Kane, Colloquium: Topological Insulators, Rev. Mod. Phys. 82 (2010) 3045-3067, and references therein.
[18] Y. Xia, D. Qian, D. Hsieh, L. Wray, A. Pal, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava, and M. Z. Hasan, Observation of a Large-Gap Topological-Insulator Class with a Single Dirac Cone on the Surface, Nat. Phys. 5, (2009) 398-402.
[19] J. G. Analytis, J. H. Chu, Y. Chen, F. Corredor, R. D. McDonald, Z. X. Shen, and Ian R. Fisher, Bulk Fermi Surface Coexistence with Dirac Surface State in Bi2Se3: A Comparison of Photoemission and Shubnikov de Haas Measurements, Phys. Rev. B 81 (2010) 204507-1 - 204507-5.
[20] H. B. Nielsen and M. Ninomiya, A No-Go Theorem for Regularizing Chiral Fermions, Physics Letters B 105 (1981) 219-223.
[21] K. S. Yee, Numerical Solution of Initial Boundary Value Problems Involving Maxwell’s Equations in Isotropic Media, IEEE Transactions on Antennas and Propagation 14 (1966) 302-307.
[22] A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, Artech House Norwood, (2005).
[23] J. P. Berenger, A Perfectly Matched Layer for the Absorption of Electromagnetic Waves, J. of Comp. Phys. 114 (1994) 185-200.
[24] S. G. Johnson, Notes on Perfectly Matched Layers (PMLs), Lecture notes, Massachusetts Institute of Technology, Massachusetts (2008).
[25] R. Hammer, W. Pötz, and A. Arnold, A Dispersion and Norm Preserving Finite Difference Scheme with Transparent Boundary Conditions for the Dirac Equation in (1+1)D, J. Comp. Phys. 256 (2014) 728-747.
[26] R. Hammer and W. Pötz, Staggered-Grid Leap-Frog Scheme for the (2+1)D Dirac Equation, Comp. Phys. Commun. 185 (2014) 40-52.
[27] R. Hammer and W. Pötz, Dynamics of Domain-Wall Dirac Fermions on a Topological Insulator: A Chiral Fermion Beam Splitter, Phys. Rev. B 88 (2013) 235119-1 - 235119-11.
[28] R. Hammer, C. Ertler, and W. Pötz, Solitonic Dirac Fermion Wave Guide Networks on Topological Insulator Surfaces, Appl. Phys. Lett. 102 (2013) 193514-1 - 193514-4.
[29] R. Stacey, Eliminating Lattice Fermion Doubling, Phys. Rev. D 26 (1982) 468-472.
[30] F. Alouges, F. De Vuyst, G. Le Coq, and E. Lorin, The Reservoir Technique: A Way to Make Godunov-Type Schemes Zero or Very Low Diffuse. Application to Colella Glaz Solver, European Journal of Mechanics - B/Fluids 27 (2008) 643-664.
[31] J. Kogut und L. Susskind, Hamilton Formulation of Wilson s Lattice Gauge Theories, Phys. Rev. D 11 (1975) 395-408.
[32] P. H. Ginsparg and K. G. Wilson, A Remnant of Chiral Symmetry on the Lattice, Phys. Rev. D 25 (1982) 2649-2657.
[33] D. B. Kaplan, A Method for Simulating Chiral Fermions on the Lattice, Phys. Lett. B 288, (1992) 342-347.
[34] W. Greiner, Relativistic Quantum Mechanics: Wave Equations, Springer Berlin (2000).
[35] A. Borzì and E. Decker, Analysis of a Leap-Frog Pseudospectral Scheme for the Schrödinger Equation, Journal of Computational and Applied Mathematics 193 (2006) 65-88.
[36] J. C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, SIAM Philadelphia (2004).
[37] R. Peierls, Zur Theorie des Diamagnetismus von Leitungselektronen, Z. Phys. 80 (1933) 763-791.
[38] M. Graf und P. Vogl, Electromagnetic Fields and Dielectric Response in Empirical Tight-Binding Theory, Phys. Rev. B 51 (1995) 4940 4949, and references therein.
[39] J. Xu, S. Shao, and H. Tang, Numerical Methods for the Non-Linear Dirac Equation, J. Comp. Phys. 245 (2013) 131-149.
[40] R. Hammer, PhD Dissertation, Dec. 2013, (unpublished).

	$\displaystyle\\|(\tilde{C},\tilde{D})\\|^{2}:=$	$\displaystyle\sum_{j}\left\|\frac{C_{j+(0,0,\frac{1}{2})}+D_{j+(\frac{1}{2},0,0)}-iD_{j+(0,\frac{1}{2},0)}}{2\sqrt{3}}+\frac{C_{j-(0,0,\frac{1}{2})}+D_{j-(\frac{1}{2},0,0)}-iD_{j-(0,\frac{1}{2},0)}}{2\sqrt{3}}\right\|^{2}$		(38)
		$\displaystyle\!\!\!\!\!+\sum_{j}\left\|\frac{D_{j+(\frac{1}{2},0,1)}-C_{j+(1,0,\frac{1}{2})}-iC_{j+(\frac{1}{2},\frac{1}{2},\frac{1}{2})}}{2\sqrt{3}}+\frac{D_{j+(\frac{1}{2},0,0)}-C_{j+(0,0,\frac{1}{2})}-iC_{j+(\frac{1}{2},-\frac{1}{2},\frac{1}{2})}}{2\sqrt{3}}\right\|^{2}~,$

	$\displaystyle E^{0}_{\mathbf{r}}=E^{j_{o}}_{\mathbf{r}}:=\left\\|u^{j_{o}+1/2}\right\\|^{2}+\left\\|v^{j_{o}+1}\right\\|^{2}$	$\displaystyle+r_{x}\Re\Big{[}(\delta_{x}u^{j_{o}+1/2},v^{j_{o}+1})\Big{]}-r_{y}\Re\Big{[}(i\delta_{y}u^{j_{o}+1/2},v^{j_{o}+1})\Big{]}$
	$\displaystyle\geq\left\\|u^{j_{o}+1/2}\right\\|^{2}+\left\\|v^{j_{o}+1}\right\\|^{2}$	$\displaystyle-r_{x}~\Big{\|}\Re\Big{[}(\delta_{x}u^{j_{o}+1/2},v^{j_{o}+1})\Big{]}\Big{\|}-r_{y}~\Big{\|}\Re\Big{[}(i\delta_{y}u^{j_{o}+1/2},v^{j_{o}+1})\Big{]}\Big{\|}$
	$\displaystyle\geq\left\\|u^{j_{o}+1/2}\right\\|^{2}+\left\\|v^{j_{o}+1}\right\\|^{2}$	$\displaystyle-\big{(}r_{x}+r_{y}\big{)}\Big{(}\left\\|u^{j_{o}+1/2}\right\\|^{2}+\left\\|v^{j_{o}+1}\right\\|^{2}\Big{)}~,$

$\displaystyle\tilde{E}:$	$\displaystyle=\frac{1}{2}\sum_{j_{x},j_{y}}\bigg{\|}\frac{u_{j_{x}+1,j_{y}}-iu_{j_{x}+1/2,j_{y}+1/2}}{\sqrt{2}}+v_{j_{x}+1/2,j_{y}}\bigg{\|}^{2}+\frac{1}{2}\sum_{j_{x},j_{y}}\bigg{\|}\frac{u_{j_{x},j_{y}}-iu_{j_{x}+1/2,j_{y}-1/2}}{\sqrt{2}}-v_{j_{x}+1/2,j_{y}}\bigg{\|}^{2}$	(74)
	$\displaystyle=\sum_{j}\|u_{j}\|^{2}+\|v_{j}\|^{2}$
	$\displaystyle~~+\frac{1}{2\sqrt{2}}\sum_{j_{x},j_{y}}\bigg{[}~(u_{j_{x}+1,j_{y}}-iu_{j_{x}+1/2,j_{y}+1/2})~v_{j_{x}+1/2,j_{y}}^{}+(u_{j_{x}+1,j_{y}}-iu_{j_{x}+1/2,j_{y}+1/2})^{}v_{j_{x}+1/2,j_{y}}$
	$\displaystyle\qquad\qquad\qquad-(u_{j_{x},j_{y}}-iu_{j_{x}+1/2,j_{y}-1/2})~v_{j_{x}+1/2,j_{y}}^{}-(u_{j_{x},j_{y}}-iu_{j_{x}+1/2,j_{y}-1/2})^{}v_{j_{x}+1/2,j_{y}}\bigg{]}$
	$\displaystyle~~+\frac{i}{4}\sum_{j_{x},j_{y}}\bigg{[}u_{j_{x}+1,j_{y}}u_{j_{x}+1/2,j_{y}+1/2}^{}-u_{j_{x}+1,j_{y}}^{}u_{j_{x}+1/2,j_{y}+1/2}$
	$\displaystyle\qquad\qquad\quad+u_{j_{x},j_{y}}u_{j_{x}+1/2,j_{y}-1/2}^{}-u_{j_{x},j_{y}}^{}u_{j_{x}+1/2,j_{y}-1/2}\bigg{]}~,$

$\displaystyle\tilde{E}:$	$\displaystyle=\frac{1}{2}\sum_{j}\bigg{\|}\frac{C_{j+(0,0,\frac{1}{2})}+D_{j+(\frac{1}{2},0,0)}-iD_{j+(0,\frac{1}{2},0)}}{\sqrt{3}}-A_{j}\bigg{\|}^{2}$	(79)
	$\displaystyle+\frac{1}{2}\sum_{j}\bigg{\|}\frac{C_{j-(0,0,\frac{1}{2})}+D_{j-(\frac{1}{2},0,0)}-iD_{j-(0,\frac{1}{2},0)}}{\sqrt{3}}+A_{j}\bigg{\|}^{2}$
	$\displaystyle+\frac{1}{2}\sum_{j}\bigg{\|}\frac{D_{j+(\frac{1}{2},0,1)}-C_{j+(1,0,\frac{1}{2})}-iC_{j+(\frac{1}{2},\frac{1}{2},\frac{1}{2})}}{\sqrt{3}}+B_{j+(\frac{1}{2},0,\frac{1}{2})}\bigg{\|}^{2}$
	$\displaystyle+\frac{1}{2}\sum_{j}\bigg{\|}\frac{D_{j+(\frac{1}{2},0,0)}-C_{j+(0,0,\frac{1}{2})}-iC_{j+(\frac{1}{2},-\frac{1}{2},\frac{1}{2})}}{\sqrt{3}}-B_{j+(\frac{1}{2},0,\frac{1}{2})}\bigg{\|}^{2}~.$