Relativistic properties and invariants of the Du Fort–Frankel scheme for the one-dimensional Schrödinger equation

The Du Fort–Frankel scheme for the one-dimensional Schrödinger equation is shown to be equivalent, under a time-dependent unitary transformation, to the Ablowitz– Kruskal–Ladik scheme for the Klein–Gordon equation. The Schrödinger equation describes a non-relativistic quantum particle, while the Klein–Gordon equation describes a relativistic particle. The conditional convergence of the Du Fort–Frankel scheme to solutions of the Schrödinger equation arises because solutions of the Klein–Gordon equation only approximate solutions of the Schrödinger equation in the non-relativistic limit. The time- dependent unitary transformation is the discrete analog of the transformation that arises from seeking a non-relativistic limit using the interaction picture of quantum mechanics to decompose the Klein–Gordon Hamiltonian into the relativistic rest energy and a remainder. The Ablowitz–Kruskal–Ladik scheme is in turn decomposed into a quantum lattice gas automaton for the one-dimensional Dirac equation, which is also the one-dimensional discrete time quantum walk. This relativistic interpretation clariﬁes the origin of the known discrete invariant of the Du Fort–Frankel scheme as expressing conservation of probability for the 2-component wavefunction in the one-dimensional Dirac equation under discrete unitary evolution. It also leads to a second invariant, the matrix element of the evolution operator, whose imaginary part gives a discrete approximation to the expectation of the non-relativistic Schrödinger Hamiltonian. abstract Schrödinger equation (2) that generates evolution through the unitary operator exp ( − it H / ¯ h ) . In one spatial dimension this is accomplished by the system takes acting on the 2-component wavefunction (cid:11) 2 = u , d ) T . Eliminating d between these two equations leads back to the Klein– Gordon equation (see Section 6). The left hand sides of (20) describe the propagation of u and d along characteristics with speeds ± c , as appropriate for a massless particle. The mass terms on the right hand sides mix these two propagation directions, leading to the characteristic “Zitterbewegung” or trembling motion of a relativistic particle [1–3]. factorisation of the one-dimensional Klein–Gordon equation into a ﬁrst-order system with the self-adjoint Hamiltonian (21) relies on the factorisation ∂ tt − c 2 ∂ zz = (∂ t + c ∂ z )(∂ t − c ∂ z ) . More generally, the three-dimensional Dirac equation requires a 4-component wavefunction (cid:11) 4 , and factorises the three-dimensional Laplacian as ( σ ·∇ ) 2 = I 2 ∇ 2 , where I 2 is the 2 × 2 identity matrix, using three Pauli spin matrices


Introduction
The Schrödinger equation for a particle of mass m in a potential V is (1) where h is the reduced Planck's constant [1][2][3]. This equation describes the evolution of a wavefunction ψ(x, t), with the interpretation that |ψ(x, t)| 2 is the probability density for the particle to be located at position x at time t.
More generally, an abstract Schrödinger equation takes the form  where H is a self-adjoint linear operator with no explicit time dependence that acts on ψ . The solution of (2) may be written formally as ψ(·, t) = exp(−i(t/h)H)ψ (·, 0).

(3)
The evolution is unitary, since exp(−i(t/h)tH) is a unitary operator because H is self-adjoint, and time-reversible since exp(i(t/h)H) is the inverse of exp(−i(t/h)H). It thus preserves the total probability ||ψ|| 2 = |ψ(x, t)| 2 d 3 x (the evolution being unitary restricts the boundary conditions to be consistent with this conservation property). A good numerical scheme for a Schrödinger equation should also generate unitary and time-reversible evolution, but in discrete time steps t. In particular, unitary evolution implies stability in the 2 norm.
The leapfrog or second order difference scheme is [4][5][6][7][8][9] ih ψ n+1 j − ψ n−1 The wavefunction ψ for a particle moving in one dimension can be approximated by its values ψ n j = ψ( j x, n t) on a discrete spatial grid indexed by j at discrete times t = n t indexed by n. A natural discrete Hamiltonian matrix H for a free particle (one with V = 0) is defined by discretising the Laplacian in (1)  x 2 ψ n j+1 − 2ψ n j + ψ n j−1 , This paper is primarily concerned with the Du Fort-Frankel scheme introduced by Wu [10] ih ψ n+1 j − ψ n−1 This scheme is obtained by replacing 2ψ n j in (6) with the average ψ n+1 j + ψ n−1 j , the same replacement made by Du Fort & Frankel to construct their scheme for the real diffusion equation [11]. The scheme thus uses the four-point stencil shown in Fig. 1(b). It omits ψ n j located at the central point of the stencil. The Du Fort-Frankel scheme (7) is unconditionally stable, unlike the leapfrog scheme (6), and one can rearrange (7) to give an explicit formula for ψ n+1 j locally at each grid point (see section 3).
These properties make the Du Fort-Frankel scheme an attractive alternative to the Crank-Nicholson scheme. The latter uses the Cayley transform [12] I + 1 2 i( t/h)H −1 I − 1 2 i( t/h)H (8) to construct a unitary and time-reversible approximation to exp(−i tH/h) for the evolution between two time levels, given concretely by This scheme uses the six point stencil shown in Fig. 1(a). It is implicit, because each timestep requires the solution of a tridiagonal linear system involving the operator I + 1 2 i( t/h)H to determine the ψ n+1 j at all grid points j. This complicates a parallel implementation. The widely used Thomas algorithm for solving tridiagonal systems is inherently serial, though parallel alternatives, such as cyclic reduction or diagonalisation by discrete Fourier transform, are available. The implicitness incurs significantly more computational complexity and expense in multiple spatial dimensions, as some iterative algorithm will typically be required to solve the linear system. Like the above Du Fort-Frankel scheme, this is the complex version of the Crank-Nicolson [13] scheme for the real diffusion equation.
The Schrödinger equation (1) contains one time derivative but two spatial derivatives. The same asymmetry appears in the stencil for the Crank-Nicholson scheme, which contains three spatial points at each of two time levels. This asymmetry between space and time makes the Schrödinger equation incompatible with special relativity, which requires an equation to be invariant under Lorentz transformations that mix the space and time coordinates.
By contrast, space and time appear symmetrically in the four-point stencil for the Du Fort-Frankel scheme shown in Fig. 1(b). This paper will explore the properties of the Du Fort-Frankel scheme as a discretisation of a relativistic quantum wave equation, not of the Schrödinger equation (1).

Relativistic wave equations
The Schrödinger equation (1) for a free particle coincides with the Newtonian expression E = |p| 2 /(2m) for the energy of a free classical point particle with momentum p and mass m when one replaces E and p by differential operators: Making the same replacements in the relativistic energy-momentum relation E 2 = c 2 |p| 2 + m 2 c 4 for a free classical point particle gives the Klein-Gordon equation [1][2][3] ∂ tt u − c 2 ∇ 2 u = −(mc 2 /h) 2 u. (11) The line element ds 2 = |dx| 2 − c 2 dt 2 , and hence the differential operator ∂ tt − c 2 ∇ 2 on the left hand side, are invariant under Lorentz transformations. The whole equation is thus Lorentz-invariant if we treat u as a Lorentz-invariant scalar field.
By construction, the Klein-Gordon equation has plane wave solutions proportional to exp(i(k · x − ωt)) that satisfy the dispersion relation This is just a rewriting of the above energy-momentum relation using the frequency ω = E/h and wave vector k = p/h. Unlike the previous Schrödinger equation, there are two branches of solutions: one with positive frequencies ω ≥ mc 2 /h, and one with negative frequencies ω ≤ −mc 2 /h. Expanding the positive branch for |k| mc/h gives 4 ). (13) This recovers the dispersion relation ω =h|k| 2 /(2m) for the Schrödinger equation, but offset by the Compton frequency mc 2 /h associated with the relativistic rest energy mc 2 . The group velocity for both signs of ω is ∇ k ω = (c 2 /ω)k. (14) This satisfies the bound |∇ k ω| < c, and coincides with the expression v = p m 2 + |p| 2 /c 2 (15) for the velocity v of a relativistic particle with momentum p.
The fundamental solution of the one-dimensional Klein-Gordon equation for a concentrated Dirac δ-function source at (ξ, τ ) is [14] u( where J 0 is the Bessel function of the first kind of order zero. The support of this solution is thus confined to the boundary and interior of the light cone |x − ξ | = c|t − τ |, as one would expect from the bound on the group velocity. However, the support occupies the whole interior of the light cone, rather than being confined to the light cone itself, as for the corresponding fundamental solution of the wave equation. By contrast, the support of the fundamental solution ψ(x, t) = (2π i th/m) −1/2 exp i x 2 /(2th/m) (17) of the one-dimensional Schrödinger equation extends over the whole domain.
The Schrödinger equation (1) implies an equation for conservation of probability, with a flux called the probability current [1][2][3]. An overbar denotes a complex conjugate. The Klein-Gordon equation (11) implies another conservation law with the same flux: However, the quantity appearing in the time derivative cannot be interpreted as a probability density. It is not positive definite, as we can specify u and u t independently as initial conditions for the second-order Klein-Gordon equation.

Dirac equation
This difficulty led Dirac to seek a relativistic wave equation involving only first derivatives in space and time. This also fits the abstract Schrödinger equation (2) that generates evolution through the unitary operator exp(−itH/h). In one spatial dimension this is accomplished by the system speeds ±c, as appropriate for a massless particle. The mass terms on the right hand sides mix these two propagation directions, leading to the characteristic "Zitterbewegung" or trembling motion of a relativistic particle [1][2][3]. This factorisation of the one-dimensional Klein-Gordon equation into a first-order system with the self-adjoint Hamiltonian (21) relies on the factorisation ∂ tt − c 2 ∂ zz = (∂ t + c∂ z )(∂ t − c∂ z ). More generally, the three-dimensional Dirac equation requires a 4-component wavefunction 4 , and factorises the three-dimensional Laplacian as (σ · ∇) 2 = I 2 ∇ 2 , where I 2 is the 2 × 2 identity matrix, using the three Pauli spin matrices [1][2][3] The three-dimensional Dirac equation reduces to two independent copies of the system (20) for wavefunctions that depend only on z and t. It is convenient to choose z as the spatial coordinate because σ z is diagonal.
More generally, given a wavefunction that solves the Dirac equation, every component individually solves the Klein-Gordon equation (which just expresses the relativistic energy-momentum relation). However, the extra structure in the Dirac equation as a first order system changes the interpretation of the solution. In particular, the system (20a,b) implies a conservation law for the positive-definite probability density |u| 2 + |d| 2 ,

Interaction picture and non-relativistic limit
The abstract Schrödinger equation (2) describes the evolution of the wavefunction ψ in time. Physically meaningful quantities are expressed as expectations of self-adjoint linear operators A, known as observables, defined by The expectation is invariant under global changes of phase: ψ → ψe iα for any constant α.
The operators A usually have no explicit time dependence. For example, the probability of the particle described by the concrete Schrödinger equation (1) being in some subset X of space is given by taking A to be the indicator function for x ∈ X . It is common to normalise the wavefunction so that ||ψ|| = 1. The expectation E(A) then coincides with the matrix element A . However, the recovery of the expectation of the Hamiltonian for the Schrödinger equation in the non-relativistic limit depends crucially on the different normalisations of the relativistic and non-relativistic theories.
This Schrödinger picture is just one approach to formulating quantum mechanics. The interaction, or Dirac, picture decomposes the Hamiltonian into H = H 0 + H I . If φ is the wavefunction in the Schrödinger picture that satisfies in the interaction picture evolves according to [3] ih∂ t ψ = e iH 0 t/h H I e −iH 0 t/h ψ. (26) We use ψ for the transformed wavefunction so that we later obtain the familiar non-relativistic Schrödinger equation for ψ by applying the interaction picture to a relativistic wave equation for φ. More generally, any time-independent operator A in the Schrödinger picture (other than the Hamiltonian) becomes in the interaction picture. The interaction picture thus interpolates between the Schrödinger picture, for which H 0 = 0, and the Heisenberg picture, for which H 0 = H and H I = 0. The interaction picture is commonly used to construct the wavefunction φ using a convergent sequence of nested integrals when H I is treated as a perturbation to H 0 [3].
To apply the interaction picture to constructing the non-relativistic limit of the one-dimensional Dirac equation above, we diagonalise the algebraic mass terms mc 2 /h in the Hamiltonian (21) by applying the unitary transformation The Hamiltonian becomes We thus isolate the rest mass Hamiltonian H 0 = mc 2 I 2 and treat H I = H − H 0 as a perturbation. As H 0 is a multiple of the identity matrix I 2 it commutes with any other operator A, so A I (t) = A is unchanged according to the prescription (27).
We will find below that, although H I is not small relative to H 0 in the operator norm, ||H I || is much smaller than ||H 0 || for wavefunctions = (ψ + , ψ − ) T describing non-relativistic solutions. The wavefunctions ψ ± evolve according to The algebraic right hand side now only affects ψ − . The slowly-varying approximation holds in the non-relativistic limit, when the wavefunction oscillates with frequencies much smaller than the Compton frequency mc 2 /h. Making this approximation in (31b) allows us to solve for ψ − in terms of the gradient of ψ + , From this relation we can estimate |ψ − |/|ψ + | ∼ |p|/(2mc). The estimate is exact for plane wave solutions, which are eigenfunctions of the momentum operator p. The ratio of the momentum p to the momentum scale mc is small for a particle moving non-relativistically.
Substituting the relation (33) into the evolution equation (31a) for ψ + leads to the one-dimensional Schrödinger equation for a free particle: As in derivations of the Navier-Stokes equations from kinetic theory, this derivation of the Schrödinger equation eliminates the fast variable ψ − to obtain a closed evolution equation for the slow variable ψ + alone [15,16]. The same idea is used in geophysical fluid dynamics to construct balanced models, such as the quasigeostrophic equations. These describe lowfrequency motions in rapidly rotating fluid flows, while eliminating high-frequency inertia-gravity waves, by formulating a closed evolution equation for the height or potential vorticity field alone [17].
The matrix element of the interaction Hamiltonian is If we insert the approximation (33) for ψ − in terms of ψ + , H I coincides with the matrix element of the Schrödinger Hamiltonian for the ψ + state, The second step requires the usual integration by parts. Another useful expression, when ψ + and ψ − are related by (33), is Only the sum ||ψ|| 2 = ||ψ + || 2 + ||ψ − || 2 is conserved by the Klein-Gordon equation. Neither term is conserved separately. However, we can use (37) to replace ||ψ − || 2 by H I /(2mc 2 ) in the non-relativistic limit, which implies We can thus recover the expectation of the Hamiltonian in the Schrödinger equation from the invariants ||ψ|| 2 and H I of the Klein-Gordon equation in the interaction picture. The transformation (30) based on the interaction picture is equivalent to the substitution made by Pauli [2] to break the symmetry between φ + and φ − by seeking solutions whose time-dependence is close to e −imc 2 t/h . The same Wentzel-Kramers-Brillouin (WKB)-like ansatz [18] u( transforms the Klein-Gordon equation (11) into the complex telegraph equation We can now motivate the Schrödinger equation as describing solutions ψ(x, t) that vary sufficiently slowly in time, with frequencies much less than the Compton frequency mc 2 /h, that we can neglect the first term relative to the second term on the left hand side. We can also motivate rewriting the Du Fort-Frankel scheme (7) as a centred finite difference discretisation of the one-dimensional complex telegraph equation: This follows Du Fort & Frankel's observation that solutions of their scheme converge to solutions of a real telegraph equation, rather than to solutions of a diffusion equation [11]. The form (42) also leads to a more general class of Du Fort-Frankel schemes in which one takes an existing scheme with a finite difference approximation to a first time derivative ∂ t ψ , as illustrated by the second term on the left hand side of (42), and adds a constant multiple of the centred finite difference approximation to ∂ tt ψ with a constant chosen to optimise the stability of the overall scheme [19,20].
The transformation (30) can also be interpreted as a special case of the general invariance of the Schrödinger and Klein-Gordon equations under the gauge transformations when appropriate scalar V = ∂ t χ and vector A = ∇χ potentials are included. This transformation is thus distinct from the earlier transformation under (24) that changed the phase of the wavefunction by a constant. By taking the gauge function χ to be an affine function of t alone, one can absorb the relativistic energy due to the rest mass into a constant scalar potential, leaving the equivalent of the interaction Hamiltonian H I .

Dispersion relations for the Du Fort-Frankel and other schemes
We first rewrite Wu's Du Fort-Frankel scheme (7) as with a single dimensionless parameter λ = m x 2 /(h t). This becomes λ = (mc 2 /h) t on putting x = c t. Thus λ is the change in phase over a timestep t for an oscillation at the Compton frequency mc 2 /h. From now on it is more convenient to adopt dimensionless, so-called natural, units in which c = 1 and h = 1, so λ = m t.
Equation (44) is formally implicit, since ψ n+1 j appears on both sides, but we can solve locally for There is no need to solve a linear system that couples different grid points, as needed for the Crank-Nicolson scheme (9). Equation (45) may be rewritten more simply as There are two real frequencies ω for each k and α, so the scheme is indeed unconditionally stable. The left-hand side is an even function of ω t + α, and for k = 0 the two roots are ω = 0 and ω t = −2α. For small k there is a branch with ω ∼ k 2 /(2m), which agrees with the dispersion relation of the Schrödinger equation.
The dispersion relation for the Crank-Nicolson scheme (9) is This equation always has a solution for a real frequency ω C N for all real k and λ, so this scheme is also unconditionally stable. The dispersion relation for the leapfrog scheme in (6) is This only has a solution for a real frequency ω L F for all real k when λ ≥ 2, as found by Harmuth [5]. Moreover, there are then two real solutions, one with ω L F → 0 as k → 0, and a second with ω L F → −π / t as k → 0. For the special case λ = 2, shown in Fig. 2(b), the group velocity dω L F /dk is discontinuous at k = ±π / x. For larger λ, the group velocity vanishes at k = ±π / x for the leapfrog scheme, just as it does for the Du Fort-Frankel and Crank-Nicolson schemes. The dispersion relations for the numerical schemes are all 2π -periodic, unlike the dispersion relations for the partial differential equations (PDEs).   the CN dispersion relation. The blue dots mark the ends of the stable region for the LF scheme, whose lower branch passes through ω t = −π at k = 0.
The right-hand legend applies to both plots.

The Du Fort-Frankel scheme as a discretisation of the Klein-Gordon equation
To move the symmetry in the dispersion relation of the Du Fort-Frankel scheme to be around ω = 0, as in the unshifted Klein-Gordon equation, we make the change of variables ψ n j = u n j e inα .

Leapfrog formulation
The further substitution u n j = w n j exp(−inπ /2) This matches the abstract leapfrog, or second order difference, scheme in (4) for the discrete Hamiltonian H defined by  (5), as required to fit the four-point stencil in Fig. 1(b). The Du Fort-Frankel scheme is thus unitarily equivalent to the leapfrog, or second order difference, scheme for the three-point Hamiltonian (5) and a constant potential. However, the mass in the Hamiltonian is rescaled by 1/(sin α), which diverges as t → 0. A related result for the leapfrog scheme using the five-point stencil is described in the appendix.

Two-time-level formulation
The one-dimensional Klein-Gordon equation can be obtained by eliminating d from the one-dimensional Dirac equation Similarly, the discrete equation (52) that determines u at time level n + 1 in terms of u at the two previous time levels n and n − 1 can be rewritten as the system u n+1 j+1 = au n j + bd n j , relating two variables u and d on the 2 time levels n and n + 1, with the real coefficients a = cos α, b = sin α.
A closely related result was used previously to relate the Du Fort-Frankel scheme for the real diffusion equation to some 2-component lattice Boltzmann schemes [23,24]. The algorithm (62a,b) is equivalent to the composition of two unitary linear operators, streaming S and collisions C, defined by for periodic boundary conditions, where (· · · ) T indicates the (non-Hermitian) transpose of a vector, and is applied pointwise to each pair (u n j , d n j ) T . Each operator exactly solves one of the decoupled systems over a timestep t, so their composition SC defines a discrete unitary evolution that approximates the continuous system (61a,b). This system has arisen before as a quantum lattice gas automaton, a "quantum lattice Boltzmann equation," a quantum generalisation of a random walk on a lattice called the discrete time quantum walk, and as a realisation of the Feynman checkerboard model [22,[25][26][27][28][29][30]. It can also be derived by integrating the system (61a,b) along its characteristics for a timestep t and applying a unitary change of variables [22,31]. This discrete system approximates the solution of the continuous system (61a,b) with a local error of O ( t 2 ), and corresponding global error of O ( t), inherent in the splitting of (61a,b) into two subsystems (66) whose evolutions are given by the non-commuting operators S and C. Moreover, the ordering of the decomposition as SC rather than CS breaks the time-reversal symmetry of the original formulation (52) for u alone across 3 time levels.
Both defects are remedied by the symmetric second-order Strang [32] splitting C 1/2 SC 1/2 , where C 1/2 is a rotation by angle α/2 analogous to (65). Writing this out in full gives u n+1 j =â âu n j−1 +bd n j−1 +b âd n j+1 −bu n j+1 , with â = cos(α/2) and b = sin(α/2). Moreover, by writing the solution after n timesteps as this symmetric second-order splitting becomes equivalent to the original splitting SC applied to the transformed variables ũ j and d j defined by [22,33] ũ n j d n where C −1/2 is a rotation by angle −α/2 applied pointwise to each pair (u n j , d n j ) T .

Wu's discrete invariant
The above two-time-level formulations generate a discrete unitary evolution of the u and d variables from time level n to time level n + 1. They therefore conserve the total probability P n = j |u n j | 2 + |d n j | 2 = j |ũ n j | 2 + |d n j | 2 .
Rewriting this expression using the ψ n j variables gives the invariant found by Wu [10]: Wu [10] called this expression "energy", following common usage for a quadratic invariant of a linear evolution equation in PDE theory. We have shown here that it corresponds to the squared 2 norm of the solution of the first order system written using u n j and d n j , which is conserved by unitary evolution. We will find an invariant related to the quantum-mechanical energy, the expectation of the Hamiltonian operator, in the next two sections.

Two-time-level schemes in ψ ± variables
Having constructed two-time-level schemes in the u and d variables, we can transform them into two-time-level schemes for the ψ ± variables. To simplify notation we write ψ + ( j x, n t) = P n j and ψ − ( j x, n t) = M n j . Starting from the scheme (61a,b) and applying the unitary transformation with a phase rotation to undo the earlier rotation that led to the u and d variables, gives P n+1 This is a unitary, but only first-order accurate, finite difference scheme for the PDE system (31a,b). There is another firstorder accurate scheme based on the reversed splitting CS.
Applying the same transformation to the scheme (67a,b) based on the Strang splitting C 1/2 SC 1/2 gives This is a unitary and second-order accurate approximation to the PDE system (31a,b).
For a massless particle with α = 0 this is the Lax-Friedrichs [34,35] scheme for the hyperbolic system P t + M x = 0, The three variants differ only in the phase of M relative to P , since the matrix C containing the algebraic terms is diagonal in these variables. Eliminating M leads back to the original Du Fort-Frankel scheme for P in the form P n+1 j = −e 2iα P n−1 j + 1 2 1 + e 2iα P n j+1 + P n j−1 .
This is a second-order accurate finite difference approximation to the complex telegraph equation All three variants generate discrete unitary evolution, and thus conserve the total probability This expression is invariant under the different phases of M relative to P in the three variants. However, unlike the u and d formulations, there is no local reconstruction of M from P over the three time levels. One would have to solve a tridiagonal system involving M n j+1 − M n j−1 , so it is more efficient to solve the pair of equations over two time levels. The advantage of this formulation is that it cleanly separates the Schrödinger-like behaviour in P from the relativistic corrections in M. Only the combined probability (79) is a discrete invariant, but the magnitude of j |M n j | 2 gives a measure of the size of the relativistic corrections.

A discrete energy invariant
We can construct a second discrete invariant related to the quantum-mechanical energy, the expectation of the Hamiltonian operator [31]. We write the system (76a,b) schematically as where U is a unitary discrete-time evolution operator that has no explicit dependence on n. It evolves the wavefunction vector n = (P n 1 , . . . , P n N , M n 1 , . . . , M n N ) T of P and M values on the spatial grid at time level n. We define a complex-valued discrete inner product with complex conjugation applied to the first entry. The sum ranges over the 2N elements of n . This corresponds to treating and as 2-component vectors, and approximating the integral inner product † dz by the trapezoidal rule. The trapezoidal rule is exponentially accurate for domains with equally spaced points and periodic boundary conditions.
The complex-valued matrix element of the evolution operator is then conserved: since U † U = I is the identity operator as U is unitary. Moreover, U approximates the exact solution operator exp(−i tH) for the PDE system (31a,b) with O ( t 3 ) local error, since (76a,b) is a globally second-order accurate discretisation. Expanding for small t gives Here H is some discrete Hamiltonian that is compatible with U = exp(−i tH) to within an O ( t 3 ) error. All of the inner products in the last expression are real, since H is self-adjoint. The imaginary part of − U n / t is therefore a discrete invariant that approximates H n with an O ( t 2 ) relative error. The discrete matrix element U n is given concretely by We can subtract the real quantity I n = z j |P n j | 2 + |M n j | 2 to isolate the imaginary part: which rearranges to give This concretely defines the discrete H whose matrix element appears in (83). The |M n j | 2 and M n j (P n j−1 − P n j+1 ) terms are recognisable discrete approximations to terms in (35) for the matrix element of the continuous interaction Hamiltonian H I . The replacement of the mass prefactor 2m by 2 sin α, and the extra term involving M n j+1 − 2M n j + M n j−1 , are needed to create a discrete invariant.
Given this approximation to H I , we can use (39) to construct a discrete invariant that approximates the expectation E(H S ), The denominator contains the discrete squared norm ||ψ|| 2 = z j {|P n j | 2 + |M n j | 2 }.
The overall cos α prefactor in (86) improves the quantitative accuracy of the approximation to H I in the non-relativistic regime, although the difference remains O ( t 2 ) at fixed mass. To see why, consider using the centred finite difference approximation for ∂ z P , and then the slowly varying approximation ∂ z P = 2imM from (33), so Im j M n j P n Combining this expression with the |M n j | 2 term in (86) gives the leading order approximation The improved accuracy follows from the observation that  cos(α)(4 tan α − 2 sin α)

Numerical experiments
The Gaussian wavepacket is an exact solution of the Schrödinger equation (1) that depends on time only through the ratio τ = t/m, as expected from the structure of (1). Its width grows with time in proportion to (1 + 4σ 2 t 2 /m 2 ) 1/2 , as shown in Fig. 3. Figure 4 shows the discrete 2 difference between numerical solutions of the centred split form (76a,b) of the du Fort-Frankel scheme and the exact solution (93) at t = 9m for increasing numbers of grid points N. Each simulation was initialised with ψ + equal to the solution (93) at t = 0, and with ψ − proportional to the spatial derivative of this solution according to (33). While the differences initially decrease with increasing N, they eventually reach a plateau. Increasing m delays the onset of this plateau to larger N. By contrast, for each fixed m the numerical solutions of the du Fort-Frankel scheme show the expected second order convergence to reference solutions of the transformed Klein-Gordon system in (31a,b). The reference solutions were obtained by expressing the initial conditions for ψ ± as a finite Fourier series. The PDE system (31a,b) then decouples into pairs of ordinary differential equations for each discrete wavenumber k, whose solutions may be expressed using exponentials of 2 × 2 matrices. In other words, the Hamiltonian becomes purely algebraic, with only 2 × 2 blocks on the diagonal. Figure 4 confirms that the du Fort-Frankel scheme converges to solutions of the relativistic Klein-Gordon equation, rather than to solutions of the Schrödinger equation. This is exactly analogous to Du Fort & Frankel's observation that solutions of their scheme converge to solutions of a telegraph equation, not to solutions of a diffusion equation [11]. This behaviour is required by the Lax equivalent theorem [35,36] which states that a consistent finite difference scheme for a well-posed linear initial value problem is convergent if and only if it is stable. The du Fort-Frankel scheme is stable, but does not converge towards solutions of the Schrödinger equation, so it cannot be a consistent finite difference scheme for the Schrödinger equation.

This allows an explicit calculation of exp(−i(t/h)H) in the formal solution (3) of the abstract Schrödinger equation.
However, the du Fort-Frankel scheme will converge to solutions of the Schrödinger equation if one rescales both time and the particle mass with increasing resolution N. A more massive particle has a smaller velocity for the same momentum p, or the same characteristic wavenumber k = p/h, and is thus less influenced by relativistic effects. The ratio of ψ − to ψ + scales as p/(mc), which becomes smaller as m increases for fixed p.
The difference between the phases of the Schrödinger and Klein-Gordon solutions at time t = mτ is This source of error becomes O ( t 2 ), the same as the truncation error of the numerical scheme, if we take m = M/ t for fixed M and τ . The parameter λ = m t is now fixed as m increases and t decreases. However, the du Fort-Frankel dispersion relation (47) gives The leading term is the Schrödinger dispersion relation for a particle of mass M, so the scheme is still second-order accurate as k x → 0, even for m t = M fixed. The change in phase over a time t = mτ is ωt = (1/2)τ k 2 c 2 + O (k 4 t 2 ), which tends to (1/2)τ k 2 as t → 0 with τ and c = x/ t fixed. Figure 5 confirms that this scaling leads to the expected second-order convergence of the numerical solutions for ψ + towards the exact solution (93). However, Fig. 6 shows that M n j converges to −i/(2m)∂ x ψ of the exact solution with only first-order accuracy. This lower accuracy can be attributed to the explicit scaling of −i/(2m)∂ x ψ with 1/m. Visually indistinguishable results were obtained from the single three-time-level equation (77) by setting P 0 j and P 1 j to the exact solution (93) at t = 0 and t = t respectively. The solution ψ n j on space-time points with j + n even is decoupled from the solution on points with j + n odd for any scheme using the four-point stencil in Fig. 1(b). In principle, the solutions on the odd and even sets of points may drift  away from each other, but no such artifacts are visible in the wavefunctions in Fig. 3, or in the convergence behaviour. This may be because the different spatial Fourier modes are themselves decoupled in the du Fort-Frankel scheme with periodic boundary conditions, and the initial amplitude of the (−1) j zig-zag mode is negligibly small. Figure 7 shows the evolution of the differences between the exact value E(H S ) = 1/(2m) and two approximations: with increasing resolution at fixed m, due to finite relativistic effects.
The initialisation using (33) eliminates these relativistic effects in E n . Figure 7(a) shows second-order convergence of E n to E(H S ) under grid refinement for m = 4. However, the convergence behaviour for m = 32 shown in Fig. 7(b) is more complicated. The convergence behaviour of E n is shown more clearly in Fig. 8, which shows the convergence of E n for the discrete initial conditions towards E(H S ) under grid refinement for particles with masses ranging from 2 to 64. The convergence behaviour is simpler without the overall cos α prefactor in (86), for which the O (N −2 ) asymptotic regime is reached for smaller N. However, the prefactor in the error in the O (N −2 ) asymptotic regime is made significantly smaller by including the cos α prefactor, as shown by comparing the two plots in Figs. 8(a) and 8(b).

Conclusion
We have shown that the Du Fort-Frankel scheme for the Schrödinger equation for a non-relativistic free particle is equivalent, under a time-dependent unitary transformation, to the Ablowitz-Kruskal-Ladik (AKL) scheme for the Klein-Gordon equation from relativistic quantum mechanics. The AKL scheme is a variational integrator with a discrete action principle [22]. It is stable for x = c t. It then has the same maximum signal propagation speed c as the Klein-Gordon equation. The conditional convergence of the Du Fort-Frankel scheme to solutions of the Schrödinger equation arises because it is not sufficient to have a converged solution of the Klein-Gordon equation. Solutions of the latter only converge to solutions of the Schrödinger equation in the non-relativistic limit.
The required unitary transformation is a discrete analog of the transformation introduced by Pauli [2] to recover the Schrödinger equation as describing slowly varying solutions of the Klein-Gordon equation in the non-relativistic limit. It is equivalent to a shift in the zero-point of the energy so that the positive energy branch of the Klein-Gordon equation has its minimum energy at zero, rather than at mc 2 . We have shown that the same transformation can be motivated using the interaction picture of quantum mechanics, by decomposing the relativistic Hamiltonian into a multiple of the relativistic rest energy, H 0 = mc 2 I 2 , and a remainder H I . Moreover, the matrix element H I coincides with the matrix element H S for the Schrödinger equation when ψ − is related to ψ + by the slowly varying relation (33).
The AKL scheme for the Klein-Gordon equation is in turn unitarily equivalent to the Feynman checkerboard model, to a one-dimensional quantum lattice gas algorithm, to a "quantum lattice Boltzmann equation," and to the one-dimensional discrete time quantum walk [22,[25][26][27][28][29][30]. These are all discretisations of the one-dimensional Dirac equation, which is a first-order hyperbolic system for a pair of variables u and d, in contrast to the second-order scalar Klein-Gordon equation.
These discretisations thus give a unitary evolution of u and d from time level n to time level n + 1, in contrast to the AKL scheme that determines u at time level n + 1 in terms of the single variable u at the two preceding times n and n − 1. This reformulation as a first-order system gives a simple interpretation of the discrete invariant found by Wu [10] as describing conservation of the total probability j {|u n j | 2 + |d n j | 2 } for the two variables in this formulation, not just of j |u n j | 2 alone. A further unitary transformation gives a decomposition of the original Du Fort-Frankel scheme into discrete ψ ± variables. These lead more easily to the non-relativistic limit, in which |ψ − |/|ψ + | ∼ |p|/mc 1, than the formulation in the u and d variables that have comparable magnitudes. Any scheme that generates unitary evolution between two time levels through an operator with no explicit time dependence has a second discrete invariant, the matrix element of the evolution operator itself [22]. Expressing this invariant in the discrete ψ ± variables, and converting from the normalisation by ||ψ + || 2 + ||ψ − || 2 of the relativistic 2-component wavefunction to the non-relativistic normalisation by ||ψ + || 2 in the Schrödinger equation, gave a discrete invariant E n that approximates the expectation of the non-relativistic Schrödinger Hamiltonian operator H S .

Conflict of interest statement
The author declares that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.