We investigate the convergence properties of a three-dimensional quantum lattice Boltzmann scheme for the Dirac equation. These schemes were constructed as discretizations of the Dirac equation based on operator splitting to separate the streaming along the three coordinate axes, but their output has previously only been compared against solutions of the Schrödinger equation. The Schrödinger equation arises as the non-relativistic limit of the Dirac equation, describing solutions that vary slowly compared with the Compton frequency. We demonstrate first-order convergence towards solutions of the Dirac equation obtained by an independent numerical method based on fast Fourier transforms and matrix exponentiation.
The lattice Boltzmann approach to computational hydrodynamics relies on the Navier–Stokes equations describing slowly varying solutions of the Boltzmann equation. The same property holds when the continuous velocity space of the standard kinetic theory is truncated to obtain a discrete Boltzmann equation, which we write as [1–3] 1.1for i=0,…,n. The distribution functions fi propagate in the directions given by the discrete velocities ξi, and relax towards their equilibrium values through the action of the collision matrix Ωij.
The quantum lattice Boltzmann (QLB) schemes introduced by Succi & Benzi  exploit the structural similarities between the discrete Boltzmann equation and the Dirac equation of quantum electrodynamics . Both are linear, symmetric, hyperbolic systems with algebraic source terms, so numerical algorithms for the Dirac equation may be obtained by applying the same techniques that lead from the discrete Boltzmann equation to lattice Boltzmann schemes for hydrodynamics (see [2,6,7]). QLB schemes may also be interpreted as cellular automata that evolve probability amplitudes (qbits) rather than bits [8,9]. Palpacelli & Succi  have surveyed recent applications of QLB schemes to simulate the Schrödinger equation of non-relativistic quantum mechanics. The Schrödinger equation describes slowly varying solutions of the Dirac equation, in close analogy with the derivation of the Navier–Stokes equations as describing slowly varying solutions of the Boltzmann equation. However, attention has not previously been paid to the properties of QLB schemes for computing solutions to the Dirac equation itself.
2. The Dirac equation
The Dirac equation offers a quantum mechanical description of an electron that is compatible with special relativity . Its standard form is 2.1where c is the light speed, is the reduced Planck’s constant and is the Compton frequency for a particle of mass m. The wave function ψ is a column vector with four components, β is a 4×4 matrix, and α=(αx,αy,αz) is a collection of three 4×4 matrices, so that α⋅∇=αx∂x+αy∂y+αz∂z. The last term couples the wave function to an applied scalar potential V via the coefficient , where q>0 is the modulus of the charge on an electron. This sign convention was used in previous QLB schemes, but it differs from that in Berestetskii et al. .
The 4×4 matrices αi and β may be written in block form as 2.2where I is the 2×2 identity matrix. The off-diagonal blocks are Pauli spin matrices, 2.3which satisfy the well-known commutation relations  2.4where ϵijk is the alternating Levi-Civita tensor. The α matrices are all Hermitian, and hence diagonalizable, but the commutators of their non-zero blocks do not vanish. There is thus no basis in which the α matrices are simultaneously diagonal.
Writing the wave function as ψ=(ϕ+1,ϕ+2,ϕ−1,ϕ−2)T, the ± superscripts correspond to components with positive and negative energies. For a free particle, these oscillate in proportion to in the long-wave limit. The subscripts 1 and 2 label the two different spin states. To derive the Schrödinger equation, one writes to transform away the rest energy of the + states, assumes ∂tΦ−1,2≪2iωcΦ−1,2, equivalent to a non-relativistic limit, and substitutes the resulting adiabatic approximations for Φ−1,2 into the evolution equations for Φ+1,2. Palpacelli & Succi  give a derivation in this notation.
3. Quantum lattice Boltzmann schemes
The Dirac and discrete Boltzmann equations are both linear, symmetric, hyperbolic systems with algebraic source terms. However, the discrete Boltzmann equation is highly unusual among multi-dimensional hyperbolic systems in possessing one-dimensional characteristic curves of the form (x,t)=(x0+sξi,t0+s), as parametrized by s. All one may expect in general is the existence of three-dimensional characteristic surfaces in four-dimensional (x,t) space .
The three α matrices are not simultaneously diagonalizable, so the characteristics of the Dirac equation are indeed three-dimensional surfaces, the light cones of special relativity. Discontinuities in the initial conditions are confined to propagating along these surfaces, but there are no one-dimensional characteristic curves to integrate along. The approach that leads from the discrete Boltzmann equation to a lattice Boltzmann scheme (as in ) thus cannot be employed.
Instead, multi-dimensional QLB schemes use operator splitting to approximate solutions to the three-dimensional Dirac equation using solutions of one-dimensional Dirac equations in which spatial dependences in the other two directions are temporarily neglected. For instance, neglecting dependence on y and z leads to 3.1The factor of 1/3 on the right-hand side arises from dividing the algebraic terms into three equal parts, one associated with each spatial direction . Each of these one-dimensional Dirac equations does have characteristic curves, since the α matrices are Hermitian and thus diagonalizable. Equation (3.1) and its analogues in y and z may thus be evolved forward in time using the one-dimensional QLB algorithm described by Succi & Benzi  and Succi . It is convenient to choose units in which c=1 and , so Δx=Δt and ωc=m.
Existing QLB schemes begin with the Majorana form of the Dirac equation , but the same approach may be applied to the standard form (2.1) in which the algebraic terms are diagonal. The only disadvantage is that αy must be diagonalized by a complex matrix. The two approaches based on the Majorana and standard forms agree to within numerical round-off errors.
The three α matrices in the standard form (2.1) may be diagonalized as 3.2where † denotes the Hermitian transpose, by the three unitary matrices 3.3Thus, the streaming step in x consists of multiplying ψ by X†, streaming in the x-directions by ±Δx, and finally multiplying by X to undo the transformation. The first QLB schemes suffered from numerical artefacts owing to applying the X and X† matrices in the wrong order . We also use operator splitting to separate the spatial derivatives from the algebraic terms, which are diagonal in the standard representation. The + components are 3.4which we discretize using the Crank–Nicolson scheme to obtain 3.5The solution, and its analogue for ϕ−1,2(t+Δt), may be expressed in the manifestly unitary form 3.6Writing the transformed post-collisional wave function as , the QLB step for colliding and streaming in x may be written as 3.7The full QLB scheme combines this step with analogous collision and streaming steps in the y- and z-directions, as diagonalized by the Y and Z matrices.
4. Free particle
We apply the QLB scheme described above to initial conditions in which the positive energy, spin-up component is a spherically symmetric Gaussian wavepacket with spread Δ0, as in Palpacelli & Succi , 4.1The other three components and are initially set to zero. We study the dispersion of the wavepacket over time, as measured by the spread Δ defined by 4.2We approximate these integrals using the trapezium rule, which is exponentially accurate for equally spaced points in a domain with periodic boundary conditions.
According to the Schrödinger equation, the spread of this wavepacket grows as 4.3Figure 1a shows the measured spread for a free particle with m=0.35 and Δ0=14 in a cube with side length ℓ=200 and periodic boundary conditions. All these quantities are expressed in natural units with c=1 and . The cube was discretized using N3 points for N=128, 256, 512, 1024. The computed spreads on the different grids all follow the general trend of the Schrödinger solution ΔS, but they are all subject to superimposed high-frequency oscillations. Similar oscillations were seen in one-dimensional computations by Valdivieso & Muñoz .
Figure 1b shows the differences between the measured spreads and the spread ΔS predicted by the Schrödinger equation. Also shown is the spread ΔD for a highly accurate numerical solution to the Dirac equation, computed as described in §4a. High-frequency oscillations are also present in ΔD. These oscillations, at frequencies comparable to the Compton frequency, are an intrinsic feature of the Dirac equation as a relativistic theory that accounts for a particle’s rest energy mc2 as well as kinetic energy owing to its motion. Figure 2a shows the second-order convergence of the measured spreads towards the spread ΔD of the reference solution. Further confirmation is given in figure 2b, showing the first-order convergence of the four fields towards the reference solution in the L2 norm.
(a) Reference solutions
In the absence of a potential (g=0), the Dirac equation becomes a linear partial differential equation with constant coefficients. It is thus amenable to Fourier transform techniques. Substituting into equation (2.1) leads to 4.4with the 4×4 matrix Hamiltonian Hk=mβ+α⋅k. The solution to equation (4.4) is given by , where 4.5and Ω=(m2+|k|2)1/2. Exponentially accurate numerical solutions in periodic domains may be computed by expressing the initial ψ as a Fourier series, and evolving each Fourier coefficient using equation (4.5). Transformations between grid-point values and Fourier coefficients may be computed efficiently using fast Fourier transforms.
5. Particle in a harmonic potential
We consider the same initial conditions, with zero initial momentum, confined by a harmonic potential under our sign convention. Setting leads to the initial spread Δ0 of the wavepacket being preserved by subsequent evolution under the Schrödinger equation. For this test case, we consider a cube with side length ℓ=200, mass m=0.1 and initial spread Δ0=14. To avoid the convergence properties being affected by discontinuities in the gradient of V at periodic boundaries, we approximate |x|2 by v(x)+v(y)+v(z), where the function v(x) is the first six terms of a Fourier cosine series approximation to x2 in [−ℓ,ℓ].
(a) Reference solutions
A spatially varying potential couples the different Fourier modes, so exact solutions cannot be computed as Fourier series. However, accurate approximations may be obtained by splitting the Hamiltonian. The free-particle Hamiltonian Hk decouples into one 4×4 block for each Fourier mode, as above. The coupling Hamiltonian Hg=−g(x)I is diagonal in physical space, with a time-evolution operator given by 5.1Thus, g is minus the dimensionless potential energy in our sign convention. Strang splitting approximates the coupled evolution to second-order accuracy by using 5.2For the reference computations, we used Forest & Ruth’s  fourth-order accurate splitting into seven steps, with an overall time step Δt=1.0 and 1283 Fourier modes. This solution was indistinguishable, for comparison purposes, from one computed with Δt=0.25 and 2563 Fourier modes.
Figure 3a shows the evolution of the spread Δ for the reference solution, and for four QLB solutions on N3 grids with N=128, 256, 512, 1024. All solutions are in good agreement for , after which the initially regular oscillations are disrupted by negative-energy states reflecting from the periodic boundaries. The solutions appear very sensitive to small errors owing to the finite resolution, as shown by the diverging trajectories. Figure 3b shows similar behaviour in the squared L2 norms of the negative-energy components , as normalized by and 2M, respectively. These two components are initially related by ||ϕ−2||2=2||ϕ−1||2, as expected for solutions in which ϕ+1 is spherically symmetric and ||ϕ+2||≪||ϕ+1||. Again, the numerical solutions diverge away from this relation when t≈100. However, at t=100, we still have the first-order convergence of all four components of the QLB solutions towards the reference solution, as shown in figure 4a. The rate of convergence deteriorates at later times, owing to the apparent sensitivity of the solution for t≈100. The measured convergence rate at t=400 is only N−1/4 over the range of resolutions used, as shown in figure 4b.
The QLB schemes are viable numerical schemes for solving the Dirac equation, not just for recovering its non-relativistic limit. Each step within the schemes is unitary, and at least second-order accurate. However, the overall schemes are only first-order accurate owing to operator splitting. We presented numerical evidence of first-order convergence towards highly accurate reference solutions for all four components of the Dirac wave function. The oscillations seen in the numerical output of QLB schemes that are not present in the Schrödinger equation are explained as relativistic effects owing to the finite mass of the particle.
We thank Silvia Palpacelli and Sauro Succi for useful conversations. P.J.D.’s research is supported by an EPSRC Advanced Research Fellowship, grant number EP/E054625/1. D.L. is supported by an attached project studentship. The computations made use of the FFTW3 library  and the Oxford Supercomputing Centre facilities.
One contribution of 26 to a Theme Issue ‘Discrete simulation of fluid dynamics: methods’.
- This journal is © 2011 The Royal Society