Formal analogy between the Dirac equation in its Majorana form and the discrete-velocity version of the Boltzmann kinetic equation
F. Fillion-Gourdeau, H.J. Herrmann, M. Mendoza, S. Palpacelli, S. Succi
aa r X i v : . [ qu a n t - ph ] O c t Formal analogy between the Dirac equation in its Majorana form and thediscrete-velocity version of the Boltzmann kinetic equation
F. Fillion-Gourdeau, ∗ H.J. Herrmann, M. Mendoza, S. Palpacelli, † and S. Succi ‡ Centre de Recherches Math´ematiques, Universit´e de Montr´eal, Montr´eal, Canada, H3T 1J4 ETH Zurich, Computational Physics for Engineering Materials,Institute for Building Materials, Schafmattstrasse 6, HIF, CH-8093 Zurich (Switzerland) Numidia srl, via G. Peroni, 130, 00131, Roma, Italy Istituto Applicazioni Calcolo, CNR, via dei Taurini 19, 00185, Roma, Italy § (Dated: October 31, 2018)We point out a formal analogy between the Dirac equation in Majorana form and the discrete-velocity version of the Boltzmann kinetic equation. By a systematic analysis based on the theoryof operator splitting, this analogy is shown to turn into a concrete and efficient computationalmethod, providing a unified treatment of relativistic and non-relativistic quantum mechanics. Thismight have potentially far-reaching implications for both classical and quantum computing, becauseit shows that, by splitting time along the three spatial directions, quantum information (Dirac-Majorana wavefunction) propagates in space-time as a classical statistical process (Boltzmann dis-tribution). PACS numbers: 02.60.Cb,03.65.Pm,03.67.Ac,11.10.Lm
BOLTZMANN AND DIRAC
Analogies between the non-relativistic Schr¨odingerequation and fluid dynamics have been noted since theearly days of quantum mechanics. In particular, back in1927, Erwin Madelung noticed that by expressing thewavefunction in eikonal form, i.e. Ψ =
R e iS/ ~ , theSchr¨odinger equation turns into the hydrodynamic equa-tion of a compressible, inviscid fluid, with number den-sity ρ = R and velocity ~u = −∇ S/m . The quantumfluid is subject to the classical potential V c ( ~x ), plus thequantum potential V q ( ~x ) = − ~ m (∆ R ) /R . Although thehydrodynamic analogy is commonly regarded as purelyformal in nature, lately, its connections with Bohm’s the-ory of hidden variables and De Broglie’s pilot wave pic-ture have known of surge of interest, mostly in connectionwith experimental investigations on the non-local natureof quantum physics [1].The quantum relativistic fluid analogy seems to havereceived comparatively less attention. Back in 1993, itwas noted that the Dirac equation can be regarded asa special form of a discrete Boltzmann kinetic equation,in which the particle velocities are confined to a handfulof discrete values [2]. The discrete components of theBoltzmann distribution, f i ( ~x ; t ) ≡ f ( ~x, ~v = ~v i ; t ), wherethe index i labels the discrete velocities, are then identi-fied with the spinor components ψ i of the Dirac equation.This opens up an interesting connection between classicalkinetic theory and relativistic quantum mechanics.Mathematically, the connection is not so surprising,since both Boltzmann and Dirac equations are hyperbolicsupersets of the Navier-Stokes and Schr¨odinger equa-tions, respectively.The interesting point, however, is that the connectionbecomes much more direct and compelling by considering the discrete-velocity version of the Boltzmann equation,in relation to the Majorana form of the Dirac equation,in which all matrices are real [3].Majorana particles have attracted significant interestin recent years, mostly in connection with the fact thatthey coincide with their own antiparticles, as beautifullydiscussed in a recent essay by F. Wilczek [4].Here, we wish to put forward a different angle of in-terest of the Majorana representation, namely the factthat not only it makes Boltzmann-Dirac analogy concep-tually more poignant, but it also turns it into a concrete unified computational scheme for the simulation of bothrelativistic and non-relativistic quantum wave equations,on both classical and quantum computers. The corre-sponding method is known as quantum lattice Boltzmann(QLB) method [2]. The QLB is based on the identifi-cation of the discrete Boltzmann distribution with thespinorial wavefunction f i ( ~x ; t ) ↔ ψ i ( ~x ; t ). Even thoughboth objects are real, they still face a mismatch of de-grees of freedom in more than one spatial dimensions,since a spinor of order s consists of 2 s +1 components, re-gardless of the number of dimensions, while the discretedistribution requires (at least) 2 d discrete componentsin d spatial dimensions. Moreover, the Dirac-Majoranamatrices cannot be simultaneously diagonalized, reflect-ing the basic fact that spinors are not ordinary vectors.As a result, in more than one spatial dimensions, it is inprinciple not possible to keep the particle velocity alignedwith its spin.Remarkably, both problems can be circumvented byresorting to operator splitting. Essentially, this amountsto splitting the spinor propagation along the three spatialdimensions into a series of three one-dimensional prop-agations, each using the diagonalized form of the corre-sponding Dirac-Majorana streaming matrix. As a result,t each propagation step the particle spin is kept alignedwith its velocity, so that the identification f i ↔ ψ i con-tinues to hold.In this Letter we show that this “heuristic stratagem”is backed up by a rigorous mathematical treatment,which leads to a unified computational approach toquantum wave mechanics. The resulting computationalscheme offers outstanding amenability to parallel com-puting on electronic computers [5] and is also suitable toprospective quantum computing simulations [6–8].To show its versatility also towards the inclusion ofnon-linear interactions, as an application, we shall solvea specific form of the non-linear Dirac equation includingdynamical-symmetry breaking term, as first proposed byNambu and Jona-Lasinio. Discrete Boltzmann and Dirac
To set up the framework, let us write down the twoequations in full display. The discrete Boltzmann equa-tion reads as follows: ∂ t f + v ai ∇ a f i = Ω ij ( f j − f ej ) (1)where f i = f ( ~x, ~v = ~v i ; t ) is the probability density offinding a particle around position ~x at time t with discretevelocity ~v i . The latin index a = x, y, z runs over spatialdimensions and Einstein summation rule is assumed. Theleft hand side represents the particle free-streaming (inthe absence of external forces, for simplicity), while theright-hand side is the collisional step steering the distri-bution function towards a local Maxwell equilibrium f ei .The (symmetric) scattering matrix Ω ij encodes the mass-momentum-energy conservation laws underpinning fluiddynamic behavior.The Dirac equation, in Majorana form, reads as follows ∂ t ψ i + S aij ∇ a ψ j = M ij ψ j (2)where S aij are the three Majorana streaming matrices and M ij is the (anti-symmetric) mass matrix, acting upon thereal spinor ψ i , i = 1 , s + 1. This clearly shows a formalanalogy with the Boltzmann equation: the lhs describesthe free streaming of the spinors, while the rhs can beregarded as a simple form of local collision between thevarious spinorial components. Note that the mass matrixhas dimensions of an inverse time scale, typically givenby the Compton frequency ω c = mc / ~ . In 1D, thisanalogy is “exact”: by choosing a representation wherethe Dirac matrix is diagonal (Majorana representation),we recover Eq. (1). In multiple dimensions however, thestory is different: the connection can be realized only byresorting to operator splitting, whereby each step can bewritten in the form of Eq. (1). This will be discussed inthe following. Quantum lattice Boltzmann
Let us consider the case of spin s = 1 / β, α a in the Dirac representation. The goal here is to find the discrete time evolution of the wave function byusing the formal analogy with the Boltzmann equation.In the QLB setting, this time evolution proceeds by asequence of streaming and collisional steps, given by (weuse natural units where c = ~ = 1): ∂ t ψ ( x ) ( t ) = − α x ∂ x ψ ( x ) ( t ) , ψ ( x ) ( t n ) = ψ ( t n ) , (3) ∂ t ψ ( y ) ( t ) = − α y ∂ y ψ ( y ) ( t ) , ψ ( y ) ( t n ) = ψ ( x ) ( t n +1 ) , (4) ∂ t ψ ( z ) ( t ) = − α z ∂ z ψ ( z ) ( t ) , ψ ( z ) ( t n ) = ψ ( y ) ( t n +1 ) , (5) ∂ t ψ ( c ) ( t ) = − iβmψ ( c ) ( t ) , ψ ( c ) ( t n ) = ψ ( z ) ( t n +1 ) , (6)and ψ ( t n +1 ) = ψ ( c ) ( t n +1 ) , (7)where the superscript labels the step of the splitting and t n = n ∆ t is the time after n iterations. In these equa-tions, the calculated solution at a given step provides aninitial condition for the next step in the sequence. Eqs.(3) to (5) correspond to streaming while the last step inEq. (6) is collisional.The streaming steps for a given coordinate a proceed asfollows. First, it should be noted that the matrix α a (for a = x, y, z ) is not diagonal and thus, the Dirac equationis not in the form of Eq. (1). However, the latter can berecovered by using the unitary transformation of spinors S a = √ ( β + α a ) . This equation allows to transform theDirac matrices to a Majorana-like representation, wherethe matrix ˜ α a = S † a α a S a = β is diagonal, with eigenval-ues ±
1. Then, by introducing the transformed spinor as˜ ψ ( a ) = S − a ψ ( a ) , the streaming steps can be turned into ∂ t ˜ ψ ( a ) ( t ) = − β∂ a ˜ ψ ( a ) ( t ) , (8)which is clearly in the form of Eq. (1) without collisionalterm. This has a solution given by˜ ψ ( a )1 , ( t n +1 , x ) = ˜ ψ ( a )1 , ( t n , x a − ∆ t ) , (9)˜ ψ ( a )3 , ( t n +1 , x ) = ˜ ψ ( a )3 , ( t n , x a + ∆ t ) . (10)where x a + v ai = x a ∓ ∆ t , i = − ,
1, is the lattice neighborpointed by the discrete speed v ai = ∓ c . This correspondsto an exact integration of the streaming operator alongthe characteristics ∆ x a = ± c ∆ t (light-cones), which istypical of the Lattice Boltzmann (LB) method.The collision step can also be integrated exactly byusing the solution ψ ( c ) ( t n +1 ) = e − iβm ∆ t ψ ( c ) ( t n ) ≡ Cψ ( c ) ( t n ) . (11)It is then possible to write C = e − M ∆ t explicitly as a4 × t . Moreover, itlooks like a classical motion of two discrete walkers, hop-ping by one lattice unit along every coordinate at eachtime-step and colliding according to the scattering ma-trix M = iβm . More complex interactions can be treatedin a similar way by including the interaction terms intothe scattering matrix. As long as the matrix is local, it2s not necessary to diagonalize S and M simultaneouslyand due to the operator splitting, the simplicity of theLB formalism is not compromised. For instance, for thecoupling to an electromagnetic field, the scattering ma-trix is given by M = iβm − ieα a A a ( x, t )+ ieV ( x, t ) where( A a , V ) is the electromagnetic potential.Symbolically, the 3D evolution of the Dirac spinorreads like a sequence of three one-dimensional streamsteps and one collisional step: ψ ( t n +1 , x ) = C ( S z P z S − z )( S y P y S − y )( S x P x S − x ) ψ ( t n , x )(12)where P a = e − ∆ tβ∂ a is a translation operator along thedirection a . The latter shifts the “1,2” and “3,4” spinorcomponents by ∓ ∆ t , respectively.Of course, this procedure is not exact: as shown in thefollowing, it corresponds to an operator splitting methodwhere the streaming and collision matrices do not com-mute. However, each step of the splitting -is- exact andthus, the only source of error comes from the splittingwhich scales like O (∆ t ) (second order accuracy). Werefer the reader to [10] for the numerical analysis ofthe scheme. Other schemes where the error scales like O (∆ t ) can also be obtained [5, 10]. Most importantly,it does not spoil the unitarity of the scheme for any valueof the timestep: this is required to conserve the proba-bility density ( L norm). Full details of the algorithmcan be found in [5] and slightly different versions are in[2, 11]. The general operator-splitting framework
The QLB was derived on heuristic grounds, based ona intuitive analogy between a genuinely quantum vari-able, the particle spin, and a discrete one, the particlemomentum in the lattice formulation of the Boltzmannequation. Since quantization is a physical concept whilediscretization is a numerical one, it might be argued thatthe analogy is somewhat artificial, hence perhaps coinci-dental and of limited applicability.In the sequel, we shall show that this is not the case:QLB can be shown to fall within the general theory ofoperator splitting, as applied to the Dirac equation.This might have potentially deep implications for bothclassical and quantum computing, because it impliesthat, by splitting time along the three spatial directions,and augmenting the stream-collide dynamics with properglobal rotations, quantum information (the Dirac wave-function) propagates in space-time as a classical statis-tical process (Boltzmann distribution). It would be ofgreat interest to explore whether such insight could beused to simulate the Dirac equation on trapped-ion ana-logue computers based on the QLB dynamics [12].The starting point of the general operator splitting the- ory is the formal solution of the Dirac equation given by ψ ( t n +1 ) = T exp (cid:20) − i Z t n +1 t n H ( t ) dt (cid:21) ψ ( t n ) , (13)= e − i ∆ t ( H ( t n )+ T ) ψ ( t n ) (14)where H ( t ) is the Dirac Hamiltonian, T is the time-ordering operator and T = i ←− ∂ t n is the “left” time-shiftingoperator. The second form of the solution was obtainedin [13] and constitutes a great starting point for deriv-ing approximation schemes. Then, the operator split-ting method consists in decomposing the Hamiltonian as H ( t ) = P Nj =1 H j ( t ) and to approximate the evolutionoperator in Eq. (14) by a sequence of exponentials in theform: ψ ( t n +1 ) ≈ N seq Y k =1 e − is ( k )0 ∆ t T N Y j =1 e − is ( k ) j ∆ tH j ( t n ) ψ ( t n )(15)where the coefficients N seq ∈ N and s ( k ) j ∈ R are cho-sen to obtain an approximation with a given order ofaccuracy. It is then straightforward to conclude thatthe QLB scheme, shown in Eq. (12) and in Eqs. (3)to (6), corresponds to a particular decomposition of theHamiltonian[14] and to a specific realization of Eq. (15).The conclusion is far reaching; the Majorana represen-tation exposes a concrete connection between the (dis-crete) Boltzmann equation and the Dirac equation inMajorana form. As a result, the information contained inthe quantum relativistic four-spinor ψ ( t, x ) can be pro-cessed on entirely classical terms, i.e free-streaming alongconstant directions and local collisions, complementedwith diagonalization steps to keep speed and spin con-stantly aligned. Remarkably, the scheme is also viablefor prospective quantum computer implementations [6–8, 15].The QLB has been applied to a variety of quantumwave problems, mostly in the non-relativistic context,[16–18]. Here we present a new application to an im-portant non-linear relativistic problem, namely the Diracequation augmented with Nambu-Jona-Lasinio dynamicsymmetry breaking terms. The NJL-Dirac equation
The Nambu-Jona-Lasinio (NJL) model was promptedout by a profound analogy between the Bardeen-Cooper-Schrieffer theory of superconductivity and chiral symme-try breaking in relativistic quantum field theories [19, 20]and it has served ever since as a model paradigm to studysymmetry-breaking phenomena in both fields.The NJL Lagrangian reads [19] L NJL = ¯ ψ ( iγ µ ∂ µ − m ) ψ + g ψψ ) − ( ¯ ψγ ψ ) ] . (16)This corresponds to the free-particle Dirac Lagrangian,plus an interaction term, driven by the coupling param-3ter g . This coupling term reflects four-fermion interac-tions, in direct analogy with the BCS theory of supercon-ductivity. By imposing the chiral symmetry, the NJL la-grangian should not present any explicit bare mass term,so we set m = 0. However, the NJL dynamics leads to theformation of a chiral condensate, corresponding to an ef-fective mass term and a spontaneous symmetry breakingof the chiral symmetry. Much of the current interest inthe NJL model is motivated by the fact that it serves asa phenomenological model of quantum chromodynamics(for a full account see [21]).The associated equation of motion reads (see Ap-pendix)( ∂ t + α a ∂ a + imβ ) ψ = igβ [( ψ † βψ ) − ( ψ † βγ ψ ) γ ] ψ. (17)A solution of this equation is required for the quantumstudy of this model, in the mean-field approximation.The space-time discretization of NJL-Dirac can be castin the standard QLB format by adding the non-linearterm into the collision step as in the case of the electro-magnetic field, by replacing C → C NJL . The collisionstep becomes ψ ( c ) ( t n +1 ) = C NJL ψ ( c ) ( t n ) , = T exp (cid:20)Z t n +1 t n dtM NJL ( t ) (cid:21) ψ ( c ) ( t n ) , (18) M NJL ( t ) ≡ − iβ ( m − gρ S ( t )) − gρ A ( t )Σ , (19)where ρ S ≡ ψ † βψ and ρ A ≡ iψ † βγ ψ depend on time,hence the time-ordering operator, and Σ ≡ βγ . Thetime-ordering can be approximated by using Eqs. (14)and (15): the ensuing ordinary exponential can be con-verted exactly to a 4 × C NJL . A similartreatment of the nonlinear term, albeit using spectralmethods, can be found in [22].
Numerical application
As an application of the QLB scheme, we simulate theemergence of a dynamic fermion mass as a result of thespontaneous breaking of the chiral symmetry of the NJLequation.For this purpose, let us consider an initial conditiongiven by the following Gaussian minimum-uncertaintywave packet ψ ( t = 0 , z ) = S y − C u e ikz + C d e − ikz C u e ikz − C d e − ikz C u e ikz + C d e − ikz C u e ikz + C d e − ikz e − z σ (2 πσ ) (20)centered about z = 0, with initial width σ . Let ω = k be the initial energy of the wave packet. The coefficients C u and C d obey the condition 2 C u + 2 C d = 1, so that ψ † ψ = | G | . Moreover, an asymmetry can be set bytuning the ratio C u /C d ≡ α = 1.We analyze our numerical results for the case of m = 0, which ensures that the axial current is conserved by thefree part of the equation, as a function of the coupling co-efficient g . For this test, the following parameter settingis used: k = 0 . σ = 48, C u = 1 .
177 and C d = 0 . ρ ( z ) = | ψ | at times t = 10, 50, 100and 200, for the case g = 0, 1 and 2 are shown in Fig. 1.This calculation requires 200 time-steps (for a meshsize of 1024 lattice sites) and about 0.01 CPU seconds ona standard PC. This amounts to a processing speed ofabout 20 MLUPS (Million lattice updates per second),which is in line with the performance of Lattice Boltz-mann schemes for classical fluids. Since the latter isknown to be very competitive, the same conclusion islikely to hold for the quantum case. A final statement inthis direction must be left to detailed head-on compari-son between QLB and state-of-the art numerical methodsfor the Dirac and Schr¨odinger equations. −400 −300 −200 −100 0 100 200 300 4000123456789x 10 −3 z ρ g=0, t=10g=1, t=10g=2, t=10 (a) −400 −300 −200 −100 0 100 200 300 40001234567x 10 −3 z ρ g=0, t=50g=1, t=50g=2, t=50 (b) −400 −300 −200 −100 0 100 200 300 4000123456x 10 −3 z ρ g=0, t=100g=1, t=100g=2, t=100 (c) −400 −300 −200 −100 0 100 200 300 4000123456x 10 −3 z ρ g=0, t=200g=1, t=200g=2, t=200 (d) FIG. 1: ρ = | ψ | at times (a) t = 10, (b) t = 50, (c) t = 100 and (d) t = 200 for g = 0, 1 and 2. The figureshows the separation of the left and right movingwavepackets in the course of the evolution. Thenon-interacting case shows no deformation of theGaussian profile, as expected, while the interacting caseleads to a slow-down and deformation of bothwavepackets.From Fig. 1, a symmetry breaking between the left andright moving wavepackets is clearly seen at increasing val-ues of g . The generation of a dynamic mass is expectedto reflect into a slowing-down of the group velocity of thewavepackets, according to v group c = k √ k + m ′ < m ′ = − gρ S is the dynamic mass of zero-rest mass parti-cles. Indeed, since the initial condition is symmetric withrespect to 1 ↔ ρ A is initiallyzero, and remains such all along the simulation.The results can be checked against the analytic solution4o Eq. (17) in 1-D and for the case of small g [23], whichgives v mean /c ≃ − . g + O ( g ) at early times. It canbe checked (not shown for space limitations) that this isconsistent with the numerical results in Fig. 1.The same phenomenon can be simulated in two dimen-sions, and the details shall be presented in a future andlenghtier publication.Extending the above work to the case of quantummany-body systems and non-linear multidimensionalquantum field theory [24], represents an outstandingchallenge for future research in the field. Appendix: NJL-Dirac equation using Pauli representa-tion
From the NJL Lagrangian of Eq. (16), the associatedequation of motion Eq. (17) is derived as follows. Varia-tion of Eq. (16) against ¯ ψ delivers( iγ µ ∂ µ − m ) ψ + g (cid:2) ( ¯ ψψ ) ψ + ( ¯ ψγ ψ ) γ ψ (cid:3) = 0 , (21)where γ ≡ iγ γ γ γ and ¯ ψ = ψ † γ .The actual definition of the gamma matrices depends onthe specific chosen representation. By using Pauli-Diracrepresentation, γ i matrices are defined as follows [25]: γ = β, γ i = βα i , with i = 1 , . . . , , (22)where β and α i are the standard Dirac matrices.Inserting these definitions into Eq. (21), yields( β∂ t + βα a ∂ a + im ) ψ = ig [( ψ † βψ ) − ( ψ † βγ ψ ) γ ] ψ. (23) ∗ fi[email protected] † [email protected] ‡ [email protected] § Also at Physics Department, Harvard University, Cam-bridge MA, 02138, USA [1] S. Haroche, J.-M. Raimond, and P. Meystre, Physics To-day , 080000 (2007).[2] S. Succi and R. Benzi, Physica D: Nonlinear Phenomena , 327 (1993).[3] C. Itzykson and J. B. Zuber, Quantum Field Theory (Mcgraw-hill, 1980).[4] F. Wilczek, Nature Physics , 614 (2009).[5] F. Fillion-Gourdeau, E. Lorin, and A. D. Bandrauk,Computer Physics Communications , 1403 (2012).[6] R. P. Feynman, International journal of theoreticalphysics , 467 (1982).[7] S. Lloyd, Science , 1073 (1996).[8] B. M. Boghosian and W. Taylor, Physica D: NonlinearPhenomena , 30 (1998).[9] Note1, it is given by C = cos( m ∆ t ) − iβ sin( m ∆ t ).[10] E. Lorin and A. Bandrauk, Nonlinear Analysis: RealWorld Applications , 190 (2011).[11] D. Lapitski and P. J. Dellar, Philosophical Transactionsof the Royal Society A: Mathematical, Physical and En-gineering Sciences , 2155 (2011).[12] R. Gerritsma, G. Kirchmair, F. Z¨ahringer, E. Solano,R. Blatt, and C. Roos, Nature , 68 (2010).[13] M. Suzuki, Proceedings of the Japan Academy. Ser. B:Physical and Biological Sciences , 161 (1993).[14] Note2, the decomposition is such that H = − iα x ∂ x , H = − iα y ∂ y , H = − iα z ∂ z and H = βm .[15] J. Yepez, Quantum Information Processing , 471 (2005).[16] S. Succi, Computer Physics Communications , 317(2002).[17] S. Palpacelli and S. Succi, Phys. Rev. E , 066704(2007).[18] S. Palpacelli and S. Succi, Phys. Rev. E , 066708(2008).[19] Y. Nambu and G. Jona-Lasinio, Phys. Rev. , 345(1961).[20] Y. Nambu and G. Jona-Lasinio, Phys. Rev. , 246(1961).[21] S. P. Klevansky, Rev. Mod. Phys. , 649 (1992).[22] J. Xu, S. Shao, and H. Tang, Journal of ComputationalPhysics , 131 (2013).[23] S. Palpacelli, P. Romatsckhe, and S. Succi, InternationalJournal of Modern Physics C (2013), accepted for publi-cation.[24] S. Succi, Journal of Physics A: Mathematical and Theo-retical , F559 (2007).[25] P. A. M. Dirac, The principles of quantum mechanics (1947).(1947).