Quantum Lattice Boltzmann is a quantum walk

Numerical methods for the 1-D Dirac equation based on operator splitting and on the quantum lattice Boltzmann (QLB) schemes are reviewed. It is shown that these discretizations fall within the class of quantum walks, i.e. discrete maps for complex fields, whose continuum limit delivers Dirac-like relativistic quantum wave equations. The correspondence between the quantum walk dynamics and these numerical schemes is given explicitly, allowing a connection between quantum computations, numerical analysis and lattice Boltzmann methods. The QLB method is then extended to the Dirac equation in curved spaces and it is demonstrated that the quantum walk structure is preserved. Finally, it is argued that the existence of this link between the discretized Dirac equation and quantum walks may be employed to simulate relativistic quantum dynamics on quantum computers.


Introduction
A quantum walk (QW) is defined as the quantum analogue of a classical random walk, where the "quantum walker" is in a superposition of states instead of being described by a probability distribution. One of the earliest realization of this concept was proposed by Feynman as a discrete version of the massive Dirac equation [1]. In recent times, there has been a surge of interest for this topic, due to the conceptual and possible practical import of QW's as discrete realizations of stochastic quantum processes and because they can solve certain problems with an exponential speedup, i.e. using exponentially less operations than classical computations [2]. Moreover, QW's are amenable to a number of experimental realizations, such as ion traps [3,4,5], liquid-state nuclear-magnetic-resonance quantum-information processor [6], photonic devices [7] and other types of optical devices [8]. As a result, they hold promise of playing an important role in many areas of modern physics and quantum technology, such as quantum computing, foundational quantum mechanics and biophysics [9].
One of the most interesting features of discrete QW's is their continuum limit, which recovers a broad variety of relativistic quantum wave equations [10,11,12]. As stated earlier, this was first discussed by Feynman and is now known as Feynman's checkerboard [13]. This was originally formulated for the free Dirac equation but extensions of these ideas, which include the coupling to external fields, have been investigated [10,11,12]. Pioneering work have been performed in which the Dirac equation is related to cellular automata [14,15]. Lately, the link between QW's and the Dirac equation have been discussed extensively [10,11,12,16]. In these studies, the starting point is a general QW formulation from which the continuum limit is evaluated and then related to the Dirac equation. arXiv:1504.03158v1 [quant-ph] 13 Apr 2015 In this article, QW's and their relations to some known numerical schemes of the Dirac equation are reviewed from a slightly different perspective: it will be demonstrated that the most general QW's are obtained from lattice discretizations of the relativistic quantum wave equation for spin-1/2 particles. More precisely, starting from the continuum Dirac equation, it is shown that QW's can be placed in one-to-one correspondence with numerical schemes based on operator splitting and the QLB scheme. These numerical methods have been developed and employed as efficient numerical tools to solve relativistic quantum mechanics problems on classical computers [17,18,19]. They have a number of interesting properties: they are easy to code, they can be easily parallellized and are very versatile. Moreover, their mathematical structure and the fact that the time discretization is realized by a set of unitary transformations makes the link to QW's possible. This connection is explored below and represents one of the main purposes of this article. Much of the material discussed here, in particular the numerical simulations, is not completely new, which is in line with the main purpose of this paper, namely an attempt to bridge the techniques utilized in numerical analysis and (quantum) Lattice Boltzmann theory to the language of quantum computing.
This article is organized as follows. In Section 2, a general formulation of QW is presented, where the transfer matrix is time and space dependent. In Section 3, the split operator method for the Dirac equation is presented, along with its exact correspondence with QW. Section 4 is devoted to the QLB method and connections with QW. Section 5 is devoted to a qualitative discussion of the link between these numerical schemes and quantum computation. In Section 6, the schemes are casted in the form of a propagation-relaxation process and the notion of quantum equilibrium is introduced. Based on the analogy between QW and QLB, a new QLB scheme for the (1+1) Dirac equation in curved space is proposed in Section 7. Finally, the generalization of these methods to many dimensions is briefly discussed and numerical results are presented in Section 8

Quantum walks
Let us consider a (1 + 1) quantum walk on the line for a pair of complex wave functions (bi-spinor ψ), obeying a discrete space-time evolution equations described by the following discrete map [10,11,12,20]: Here, the indices j, m ∈ N ⊗ Z label points on a discretization of space and time, respectively. The object B j,n is a two-by-two matrix with components [12] B j,n := e −iξj,n e iαj,n cos θ j,n e iβj,n sin θ j,n −e −iβj,n sin θ j,n e −iαj,n cos θ j,n .
This matrix is a U (2) operator (B ∈ U (2)) parametrized by the three space-time dependent Euler angles θ j,n , α j,n , β j,n and a space-time dependent phase ξ j,n . The latter is relevant when it depends on time and space, i.e. when it is local. If it is global (ξ j,n = ξ), it disappears from any observables and becomes unimportant. This occurs because when the phase is space-time dependent, it does not commute with the time and space translation operators. The matrix obeys B ∈ SU (2) only when ξ j,n = kπ for all j, n, with k ∈ Z. Thus, this formulation is slightly more general than QW's considered in [10,11,20] where B ∈ SU (2) is studied. As seen in the next section, the choice B ∈ U (2) will be important to have a general connection between mass terms of the Dirac equation and QW's. Finally, the U (2) QW can also be implemented on quantum computers because the matrix B is a unitary transformation: it represents the most general QW consistent with quantum computations.
In the above, the amplitudes ψ 1,2 code for the probability of the quantum walker to move up (down) along the lattice site j ∈ Z at the time step n ∈ N. This is a very rich structure, which has been shown to recover a variety of important quantum wave equations, as soon as the Euler angles are allowed to acquire a spacetime dependence [12]. In addition, it provides a wealth of potential algorithms for quantum computing. This was studied extensively in [10,11,12] by analysing the continuum limit of these QW, yielding different versions of the Dirac equation. In this work, the opposite path is taken: it is shown that specific discretizations of the Dirac equation, using either a split-operator approach or the lattice Boltzmann technique, naturally lead to a QW formulation.

Split-operator and quantum walks
The starting point of this discussion is the 1-D Dirac equation in Majorana representation written as (in units where c = = 1): with the bi-spinor ψ ∈ L 2 (R, C 2 ). The generalized "mass" matrix M is space and time dependent and may include contributions from the physical mass, the coupling to an electromagnetic potential or any other type of coupling. One requirement, however, is that M is a Hermitian local operator without any derivatives. Generally, it can be written as where I 2 is the two-by-two identity matrix, (σ i ) i=x,y,z are Pauli matrices and the coefficients M 0,x,y,z represent the time-and space-dependent external fields which couple to the spinor. The formal solution of Eq. (3) is given by whereT is the time-ordering operator, t 0 is the initial time and ∆t = t − t 0 . Using an operator splitting technique, the solution in Eq. (6) can be approximated by [19] The first exponential is a translation (streaming) operator which shifts the spinor components according to: This suggests to use a spatial discretization where ∆z = ∆t (this corresponds to a Courant-Friedrichs-Lewy (CFL) condition C = c∆t/∆z = 1, c being the speed of light) such that the translation is exact on the lattice. Eq. (7) can then be written as: where we defined the following quantities on the lattice: M j,n := M (n∆t, j∆z) and ψ n j := ψ(n∆t, j∆z). This last equation yields a numerical scheme to solve the Dirac equation. This numerical scheme has interesting properties, as discussed extensively in [19,17,21]: it can be extended to higher order accuracy, it can be easily parallellized and it can be easily coded on a computer. In the following, it is also demonstrated that it is completely equivalent to the U (2) QW described in the last section.
Eq. (9) is in the form of Eq. (1). Moreover, the exponential is also a unitary matrix: and thus, there clearly exists a connection between B and B . They are expressed in different representation: B uses the Euler angle parametrization while B is expressed in the canonical representation obtained by the exponential mapping of the Lie algebra. The latter is given explicitly by where |M j,n | = M 2 x,j,n + M 2 y,j,n + M 2 z,j,n .
It is well-known that parametrizations of U (2) matrices are related to each other [22] and it can be determined that the mapping between both representations is given by tan(θ j,n ) = tan(|M j,n |∆t) .
When M 0 = 0, one recovers the SU (2) QW because then,B := B | M0=0 ∈ SU (2). These last equations give a one-to-one correspondence between the QW formulation, characterized by the parameters α j,n , β j,n , θ j,n , ξ j,n , and the discretization of the Dirac equation with a generalized mass term M j,n . Therefore, the last results show that in the continuum limit, every space-time dependent QW on the line becomes a time-dependent 1D Dirac equation with a specific mass matrix. A given QW can thus be fully characterized by the relativistic dynamics of an electron coupled to a space-time dependent external field. This occurs because there is an equivalence between the discretization of the Dirac equation, based on operator splitting, and the QW formulation. Moreover, this connection may be the base for the implementation of a quantum algorithm that solves the Dirac equation on quantum computers. The operator splitting technique presented here also bears a close relationship with the QLB technique, to which we now turn.
where a, b are spinor indices and the "microscopic velocities" are given by v 1,2 = ±1. This is clearly in a "Boltzmann-like" form with two discrete velocities and a collision term M . When QLB is employed in fluid mechanics to solve the continuum equations of motion, there is a family of possible numbers of discrete velocities for a given lattice [23] and each choice yields a different numerical scheme (for instance, the 9 velocities scheme in 2-D and the 27 velocities scheme in 3-D are popular choices on square lattices [24]). For the Dirac equation, this choice is dictated by the mathematical structure of the equation.
The lattice Boltzmann equation is then written as [25,26] As for the splitting method described in the last section, this suggests to use a space discretization where ∆z = ∆t. Using a "naive" approach in line with LB methods in fluid mechanics, the last equation would become Then, the matrix −iM on the right acts as collision operator while the streaming is executed by the left part of the equation. This scheme is derived by using the formal analogy between the Dirac equation, the Boltzmann equation and the LB technique. However, the resulting numerical method is unstable and the L 2 norm is not preserved [27]. It is possible, however, to recover a stable and norm-preserving method. Instead of using the "naive" LB equation discretization given in Eq. (20), the mass term on the right-hand side of Eq. (19) is discretized by using an implicit Crank-Nicolson average: Reporting this into Eq. (19), one obtains the second order accurate QLB scheme: where the transfer matrix is given by To link these results with the ones of the last sections, it is convenient to shift the spinor components by one lattice point (we let ψ 1,j+1 → ψ 1,j and ψ 2,j−1 → ψ 2,j ). Then, the QLB scheme becomes which is in the same form as the QW and the splitting method presented previously. The transfer matrix can be evaluated explicitly. It is given by where C j,n = 1 + i∆tM 0,j,n − ∆t 2 4 M 2 j,n and M 2 j,n = M 2 0,j,n − M 2 j,n . It can be readily verified that T j,n T † j,n = I 2 . As a consequence, the transfer matrix is unitary T j,n ∈ U (2), for any size of the time step and thus, conserves the L 2 norm. Moreover, just like the splitting method, there exists a correspondence with QW's because both have U (2) collision matrices. The identification with the general QW in Eq. (1) yields These relations give a correspondence between the QLB technique and the QW.
They are similar to the ones for the splitting method displayed in Eqs. (14) to (17) and actually serve the same purpose: they allow to map a numerical scheme to the QW. The differences are obviously due to the fact that QLB is based on a LB formulation combined with a Crank-Nicolson average to insure stability while the splitting scheme separates the Dirac equation into different operators which can be integrated exactly. However, both methods share the same general structure where a streaming step is followed by a collision step.

An explicit example: the free case
The Dirac equation for the massive and free case (no external field coupling) is where M y = m is a physical fermion mass. This representation of the Dirac equation yields real spinor components, as readily seen by writing the equation componentwise: The QLB scheme, for the massive and free case, reads as follows [18]: where the transfer matrix is can be obtained from Eq. (25) by setting M 0 = M x = M z = 0 and M y = m: It is readily seen that a(m) and b(m) are second-order Pade'-like approximants of cos(m) and sin(m), respectively. This is the natural consequence of the implicit time-marching (Crank-Nicolson) scheme as applied to the (1+1) Dirac equation in Majorana form. For a free particle, the correspondence to QW is particularly simple and is given by This mapping of the QLB to the QW is exact for any value of m.

Prospects for quantum simulation
The numerical methods described in Sections 3-4 can be straightforwardly implemented on classical computers. However, the stream-collide structure of these numerical scheme makes them suitable for an efficient implementation on quantum computers as well. In particular, they can be written as: where S j is the shift operator, defined as: This shift operator is a unitary operation and it can be realized experimentally by using fundamental quantum gates [6,20]. The rotation operator B j,n belongs to U (2) and therefore can also be realized by these quantum gates, as any other unitary transformations [28]. Therefore, it is possible to map the numerical method to quantum walk, which can be implemented efficiently on quantum computers. This mapping is possible because each step of the scheme is a unitary transformation: this makes these schemes norm-preserving and sets the link with quantum computations.
The latter would be particularly useful for the study of relativistic quantum systems where a time-dependent solution of the Dirac equation is required, such as in very high intensity laser physics [29] or graphene physics [30].
Another subject of major interest for future research is the extension of the QLB methodology to quantum many body systems and quantum field theory, two paramount sectors of modern physics which are particularly exposed to the limitations of classical (non-quantum) electronic computing. Progress in this direction depends on the ability to replace the quantum wavefunction by the corresponding second-quantized quantum operators, and show that the dynamics of the secondquantized QLB scheme still preserves the appropriate equal-time commutation relations. Preliminary efforts along this line have been developed in [31] in 1 + 1 dimensions. Extensions to strongly non-linear field theories in d > 1 remain to be explored. As to quantum-many body problems, LB-like methods have been recently adapted to electronic structure simulations [32] In this work, a classical LB scheme is employed to solve the Kohn-Sham equations of density functional theory in the form diffusion-reaction equations in imaginary time. Allied QLB schemes could prove very useful to solve the corresponding real-time quantum many-body transport problems within the framework of time-dependent density functional theory.
Finally, we wish to point out the intriguing possibility of realizing both quantum and classical LB schemes on quantum analogue simulators, as recently explored in [33].

Quantum equilibria
In view of quantum computing implementations, it is of interest to cast the Dirac equation in the form of a propagation-relaxation process, where the collision matrix is now interpreted as a scattering process, relaxing the spinor component around a local quantum equilibrium.
For this purpose, it is useful to reconsider the 1-D Dirac equation in the form: Then, it is formally possible to define a local equilibrium as: where U is a unitary matrix that depends on M . This transformation is chosen to recast the Dirac equation in relaxation form: where Ω = iM [I + U ] −1 . The explicit value of U, Ω is not unique but a convenient choice is U = e iM τ . In this vests, the Dirac equation looks formally like a linear Boltzmann equation for two-component models in the single relaxation time approximation [25]. Therefore, it can be interpreted as a propagation-relaxation process in imaginary time, whereby collisions, implemented by the scattering operator Ω, drive oscillations around the equilibrium distribution ψ eq . By defining a post-collision wavefunction as: the Dirac equation takes the most compact form This is particularly suitable to quantum computing implementations in the form of a classical stream-collide dynamics. The collision (relaxation) gate transforms the pre-collisional spinor ψ into the post-collisional state ψ , and the streaming gate moves the post-collisional spinor to its destination location z ± ∆z. Both operations are unitary and can be encoded in logical gates for quantum computing purposes [34,35].
The expression (40) shows that the local equilibria is a linear function of the actual wavefunction ψ, hence itself a function of space and time.
The question is: will the actual wavefunctions ever reach this moving target, i.e, ψ = ψ eq ? Based on its definition, this can only occur once ψ eq lies in the null-space of the scattering matrix M , namely: For the case of a free massive particle, it can be checked explicitly that the only solution is the trivial vacuum ψ eq ≡ 0. This means that the spinorial wavefunction is a superposition of (both slow and fast) zero-average oscillations around a local equilibrium which, consistently with the reversible nature of quantum mechanics, is actually never attained.
Although this remains to be checked in detail, we conjecture that the same holds true for the case of a massive particle in an external potential, because in this case the Dirac equation is still linear.
Based on (44), the condition for the local equilibrium to depart from a trivial vacuum is that the matrix M be singular, i.e. the local equilibrium is a zero-mode of the scattering matrix.
Non-trivial quantum zero-modes, ψ eq = 0, may indeed arise for non-linear quantum wave equations, such as the Gross-Pitaevski or the mean-field version of the Nambu-Jona-Lasinio model, to be dicussed shortly. A non-trivial local equilibrium would then signal a spontaneously broken symmetry, which is indeed the distinctive trait of the aforementioned non-linear quantum wave equations.
Even though the notion of quantum equilibrium remains purely formal in nature, it is argued that it might nonetheless facilitate quantum computing implementations based on the compact expressions (42) and (43). This stands out a very interesting topic for future research.

QLB in curved space-time
Quantum walks have been shown to map into Dirac-like equations in curved space as well by evaluating the continuum limit of certain QW's [36,12]. Here, in the same spirit as other sections, a QW structure is obtained by discretizing the Dirac equation in curved space time using a QLB-like approach. However, because the wave function propagates on a curved manifold, the structure of the resulting QW is different from Eq. (1) and should include a residency matrix that corrects the streaming step, which is strictly valid only in flat space. This is different from the result obtained in [36,12].
The Dirac equation in a static (1+1) curved space writes as: where a = 0, 1 and µ = t, z. In the above e µ a is the two-dimensional vierbein (zweibein) relating the components of the locally tangent Minkowski space basis (e 0 , e 1 ), described by the metric η ab , to the global space basis (e t , e z ), described by the general metric g µν . Also, g := det(g µν ) is the determinant of the metric tensor.
The general form of the corresponding partial differential equation is with A(z) ∈ R a function of space and Q(z, t) the two-by-two gravitational collision matrix associated with Eq. (45). This becomes a hyperbolic system of equations where the advection speed A(z) should not be confused with the vector electrodynamic potential. Since the advection term is heterogeneous, a strict QLB structure, i.e. streaming along constant lightcones ∆z = ∓c∆t, is no longer viable. Here, the situation is very similar to classical LB schemes on non-uniform grids (unsurprisingly, since in dimension D = 1, gravity is basically a stretching of the metric). For this purpose, several finite-volume LB (FVLB) schemes have been formulated, whose main outcome is that the streaming operator is no longer a diagonal matrix with eigenvalues ±1, as required by the formal QLB structure [37,38]. In this respect, a finite-volume QLB scheme for the Dirac equation in curved space can be obtained by writing the Dirac equation as with Q (z, t) := Q(z, t) + ∂ z A(z)σ z . Then, this equation is integrated over a control volume V (enclosed by a surface S), extending from (j − 1 2 )∆z to (j + 1 2 )∆z. Finally, applying a Crank-Nicolson time-marching and combining with upwind finitedifferences, the discretized equation is given by where c := ∆z/∆t is the uniform lattice light speed (or CFL condition). The boundary integrals are given by an thus, require an interpolation from the cell centers j, to the north and south boundaries at j ± 1/2, respectively. Following a common practice in finite-volume formulations of hyperbolic problems, the flux terms are approximated by [39]: As a result: It is readily seen that this reduces to a standard QLB in the limit of a uniform grid, when A j /c = 1. In this case the streaming is diagonal with speed ±c and the spinors at (j, n + 1) are connected to the corresponding spinors at (j ± 1, n) by a local 2 × 2 matrix, which can be readily inverted to deliver a fully explicit map. However, when A j /c = 1, the spinors at (j, n) also enter the map, so that local inversion delivers a slightly more elaborated structure, namely: which can be written more explicitly as In the above T is the local 2 × 2 transfer matrix including collisions, S is the streaming operator and R the local residency matrix, expressing the fraction of spinors which are left in the cell centered about z as the quantum system advances from t to t + ∆t. This is depicted in Fig. 1. Clearly, the residency matrix vanishes in the case of a uniform mesh, i.e. no gravity. The mapping (58) represents the "gravitational" QLB. The detailed expressions of the streaming and residency matrices depend on the specific form of the metric tensor and associated vierbeins. Moreover, this analysis concentrates on the mathematical structure (streaming and collision steps) of the resulting scheme rather than on its numerical properties (convergence, stability, etc). These topics shall make the object of a future publication.

Multi-dimensions
The discretization presented in this work extends to the D + 1 dimensional case by applying the notion of operator splitting. This implies the inclusion of a new dynamic step which is entirely quantum: namely a "rotation", designed so as to keep the spin aligned with the momentum along each of the three spatial directions. Schemes using this strategy can be found in [18,19] and in [17] for higher order splittings.
It might be that such rotation is not needed by formulating the Dirac equation as a random walk on other lattice with more natural topologies (the diamond) lattice. The QLB-QW equivalence in multi-dimensions will be discussed in a future publication. However, to demonstrate the strength of the numerical schemes presented here and to show some possible applications for quantum computing, numerical results in 2-D are presented in the following.

Numerical results
As an example of possible applications of QLB scheme, we present two representative simulations: Klein tunnelling in the presence of random impurities and Dirac equation with Nambu-Jona-Lasinio (NJL) interactions in 2 + 1 space-time dimensions. Details on the numerical methods used to obtain these results are given in [18,19,40]. Also, these results are not completely new as similar systems have been studied in [19,41].

Graphene with random impurities
In the first numerical test, the propagation of a Gaussian wave packet through a graphene sample with randomly distributed impurities is simulated [41]. In Ref. [41], simulations are performed for different values of the impurity concentration and the potential barrier, in order to provide an estimate of the effect of impurity concentration on the conductivity of the graphene sample. In Fig. 2, we report some representative snapshots of the first 1800 time steps of the simulation, at an impurity percentage= 0.5% and V = 50 MeV. A lattice of size 2048 × 512 cells is used and the cell size is chosen to be ∆z = 0.96 nm, while the spreading of the initial Gaussian wave packet is σ = 48 (in numerical units), leading to a Fermi frequency k F = 0.117 (80 MeV in physical units). In this simulation, a fully relativistic particle (m = 0) is considered.
From Fig. 2, we can see that the wave packet is scattered by the impurities, giving rise to a plane front out of the initial Gaussian configuration. As a consequence of the randomness induced in the wave function by the disordered media, there is a momentum loss and therefore the motion of the wave packet is found to experience a corresponding slow down. It is also found that the wave packet takes more time to regroup as the impurity concentration and impurity potential are increased.
A grid size of N z × N y = 1024 2 elements is used and the initial wave packet spread is set at σ = 48, a fully relativistic particle (m = 0) is considered. In these simulations, we impose g = 0 and g = 1000 and vary the initial energy of the wave packet k ≡ k z = k y in order to inspect the effect of this parameter on the wave packet separation, which, in turn, informs on the effective mass acquired by the up and down propagating modes.
In Fig. 3, the wave function density at time t = 200 for k = 0.004, 0.04 and 0.4 is shown for g = 0 and g = 1000, respectively. The figure shows that sufficient energy, k > 0.004, is needed to observe the splitting of the wavepacket. The effects of nonlinear interactions, fringes and distortions, are also well visible in the right column of Figure 3, A quantitative analysis in the one-dimensional case led to satisfactory agreement with asymptotic solutions for the dynamic mass as a function of the interaction strength g. A similar analysis in two spatial dimensions remains to be developed.

Summary and outlook
Summarizing, we have reviewed discretizations of the Dirac equation and described their mapping into QW's. These relations may allow the solution of the Dirac equation on quantum computers. In the first part, a general argument is given, using the operator splitting method. Then, the QLB scheme is studied within the same g = 0 g = 1000 Figure 3 Wave packet density ρ = |ψ 1 | 2 + |ψ 2 | 2 + |ψ 3 | 2 + |ψ 4 | 2 for the scheme solving Dirac equation with NJL interaction. These snapshots are taken at time t = 200 for k = 0.004, 0.04 and 0.4 (from top to bottom). The figure shows how the initial energy affects the separation phenomenon. In particular, at low energy, k = 0.004, no splitting of the wavepackest is observed.
perspective and a similar relation is found. We have also shown that a similar structure remains in curved space, using a scheme based on a finite volume formulation, with the important caveat that the exact nature of the streaming operator, typical of QLB, is no longer preserved. Rather, one sees the appearance of the residency matrix, which characterizes the fraction of spinor which is left in the cell after one step in the time evolution. This scheme, along with its generalization to many dimensions, will be studied in future work.