A relativistic discrete spacetime formulation of 3+1 QED

This work provides a relativistic, digital quantum simulation scheme for both $2+1$ and $3+1$ dimensional quantum electrodynamics (QED), based on a discrete spacetime formulation of theory. It takes the form of a quantum circuit, infinitely repeating across space and time, parametrised by the discretization step $\Delta_t=\Delta_x$. Strict causality at each step is ensured as circuit wires coincide with the lightlike worldlines of QED; simulation time under decoherence is optimized. The construction replays the logic that leads to the QED Lagrangian. Namely, it starts from the Dirac quantum walk, well-known to converge towards free relativistic fermions. It then extends the quantum walk into a multi-particle sector quantum cellular automata in a way which respects the fermionic anti-commutation relations and the discrete gauge invariance symmetry. Both requirements can only be achieved at cost of introducing the gauge field. Lastly the gauge field is given its own electromagnetic dynamics, which can be formulated as a quantum walk at each plaquette.


Introduction
The simulation of quantum phenomena, in order to scale in size, requires the use of a quantum mechanical device.A short term application is to find ground states of Hamiltonians, that hold the key to certain molecular structures or condensed matter properties [1].In the longer term, however, it may even be used to simulate the dynamics of these from first principles, based on their constituent fundamental particles' dynamics.This has motivated a strand of works on the quantum simulation of quantum field theories (QFT) [2,3,4].All of them rely on a prior spatial discretization, but some are actually based on a spacetime discretization, allowing for a natively discrete account of both the relativistic and gauge symmetries.This work provides a quantum simulation scheme based on a discrete spacetime formulation of 3 + 1-quantum electrodynamics.
Lattice QFT.Lattice gauge theories appear for instance in quantum error correction theory (e.g.Kitaev's toric code [5,6,7]) and condensed matter (e.g. through their applications in spin liquids).But they originate from particle physics, where lattice quantum chromodynamics [8] has extensively been used in order to obtain numerical values, to then be compared against experimental values from particle accelerators.This procedure is paradigmatic of the way new physics is discovered, making simulation take a central role.However, these techniques are computationally heavy.Finding a way to simulate lattice QFT efficiently and accurately through a quantum device would be a game changer.
Non-relativistic, analogue simulation of QFT.The standard ways to quantum simulate QFT are fundamentally non-relativistic, as they begin asymmetrically by discretizing space but not time, by means of a Kogut-Susskind Hamiltonian [9,10].Next the matter (fermions) and the gauge field (bosons) degrees of freedom are encoded as quantum systems on the simulating device, whose interactions will mimic those of the Hamiltonian [3].These interactions are sometimes implemented as discrete-time unitaries, but even then these are short-time approximations of the Hamiltonian, as obtained by the Trotter formula under the non-relativistic ∆ t ≪ ∆ x assumption.This approach was recently realized experimentally using an ion trap architecture [11].Recent, classical but quantum-inspired tensor networks techniques, come to complement this standard approach [12].These use compact, approximate description of quantum states [13] such as the multiscale entanglement renormalization ansatz (MERA) [14,15], discarding hopefully unwanted information about the states as they evolve, so that the description keeps a manageable size, whilst attempting to keep track of interesting ingredients, including entanglement.Tensor network techniques, however, mainly focus on finding low energy states and will inevitably hit a scalability and precision barrier when dealing with many-body states and their dynamics.
Relativistic, digital simulation of QFT.In order to quantum simulate QFT in a relativistic manner, we must place space and time on an equal footing, discretizing both simultaneously, with parameter ∆ t = ∆ x .This leads of an infinitely repeating quantum circuit, across space and time, namely a quantum cellular automata (QCA).The speed of light in the simulated QFT will then, at each step, strictly correspond with the 'circuit speed', i.e. the maximal speed allowed by the wires.And whilst it is true that after multiple steps the circuit speed could potentially produce square lightcones instead of round lightcones, a robust continuum limit argument shows that this does not happen for Dirac QCA [16].In fact such schemes have been shown to recover discrete counterparts to Lorentz covariance in three different ways [17,18,19].This is in stark contrast with the earlier mentioned analogue simulation paradigm, where the ∆ t ≪ ∆ x assumption kills off Lorentz covariance and yields a non-strict and much lower speed of light for the simulated QFT, hopefully matching the Lieb-Robinson bound [20]-a fragile process however [21,22].Fig. 1 illustrates the circuit and light speed under both simulation paradigms.From a theoretical standpoint, relativistic, digital quantum simulation is therefore advantageous: (i) strict causality at each step is ensured as the circuit wires match the lightlike worldlines of the simulated QFT ; (ii) space and time are treated on an equal footing as demanded by relativity.From a practical standpoint these advantages are expected to translate into longer simulation times.For instance, where a quantum simulation device suffers from a given typical decoherence time τ , we could still simulate the QFT over a period of logical time of about τ .
A natively discrete approach to QFT.In the continuous, relativistic settings, the standard way to express a QFT is by means of a Lagrangian, i.e. a 'local cost function', which integrated over a possible history provides the action.The action is to be minimized in the classical theory, or to serve as phase in the quantum theory.The use of a particular Lagrangian is justified by means of relativistic and gauge symmetries.But the way it is then brought to quantum theory is only through regularization, typically through an energy cut-off, which is essentially a discretization of spacetime.Then a continuum limit is worked out on a per-case basis through renormalization, a process which is heavily dependent on the cut-off parameters.Hence, even though the Lagrangian approach starts as continuous, it does not actually solve the problem of obtaining a well-defined quantum theory in the continuum-that is except in the case of non-interactive theories.In general, given some Lagrangian, an elaborate continuum limiting process still needs to be worked out.Our aim is to acknowledge this fact, and instead express the QFT directly as a family of infinitely repeating circuits of local quantum gates, parametrised by the discretization step.In order to justify the use of a particular QCA we must then, just like in the Lagrangian approach, begin with a quantum walk (QW) accounting for free relativistic fermions and then extend it to the multi-particle sector QCA by imposing the fermionic anti-commutation relations as well as discrete gauge invariance, thereby deducing the need for a gauge field.We can then 'turn on' the interaction by providing the gauge field with a simple dynamics.That is, we must transpose the logics of construction that leads to a particular QFT in the Lagrangian approach, to a natively discrete setting, whose discretization parameter can then be made arbitrary small.This ought to provide a rigorous, natively discrete formulation of QFT.Closest work. 1 + 1 quantum electrodynamics (QED), also known as the Schwinger model [23], has been recovered under the relativistic, digital quantum simulation paradigm, by discretizing through ∆ t = ∆ x and following gauge theoretical justifications in [24].Next, this was generalized by [25,26,27] in order to allow for arbitrary ∆ t ≤ ∆ x , so that both the continuous spacetime limits (which exists when the interaction is turned off) and the continuous time discrete space limits (which always exists) may be taken, the latter coinciding with the Kogut-Susskind Hamiltonian.The aim of the present paper is to take the work of [24] to two and three spatial dimensions, yielding the first natively discrete spacetime formulation of a 'real-life' QFT, namely 3 + 1 QED.
Implementing the anti-commutation of fermions.Digital quantum simulation has been very successful at describing relativistic particles in different fields [28,29,30], but only a handful of works deal with interacting QFT with more than one particle [31,32,24,27].One of the difficulties is that in order to encode multiple fermions as qubits, one must enforce the anti-commutation of their creation/annihilation operators, e.g. through the Jordan-Wigner transformation.However, this method has all the looks of breaking locality, especially as soon as one considers more than one dimension of space.This was even formulated as a no-go result [33], stating that any QCA implementing the fermionic anti-commutation relations in two spatial dimensions would have very high internal space dimension, as in [34,35].In the tensor network community, anti-commutation, locality and low internal space dimensions do coexist, at the cost of introducing a cut-off for the gauge field and two extra fermions per links, called rishons.Moreover, in lattice gauge theory, a solution where the fermionic degrees of freedom are replaced by bosonic ones, at the cost of introducing two new fermionic degrees of freedom, has been developed [36].A similar idea, where just the parity of the gauge field is treated as a fermion, was sketched in Farrelly's PhD thesis [37].The first main contribution of this paper is to combine these ideas and formalize them in the discrete spacetime setting.We introduce no extra field, but replicate the gauge field information once for each direction.For each direction, its parity provides a rishon.Then the Jordan-Wigner transform needs only be implemented locally, at the level of each site.This does allow for a QCA of low internal space dimension, while enforcing fermionic anti-commutation.Ultimately, the reason why the no-go result [33] is circumvented is the presence of the gauge field, as well as our focus on expressing the actual dynamics in local manner-rather than the creation/annihilation operators.
Fully discrete magnetic contribution.In QED's gauge invariant states, the fermions are the sources of the gauge field lines.In one spatial dimension there is no magnetic term, lines are confined to the unique dimension, and thus they have no dynamics [24].But in two and three spatial dimension, the Hamiltonian has an added magnetic term, a.k.a. the plaquette term.The second main contribution of this paper is to introduce two possible discrete spacetime counterparts to the plaquette term.The first proposition works by simply integrating the plaquette term in the Fourier basis.But it requires a prior cut-off in the gauge field degrees of freedom, and allows for arbitrary changes in values within that cut-off, even in one time step.The second proposition takes the form of a local quantum walk (QW)-like evolution in the local gauge field degrees of freedom of each plaquette.It does not require a prior cut-off and ensures that gauge field values only change one step at a time.Both constructions agree in the continuum limit.
Further spin-dimensions.In one and two spatial dimensions, the Dirac equation is a PDE on a wave function having two complex amplitudes at each site, corresponding to the presence of a particle or an antiparticle.In the multi-particle settings, and because there can be no more than one particle in a given state and site, the four occupation numbers of a site are thus: no fermion, one particle, one antiparticle, and, both a particle and an antiparticle.This can be encoded as 2 qubits in the 2 + 1 Dirac QCA.Moving on to three spatial dimensions, the Dirac equation is a PDE having four complex amplitudes at each site, corresponding to the presence of a particle spin up, a particle spin down, an antiparticle spin down or an antiparticle spin up.In the multi-particle settings, the number of qubits per site has then to be increased from 2 to 4 qubits, which makes the 3 + 1 QCA more involved.The third main contribution is therefore to provide a construction of the 3 + 1 QED QCA in spite of this added complexity.
The paper is organized as follows.In Sec. 2 we set the conventions and show how to enforce the fermionic anti-commutation relations while allowing for a local definition of the gauge invariant operators that govern the dynamics of the theory.In Sec. 3, we gradually derive the gauge invariant dynamics of the 2 + 1 QED QCA starting from the 2 + 1 Dirac QCA and adding the electric and magnetic contributions-the simpler, two spatial dimensional case makes the argument clearer.In Sec. 4, we reach the 3 + 1 QED QCA.Finally, we provide some perspectives.

Enforcing anti-commutation, locality and gauge invariance
In QFT, fermionic particles are represented by means of operators that annihilate them or create them at position x.These are denoted a x and a † x respectively.Applying the annihilator on quantum states takes occupation number |1⟩ x to |0⟩ x and produces the null vector otherwise.The creator takes |0⟩ x to |1⟩ x and produces the null vector otherwise.Moreover, these operators are required to have the specific anti-commutation relations {a x , a † y } = δ x,y where {•, •} denotes the anti-commutator and δ the Kronecker delta.In order to obtain a quantum numerical scheme for a QFT, to be run on a generic quantum computer or some specific-purpose quantum simulation device, we need to encode the QFT degrees of freedom as a lattice of qubits.The natural point of departure is to interpret the occupation number degrees of freedom at each site (i.e.|1⟩ x vs |0⟩ x ), as qubits, thereby obtaining a lattice of qubits.But enforcing the fermionic anti-commutation whilst remaining local is non-trivial.
In Subsec.2.2, the impact of the fermionic anti-commutation relations upon the local operators that are needed to express the discrete time dynamics, will be carefully worked out.In order to obtain them, however, we crucially rely on there being a gauge field, as demanded by gauge invariance.

Introducing the gauge field
The QED Lagrangian is built by considering the Dirac Lagrangian for free fermions, and then demanding that it be gauge invariant, under U (1) gauge-transformations.This is impossible without introducing a new field, the gauge field, which in the case of QED turns out to be the electromagnetic field.
We will proceed in the same manner in the discrete.One reason for that is that in numerical analysis, the fact that a numerical scheme conserves the original symmetries is desirable and seen as a good sign of numerical stability.The other reason is more fundamental, as we aim to show that a natively discrete spacetime formulation of QED is just as legitimate at the Lagrangian formulation, in terms of its justification through symmetries.
Discrete gauge invariance in the context of classical cellular automata has been formalized in [38,39,40].Together with the treatment of gauge invariance upon quantum walks [41] and in lattice gauge theories in general [3], it inspired [24] to formulate discrete gauge transformations in the context of QED QCA in the following manner.
Let φ : Z → R. We call discrete gauge transformation the global operator g φ = x g(φ(x)) which associates, to each position x, an on-site gauge transformation g(φ(x)) parametrized by φ(x) and defined as Thus, a discrete gauge transformation g φ is essentially a space-dependent phase, exempt of any regularity requirement, applied at every point of the lattice in accordance to the occupation number at that point.To be gauge invariant, the evolution of a QCA must commute with every possible gauge transformation.The Dirac QCA, which solely describes moving fermions, is not gauge invariant unless we introduce the gauge field.The argument boils down to the elementary fact that, as a particle moves from position x to the adjacent position x + η, the discrete gauge transformation will trigger a phase φ(x) applied beforehand, or a phase φ(x + η) applied afterwards.It follows that fermionic transport does not commute with gauge transformations.
In order to fix this, one needs to introduce a gauge field on the links.Consider a link between sites x and x + η.A gauge field is much like a doorman/bouncer counter standing on the link.It adds one whenever a particle crosses the link.Actually, in the present paper, we will not place just one bouncer per link, but two-one at each end, where the ends are denoted x : η and x + η : −η respectively.The bouncer at x : η counts positively the number of fermions leaving x towards the link, and negatively those entering x from the link.The bouncer at (x + η) : −η counts negatively the number entering x + η from the link, and positively those leaving x + η towards the link.Now fermionic transport acts as and the discrete gauge transformation triggers a phase (l + 1)φ(x) − lφ(x + η) regardless of whether it is applied before of after the move.The restriction that the two gauge fields of a link be opposite of signs, as in |l⟩ x:η and |−l⟩ (x+η):−η , is quite natural if we think of them as holding the total number of fermions that went through the link, ever, and nothing else.Still, this restriction is not a necessity for gauge invariance.However, lifting it would allow for unnecessary degrees of freedom, that are not demanded by gauge invariance and seem to be unphysical in practice.So, we must not just ask for gauge invariance, but for 'minimal' gauge invariance, i.e. the demand the gauge field be obtained as a relative gauge extension [40], making this restriction a necessity.We therefore impose this restriction across the grid, except of course at the boundaries, where it becomes vacuous.
Placing two gauge fields per link is non-standard but has several advantages: (i) each gauge field is well-localised on a site, (ii) it gauge transforms in the same way as the fermions on that site, (iii) this is number conserving, (iv) it will help to implement the fermionic anti-commutation relations in a local manner as we will now see.

(Anti-)commuting annihilation and lowering operators
Consider the lattice generated by unit vectors describing the space directions.These vectors are denoted µ, ν and κ, in two spatial dimensions, only µ and ν are used.
At each lattice site x lies a group of d qubits, each stating whether a fermion in mode j ∈ 0 . . .(d − 1) is present at the site.The possible modes correspond to the number of internal degrees of freedom (e.g. the spin) of the fermions.This encoding captures the Pauli exclusion principle as there cannot be two fermions in the same mode at the same site.
Each link between x and x + η (where η ∈ {±µ, ±ν, ±κ}) has a gauge field attached at both ends: one at x : η and the other a (x + η) : −η.Each gauge field lives in the Hilbert space of integers H Z .
The local electric counting operator, denoted E x:η , is the local observable of the gauge fields, i.e. it acts as and as the identity elsewhere.Following the restriction that the two gauge fields of a link be of opposite signs, one has E x:η = −E x+η:−η .
The lowering operator of the gauge fields is r x:η .It acts as r x:η |l⟩ x:η = |l − 1⟩ x:η and as the identity elsewhere.Most often we need to act on both the gauge fields of a link with One has U x:η = U † x+η:−η -i.e.lowering the gauge field value attached to one site corresponds to raising on the other site.
The local parity observable, denoted Z x , is defined by acting with the on-site operator Z at site x, and as the identity elsewhere.On each qubit at the site, Z acts like the third Pauli matrix σ 3 .On each gauge field at the site, it acts as Z |2l⟩ x:η = |2l⟩ x:η and Z |2l + 1⟩ x:η = − |2l + 1⟩ x:η .
All the operators described so far were 'local' operators, by which we mean that they were of the form L x = M x ⊗ I with M x acting just at site x (and sometimes direct neighbours of x).We will now represent the fermionic annihilation a x,j and gauge lowering operators V x:η of the QFT as products of these local operators.These implementations will not be local, as we must meet the desired (anti-)commutation relations: •] the commutator and δ x,y corresponding to Kronecker delta.These commutation relations are commonly implemented by means of the Jordan-Wigner (JW) transform.This will be the basis for the redefined version of the operators.However, in 2 + 1 and 3 + 1 dimensions, it leads to non-locality of the operators expressing dynamics if used as is [33].In order to fix this we use the idea hinted in [37] and treat the parity of the gauge fields as a fermion.Moreover, we use two gauge fields per link, so that this parity plays the role of the rishons of lattice gauge theories.
Let ≺ denote a so-called JW order between all the fermionic occupation number qubits and gauge fields present in the model.The sole requirement to make this work is that the degrees of freedom at any given site be contiguous in the JW order, i.e. that they follow each other.In what follows we will use, quite arbitrarily, that for a given site x, the fermionic modes and the gauge fields are ordered as We define fermionic creation and gauge lowering operators based on this order.Let x be a position, j a mode and η a direction: where a, b is short for [min{a, b}, max{a, b}[.When x : η ≺ x + η : −η, this interval does not contain x + η : −η and r x+η:−η commutes the Z y s leading to the simplification: The following fermionic anti-commutation relations are ensured by the JW transform: Since V x:η is composed of two fermion-like operators, Eqs. ( 7) and (8) yield: Hence, the (anti-)commutation of Eqs.(3), ( 4) and (5) are ensured.

Local evolution operators
The annihilator operators a † x,j are not local since they have an infinite trail of Z on-site operators.The lowering operators V x:η can also be non-local as soon as both ends of the link are far apart in the JW order.However, these operators are never used by themselves in the evolution.They may perhaps be used by themselves for some initial state preparation.But Physics does not require that initial state preparation be local.Moreover, in the context of initial state preparation, a † x,j could be interpreted as a one-step implementation of having the particle come from infinity (or from the lattice border) to its current position through a product of local steps.
What matters physically is the evolution be local.One can then consider the operators allowing to express that: (i) the movement of fermions be that of free fermionic QCA, (ii) the gauge field induce an interaction between fermions-through the electric contribution, (iii) the gauge field vibrates-through the magnetic contribution.From these requirements, we can define the simplest local gauge invariant operators (checking for gauge invariance is postponed till Subsec.2.4): E 2 x:η (11) with η, ζ ∈ {µ, ν, κ} and η ̸ = ζ.Eq. ( 9) represents a mass term (i.e. a fermionic only operation, local to a site, changing the mode), Eq. (10) corresponds to a fermion hopping term (i.e. a transport), Eq. ( 11) is the squared electric operator defined previously, and Eq. ( 12) is a plaquette term (i.e. a local vibration of the gauge field).The QED QCA evolutions will be defined based upon these local operators.The mass term is indeed local because it involves a pair of fermionic operators acting on the same site, cancelling out the trail of Z outside this site.The electric term is local by definition.
For the locality of the hopping term, developing the gauge field operators V x:η as pairs s x:η and s † x+η:−η on each side of the link, allows for the pairing of an annihilator a with an s on each site, thus cancelling out the Z outside the sites: The remaining Zs are local to sites x and x + η.The locality of the pair s † x:η a x,k , which is the local brick from which the hopping term is built, is illustrated in Fig. 2. x:ν a x .The coloured dots and lines corresponds to Z operators acting on fermions and gauge fields respectively.Each operator Z y is applied exactly twice which cancels them out except on site x.Here, only one fermionic mode per site is represented, for clarity.
For the locality of plaquette term-i.e.Eq. ( 12)-notice that it forms a small loop of four gauge field links.Developing each gauge operator V x:η as a pair s x:η and s † x+η:−η , and reordering the resulting product of s, one gets two anti-commuting operators per site.Hence, the string of Z cancels out outside the site they act on: where the minus sign that appears between the second and third line comes from the anti-commutation of s x:η with the other 7 operators.In order to better understand the structure of the local, minus signs in this term, it is helpful to break down the plaquette term into constituent, smaller local operators.Indeed, let us define the corner operators c x:η,ζ as These are local to x since the remainder Z is on-site x.A corner operator would be the result of a fermion having come from direction η, passed through site x, and left in the direction ζ.The plaquette operator can then be redefined as: Thus, Eqs. ( 13) and (15) allow the fermionic and bosonic dynamics terms (left-hand-side of the equations) to be expressed into simpler local dynamics (right-hand-side).
Notice that the only dependency of the definition of these local evolution operators w.r.t the JW order, is per-site.In fact the per-site JW order could even have been chosen different from one site to the other.

Gauge invariance
Gauge invariance was introduced in Sec.2.1 to motivate the need for a gauge field.It remains to be checked that above-defined local evolution operators are gauge invariant.

Gauge invariant operators
Gauge invariance is the commutation with gauge transformations, i.e. space-dependent phases proportional to the occupation number.Because the annihilation and lowering operators act on the occupation number, they individually are not gauge invariant.However, the local evolution operators combine multiple annihilation and lowering operators so that the total occupation number of a site is unchanged (only the distribution inside the site is modified)-e.g. a fermion leaving site x will reduce the number of fermions by one, but increase one of the gauge fields at x by one.
In order to formally show this, let us first write the commutation relations for the individual annihilation operators a x,j a x,j g φ = e iφ(x) g φ a x,j a † x,j g φ = e −iφ(x) g φ a † x,j , the gauge field operators on half-links s x,η and the gauge field operators V x:η Based on these commutation relations, we derive those of the local evolution operators.For the mass term: Hence, the mass term is gauge invariant.For the hopping term: Hence, the hopping term is gauge invariant.The electric term does not change the occupation number, hence it is directly gauge invariant: As for the plaquette term, we start by showing the gauge invariance of the corner operator: Since the plaquette term is a product of four corner operators, its gauge invariance is ensured.
The QCA operators that will be defined in the next sections are linear combination or exponentials of the hopping, electric and magnetic terms.Therefore, their gauge invariance is ensured by the individual gauge invariance of those terms.

Gauge invariant states
In a gauge theory, gauge transformations should not be observable.Let O be an observable, ρ a density matrix representing a state and g φ a gauge transformation.That gauge transformations should not be observable amounts to asking that Tr Og φ ρg † φ = Tr(Oρ).There are two common ways to enforce this: namely to restrict the set of observables to the gauge invariant ones [O, g φ ] = 0, so that Tr Og φ ρg † φ = Tr g φ Oρg † φ = Tr g † φ g φ Oρ = Tr(Oρ) or to restrict the set of states to the gauge invariant ones, so that Tr Og φ ρg † φ = Tr Oρg φ g † φ = Tr(Oρ).We opt for the second, demanding that for every gauge transformation, Let us draw the consequences of this demand that the states be gauge-invariant.Any density matrix ρ is a convex linear distribution over pure states.In the case of a pure state |ψ⟩ ⟨ψ|, this commutation relation amounts to forbidding superposition across any two eigenspaces of a gauge transformation, so that g φ |ψ⟩ ⟨ψ| = e iλ |ψ⟩ ⟨ψ| = |ψ⟩ ⟨ψ| e iλ = |ψ⟩ ⟨ψ| g φ with e iλ the eigenvalue of the eigenspace of g φ to which |ψ⟩ pertains.In other words, |ψ⟩ = c∈S α c |c⟩ is a superposition of particular basis states |c⟩, i.e. taken in some subset S such that for all φ there exists λ, such that for all c ∈ S, g φ |c⟩ = e iλ |c⟩.Let f (x) denote the sum of the occupation numbers, for both the fermions and the gauge fields, at each site x of c.Because the λ = x φ(x)f (x) and φ(x) is arbitrary, the gauge-invariance therefore imposes that f (x) be the same for c ∈ S. Thus, gauge-invariance demands that there exists some fixed, classical occupation number function f (x), and that pure states be considered physical if and only if, at each position x, the occupation number operator yields f (x): The physical pure states then form a subspace, call it the f (x)-subspace.The choice of a particular f (x) can be interpreted physically as the choice of a classical fixed, external electromagnetic field.
Having fixed f (x) and thus S, one may wonder about the operators which allow us to prepare some basic state c ′ ∈ S given initial basic state c ∈ S, e.g. by creating a fermion.As seen previously, the fermionic annihilators a x,j is not gauge invariant; it does not preserve occupation numbers and will take us out of S. Following ideas from [42], each fermion creation operator a † x,j could be turned into a gauge invariant state preparation by accompanying it with a gauge field lowering operator V x:η , but this changes the occupation number at position x + η, which in turn has to be compensated by a V x+η:ζ , and so on until a boundary (possibly infinitely far) is reached.The following defined a gauge invariant state preparation: a † x,j,p = a † x,j y:η∈p where p is a gauge field path from position x to the space boundary (possibly infinite).
If the lattice is finite, the last operator is understood to be a half-link transformation s, instead of V , since there would be no end to the last link.(Notice that a † x,j s x:η is also gauge invariant operator, but disallowed as it breaks the restriction that E x:η = −E x+η:−η , except at the boundary.) The state preparation a † x,j,p follows the prescribed anti-commutation relations.Indeed, recalling Eq. (5), we have that a x,j commutes with any V y:η , hence a † x,j,p also does.Moreover, within two operators a x,j,p and a y,k,q , the fermionic parts a x,j and a y,k anti-commute, while every other pair of operators commutes, enforcing the fermionic anti-commutation relation (3).
Another gauge invariant state preparation is the creation of a gauge field loop, i.e.
y:η∈p V † y:η where p is a cyclic gauge field path, or one that begins and ends at a boundary.
Pair creation.Our construction is completely general and offers the possibility of working in specific subspaces to mimic physically-relevant phenomena, such the pair production, i.e., the spontaneous creation of electrons and positrons from the QED vacuum [3,43].For this purpose, let us consider the sector in which for all site x, f (x) = d/2.In this subspace, a natural choice for representing the canonical vacuum is the state in which all the antiparticle degrees of freedom are populated, whereas no particle and gauge fields are present (this state can be viewed as the reference "Dirac sea").For instance in 2 + 1 dimensions, when d = 2, the canonical vacuum would have |10⟩ (x,1)(x,0) at each site x, with all gauge fields set to |0⟩ x:η .Acting on such a state, a term like a x,1 a † x,0 yields |01⟩ (x,1)(x,0) .This can be understood as 'creating particle, namely an electron, at (x, 0)' at the cost of 'creating an antiparticle hole, namely a positron, at (x, 0)'.This evolution describes the process of electron-positron pair creation from the reference vacuum.More precisely it models the very moment when, at site x, the two particles are created : it is only once they move away from each other that the gauge fields will get updated, so as to account for a gauge field line between them two.The reverse evolution describes electron-positron pair annihilation.

2 + 1 Quantum Cellular Automaton
Using the previously defined local, gauge invariant operators, it is possible to define a QCA that accounts for QED in 2 + 1 dimensions.We proceed in three steps, following the same principles as that leading to the QED Lagrangian.First, a gauge invariant free dynamics for the fermions is defined through a generalization of a Dirac quantum walk (QW) to multiple walkers.Second, the electric contribution is defined as one of the simplest gauge invariant operator acting on the gauge field based on the electric operator, following ideas from [24].It is found to match the Trotterization of the electric part of the Kogut-Susskind Hamiltonian [9].Third, the magnetic part is added as one the simplest gauge invariant operator acting on the gauge field with the lowering and raising operators.In fact, two constructions are provided for this magnetic term, both agreeing in the limit and matching a Trotterization of the magnetic Kogut-Susskind Hamiltonian.

Fermionic dynamics
Let us first recall the 2 + 1 dynamics of the fermionic field, without any electromagnetic contribution, without any gauge field even, and restricting to the one particle sector.This is the well-known Dirac quantum walk (QW) [44,16].
At each lattice site x lies a group of 2 qubits, each stating whether a fermion in mode j ∈ 0 . . . 1 is present at the site (remember that in 2 + 1 dimensions the Dirac Eq. is a PDE on a wave function having two complex amplitudes).Again, this encoding captures the Pauli exclusion principle as there cannot be two fermions in the same mode at the same site.Moreover, as we focus on the one particle sector first, we temporarily look at the case where only one of the qubits can be in the state |1⟩.
One time-step is divided in three sub steps: a vertical translation, a horizontal translation and a mass term.Moreover, each translation is decomposed into two swaps.Each of these terms thus acts on a pair of qubits: on a single site for the mass and the first swap (S), and across two neighbouring sites for the second swaps (T ).In the one particle sector, the possible input states are |00⟩, |01⟩ and |10⟩.Number conservation forces |00⟩ to be mapped to itself.Without loss of generality we can assume that this absence of particle triggers no phase.Hence, our on-site (or across-two-sites), two qubit operators for the Dirac QW are of the form W = 1 ⊕ M where 1 leaves |00⟩ unchanged and M is the unitary acting on the subspace spanned by |01⟩ and |10⟩.The Dirac QW is a global operator expressed in terms of a large product of such on-sites operators (which individually are independent of x): where C ϵ = 1 ⊕ e −imϵY is the mass term, H µ = 1 ⊕ H with H the Hadamard operator (such that HZH = X), S swaps qubits (x, 0) and (x, 1), and T η swaps qubits (x, 1) and (x+η, 0).It follows that (x,1),(x+η,0) T η x S = 1⊕e ϵZ∂η is the displacement operators in the η direction by a factor ϵ, moving the first qubit in the positive direction and the second in the negative direction.The X, Y, Z are Pauli matrices and ∂ η the partial derivatives.Convergence towards the Dirac Eq. is rigorously proven in [44,16].Again this was without any electromagnetic contribution, without any gauge field even, and restricting to the one particle sector.
Moving on to multiple walkers, the Dirac QW becomes the Dirac QCA.The dynamics needs to be extended to take into account the case where multiple qubits are in state |1⟩.Since the QW on-sites operator acted on at most two qubits, and due to number conservation, only the input state |11⟩ requires our attention, and it can only be sent to itself, up to a phase.To find out exactly which phase has to be applied, let us move to the Heisenberg picture.
The Heisenberg picture tells about the future impact of our past actions.For instance, say that the overall evolution from time t to time t + 1 is governed by a unitary operator W, e.g.mapping |ψ⟩ to |ψ ′ ⟩.Then, the past action of creating a fermion at (x, 0) at time t, as implemented by a † x,0 , e.g.mapping |ψ⟩ to a † x,0 |ψ⟩, will have future impact Wa † x,0 W † at time t + 1, i.e. |ψ ′ ⟩ to Wa † x,0 W † |ψ ′ ⟩.More specifically consider S the multi-particle sector extension of S which is such that S x a † x,0 S † x = a † x,1 and S x a † x,1 S † x = a † x,0 , with S x acting as S on site x and as the identity elsewhere.Such an S is called 'fermionic swap'.Notice there exists S ′ which coincide with S (and thus with S) in the one-particle sector, but that do not obey these two equations.We discard them on the basis that S is the 'rightful non-interacting extension' of S, as the future impact of a † x,0 ought to be a † x,1 regardless of there being other particles or not.
We see that this phase got determined as a consequence of the (anti-)commutation relations hypothesis discussed in Sec. 2, as well as the way we chose to implement the annihilation operators, under which Jordan-Wigner transform etc., as defined in that same section.This requires the presence of gauge field, which we had temporarily ignored for describing the Dirac QW, but we now restore for describing the Dirac QCA.We use bold fonts to denote the Dirac QCA counterparts of the Dirac QW operators.
The transport sub-step takes a right-moving (resp.left-moving) qubit at position x (resp.x + η) and maps it to the right-moving (resp.left-moving) qubit at position x + η (resp.x).It does so by first swapping the qubits on each site using the on-site operator S, and then moving them through an across-two-sites operator T η , whilst changing the gauge field accordingly.Here is S as acting on some |ψ⟩ (x,1)(x,0) , Thus, Here is T η as acting on some |ψ⟩ (x+η,0)(x, 1) , where K η is with j = 1, k = 0, and U x:η the action of the gauge field lowering operator on (x : η)(x + η : −η).Let T x,η be the local operator which acts as T η across sites x and x + η and as the identity elsewhere.We have The equations using fermionic annihilators and creators are based on the local evolution operators from Eqs. ( 9) and (10).The product of Z y comes from the hopping term defined in Eq. (13), with x : η the link along which the swap takes place.The minus one, when the input qubits are in state |11⟩ (x+η,0)(x,1) , is the exchange phase for crossing fermions.Again full justification is given in Appendix A.
The transport sub-step of the Dirac QCA is illustrated in Fig. 3.The basis change H is similar to the mass term in that it is an on-site operator that can be written as a direct sum for the case with 0, 1 or 2 particles as follows (as acting on some |ψ⟩ (x,1)(x,0) ): Let H x be the local operator acting as H on site x and as the identity elsewhere.We have ) .

This last equation corresponds to the basis change expressed using the local evolution operator (9);
The minus one is justified in Appendix A. Since H † = H, we have H † = H.The complete Dirac QCA is: It is represented in Fig.  Gauge invariance.Each of the operators constitutive of the Dirac QCA has been expressed as a linear combination of the local evolution operators given in Subsec.2.3, which were proven to be gauge invariant in Subsec.2.4.Therefore, the Dirac QCA is gauge invariant.

Electric contribution
Let us now define the electric contribution.To do so, we follow the same idea as in the Lagrangian formalism, which is to take some of the simplest gauge invariant electric term.
We then check that it matches the Trotterization of the Kogut-Susskind Hamiltonian.The construction proposed here is highly inspired by [24].
A simple electric contribution ought to act on the gauge field according to a functional of the local electric counting operator E x:η for η ∈ {µ, ν}, as defined in (1).We could also demand invariance under E x+η:−η instead.However, the electric operator itself does not match this requirement because E x:η = −E x+η:−η .The squared electric term E 2 x:η , does.The squared electric term itself is not unitary, but its exponential e iE 2 x:η , is.Notice also that the spacetime discretization should impact the amplitude of the phase: when any of the space or the time discretization parameter goes to zero, the exponential should tend towards the identity.Taking these considerations into account, one obtains: where the ϵ 2 comes from the simultaneous discretization in space and time with ∆ x = ∆ t = ϵ, and g E is the coupling constant as determined experimentally.
Without the electric contribution, the gauge field just keeps track of fermions passing by, but has no influence upon them.With the electric contribution, the gauge-field dependent phase, mediates the interaction.
Let us compare the electric contribution with the electric part of the Kogut-Susskind Hamiltonian [9] for QED: Integration over a ∆ t period of time-i.e.computing e i∆tH E -and Trotterizing to separate the spatial sum into a product of exponentials, yield D E .
Remark 1 (Truncation of the electric field).If 1 2 ϵ 2 g 2 E = 2π k exactly, with k an integer, then the phase of the electric contribution will wrap up around 2π as soon as E 2 reaches k.This happens when ϵ = 4π kg 2

E
. If we restrict ourselves to such values of ϵ, i.e. decreasing it by augmenting k, then, as far as the electric contribution is concerned, the gauge field can equally well be represented with a k-dimensional Hilbert space H k .Indeed, when E outputs k + l, its square gives k 2 + 2kl + l 2 , and the phase 1  2 Notice that this k still goes to infinity when taking ϵ to zero, augmenting proportionally to 1/ϵ 2 .This idea of restricting the gauge field to finite-dimensions labelling roots of unity was suggested in [24].Truncations in lattice gauge theory simulations were evaluated in [45].
Gauge invariance.The electric contribution is gauge invariant as an exponential of a gauge invariant operator.

Magnetic contribution
In order to define the magnetic contribution, we follow the same path as for the electric contribution, that is to say define some of the simplest gauge field-only term that is local, unitary and gauge invariant.

Formulation of the magnetic contribution. The magnetic contribution acts by lowering or raising the gauge field.
As discussed in Sec. 2 operators that are defined solely using the lowering and raising operators need to form a loop (spatially) in order to both be local and ensure gauge invariance.Indeed, raising the gauge field at one end of a link, implies lowering it at the other, hence imposing that the magnetic contribution be a loop; the smallest loop is realized by the plaquette local evolution operator P x:η,ζ of Eq. (15).
Since the loop is oriented, one may wish to symmetrize.This is done through a sum P x:η,ζ + P † x:η,ζ , but this breaks unitarity.Just like for the electric field, unitarity is restored through exponentiation, and the space and time discretization parameters need be used in the exponential to ensure that it goes to the identity when these go to zero.Taking these considerations into account, one obtains: with ϵ = ∆ t = ∆ x and g 2 M /2 a constant.Taking the limit of D M when ϵ goes to zero yields: Let us compare the obtained magnetic contribution with magnetic part of the Kogut-Susskind Hamiltonian: where g M = 1 ∆xg E .Integrating this Hamiltonian over a ∆ t period of time-i.e.computing e i∆tH M -and Trotterizing to separate the spatial sum leads to the same local operator Two formulations in terms of gates.In order to define D M in terms of quantum gates, one needs to explicitly compute or closely approximate the exponential.We provide two constructions in order to do so.The first construction is to diagonalize the plaquette term so that taking its exponential simply amounts to exponentiating the eigenvalues, that is to say, we go from an electric basis to a magnetic one [46].The quantum Fourier transform is used in the diagonalization.The second construction consist in a more subtle approach were the term P x:µ,ν + P † x:µ,ν is written as the sum of two terms P x:µ,ν and Q x:µ,ν whose exponentials are simple to compute.This second construction yields a gate formalism reminiscent to that of a QW.
Ordering the operators.Two neighbouring plaquette local evolution operators P x:µ,ν and P x+µ:µ,ν both act on the two gauge fields of the link between sites x + µ and x + µ + ν.Therefore, the order of the operations could have been relevant.However, the plaquette local evolution operators actually commute.If one insists on not acting simultaneously with two different operators on a same system, then any arbitrary ordering can be chosen, such as applying the plaquette local evolution operators at even positions (x 1 + x 2 + x 3 mod 2 = 0) first, and then at odd positions.
Redefining the states.Before diving into the two constructions, let us introduce a new way of representing the gauge field around a plaquette.Forgetting to denote the second gauge field of each link, which is just the opposite of the first, the gauge fields of the four links of the plaquette are Let us introduce 'plaquette-like' states: and see how these may be affected by a loop of lowering operators, first using U operators of (2) instead of the V , i.e.
This will just act as a shift on the fourth number |abcn⟩ → |abc, n + 1⟩ (24) and as the identity elsewhere.The plaquette local evolution operator of ( 6) and (12), is a loop of V however, i.e. a loop of U combined with Z operators.These Z operators only add a plus or minus sign to the state |abcn⟩.This sign may be influenced by the other gauge fields and qubits present at the four sites where the plaquette term is applied, but not beyond, cf. the locality stated in (15).In order to be able to determine this sign, we extend the plaquette-like states to encompass the states of all the qubits and gauge fields of the plaquette that were not already accounted for in |abcn⟩.More specifically, we now take as plaquette states |abcno⟩ where |o⟩ is a canonical basis state of the systems We can then let φ abcno ∈ {0, 1} be such that where P is the on-plaquette operator such that P x:µ,ν acts as P on the four sites of the x : µ, ν plaquette, and as the identity elsewhere.Starting from state |abc, 0, o⟩, one reaches state |abc, n, o⟩ by applying the plaquette operator n times.We can then let ψ abcno ∈ {0, 1} be such (−1) ψ abcno |abcno⟩ = P n |abc, 0, o⟩ .
The ψ abcno can be given explicitly in terms of φ abcno , so that ψ abcno + φ abcno = ψ abc,n+1,o : This leads to the definition of 'plaquette states' | abcno⟩ = (−1) ψ abcno |abcno⟩ which verify the same relation as that of Eq. ( 24), with V operators instead of U : Again, these plaquette states help understand the plaquette local evolution operator as a shift, which will be useful for both decompositions of the magnetic contribution in terms of quantum gates.The state and plaquette operator are illustrated in Fig. 5.

Gauge field truncation.
If truncating the gauge field to H k , one needs | abc, k, o⟩ = | abc, 0, o⟩ so that any sign discrepancy in Eq. ( 25) is avoided.This paragraph is just to overcome this technicality.In terms of phase the condition is that ψ abcko = ψ abc0o = 0, i.e. that ψ abcko = 0≤n<k φ abcno be even.Each term φ abcno is actually a sum φ abcn + φ o , where φ abcn depends on the four links under transformation, whereas φ o depends on the other gauge fields and qubits at the four sites, which are not modified by the plaquette operator.Both of these depend on the JW order chosen.Since φ o is constant when acting only using the plaquette term, for an even k its contribution to the sum 0≤n<k φ o is even.
(x, µ) Additionally, for any n, we have that φ abc,n = φ abc,n+2 as all the parities are equal.As a consequence φ abc,n + φ abc,n+1 + φ abc,n+2 + φ abc,n+3 is even.It follows that for k a multiple of 4, the sum 0≤n<k φ abcn is even.Hence, for k = 4q, with q an integer, the truncation is valid.Relying on a specific JW order actually allows for any even k.Indeed, we can enforce that φ abcn always be equal to 0: this is the case when (x : −η) and x : η are both inferior (or both superior) to (x : −ζ) and (x : ζ) for η and ζ two distinct directions.In this case, the four corner operators that form a plaquette defined in Eq. ( 14) will induce a phase (−1) φ abcn through the following Z operators (disregarding the part which contributes to φ o ): Notice that each link appear twice, hence they do not induce any phase, i.e. we have that P | abcno⟩ = (−1) φo | abc, n + 1, o⟩.Thus, the truncation is well-defined for k even for that specific JW orders.

Derivation through diagonalization
The first construction moves from the electric to the magnetic basis [46].In other words, it works by computing the eigenvectors and eigenvalues of the plaquette operator, so that taking the exponential of the operator simply amounts to exponentiating the eigenvalues.
A single plaquette local evolution can be seen as a shift operator upon plaquette states as illustrated Fig. 5b.That shift amounts to a phase in the Fourier basis, therefore the Fourier transform diagonalizes the plaquette.To define this Fourier transform, the state space of the gauge field needs be truncated to Z k , for instance according to Remark 1, and in accordance with the restriction of Subsec.(3.3).Then the eigenvectors of a plaquette local evolution are for a, b, c and p integers in 0, ..., k − 1.These are the Fourier transform of the | abc, n⟩ states, that can be obtained through a quantum Fourier transform (FT) whose circuit representation is well-known [47].The corresponding eigenvalues are The eigenvectors of the hermitian conjugate of the plaquette operator are the same but with eigenvalues λ −p .For the sum of the plaquette and its hermitian conjugate, the eigenvalues are: e 2πip/k + e −2πip/k = 2 cos(2πp/k).
Having found an eigenbasis of the plaquette term, it is now easy to exponentiate it.Doing so one gets the eigenvalue equation: This defines a diagonal operator Diag that contains the eigenvalues exp iϵ 2 g 2 M cos(2πp/k) .The magnetic evolution for the QCA can thus be written as: where FT is the quantum Fourier transform.Hence, Eq. ( 26) thus defines a circuit of quantum gates for the magnetic term.In this construction, no approximation has been done, hence the limit is exactly the one given in Eq. (23).

Quantum walk-like derivation
The second construction uses a reformulation of the plaquette operators to make it look like a quantum walk.Here we sometime abuse notations and write | abcn⟩ or even just | n⟩ to talk about a state | abcno⟩.Indeed, the plaquette term can be divided into two operators, one which acts as | 2n⟩ ⟨ 2n + 1|+| 2n + 1⟩ ⟨ 2n| and the other as | 2n + 1⟩ ⟨ 2n + 2|+ | 2n + 2⟩ ⟨ 2n + 1|-i.e. the two operators act as shifted swaps between pairs of sites.Let P and Q denote those operators: We have that P + Q = P + P † .
P and Q are hermitian, unitary and their eigenvectors are with n even for P and odd for Q.The corresponding eigenvalues are plus and minus one.Now taking the exponential of P , we have, with n even: These equations are identical for exp iϵ 2 g 2 M 2 Q when taking n odd.This is reminiscent of a one dimensional quantum walk with coin cos(θ) i sin(θ) i sin(θ) cos(θ) .
Such QW is homogeneous in the 'tilde' basis, but at first glance, it seems inhomogeneous in the canonical basis because of the operator which induces a phase φ abcno , dependent upon n, in the canonical basis.However, as explained in Subsec.3.3, this phase can be made independent of n through a specific choice of JW order, recovering a homogeneous QW even in the canonical basis.Consider A way to implement the plaquette term is to apply those successively in one time step: D Indeed, taking the limit when ϵ goes to zero after applying both the exponentials of P and Q to | n⟩ (where n is even) gives: For n odd, the exponential of P would yield the state | n − 1⟩ instead of | n + 1⟩ (in the second line of the equation) and the exponential of Q the state | n + 1⟩ instead of | n − 1⟩ (in the third line).Therefore, the same limit would be reached.
In D M this is done at every plaquette, matching Eq. (23).Thus, both quantum gate implementations agree in the limit as Gauge invariance.Gauge invariance of the plaquette term has been verified in Sec.2.4.1.Its exponentiation is a linear combination of them, thus still gauge invariant.

Complete dynamics
Combining the fermionic dynamics, the electric and magnetic contribution, one obtains the following complete dynamics for the 2 + 1 QED QCA: where D F (fermionic term) refers to Eq. (20), D E (electric term) refers to Eq. ( 21) and D M (magnetic term) refers to Eq. ( 22) or (27) depending on the quantum gate implementation chosen.

3+1 Quantum Cellular Automaton
This section extends the previous 2 + 1 QED QCA construction, to reach a 3 + 1 QED QCA.In 3 + 1 dimensions, the Dirac Eq. is a PDE on a wave function having four complex amplitudes, i.e. there are four fermionic modes instead of two.As for the gauge field, no additional degree of freedom is required, but the third dimension needs to be taken into account when considering the electronic and plaquette terms.
In the following, to avoid mixing notations, every operator referring to the 3 + 1 case will be overlined, e.g. the on-site operator modelling the mass in the 3 + 1 Dirac QW will be denoted C ϵ .Again the operators of the multi-particle sector QCA are in bold, e.g.C ϵ .

Fermionic dynamics
The procedure is the same here as in the 2 + 1 case: (i) define the 3 + 1 Dirac QW, (ii) extend each on-site operator using the Heisenberg picture so that it acts on the full state space and not just in the one-particle sector-cf.Appendix A.
Let γ j be the following generalized Pauli matrix: They respect the anti-commutation relation with δ j,k the Kronecker delta.The state space for the Dirac QW is that of 4 qubits restricted to the one-particle sector, i.e. there are 5 possible states: |0000⟩, |0001⟩, |0010⟩, |0100⟩ and |1000⟩.As in the 2 + 1 case, state |0000⟩ is mapped to itself because of number conservation.Hence, every on-site operator is of the form 1 ⊕ U where U acts on the remaining four dimensional subspace.The 3 + 1 Dirac QW has the same structure as that in 2 + 1: transport sub-steps in each direction, each of these being surrounded by a basis change, and lastly a mass sub-step [44,16,48]: where S swaps qubits (x, 0) and (x, 1) with (x, 2) and (x, 3), and T η swaps qubits (x, 2) and (x, 3) with (x + η, 0) and (x + η, 1) so that x T η x S = 1 ⊕ e ϵ(Z⊗I)∂η , transporting the first two qubits in the positive direction and the last two in the opposite one.
The operators H η are defined so that, as acting on some |ψ⟩ (x,3)(x,2)(x,1)(x,0) : One possible choice for them is as follows, as acting on some |ψ⟩ (x,3)(x,2)(x,1)(x,0) : where H is the Hadamard matrix such that HZH = X and is a matrix such that F Y F = Z.These last equality can be checked easily when noticing that Then The same line of reasoning applies to H = X+Z √ 2 .Note that all three H η are Hermitian.The on-site operator for the mass in 3 + 1 dimensions is quite similar to that of the 2 + 1 dimensional case.We have, as acting on some |ψ⟩ (x,3)(x,2)(x,1)(x,0) : with c = cos(mϵ) and s = sin(mϵ).
The derivation of the operators making up the 3 + 1 Dirac QCA requires us to extend the on-site operators any configuration for the 4 qubits, hence possibly more than two fermions crossing.In general this could be complicated to analyse.Fortunately, all of the above 4-qubit on-site operators can be reexpressed as products of 2-qubit on-site gates, with the 2 qubits being adjacent in the JW order.Thanks to this, the 2 + 1 methodology and results readily apply, see Appendix A.
The mass sub-step of the 3 + 1 Dirac QCA is reached by: 1/ Seeking to represent the 3 + 1 Dirac QW on-site operator C ϵ as a circuit of one-particle-sector 2-qubit gates, as done in Eq. (31).2/ Replacing each 2-qubit gate of this circuit by its multi-particle sector extension, so as to obtain the 3 + 1 Dirac QCA on-site operator C ϵ .3/ Repeating C ϵ across space.That is because the QCA should act as the QW in the one-particle sector.This is represented in Fig. 6.We get, as acting on some |ψ⟩ (x,3)(x,2)(x,1)(x,0) : where S is the swap defined in Eq. ( 17) and C ϵ is the 2+1D mass term defined in Eq. ( 16).The transport sub-step of the 3+1 Dirac QCA is recovered from that of the 3+1 Dirac QW in a similar fashion.It moves the first two qubits in the η direction while moving the last two in the direction −η.This is done by first swapping the first two qubits with the last two using a on-site operator S, then applying a transport operator which will swap these pairs of qubits with neighbouring ones in the η direction using an across-two-sites operator T η .Again each of these is obtained by decomposing the corresponding 3 + 1 Dirac QW on-sites operators into two qubit-gates, and replacing them by their multi-particle sector extensions, as illustrated in Fig. 7.
For S we get, as acting on some |ψ⟩ (x,3)(x,2)(x,1)(x,0) : For T η we get, as acting on some |ψ⟩ (x+η,1)(x+η,0)(x,3)(x,2) : where S is defined as in Eq. ( 17), and T η is defined as in (18) but with j = 3.The basis changes are of two kinds in the 3+1 dimensional case.The first kind is that of H κ from (30) which exchanges the second and last qubits of the QW.This can be realized through the operator H κ represented as a circuit in Fig. 8 and defined, as acting on some |ψ⟩ (x,3)(x,2)(x,1)(x,0) , by: The other two basis change-H µ and H ν -share a common formalism, the QW counterpart is simply the operator H (resp. F ) applied independently to the first and last two qubits, followed by the H κ operator.This can be immediately translated into QCA operators: where H is defined in Eq. (19) and F is the generalization of F using the procedure in appendix A, i.e., as acting on some |ψ⟩ (x,3)(x,2)(x,1)(x,0) : The full fermionic dynamics is therefore:

Electric and magnetic contribution
The electric contribution works in the same way as in the 2 + 1 dimensional case, i.e. the exponentiated squared electric operator is applied at every link: The magnetic contribution needs to be generalized to take into account the three dimensions, and thus the three possible directions for the plaquettes.Therefore, the magnetic contribution of the 3 + 1 QED QCA is the same as in the 2 + 1 case, but it is applied three times: one for each pair of directions.Let D M,η,ζ from Eq. ( 22) or (27) (depending on the formulation one chooses) denote the magnetic contribution along the two spatial dimension η and ζ.We now have: Again this evolution coincides, in the limit, with the magnetic part of the Kogut-Susskind Hamiltonian: x η,ζ∈{µ,ν,κ} η̸ =ζ

Complete dynamics
Combining the 3 + 1 Dirac QCA of (32) with the electric (33) and magnetic (34) contributions, one obtains the 3 + 1 QED QCA: Its gauge invariance is ensured by the same arguments as in the 2 + 1 case.

Conclusion
Summary of contributions.In this paper we constructed a quantum cellular automata (QCA) accounting for QED in 2+1 and 3+1 dimensions.The construction follows the same principles used to build the Lagrangian formulation of QED-i.e. free anti-commuting fermions, gauge invariant, simplest electric and magnetic term.But here spacetime is discrete, and space and time are treated on equal footing.The evolution is described in terms of local quantum gates, whose wirings coincide exactly with the speed of light of the QED.To reach our goal, we needed three contributions.The first contribution was the formulation of gauge invariant local evolution operators-Eqs.( 13) and (15)-that meet the specifications imposed by the (anti-)commutation relations (3), ( 4) and (5) of the fermionic and bosonic annihilators/lowering operators they are made of.It was, to us, a surprise that this could be achieved since it is in apparent contradiction with the no-go result of [33].But, following ideas of [37,12] the gauge field came to the rescue.This was used to obtain a 2 + 1 Dirac QCA, i.e. a generalization of the 2 + 1 Dirac QW to multiple walkers, recovering the free fermionic dynamics D F (20) which is represented in Fig. 4.
The second contribution was to derive the electric D E (21) and magnetic D M (22) contributions, leading to the 2+1 QED QCA in Eq. (28).The magnetic contributions was the most challenging, but in the end two possible quantum circuit representations were found.The first uses a diagonalization, through a Fourier transform, in order to exponentiate plaquette terms (26), albeit up to truncation.The second approach formulates the magnetic contribution in terms of a quantum walk over the gauge field Hilbert space H Z (27).Both approaches lead to the same continuum limit (23) and were checked to correspond to integration of the magnetic part of the Kogut-Susskind Hamiltonian.This provides a first-to the authors' knowledge-discrete spacetime formulation of the magnetic term.
The third contribution was the extension of the QCA to a 3 + 1 QED QCA (35) which required raising the state space from 2 to 4 qubits per sites to represent fermions.On paper, this meant working out the consequences of the anti-commutation relations for up to four fermions, under sophisticated changes of basis.Fortunately, these 4-qubit gates were decomposable as 2-qubit gates.The electric and magnetic contribution were straightforwardly extended from the 2 + 1 case.Altogether, the 3 + 1 QED QCA provides a first, relativistic discrete spacetime formulation of a real-life quantum field theory.
It is not clear whether the proposed 2 + 1 QED QCA and 3 + 1 QED QCA hereby constructed are 'staggered' in the traditional sense, because particles and antiparticles live at the same positions.But they are at least staggered-like in the sense of the 1 + 1 QED QCA of [24].In 1 + 1, this suffices to escape the fermion-doubling problem, whereby high momentum wavepackets result in low energy states.We leave it as an open question whether fermion-doubling is escaped in a similar way in the proposed 2 + 1 and 3 + 1 models.
Perspectives.These two and three-dimensional QED QCA are quantum circuits.D F , D E and D M can be expressed in terms of standard universal gates such as CNot, Phase, Hadamard.Thus, the QCA is directly interpretable as a digital quantum simulation algorithm, to be run on a Quantum Computer.A first perspective is the implementation of this QCA on quantum computers.This simulation scheme is efficient, in that it requires O(s d /∆ d x t/∆ t ) gates in order to simulate a chunk of space of size s, over t time steps with ∆ x the space resolution, ∆ t the time resolution, and d the space dimension.The output of this is a quantum state and the simulation may need to be run multiple times in order to obtain meaningful statistics.However, classically just the state space itself is an O(e s d /∆ d x ) as it grows exponentially with the number of quantum systems to be simulated.The classical time complexity is thus O(e s d /∆ d x t/∆ t ).The exponential gain here clearly comes from the fact that the scheme simulates a multi-particle systems, just like in Hamiltonian-based multi-particle quantum simulation schemes.Quantum walkbased simulation schemes on the other hand, are by definition in the one-particle sector, and thus can only yield polynomial gains.
An immediate continuation of this work would be to further parametrize the QED QCA, so as to make it 'plastic' enough so that we may be able to take a discrete space continuous time limit of the model, and prove that one recovers the Kogut-Susskind Hamiltonian, as was done for 1 + 1-QED in [27].In this paper, the QFT considered is an abelian gauge theory.Another extension of this work is to consider non-abelian gauge transformations so as to recover the quantum chromodynamics (QCD).Looking at different geometries for the underlying space, another extension would be to define the QCA over a triangular or tetrahedral spatial discretization of space, as was done for quantum walks [49].Finally, notice that the distinction between fermions and interacting-hardcore-bosons [50, 51] is wearing thin with this local model, and yet seems to persist as embodied by the Z terms of Eqs. ( 13) and ( 14).This is intriguing, and one truly wonders whether the distinction is physically observable, or can be proven otherwise.
[47] Peter W Shor. "Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer".SIAM review

A From quantum walk to quantum cellular automata operators
We explain the process of extending a one-particle sector on-site operator making up a quantum walk, into a multi-particle sector on-site operator making up a quantum cellular automata.We start with the 2 + 1 dimensional case.
The QCA on-site operators act on 2 qubits and are number conserving, which constrains them.Indeed, the evolution for the input |00⟩ (x,1)(x,0) will always be the identity without loss of generality, and the evolution for the input |11⟩ (x,1)(x,0) is only a phase since it is the only state with occupation number equal to 2. Therefore, any QCA on-site operator W can be written as a direct sum W = W ⊕ e iφ with W = 1 ⊕ M the quantum walk on-site operator for the one-particle sector and φ a phase to be determined.To find out this phase exactly, we use the Heisenberg picture.
The Heisenberg picture describes the future impact of our past actions.Consider the past action a † x,1 at t: Its future impact at time t + 1 is Suppose that the φ we seek to determine, is equal to θ.With this supposition, the above simplifies and we have that the future impact of a † x,1 is just Wa † x,1 W † = M 11 a † 1 + M 01 a † x,0 .A contrario, with any other choice of φ we would be constructing some W ′ which, albeit coinciding with W and thus W in the one-particle sector, would for instance leave a phase between a † x,1 a † x,0 a x,0 and a † x,1 a x,0 a † x,0 in the future impact of a † x,1 .The future impact of a † x,1 would again yield a superposition of a particle at (x, 0) or at (x, 1), but with a phase depending on another particle being there or not.Such a W ′ would not, therefore, be the 'rightful non-interacting extensions of W ' to the multi-particle sector.Only W is. In other words, we were seeking to determine the φ of the last entry of W; setting it to θ fixes it to its non-interactive value.
The same process can be applied to the evolution a † x,0 leading to Wa † x,0 W † = M 00 a † x,0 + M 10 a † x,1 Considering as past action the product of these two operators a † x,0 a † x,1 .Its future impact is
In other words we have that for any QW on-site operator W = 1 ⊕ M acting on qubits, the corresponding non-interactive multi-particle extension QCA on-site operator is Let us now use this to define the on-site and transport gates of the QCA.
On-site QW gates that are mere permutations of qubits in the one-particle sector, can then be extended to the multi-particle sector directly by means of products of S, in the same way that any permutation can be obtained from local transpositions.
For the mass term and the basis changes, this is not the case, but Eq. (36) again readily applies.As an example, the mass term C ϵ = 1 ⊕ C has C 00 C 11 − C 01 C 10 = c 2 + s 2 = 1, which results in the following mass on-site operator for the QCA, expressed as acting over some |ψ⟩ (x,0),(x,1) :

A.2 Transport operators
In this paper the transport is implemented in two steps: first swap the qubits on-site so that the right-moving qubit is on the right-hand side of the site (S), then hop qubits across adjacent sites (T η ) whilst changing the gauge field accordingly.In the one particle sector, the across-two-sites operator T η , as acting on some |ψ⟩ (x+η,0)(x,1) , is given by:

A.3 3+1 dimensions
In 3 + 1 dimensions, there are four qubits per site instead of two.A QW operators W that acted in the one particle sector |0001⟩ (x,3)(x,2)(x,1)(x,0) , |0010⟩ (x,3)(x,2)(x,1)(x,0) , . .., now needs to be extended to a QCA operator W that can handle the multi-particle sector |0011⟩ (x,3)(x,2)(x,1)(x,0) , |1110⟩ (x,3)(x,2)(x,1)(x,0) , . .., which seems to leave open many more possibilities than in the 2 + 1 case.However, the 3 + 1 QW operators we are dealing with can easily be decomposed as circuits of 2-qubit gates, with the qubits following each other in the JW order.So, we extend these two qubit gates instead, through the same process as in 2 + 1 dimensions, and then recombine them to form the extension of the 3 + 1 QW operator.The 2-qubit QW operators used will turn out to be same as those of the 2 + 1 case, except for 1 ⊕ F which was not defined previously.Its QCA version F is given through Eq. (36), as acting on some |ψ⟩ (x,1)(x,0) : Let us first decompose the 4-qubit operators for the 3 + 1 QW that was given in Eq. (29), as circuits of 2-qubit gates.
One notices that it acts on the zeroth and second qubits on the one hand, and on the first and third qubits on the other.In order to obtain a circuit of adjacent gates in the JW order, The first and second qubits are swapped, using the operation 1 ⊕ X ⊕ 1, allowing for the application of the 2-qubit mass term on the leading and last two qubits separately through C ϵ ⊕ C ϵ .The swap is then reapplied such that the initial order is recovered.These swaps can be understood as crossing of wires in a circuit so that the 2-qubit gates have the correct inputs and outputs.Such swaps correspond, in the multi-particle sector, to the QCA operator S. Hence, we can extend 1 ⊕ X ⊕ 1 into 1 ⊗ S ⊗ 1, as this again swaps the two middle qubits while leaving the rest unchanged.Notice that the operator ⊗ is used instead of ⊕ since we are no longer in the one particle sector.Similarly, the second step of the circuit for the mass term is extended into C ϵ ⊗ C ϵ .In the end, one obtains the following multi-particle sector, QCA operator for the mass term, as acting on some |ψ⟩ (x,3)(x,2)(x,1)(x,0) , (i.e.w.r.t the canonical basis { |0000⟩ (x,3)(x,2)(x,1)(x,0) , |0001⟩ (x,3)(x,2)(x,1)(x,0) , |0010⟩ (x,3)(x,2)(x,1)(x,0) , |0011⟩ (x,3)(x,2)(x,1)(x,0) , . ..}

Figure 1 :
Figure 1: (a) In relativistic, digital quantum simulation, the light-like worldlines of the simulated theory coincides with circuit wires, yielding strict causality.(b) In non-relativistic, Trotterized analogue quantum simulation, light-like worldlines are approximately recovered through a Lieb-Robinson bound, and are slower.Thus, the simulation is running slower.As typical decoherence times match the depth of the circuit, the QFT is simulated over a shorter period.

Figure 2 :
Figure 2: Visualization of locality for s †x:ν a x .The coloured dots and lines corresponds to Z operators acting on fermions and gauge fields respectively.Each operator Z y is applied exactly twice which cancels them out except on site x.Here, only one fermionic mode per site is represented, for clarity.

Figure 3 :
Figure 3: Transport with qubit (x, 0) moving right (and (x, 1) moving left) and T updating the gauge fields accordingly (as represented as a single wiggly line for conciseness).

Figure 4 :
Figure 4: The 3 steps of the free evolution in the multi-particle sector.The gauge field is omitted for clarity.

Figure 5 :
Figure 5: Plaquette states and operator represented in the subspace Z 4 containing four gauge field values.

Figure 7 :
Figure 7: Representation of the transport