Quantum Electrodynamics from Quantum Cellular Automata, and the Tension Between Symmetry, Locality and Positive Energy

1. Introduction

Quantum walks (QWs) [1,2,3,4] connect information theory with physical law by providing a discrete model of particle dynamics. In such systems, a quantum walk is based on a simple set of principles and symmetries, with a particle located at any of the vertices of a graph or lattice. In each time step, the particle moves to one of the neighboring vertices (connected by an edge). The particle’s time evolution is given by a unitary transformation that, in the continuum limit, can lead to a relativistic wave equation such as the Weyl or Dirac equations, or to Maxwell’s equations [5,6,7,8,9,10,11,12,13,14]. In recent years this correspondence has been extended to quantum field theories (QFTs) [15,16,17,18,19,20,21]. The idea is that one can construct a quantum cellular automaton (QCA) from the quantum walk and then recover the desired Lorentz invariant QFT in the limit of continuous time and energies.

Quantum information processing comprises the action of a sequence of unitary operations (quantum gates) on an initial state of a set of quantum systems (e.g., qubits, which are the fundamental units of quantum information). In QFT the time development of a quantum field is given by the action of a unitary operator on a state describing quantum particles or the creation and annihilation operators corresponding to them. The QCA/QFT correspondence described here suggests that systems of the latter type can arise from systems of the former type. This reveals a deep and fundamental connection between those two seemingly unrelated systems.

For example, quantum field theories are traditionally constructed by starting with a classical field and then quantizing it, while the QCA/QFT correspondence suggests that what is fundamental is not a classical field but an underlying quantum walk. This is reminiscent of Wheeler’s suggestion “It from Bit,” in which he proposed that information plays a significant role in the foundation of physics [22]. The universe, he suggested, may be fundamentally an information processing system from which the apparent reality of matter somehow emerges. The correspondence may also have more practical implications—providing a new calculational approach to simulation on a quantum computer [23,24,25].

But there is another reason to consider the relation of discrete models to quantum field theory: the possibility that space may indeed be discrete, or at least that there is different physics at the Planck scale. Though continuous spacetime is an almost universal assumption in physics, even in classical physics the idea of a point particle such as an electron inevitably leads to divergences. Gravity, too, does not allow an arbitrarily high concentration of gravitational charge (i.e., energy) in a small region of spacetime, because if too concentrated an object will collapse to a black hole. Yet to make observations at shorter and shorter distances, one must focus ever larger energies in a small volume. In both QFT and general relativity (GR) these short-range problems become serious—both theories are plagued by singularities and infinities at short distances. In QFT the usual remedy is renormalization but the methods are mathematically dubious. In GR, it is broadly accepted that at small distance scales the theory breaks down and is no longer applicable. As a result, some have speculated that a quantum theory of gravity will naturally lead to the idea that spacetime is not continuous, with the likely values of a fundamental unit of space and time being (at least roughly) the Planck length $\ell \approx {10}^{-33}$ cm and the Planck time $t\approx {10}^{-43}$ sec. In this view, spacetime is composed of discrete, tiny cells, as is the case in the QCA models of field theory.

An obvious objection to discrete spacetime is that we would lose Lorentz invariance. But though the lattice QCA models do not admit Lorentz invariance, it is recovered in the limiting case of small lattice spacing (or long wavelength). Lorentz symmetry, then, could be an approximation, and effective symmetry valid at length scales long compared to the fundamental length. Probing shorter lengths would uncover small deviations from the Dirac evolution. In particular, a lattice structure would be reflected in a modification of the free-particle dispersion relation. In a discrete universe with a fundamental length and time on the Planck scale, current technology could not distinguish between full Lorentz invariance and the approximate invariance of the discrete system. This situation may soon change and experiments involving neutron diffraction have been proposed to detect the Lorentz non-invariance that would result if a lattice-based theory represented the true laws of nature, with the Dirac theory being only a continuous approximation [26].

In this paper we will show how free Quantum Electrodynamics (QED) arises from quantum walks on a lattice and their corresponding QCAs. We will detail how imposing a handful of basic principles leads to the Fermi (Dirac) and Bose (Maxwell) theories and discuss how they can be coupled. We show that a difficulty in going to an interacting theory is that the symmetries of these QCAs—particularly translation and time-reversal symmetry—lead to theories with both positive- and negative-energy solutions. These solutions are highly delocalized, which means that localized couplings will in general lead to the production of negative-energy particles. We discuss this problem, analyze a one-dimensional example, and look at ways to overcome it.

In Section 2 we begin by defining appropriate discrete-time quantum walks for fermions and bosons on a lattice [14,21]. The walk represents the one-particle sector of the theory. Once this is defined, it can be promoted to the multi-particle case in much the same way that one defines a Fock space from a single-particle theory. In Section 3 we impose the condition that the walk not have a preferred lattice axis. We’ll see this implies that although the dynamics will not be Lorentz invariant for any finite lattice spacing, Lorentz invariance will be recovered in the long wavelength limit. In Section 4 we diagonalize the physical space (but not the internal space) factor in the evolution operator by changing to the momentum representation of the wave function. For low-momentum fermions with mass, and bosons that are massless, we take the interval $\Delta t$ to be small and define an approximate time derivative operator, allowing us to obtain the momentum space representation of the Dirac and Maxwell equations. In Section 5 we show it is possible to promote the single-particle fermionic and bosonic theories to multi-particle QCAs that reduce to the usual Dirac and Maxwell field theories in the long-wavelength limit. In Section 6 we discuss symmetries of the lattice theories, and their roles, especially the relation of time reversal symmetry to negative energy states. In Section 7 we discuss issues in the inclusion of interactions between the Fermi and Bose theories. Section 8 presents a short summary of our results.

2. Structure of the Walks

In deriving quantum field theories of fermions and bosons from quantum cellular automata, the first step is to define an appropriate discrete-time quantum walk on a lattice [5,9,10,11,12,13,14,21]. The walk represents the one-particle sector of the theory. One this is defined, it can be promoted to the multi-particle case in much the same way that one defines a Fock space from a single-particle theory.

The quantum walk is a unitary analog of a classical random walk. In a classical random walk, in each time step the particle has some probability of moving to the neighboring vertices. In a quantum walk, the evolution is unitary, and the particle can move to a superposition of positions.

We can think of the possible particle positions as being the vertices of a graph, with edges between neighboring vertices. For the quantum walks in this paper, the graph is always regular, with every vertex having the same number of neighbors. The number of neighbors is the degree of the graph, labeled by d. We can label each of the d outgoing edges from a vertex $1,\dots ,d$, with that label describing a direction the particle can move. To start with, we will assume that this graph is finite, but we will later consider the infinite limit.

To maintain both unitarity and nontrivial dynamics, the particle generally must have an internal degree of freedom, a “coin” space. The Hilbert space has the form ${\mathcal{H}}_{\mathrm{QW}}={\mathcal{H}}_X\otimes {\mathcal{H}}_C$, where ${\mathcal{H}}_X$ is the Hilbert space of the particle position, and ${\mathcal{H}}_C$ is the coin space. As we’ll see, for fermions the coin space has even dimension and the coin space operators that appear in the time evolution are operators of half-integral spin. For bosons the corresponding operators are rotation group generators with integral spin.

The lattice we consider is the bcc lattice, consisting of a set of vertices $\mathbf{x}\equiv (x,y,z)=(i\Delta x,j\Delta x,k\Delta x)$, where $i,j,k$ are integers and $\Delta x$ is a fixed lattice spacing. Thus two vertices are neighbors if one is (x,y,z) and the other is $(x\pm \Delta x,y\pm \Delta x,z\pm \Delta x)$. A discrete-time walk on this lattice involves moving one unit $\Delta x$ in each of the three cardinal directions at each time step. In the standard formulation of a quantum walk in three spatial dimensions, the evolution from one time to the next is given by a unitary evolution operator $U_{\mathrm{QW}}$ of the form

$$|{\psi }_{t+\Delta t}\rangle =U_{\mathrm{QW}}|{\psi }_t\rangle ,$$

where

$$U_{\mathrm{QW}}|{\psi }_t\rangle =\left(\sum _j{S}_j\otimes A_j\right)\left(I\otimes C\right)|{\psi }_t\rangle .$$

(1)

The $\left\{S_j\right\}$ are shift operators that move the particle from its current position to its neighbor in the direction j and the operators $A_j$ and the unitary C act on an internal space. C scrambles the faces, so that one does not constantly move in the same direction. The overall unitarity of the walk is not guaranteed: unitarity requires that the operators $\left\{A_j\right\}$, which act on an internal space, satisfy the condition

$$U{U}^{†}=\sum _{j,j^{\prime }}{S}_j{S}_{j^{\prime }}^{†}\otimes A_j{A}_{j^{\prime }}^{†}=I.$$

(2)

In this paper we use a different formulation [14,27,28]. We decompose each time step into three successive moves, one in each cardinal direction:

$$U_{\mathrm{QW}}=e^{i\theta Q}\left(S_X{P}_X^{+}+S_X^{†}{P}_X^{-}\right)\left(S_Y{P}_Y^{+}+S_Y^{†}{P}_Y^{-}\right)\left(S_Z{P}_Z^{+}+S_Z^{†}{P}_Z^{-}\right)$$

(3)

for fermions, and

$$U_{\mathrm{QW}}=e^{i\theta Q}\left(S_X{P}_X^{+}+P_X^0+S_X^{†}{P}_X^{-}\right)\left(S_Y{P}_Y^{+}+P_Y^0+S_Y^{†}{P}_Y^{-}\right)\left(S_Z{P}_Z^{+}+P_Z^0+S_Z^{†}{P}_Z^{-}\right)$$

(4)

for bosons.

In the above, if the particle is a fermion, at each time step it will move to a neighboring position in each direction. If it is a boson, there is the additional possibility that it remains in place when one or two of the direction factors is applied, but not all three (which would correspond to nothing happening, and thus can be excised from our description). Thus, for fermions the $P_j^{\pm }$ operators are three pairs of orthogonal projectors acting on the coin space, satisfying $P_j^{+}+P_j^{-}=I$. For bosons the $P_j^{\pm }$ and $P_j^0$ operators are three triplets of orthogonal projectors acting on the coin space, satisfying $P_j^{+}+P_j^{-}+P_j^0=I$.

The unitary coin flip operator is $C\left(\theta \right)=e^{-i\theta Q}$ where Q is a Hermitian operator that acts on the internal space. Q flips the forward and backward directions, i.e., $Q{P}_j^{+}Q=P_j^{-}$. For bosons, we also have $Q{P}_j^{0}Q=P_j^0$. We can, without loss of generality, assume that $Q^2=I$ and $\mathrm{Tr}\left\{Q\right\}=0$. Note that the equations we require for Q are obeyed by the parity transformation, P: $P|\mathbf{x}\rangle =|-\mathbf{x}\rangle$ implies that $P{P}_j^{+}P=P_j^{-}$, $P{P}_j^{-}P=P_j^{+}$.

The idea is that the quantum walk proceeds in a manner analogous to the classical walk, through a series of coin flips. The projectors $P_j^{\pm }$ (and $P_j^0$ in the case of bosons) correspond to different faces of the coin, which indicate in which direction to move; the unitary C scrambles the faces, so that one does not constantly move in the same direction. Unlike classical random walks, this evolution is always invertible: interference effects play a profound role in the time development of the system.

This novel formulation is inherently unitary. For fermions the $\gamma$ matrices (and the Dirac equation as the long wavelength limit) arise from the condition that the different walk axes and directions be uncorrelated (as we explain below). We get an analogous result for bosons, leading to the Maxwell equations.

3. Internal Space Operators from Symmetry Requirements

The Dirac and Maxwell equations are Lorentz invariant but one cannot demand rotational invariance of a walk on a discrete lattice. Here we impose the condition that the walk not have a preferred lattice axis. We’ll see that although the dynamics will not be Lorentz invariant for any finite lattice spacing, Lorentz invariance will be recovered in the long wavelength limit.

By “not having a preferred axis” we mean that it should be equally likely for the walker to move in any of the eight allowed directions. This is equivalent to demanding that movement along the different axes is, generically, not correlated. A move in the positive X direction, for example, should not bias the move in the positive or negative Y or Z directions. To achieve this, we impose a noncorrelation condition on the 3D walk, which we call the equal norm condition.

To formulate the equal norm condition for fermions, we begin by noting that shifts in the positive or negative direction along an axis, j, are determined by the projectors $P_j^{+}$ or $P_j^{-}$, respectively [1,2,3,4,14]. If we require $\mathrm{Tr}\left\{P_j^{+}\right\}=\mathrm{Tr}\left\{P_j^{-}\right\}$, the walk will not be biased in the positive or negative direction. We also do not want the shifts along different axes to be correlated. We achieve that by requiring that any eigenstate of $P_i^{\pm }$ have equal amplitude in the $P_j^{+}$ and $P_j^{-}$ subspaces, for all $i\ne j=X,Y,Z$. This requirement is equivalent to the following condition on the projectors:

$$P_i^k{P}_j^{+}{P}_i^k=P_i^k{P}_j^{-}{P}_i^k={\displaystyle \frac{1}{2}}{P}_i^k,$$

(5)

where $k=\pm$, $i,j=X,Y,Z$, and $i\ne j$.

If we introduce operators $\Delta P_j=P_j^{+}-P_j^{-}$, we can see that this equal norm condition is equivalent to requiring that the operators, $\Delta P_j$ all anti-commute with each other. Given that, it is easy to see that $\{Q,\Delta P_j\}=0$ as well. Four mutually anticommuting operators require a coin space of dimension $D\ge 4$, leading us to identify

$$Q={\gamma }_0,\Delta P_X={\gamma }_0{\gamma }_1,\Delta P_Y={\gamma }_0{\gamma }_2,\Delta P_Z={\gamma }_0{\gamma }_3,$$

(6)

where the $\left\{{\gamma }_i\right\}$ are the gamma matrices from the Dirac equation. In other words, we have

$$P_j^{+}={\displaystyle \frac{1}{2}}\left(I+{\gamma }_0{\gamma }_j\right),P_j^{-}={\displaystyle \frac{1}{2}}\left(I-{\gamma }_0{\gamma }_j\right).$$

(7)

The equal norm condition has a simple geometric interpretation. There are two projection operators associated with each dimension, corresponding to the forward and backward directions, and they project onto orthogonal spaces. If the internal space is four-dimensional, for example, each of those spaces is two-dimensional, so they project onto orthogonal planes within the four-dimensional internal space. Similarly, the pairs associated with the other two space dimensions project onto orthogonal planes. However, the planes associated with X-dimension projection operators are not orthogonal to the planes associated with Y- and Z-dimension projectors. The equal norm condition requires that all those nonorthogonal planes be oriented at equal angles to each other. It implies that the form of the quantum walk is preserved by rotations that map the lattice axes into each other. Also, if we change the order in which we apply the steps in the X, Y, and Z directions, we will get another walk of the same form, with the same properties, and this walk will have the same continuum limit.

Because bosons have a different internal space, allowing for the possibility that the walker remain in place, the equal norm condition for bosons is different [21]. As before,

$$P_i^k{P}_j^{+}{P}_i^k=P_i^k{P}_j^{-}{P}_i^k=c{P}_i^k,$$

(8)

where $k=\pm$, $i,j=X,Y,Z$, $i\ne j$, and c is some positive real constant.

With regard to the $P_j^0$ operators, we also want the amplitude for the particle to remain in its initial location without moving to be zero. This condition requires that $P_X^0{P}_Y^0{P}_Z^0=0$, which implies that either $P_X^0{P}_Y^0=0$ or $P_Y^0{P}_Z^0=0$ or both. Because we want to treat motion along all three axes similarly and not have the form of the evolution operator be highly dependent on the order of the shifts, we choose the symmetric condition

$$P_X^0{P}_Y^0=P_Y^0{P}_Z^0=P_Z^0{P}_X^0=0.$$

(9)

Finally, we require an equal norm condition for the $P_j^0$ that is analogous to that above, for the $P_j^{\pm }$:

$$P_i^k{P}_j^0{P}_i^k=c^{\prime }{P}_i^k,P_i^0{P}_j^k{P}_i^0=c^{\prime }{P}_i^0,$$

(10)

where again $k=\pm$, $i,j=X,Y,Z$, $i\ne j$, and $c^{\prime }$ is some positive real constant.

Again, there is a simple geometric interpretation, analogous to that in the fermion case. Here, the $P_j^0$, $P_j^{+}$ and $P_j^{-}$ are one-dimensional projectors; each projects onto some vector ${|k\rangle }_j$, where $k=0,+,-$ and $j=X,Y,Z$. The above conditions tell us:

• For fixed j = X, Y, Z, the three vectors ${|k\rangle }_j$ are orthogonal.
• The vectors ${|0\rangle }_j$ and ${|0\rangle }_{j^{\prime }}$ with $j\ne j^{\prime }$ are orthogonal.
• The inner products ${}_j\langle k|k^{\prime }\rangle _{j^{\prime }}$ for $k,k^{\prime }=+,-$ and $j^{\prime }\ne j$ all have equal magnitude.
• The inner products ${}_j\langle 0|k\rangle _{j^{\prime }}$ for $k=+,-$ and $j^{\prime }\ne j$ all have equal magnitude.

These conditions are met if we choose

$$P_{X,Y,Z}^{\pm }={\displaystyle \frac{1}{2}}\left(J_{X,Y,Z}^2\pm J_{X,Y,Z}\right),\mathrm{and}{P}_{X,Y,Z}^0=I-J_{X,Y,Z}^2,$$

(11)

where the J matrices form a spin-1 representation of $O\left(3\right)$:

$$J_X=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& -i\\ 0& i& 0\end{array}\right),J_Y=\left(\begin{array}{ccc}0& 0& i\\ 0& 0& 0\\ -i& 0& 0\end{array}\right),J_Z=\left(\begin{array}{ccc}0& -i& 0\\ i& 0& 0\\ 0& 0& 0\end{array}\right).$$

(12)

This gives us $\Delta P_j=J_j$.

Note that $J_i$ plays the role of ${\gamma }_0{\gamma }_i$ in the above fermion expression for $\Delta P_{X,Y,Z}$. For fermions, in the 1D and 2D cases as well as the massless 3D case, only a two-dimensional internal space is necessary, with $\Delta P_i={\sigma }_i$. For the 3D massive fermion case, the dimension of the internal space must be doubled, with $\Delta P_i={\sigma }_i\to \Delta P_i=\sigma \otimes {\sigma }_i={\gamma }_0{\gamma }_i$, where the choice of the first $\sigma$ depends on the representation (it is ${\sigma }_X$ for the Dirac representation, with ${\gamma }_0={\sigma }_Z\otimes 1$). The doubling is due to the number of anticommuting operators required and leads to a theory accommodating two spin states as well as particles of both matter and antimatter.

In our theory the internal space of the boson theory must also be doubled. As we’ll see when we calculate the dynamics, the model just described already includes anti-particle/negative energy states (as well as longitudinal zero-energy states). Since the photon is its own antiparticle the negative energy sector does not carry additional information, and we will restrict ourselves to the positive energy sector of the free theory. We will see however that the dynamics of the above theory ties the polarization to the sign of the energy: right-handed photons come with a positive energy and left-handed photons with a negative energy, and there is no parity operator. To see this last, note that the parity operator $\mathbb{P}$ would satisfy $\mathbb{P}\Delta P_i\mathbb{P}=-\Delta P_i$, which is not possible for all three operators, i = 1, 2, 3 (see Equation (12). A doubling of the internal space is thus necessary to allow the inclusion of both transverse polarizations, and a parity operator. We will therefore work with a six-dimensional internal space. Versions of Maxwell’s equations with six-dimensional field vectors have been studied in other contexts (see for example Ref. [29]).

To double the internal space we proceed along the lines of the fermion case:

$$J_i\to \sigma \otimes J_i,$$

where we will choose $\sigma ={\sigma }_Z$ to give a Dirac-like representation. With this substitution, the projection operators become

$$\begin{array}{ccc}\hfill P_j^0& =& 1\otimes (1-J_j^2),\hfill \\ \hfill P_j^{+}& =& {\displaystyle \frac{1}{2}}(1\otimes J_j^2+{\sigma }_Z\otimes J_j),\hfill \\ \hfill P_j^{-}& =& {\displaystyle \frac{1}{2}}(1\otimes J_j^2-{\sigma }_Z\otimes J_j),\hfill \end{array}$$

(13)

and $\Delta P_j={\sigma }_Z\otimes J_j$.

We can now define a parity operator $\mathbb{P}={\sigma }_X\otimes 1$, which takes $P_j^{+}$ into $P_j^{-}$ and vice versa. This could be taken as the coin operator $Q=\mathbb{P}$, as we did in the fermion case. However, as we’ll see, the coin generator endows the particle with mass; so, since we are interested in (massless) photons here, we will set $Q=0$.

To preserve a common structure for both fermions and bosons it is useful to define boson gamma matrices in analogy with the usual fermion ones:

$$\begin{array}{ccc}\hfill {\gamma }_0^{\mathrm{bose}}& \equiv & {\sigma }_X\otimes 1=\mathbb{P},\hfill \\ \hfill {\gamma }_i^{\mathrm{bose}}& \equiv & -i{\sigma }_Y\otimes J_i.\hfill \end{array}$$

(14)

With these choices, $\Delta P_i={\sigma }_Z\otimes J_i={\gamma }_0^{\mathrm{bose}}{\gamma }_i^{\mathrm{bose}}$.

4. Momentum Space Picture and the Dirac and Maxwell Equations

We can diagonalize the physical space factor in the evolution operator (though not the factor that operates on the internal space) by changing to the momentum representation of the wave function. If the position states are $|x,y,z\rangle$, where $x=j\Delta x$, $y=k\Delta x$, and $z=l\Delta x$ for some integers $j,k,l$, then the momentum states are

$$|k_X,k_Y{k}_Z\rangle ={\displaystyle \frac{1}{\mathcal{N}}}\sum _{j,k,l}{e}^{-i\left(k_{X}j\Delta x+k_{Y}k\Delta x+k_{Z}l\Delta x\right)}|j\Delta x,k\Delta x,l\Delta x\rangle ,$$

(15)

where $-\pi /\Delta x<k_{X,Y,Z}\le \pi /\Delta x$. For mathematical simplicity, we assume that the overall size of the lattice is finite, albeit large, with periodic boundary conditions; this choice makes the momentum discrete, and the momentum states renormalizable. If the cubic lattice is $N\times N\times N$ then the normalization factor is $\mathcal{N}=N^{3/2}$. If necessary we can then go to a limit where $N\to \infty$. In the finite case, the momentum values $k_{X,Y,Z}$ are integer multiples of $2\pi /N\Delta x$.

The states $|k_X,k_Y{k}_Z\rangle$ are eigenstates of the shift operators $S_{X,Y,Z}$ with eigenvalues $e^{i{k}_{X,Y,Z}\Delta x}$. In this representation, the evolution operator (with $Q={\gamma }_0$ for fermions and $Q=0$ for photons) is

$$\begin{array}{ccc}\hfill U& =& e^{-i\theta Q}{e}^{i{K}_X\Delta x\Delta P_X}{e}^{i{K}_Y\Delta x\Delta P_Y}{e}^{i{K}_Z\Delta x\Delta P_Z}\hfill \\ & =& \sum _{\mathbf{k}}|\mathbf{k}\rangle \langle \mathbf{k}|\otimes U_{\mathbf{k}}\equiv \sum _{\mathbf{k}}|\mathbf{k}\rangle \langle \mathbf{k}|\otimes V_{\mathbf{k}}{\mathsf{\Lambda }}_{\mathbf{k}}{V}_{\mathbf{k}}^{†},\hfill \end{array}$$

(16)

where the $K_{X,Y,Z}$ are the operators corresponding to the components of the momentum and the $\left\{{\mathsf{\Lambda }}_{\mathbf{k}}\right\}$ are the diagonalized forms of the $\left\{U_{\mathbf{k}}\right\}$. Using the choices above for $\Delta P_{X,Y,Z}$, we have

$$U_{\mathbf{k}}=e^{-i\theta {\gamma }_0}{e}^{i{k}_X\Delta x{\gamma }_0{\gamma }_1}{e}^{i{k}_Y\Delta x{\gamma }_0{\gamma }_2}{e}^{i{k}_Z\Delta x{\gamma }_0{\gamma }_3},$$

(17)

where for fermions these are the usual gamma matrices, and for bosons the gamma matrices are the boson versions defined, and the first term is absent, i.e., $\theta =0$. (As we will see, this corresponds to the massless case.) In this representation we see that the Hilbert space decomposes into D-dimensional blocks ($D=4$ for fermions, $D=6$ for bosons) labeled by the momentum $\mathbf{k}$. The energy eigenstates $|\mathbf{k},{\lambda }_{\mathbf{k},j}\rangle$ of the quantum walk can be found by diagonalizing the $D\times D$ matrices $U_{\mathbf{k}}=V_{\mathbf{k}}{\mathsf{\Lambda }}_{\mathbf{k}}{V}_{\mathbf{k}}^{†}$. Since $U_{\mathbf{k}}$ is unitary, the matrices ${\mathsf{\Lambda }}_{\mathbf{k}}$ have eigenvalues of the form ${\lambda }_{\mathbf{k},j}=e^{i{\varphi }_{\mathbf{k},j}}$.

For fermions, if we go to the long-wavelength/low-mass limit $|\mathbf{k}\Delta x|\ll 1$ and $\left|\theta \right|\ll 1$, we get

$$U_{\mathbf{k}}\approx 1+i\Delta x\left(k_X\Delta P_X+k_Y\Delta P_Y+k_Z\Delta P_Z\right)-i\theta Q.$$

(18)

For low-momentum particles with mass, we can take the interval $\Delta t$ to be small and define an approximate time derivative operator:

$${\partial }_t\equiv {\displaystyle \frac{U-1}{\Delta t}}.$$

(19)

In terms of this we can get a Dirac equation, which holds in the long-wavelength/low-mass limit:

$${\partial }_t|\psi \rangle =i{\displaystyle \frac{\Delta x}{\Delta t}}\left(k_X\Delta P_X+k_Y\Delta P_Y+k_Z\Delta P_Z\right)|\psi \rangle -i{\displaystyle \frac{\theta }{\Delta t}}Q|\psi \rangle .$$

(20)

We can relate the physical space interval $\Delta x$ and the physical time interval $\Delta t$ in terms of the speed of light c: $c\equiv \Delta x/\Delta t$. Then, defining the particle mass $\theta /\Delta t\equiv m_e{c}^2$, and using $\Delta P_i={\gamma }_0{\gamma }_i$, we obtain the momentum-space representation of the Dirac equation:

$${\gamma }_0{\partial }_t|\psi \rangle =ic\left(k_X{\gamma }_1+k_Y{\gamma }_2+k_Z{\gamma }_3-m_{e}c{\gamma }_0\right)|\psi \rangle .$$

(21)

For bosons, we get the analogous equation, except with $\Delta P_i={\gamma }_0^{\mathrm{bose}}{\gamma }_i^{\mathrm{bose}}$. This gives, in the same small $\mathbf{k}$ limit,

$$U_{\mathbf{k}}\approx 1+i\Delta x\left(k_X{\gamma }_0^{\mathrm{bose}}{\gamma }_1^{\mathrm{bose}}+k_Y{\gamma }_0^{\mathrm{bose}}{\gamma }_2^{\mathrm{bose}}+k_Z{\gamma }_0^{\mathrm{bose}}{\gamma }_3^{\mathrm{bose}}\right).$$

(22)

We thus get

$$\begin{array}{ccc}\hfill {\partial }_t& \equiv & {\displaystyle \frac{U-1}{\Delta t}}\hfill \\ & =& ic\left(k_X{\gamma }_0^{\mathrm{bose}}{\gamma }_1^{\mathrm{bose}}+k_Y{\gamma }_0^{\mathrm{bose}}{\gamma }_2^{\mathrm{bose}}+k_Z{\gamma }_0^{\mathrm{bose}}{\gamma }_3^{\mathrm{bose}}\right)\hfill \\ & =& ic{\gamma }_0^{\mathrm{bose}}\left(\mathbf{k}·{\gamma }^{\mathrm{bose}}\right)=ic{\sigma }_Z\otimes \left(\mathbf{k}·\mathbf{J}\right).\hfill \end{array}$$

(23)

We can define an effective Hamiltonian

$$H\equiv c{\sigma }_Z\otimes \left(\mathbf{k}·\mathbf{J}\right)=ic{\sigma }_Z\otimes \left(\begin{array}{ccc}0& -k_Z& k_Y\\ k_Z& 0& -k_X\\ -k_Y& k_X& 0\end{array}\right).$$

(24)

Using the identity

$$\left(\mathbf{k}·\mathbf{J}\right)\mathbf{v}=i\mathbf{k}\times \mathbf{v},$$

(25)

we see that the eigenvectors of $\left(\mathbf{k}·\mathbf{J}\right)$ are

$${\mathbf{v}}_0={\displaystyle \frac{1}{k}}\left(\begin{array}{c}{k}_X\\ k_Y\\ k_Z\end{array}\right),{\mathbf{v}}_{\pm }={\displaystyle \frac{1}{\sqrt{2k^2\left(k_X^2+k_Z^2\right)}}}\left(\begin{array}{c}-ik{k}_Z\pm k_X{k}_Y\\ \mp \left(k_X^2+k_Z^2\right)\\ ik{k}_X\pm k_Y{k}_Z\end{array}\right),$$

(26)

with eigenvalues 0 and $\pm ik$, respectively, where $k=\left|\mathbf{k}\right|=\sqrt{k_X^2+k_Y^2+k_Z^2}$. These correspond to a longitudinal mode, and right- and left-polarized modes, respectively. The Hamiltonian thus has two zero-energy longitudinal eigenstates, and a positive and negative energy eigenstate for each circular polarization. In the free theory these do not mix, so if we consider positive-energy initial states, the dynamics do not take us out of the positive-energy subspace. That is not generally true in the presence of interactions, which we will discuss later.

When we promote the single-particle theory to allow for the creation and annihilation of photons, we’ll find that the field operators, expressed in the usual manner in terms of the creation and annihilation operators and the above eigenvectors, will satisfy the Maxwell equations in the long wavelength limit. To see that here, for the one-particle sector, we have to consider the photon wave function (see Refs. [30,31]):

$$\mathsf{\Psi }=\left(\begin{array}{c}{\mathsf{\Psi }}_{+}\\ {\mathsf{\Psi }}_{-}\end{array}\right),$$

(27)

where ${\mathsf{\Psi }}_{+}=(\mathbf{E}+i\mathbf{B})$, and ${\mathsf{\Psi }}_{-}={{\mathsf{\Psi }}_{+}}^{*}$. Using the Schrödinger equation

$$i{\partial }_t\mathsf{\Psi }=H\mathsf{\Psi },$$

(28)

and the vector identity in Equation (25), we get

$$i{\partial }_t\left(\begin{array}{c}{\mathsf{\Psi }}_{+}\\ {\mathsf{\Psi }}_{-}\end{array}\right)=\left(\begin{array}{c}c\mathbf{k}\times {\mathsf{\Psi }}_{+}\\ -c\mathbf{k}\times {\mathsf{\Psi }}_{-}\end{array}\right).$$

(29)

Substituting $-i\nabla$ for $c\mathbf{k}$ and considering the real and imaginary parts separately, we get the two time-dependent Maxwell equations:

$$\begin{array}{ccc}\hfill {\partial }_t\mathbf{E}& =& \nabla \times \mathbf{B},\hfill \\ \hfill {\partial }_t\mathbf{B}& =& -\nabla \times \mathbf{E}.\hfill \end{array}$$

(30)

The other two Maxwell equations (in the absence of charge and current) arise from the transversality of the positive energy eigenvectors.

5. Promotion to Quantum Cellular Automata

5.1. Fermions

As was shown in [21], it is possible to promote the single-particle fermionic theory to a multi-particle QCA that reduces to the usual Dirac field theory in the long-wavelength limit. To accomplish this, one begins with a maximum number of particles, $N_{\max }$, and then let $N_{\max }$ get arbitrarily large. In other words, one begins with the space

$${\mathcal{H}}_{\mathrm{total}}={\mathcal{H}}^{\left(1\right)}\otimes \cdots \otimes {\mathcal{H}}^{\left(N_{\max }\right)},$$

(31)

where ${\mathcal{H}}^{\left(j\right)}={\mathcal{H}}_{\mathrm{QW}}\oplus \mathrm{span}\left(|\omega \rangle \right)$ is the Hilbert space for particle type j, which contains either one quantum walk particle or no particle (the vacuum state $|\omega \rangle$). The system then evolves by the unitary

$$U_{\mathrm{total}}=U^{\left(1\right)}\otimes \cdots \otimes U^{\left(N_{\max }\right)},$$

(32)

where

$$U^{\left(j\right)}=U_{\mathrm{QW}}\oplus |\omega \rangle \langle \omega |.$$

(33)

We then proceed with essentially the usual Fock space construction, restricting ourselves to the totally antisymmetric subspace of ${\mathcal{H}}_{\mathrm{total}}$. We call this the physical subspace, which we define as

$${\mathcal{H}}_{\mathrm{phys}}={\mathcal{H}}_0\oplus {\mathcal{A}}_1\oplus \cdots \oplus {\mathcal{A}}_{N_{\max }}.$$

(34)

The subspace ${\mathcal{A}}_n$ is the totally antisymmetric subspace with n particles and ${\mathcal{H}}_0$ is the one-dimensional vacuum space $\mathrm{span}\left(|\mathsf{\Omega }\rangle \right)$, where $|\mathsf{\Omega }\rangle ={|\omega \rangle }^{\otimes N_{\max }}$. We can define basis vectors for the subspaces ${\mathcal{A}}_n$ in terms of the “energy eigenstates” of the quantum walk:

$$\begin{array}{ccc}\hfill |\mathsf{\Omega }\rangle & =& {|\omega \rangle }^{\otimes N_{\max }},\hfill \\ \hfill |\mathbf{k},j\rangle & =& |\mathbf{k},{\lambda }_{\mathbf{k},j}\rangle \otimes {|\omega \rangle }^{\otimes N_{\max }-1},\hfill \\ \hfill |{\mathbf{k}}_1,j_1;{\mathbf{k}}_2,j_2\rangle & =& {\displaystyle \frac{1}{\sqrt{2}}}(|{\mathbf{k}}_1,{\lambda }_{{\mathbf{k}}_1,j_1}\rangle \otimes |{\mathbf{k}}_2,{\lambda }_{{\mathbf{k}}_2,j_2}\rangle \hfill \\ & & -|{\mathbf{k}}_2,{\lambda }_{{\mathbf{k}}_2,j_2}\rangle \otimes |{\mathbf{k}}_1,{\lambda }_{{\mathbf{k}}_1,j_1}\rangle )\otimes {|\omega \rangle }^{\otimes N_{\max }-1},\hfill \end{array}$$

(35)

and so forth, defining basis states $|{\mathbf{k}}_1,j_1;\dots ;{\mathbf{k}}_n,j_n\rangle$ for particle numbers n up to $N_{\max }$. The state $|{\mathbf{k}}_1,j_1;\dots ;{\mathbf{k}}_n,j_n\rangle$ is an eigenstate of the evolution operator $U_{\mathrm{total}}$ with eigenvalue

$$\lambda =exp\left(i\sum _{\ell =1}^n{\varphi }_{{\mathbf{k}}_{\ell },j_{\ell }}\right).$$

The n-particle antisymmetrized subspace is thus preserved by the time evolution, and hence so is the full physical subspace ${\mathcal{H}}_{\mathrm{phys}}$: an antisymmetric initial state will evolve to be antisymmetric at all later times. Given this basis, we can define a set of creation and annihilation operators that transform the energy basis vectors into each other:

$$|{\mathbf{k}}_1,j_1;\dots ;{\mathbf{k}}_n,j_n\rangle \equiv a_{{\mathbf{k}}_1,j_1}^{†}\cdots a_{{\mathbf{k}}_n,j_n}^{†}|\mathsf{\Omega }\rangle ,$$

(36)

and that obey the usual anticommutation relations. To avoid redundant counting of the basis vectors, we should establish a standard ordering for the states ${\mathbf{k}}_i,j_i$. The choice of this ordering is unimportant, so long as it is unambiguous.

We can write the action of the local evolution U in terms of these operators for any $|\mathsf{\Psi }\rangle \in {\mathcal{H}}_{\mathrm{phys}}$:

$$U_{\mathrm{total}}|\mathsf{\Psi }\rangle =exp\left(i\sum _{\mathbf{k},j}{\varphi }_{\mathbf{k},j}{a}_{\mathbf{k},j}^{†}{a}_{\mathbf{k},j}\right)|\mathsf{\Psi }\rangle .$$

(37)

A very important observation should be made here. These creation and annihilation operators $a_{\mathbf{k},j}^{†}$ and $a_{\mathbf{k},j}$ are clearly nonlocal operators; so the unitary operator on the right-hand side of Equation (37) is nonlocal. But the unitary operator $U_{\mathrm{total}}$ is itself purely local. It is a collection of local quantum walk unitaries, and can easily be embedded in a quantum cellular automaton. These two sides are equal only on states in the physical subspace. Indeed, the right-hand side has not been defined for states outside that subspace.

Having defined creation and annihilation operators for the energy eigenstates, we can then define creation and annihilation operators ${\tilde{a}}_{\mathbf{k},j}^{†}$ and ${\tilde{a}}_{\mathbf{k},j}$ in the momentum representation. We can write these as operator-valued column vectors ${\tilde{\mathbf{a}}}_{\mathbf{k},j}^{†}$ and ${\tilde{\mathbf{a}}}_{\mathbf{k},j}$, where

$${\tilde{\mathbf{a}}}_{\mathbf{k},j}^{†}=\left(\begin{array}{c}{\tilde{a}}_{\mathbf{k},1}^{†}\\ ⋮\\ {\tilde{a}}_{\mathbf{k},D}^{†}\end{array}\right)=V_{\mathbf{k}}^{†}\left(\begin{array}{c}{a}_{\mathbf{k},1}^{†}\\ ⋮\\ a_{\mathbf{k},D}^{†}\end{array}\right),$$

(38)

using the matrices $V_{\mathbf{k}}^{†}$ from Equation (16). These vectors have an effective time evolution

$$U_{\mathrm{total}}{\tilde{\mathbf{a}}}_{\mathbf{k},j}^{†}{U}_{\mathrm{total}}^{†}=U_{\mathbf{k}}{\tilde{\mathbf{a}}}_{\mathbf{k},j}^{†},U_{\mathrm{total}}{\tilde{\mathbf{a}}}_{\mathbf{k},j}{U}_{\mathrm{total}}^{†}=U_{\mathbf{k}}^{*}{\tilde{\mathbf{a}}}_{\mathbf{k},j},$$

(39)

where $U_{\mathbf{k}}$ is given by Equation (16).

What about the maximum number of particles $N_{\max }$? This may seem like an unphysical feature of the theory. In the free theory, particle number is conserved, so a solution will generally be confined to a single antisymmetric subspace ${\mathcal{A}}_n$. However, in the interacting case this conservation need not hold.

One way to handle this is to start with a vacuum corresponding to a Dirac sea, in which all the negative-energy states are occupied while all positive-energy states are vacant. Interactions do not create or destroy particles, but may move a particle from a negative- to a positive-energy state, leaving a hole (corresponding to particle-antiparticle pair creation), or move a positive-energy state to a negative-energy hole (corresponding to particle-antiparticle annihilation). In such a construction we would have $N_{\max }=2N^3$, where N is the lattice size. This means that the Hilbert spaces of the local subsystems at each lattice site are enormously high-dimensional (albeit finite), but the states of the theory remain restricted to the much lower-dimensional physical subspace. In the full space the dynamics of the free theory are fully local; the effective dynamics in the physical subspace obey Fermi statistics.

5.2. Bosons

We can make a similar construction for free bosons, with $N_{\max }$ copies of a bosonic QW, representing distinguishable particles, and then confining ourselves to the completely symmetric subspace. So long as the total number of bosonic particles N is less that $N_{\max }$, this system will behave like a bosonic field theory with N excitations, and we can introduce creation and annihilation operators and rewrite the evolution unitary in terms of them, just as we did in the fermionic case.

This construction is sufficient for the free theory, where particle number is conserved. But it is unsatisfactory for an interacting theory, in which bosons may be created or destroyed. This issue is less critical for fermions, where no more than one particle may occupy any given state and we can introduce a Dirac Sea as the vacuum state. There is no analogous construction for bosons, since any number of bosons may occupy any state.

However, the restrictions on bosonic QCAs are also much weaker than for fermions: since bosonic creation and annihilation operators for different modes commute, these operators can be made local, and there is not a “no-go” theorem as there is for fermionic QCAs [18]. It is therefore possible to construct bosonic QCAs in a different way, which avoids these locality issues, by using excitations of local modes as our particles. We will lay out this construction of a free bosonic QCA theory in this subsection. The key idea is quite simple: each possible basis state (comprising a position and an internal state $|\mathbf{x},\epsilon \rangle$) of the single particle in the quantum walk becomes a harmonic oscillator mode in the QCA. The unitary evolution operators will then trade excitations between neighboring modes.

Recall that the evolution unitary for the massless bosonic quantum walk was a product of terms of the form

$$\left(S_j{P}_j^{+}+P_j^0+S_j^{†}{P}_j^{-}\right),$$

(40)

where $j=X,Y,Z$. As discussed in Section 4 we can recover a Maxwell field if we use a doubled space, where each of the projectors is two-dimensional, for an internal space dimension of 6. To simplify the derivation of the corresponding QCA, we will first consider a non-doubled space, with internal dimension 3, so we can write the projectors as outer products of the eigenstates of the $J_j$ matrix:

$$P_j^{+}={|+\rangle }_j\langle +|,P_j^{-}={|-\rangle }_j\langle -|,P_j^0={|0\rangle }_j\langle 0|.$$

Once we have derived the form of the QCA in this case, we can double the internal space, which is straightforward.

Instead of using shift operators, we can rewrite this unitary factor as a product of swaps between states at neighboring sites. Let $\widehat{\mathbf{j}}$ be the unit vector in the $j=X,Y,Z$ direction. The term in the evolution unitary can be written as

$$\begin{array}{ccc}\hfill \left(S_j{P}_j^{+}+P_j^0+S_j^{†}{P}_j^{-}\right)& =& \sum _{\mathbf{x}}(|\mathbf{x}+\Delta x\widehat{\mathbf{j}}\rangle \langle \mathbf{x}|\otimes {|+\rangle }_j\langle +|\hfill \\ & & +|\mathbf{x}-\Delta x\widehat{\mathbf{j}}\rangle \langle \mathbf{x}|\otimes {|-\rangle }_j\langle -|+|\mathbf{x}\rangle \langle \mathbf{x}|\otimes {|0\rangle }_j\langle 0|)\hfill \\ & =& \sum _{\mathbf{x}}|\mathbf{x}\rangle \langle \mathbf{x}|\otimes \left({|+\rangle }_j\langle -|+{|-\rangle }_j\langle +|+{|0\rangle }_j\langle 0|\right)\hfill \\ & & \times \sum _{\mathbf{x}}({|\mathbf{x}+\Delta x\widehat{\mathbf{j}},-\rangle }_j\langle \mathbf{x},+|\hfill \\ & & +{|\mathbf{x},+\rangle }_j\langle \mathbf{x}+\Delta x\widehat{\mathbf{j}},-|+{|\mathbf{x},0\rangle }_j\langle \mathbf{x},0|).\hfill \end{array}$$

(41)

In other words, the shifts are done in two steps: the + and − states at neighboring sites are interchanged, and then at each site + and − are switched back. (The 0 states remain unchanged throughout.)

The advantage of breaking this down into two steps is that each step is a purely local operation that only affects two neighboring sites. This makes it very easy to turn into a QCA that is definitely local. In constructing the QCA, each basis state $|\mathbf{x},\epsilon \rangle$ of the QW corresponds to a harmonic oscillator mode with creation and annihilation operators $a_{\mathbf{x},\epsilon }^{†},a_{\mathbf{x},\epsilon }$. Then the evolution unitary becomes

$$\begin{array}{ccc}\hfill {\widehat{U}}_j& =& \left(\underset{\mathbf{x}}{⨂}exp\left\{-{\displaystyle \frac{i\pi }{2}}\left(a_{\mathbf{x},+}^{†}{a}_{\mathbf{x},-}+a_{\mathbf{x},-}^{†}{a}_{\mathbf{x},+}\right)\right\}\right)\hfill \\ & & \times \left(\underset{\mathbf{x}}{⨂}exp\left\{-{\displaystyle \frac{i\pi }{2}}\left(a_{\mathbf{x}+\Delta x\widehat{\mathbf{j}},-}^{†}{a}_{\mathbf{x},+}+a_{\mathbf{x},+}^{†}{a}_{\mathbf{x}+\Delta x\widehat{\mathbf{j}},-}\right)\right\}\right),\hfill \\ & =& exp\left\{-{\displaystyle \frac{i\pi }{2}}\sum _{\mathbf{x}}\left(a_{\mathbf{x},+}^{†}{a}_{\mathbf{x},-}+a_{\mathbf{x},-}^{†}{a}_{\mathbf{x},+}\right)\right\}\hfill \\ & & \times exp\left\{-{\displaystyle \frac{i\pi }{2}}\sum _{\mathbf{x}}\left(a_{\mathbf{x}+\Delta x\widehat{\mathbf{j}},-}^{†}{a}_{\mathbf{x},+}+a_{\mathbf{x},+}^{†}{a}_{\mathbf{x}+\Delta x\widehat{\mathbf{j}},-}\right)\right\}.\hfill \end{array}$$

(42)

We can see that ${\widehat{U}}_j$ is a product of tensor products of unitaries, which act on the local subsystems $\mathbf{x}$ or on non-overlapping sets of modes at neighboring subsets $\mathbf{x}$ and $\mathbf{x}+\Delta x\widehat{\mathbf{j}}$. So ${\widehat{U}}_j$ is local, and defines a QCA. Since operators on different modes commute, we can rewrite each tensor product of exponentials as the exponential of a sum of terms (where now we understand the creation and annihilation operators as acting as the identity on other modes). Further, we see that the defined QCA preserves particle number (or the number of excitations); and on the single-particle sector, it exactly reproduces a quantum walk.

To reproduce the full QCA, of course, we need not just a single operator ${\widehat{U}}_j$, but the product of operators producing shifts in all three directions:

$${\widehat{U}}_{\mathrm{QCA}}={\widehat{U}}_X{\widehat{U}}_Y{\widehat{U}}_Z.$$

(43)

Each of the operators ${\widehat{U}}_{X,Y,Z}$ has the same form as in Equation (42). In the quantum walk, we get nontrivial evolution because the projectors $P_{X,Y,Z}^{\pm }$ for the shifts along each axis project onto subspaces that are non-orthogonal to those of the other axes. For the QCA, this means that the creation and annihilation operators in Equation (42) for one axis j are linear combinations of those for a different axis $j^{\prime }$.

Since all the expressions inside the exponentials are quadratic in creation and annihilation operators, we can rewrite Equation (42) in terms of a standard set of modes with operators $a_{\mathbf{x},\ell }^{†}$ and $a_{\mathbf{x},\ell }$:

$$\begin{array}{ccc}\hfill {\widehat{U}}_j& =& exp\left\{-i\sum _{\mathbf{x},\ell ,{\ell }^{\prime }}(m_{\ell {\ell }^{\prime }}^j+{\left(m_{{\ell }^{\prime }\ell }^j\right)}^{*})\left(a_{\mathbf{x},\ell }^{†}{a}_{\mathbf{x},{\ell }^{\prime }}\right)\right\}\hfill \\ & & \times exp\left\{-i\sum _{\mathbf{x},\ell ,{\ell }^{\prime }}\left(m_{\ell {\ell }^{\prime }}^j{a}_{\mathbf{x}+\Delta x\widehat{\mathbf{j}},\ell }^{†}{a}_{\mathbf{x},{\ell }^{\prime }}+{\left(m_{\ell {\ell }^{\prime }}^j\right)}^{*}{a}_{\mathbf{x},{\ell }^{\prime }}^{†}{a}_{\mathbf{x}+\Delta x\widehat{\mathbf{j}},\ell }\right)\right\},\hfill \end{array}$$

(44)

where $\left[m_{\ell {\ell }^{\prime }}^j\right]$ are the matrix elements of the matrix $(\pi /2){|-\rangle }_j\langle +|$ in the standard basis.

One can go to the momentum picture by introducing a new set of modes corresponding to different momentum states:

$$a_{\mathbf{x},\ell }={\displaystyle \frac{1}{N^{3/2}}}\sum _{\mathbf{k}}{e}^{-i\mathbf{x}·\mathbf{k}}{a}_{\mathbf{k},\ell },a_{\mathbf{x},\ell }^{†}={\displaystyle \frac{1}{N^{3/2}}}\sum _{\mathbf{k}}{e}^{i\mathbf{x}·\mathbf{k}}{a}_{\mathbf{k},\ell }^{†}.$$

(45)

It is easy to check that creation and annihilation operators for different values of $\mathbf{k}$ commute with each other. If we substitute this definition into Equation (44) the unitary decomposes into a product of factors acting on different momentum modes:

$$\begin{array}{ccc}\hfill {\widehat{U}}_j& =& \prod _{\mathbf{k}}exp\left\{-i\sum _{\ell ,{\ell }^{\prime }}(m_{\ell {\ell }^{\prime }}^j+{\left(m_{{\ell }^{\prime }\ell }^j\right)}^{*})\left(a_{\mathbf{k},\ell }^{†}{a}_{\mathbf{k},{\ell }^{\prime }}\right)\right\}\hfill \\ & & \times exp\left\{-i\sum _{\ell ,{\ell }^{\prime }}\left(e^{i{k}_j\Delta x}{m}_{\ell {\ell }^{\prime }}^j{a}_{\mathbf{k},\ell }^{†}{a}_{\mathbf{k},{\ell }^{\prime }}+e^{-i{k}_j\Delta x}{\left(m_{\ell {\ell }^{\prime }}^j\right)}^{*}{a}_{\mathbf{k},{\ell }^{\prime }}^{†}{a}_{\mathbf{k},\ell }\right)\right\},\hfill \end{array}$$

(46)

The coupling between neighboring sites in the position picture becomes a phase of $e^{i{k}_j\Delta x}$ in the momentum picture.

A product of unitaries produced by quadratic Hamiltonians is equivalent to a single unitary produced by a quadratic Hamiltonian. That is, given a set of field modes with annihilation operators $a_1,\dots ,a_n$ and two $n\times n$ Hermitian matrices A and B, then there exists an $n\times n$ Hermitian matrix C such that

$$exp\left\{-i{\mathbf{a}}^{†}·C\mathbf{a}\right\}=exp\left\{-i{\mathbf{a}}^{†}·B\mathbf{a}\right\}\times exp\left\{-i{\mathbf{a}}^{†}·A\mathbf{a}\right\},$$

(47)

where

$$\mathbf{a}=\left(\begin{array}{c}{a}_1\\ ⋮\\ a_n\end{array}\right),{\mathbf{a}}^{†}=\left(\begin{array}{c}{a}_1^{†}\\ ⋮\\ a_n^{†}\end{array}\right).$$

(48)

One can calculate the matrix C from the matrices A and B by solving the equation

$$exp\left\{i{C}^T\right\}=exp\left\{i{B}^T\right\}exp\left\{i{A}^T\right\}.$$

(49)

This is most readily done by diagonalizing A and B, exponentiating them, and then multiplying them. One can find C by diagonalizing the resulting matrix and taking its log. In the case of the unitary factor $U_j$, the matrices $A_j$ and $B_j$ have elements

$$\begin{array}{ccc}\hfill a_{\ell {\ell }^{\prime }}& =& e^{i{k}_j\Delta x}{m}_{\ell {\ell }^{\prime }}^j+e^{-i{k}_j\Delta x}{\left(m_{{\ell }^{\prime }\ell }^j\right)}^{*},\hfill \\ \hfill b_{\ell {\ell }^{\prime }}& =& m_{\ell {\ell }^{\prime }}^j+{\left(m_{{\ell }^{\prime }\ell }^j\right)}^{*}.\hfill \end{array}$$

(50)

The result of this series of three unitaries is a quadratic evolution with a matrix C, where

$$exp\left\{i{C}^T\right\}=exp\left\{i{B}_x^T\right\}exp\left\{i{A}_x^T\right\}exp\left\{i{B}_y^T\right\}exp\left\{i{A}_y^T\right\}exp\left\{i{B}_z^T\right\}exp\left\{i{A}_z^T\right\}.$$

(51)

These matrices are $3\times 3$, so in general it is not difficult to diagonalize them. For a bosonic QW, we suggested in Sec. 3 that the basis vectors should be eigenvectors of the three spin-1 J matrices given in Equation (12). This choice for the bosonic QCA would lead to three pairs of A and B matrices:

$$\begin{array}{ccc}\hfill e^{i{A}_x^T}=\left(\begin{array}{ccc}1& 0& 0\\ 0& -icos(k_X\Delta x)& isin(k_X\Delta x)\\ 0& isin(k_X\Delta x)& icos(k_X\Delta x)\end{array}\right),& & e^{i{B}_x^T}=\left(\begin{array}{ccc}1& 0& 0\\ 0& -i& 0\\ 0& 0& i\end{array}\right),\hfill \\ \hfill e^{i{A}_y^T}=\left(\begin{array}{ccc}icos(k_Y\Delta x)& 0& -isin(k_Y\Delta x)\\ 0& 1& 0\\ -isin(k_Y\Delta x)& 0& -icos(k_Y\Delta x)\end{array}\right),& & e^{i{B}_y^T}=\left(\begin{array}{ccc}i& 0& 0\\ 0& 1& 0\\ 0& 0& -i\end{array}\right),\hfill \\ \hfill e^{i{A}_z^T}=\left(\begin{array}{ccc}-icos(k_Z\Delta x)& isin(k_Z\Delta x)& 0\\ isin(k_Z\Delta x)& icos(k_Z\Delta x)& 0\\ 0& 0& 1\end{array}\right),& & e^{i{B}_z^T}=\left(\begin{array}{ccc}-i& 0& 0\\ 0& i& 0\\ 0& 0& 1\end{array}\right).\hfill \end{array}$$

(52)

Multiplying out these matrices gives us

$$exp\left\{i{C}^T\right\}=\left(\begin{array}{ccc}{c}_y{c}_z& -c_y{s}_z& s_y\\ c_x{s}_z+s_x{s}_y{c}_z& c_x{c}_z-s_x{s}_y{s}_z& -s_x{c}_y\\ s_x{s}_z-c_x{s}_y{c}_z& s_x{c}_z+c_x{s}_y{s}_z& c_x{c}_y\end{array}\right),$$

(53)

where we have used the compact notation $c_{x,y,z}=cos(k_{x,y,z}\Delta x)$ and $s_{x,y,z}=sin(k_{x,y,z}\Delta x)$. This matrix can, in principle, be diagonalized in closed form, but the expression is dauntingly complicated. Its eigenvalues can also be calculated, though they are also somewhat complicated:

$$\begin{array}{ccc}\hfill {\lambda }_0& =& 1=e^0,\hfill \\ \hfill {\lambda }_{+}& =& {\displaystyle \frac{1}{2}}\left(c_x{c}_y+c_y{c}_z+c_z{c}_x-s_x{s}_y{s}_z-1\right)\hfill \\ & & -i\sqrt{1-{\displaystyle \frac{1}{4}}{\left(c_x{c}_y+c_y{c}_z+c_z{c}_x-s_x{s}_y{s}_z-1\right)}^2}\equiv e^{-i{\varphi }_{+}},\hfill \\ \hfill {\lambda }_{-}& =& {\displaystyle \frac{1}{2}}\left(c_x{c}_y+c_y{c}_z+c_z{c}_x-s_x{s}_y{s}_z-1\right)\hfill \\ & & +i\sqrt{1-{\displaystyle \frac{1}{4}}{\left(c_x{c}_y+c_y{c}_z+c_z{c}_x-s_x{s}_y{s}_z-1\right)}^2}\equiv e^{-i{\varphi }_{-}}.\hfill \end{array}$$

(54)

Simple observation shows that the + and − eigenvalues have opposite phases: ${\varphi }_{-}=-{\varphi }_{+}$. These can be thought of as negative and positive energy modes. We can recover more familiar expressions by going to the long-wavelength limit, where $|\mathbf{k}\Delta x|\ll 1$. In this limit

$$exp\left\{i{C}^T\right\}\approx I-i(k_X{J}_X+k_Y{J}_Y+k_Z{J}_Z)\Delta x\approx exp\{-i\mathbf{k}·\mathbf{J}\Delta x\}.$$

(55)

In this limit, the phases corresponding to the eigenvalues are approximately ${\varphi }_0=0$, ${\varphi }_{+}\approx k\Delta x=\Delta x\sqrt{k_X^2+k_Y^2+k_Z^2}$, ${\varphi }_{-}\approx -k\Delta x$. These correspond to a zero-energy, positive-energy, and negative-energy mode.

As discussed in Section 4, we can recover the Maxwell field theory by doubling the internal space. In place of the spin matrices $J_j$ we use the $6\times 6$ matrices ${\sigma }_Z\otimes J_j$. The derivation above is essentially doubled, giving two zero-energy, two negative-energy and two positive-energy modes for each allowed momentum $\mathbf{k}$. If we confine ourselves to the positive-energy subspace, these two modes correspond to the two allowed helicity states of the electromagnetic field. The two negative-energy modes are like retarded versions of these two modes. The zero-energy modes are static and unaffected by the other modes.

We have seen that the free theories of both fermionic and bosonic QCAs include matching positive and negative energy modes. (And zero-energy modes as well, in the bosonic case.) This turns out to be a general property of QCAs based on their overall symmetry, as discussed in the next section.

6. Symmetries of QCAs

Symmetries of QCAs—especially translation and parity symmetry—have been discussed extensively elsewhere [6,7,8], so here we will just briefly summarize their effects on the types of solutions that exist for the systems we consider in this paper.

6.1. Translation Symmetry

Translation symmetry implies that the evolution operator is invariant under shifts of the spatial lattice. In a 3D cubic lattice of size N with periodic boundaries, we can define shift operators $S_{X,Y,Z}$ that move the states of the local subsystems to the site $\Delta x$ away in the positive X, Y or Z direction, respectively, wrapping around at the periodic boundaries. $S_{X,Y,Z}^{†}$ obviously represents a shift in the negative direction. Translation symmetry was extensively studied in the early QCA literature [6,7,8]. On the cubic lattice these shifts all commute.

If S is a product of any number of these unit shifts $S_{X,Y,Z}$ or $S_{X,Y,Z}^{†}$ then translation symmetry of the QCA evolution unitary $\widehat{U}$ implies

$$S\widehat{U}{S}^{†}=\widehat{U}\Rightarrow S\widehat{U}=\widehat{U}S\Rightarrow [\widehat{U},S]=0.$$

(56)

Since $\widehat{U}$ and S commute, we can find a basis of eigenstates of $\widehat{U}$ that are eigenstates S. This is what allows the momentum picture that we used both for QWs and QCAs in this paper.

6.2. Rotation

By contrast, the discrete lattice of the QCA means that it does not have rotation symmetry: waves propagating along the axes may behave differently than waves between axes. However, QCAs derived using the equal norm condition described in Section 3 recover rotation symmetry (and indeed, Lorentz invariance) in the long-wavelength limit [14]. So these continuous symmetries may arise effectively at long wavelengths.

6.3. Parity

Parity symmetry implies that the evolution unitary is invariant under spatial inversion, $\mathbf{x}=(x,y,z)\to -\mathbf{x}=(-x,-y,-z)$. For this to be possible, we must define the parity operator $\mathcal{P}$ to also act appropriately on the internal space. This is most straightforward to see in the QW case [14]. If we have $\mathcal{P}U\mathcal{P}=U$, then

$$\begin{array}{ccc}\hfill U_{QW}& =& e^{i\theta Q}\prod _{j=X,Y,Z}\left(S_j{P}_j^{+}+S_j^{†}{P}_j^{-}\right)\hfill \\ & =& \mathcal{P}{e}^{i\theta Q}\prod _{j=X,Y,Z}\left(S_j{P}_j^{+}+S_j^{†}{P}_j^{-}\right)\mathcal{P}\hfill \\ & =& e^{i\theta \mathcal{P}Q\mathcal{P}}\prod _{j=X,Y,Z}\left(\mathcal{P}{S}_j\mathcal{P}\mathcal{P}{P}_j^{+}\mathcal{P}+\mathcal{P}{S}_j^{†}\mathcal{P}\mathcal{P}{P}_j^{-}\mathcal{P}\right).\hfill \end{array}$$

(57)

Spatial inversion means that $\mathcal{P}{S}_j\mathcal{P}=S_j^{†}$, which implies that

$$\mathcal{P}Q\mathcal{P}=Q,\mathcal{P}{P}_j^{\pm }\mathcal{P}=P^{\mp }.$$

(58)

The simplest way of satisfying these conditions is if $\mathcal{P}$ acts on the internal space like Q, and if Q anticommutes with the operators $\Delta P_j=P_j^{+}-P_j^{-}$. This anticommutation plays an important role in the emergence of the Dirac equation in the long-wavelength limit of QWs [14,18,19,21].

6.4. Time-Reversal Symmetry and Negative Energies

Time-reversal symmetry captures the idea that time evolution looks the same forward and backward. Classically, this would be illustrated by reversing the momenta of all particles while leaving their positions unchanged; this is like “running the movie backwards.” In quantum mechanics we also need to flip the spins of particles. This symmetry is very important to the tension between locality of interactions, one the one hand, and positive energy on the other, as we shall see in Section 7.

We can write this symmetry in terms of a time-reversal operation $\mathcal{T}$:

$$\widehat{U}=\mathcal{T}\widehat{U}{\mathcal{T}}^{-1}.$$

(59)

Unlike the previous two examples, this operation $\mathcal{T}$ does not correspond to a unitary transformation, but to an antiunitary transformation, which includes complex conjugation of the state:

$$\mathcal{T}=VK,{\mathcal{T}}^{-1}=K{V}^{†},$$

(60)

where K is complex conjugation ($K|\psi \rangle ={|\psi \rangle }^{*}$) and V is a unitary operator. Plugging this form into Equation (59) we get

$$\widehat{U}=\mathcal{T}\widehat{U}{\mathcal{T}}^{-1}=VK\widehat{U}K{V}^{†}\Rightarrow V^{†}\widehat{U}V={\widehat{U}}^{*}.$$

(61)

Suppose that the evolution unitary $\widehat{U}$ has eigenstates $\left\{|{\varphi }_m\rangle \right\}$ with eigenvalues $\left\{e^{i{\varphi }_m}\right\}$. Then ${\widehat{U}}^{*}$ is also unitary, with eigenstates $\left\{{|{\varphi }_m\rangle }^{*}\right\}$ and eigenvalues $\left\{e^{-i{\varphi }_m}\right\}$. On the other hand, conjugating a unitary leaves its eigenvalues unchanged, so $V^{†}\widehat{U}V$ has the same eigenvalues $\left\{e^{i{\varphi }_m}\right\}$ as $\widehat{U}$. The only way a unitary $\widehat{U}$ can satisfy Equation (61), therefore, is if for every eigenvalue $e^{i{\varphi }_m}$ it must also have a corresponding eigenvalue $e^{-i{\varphi }_m}$. That is, every “positive energy” state must have a corresponding “negative energy” solution. (Zero-energy solutions, if any, can be taken to themselves or each other.)

We can put this condition in a more intuitive form. The inverse operator ${\widehat{U}}^{†}$ has the same eigenvectors $\left\{|{\varphi }_m\rangle \right\}$ as $\widehat{U}$, but the same eigenvalues $\left\{e^{-i{\varphi }_m}\right\}$ as ${\widehat{U}}^{*}$. We can define a new unitary operator $\tau$

$$\tau =\sum _m|{\varphi }_m\rangle \left({\langle {\varphi }_m|}^{*}\right)V^{†}.$$

(62)

This operator satisfies the symmetry relation

$$\tau \widehat{U}{\tau }^{†}={\widehat{U}}^{†}.$$

(63)

The existence of such a unitary $\tau$ in Equation (63) implies the time-reversal symmetry defined in Equation (59). Intuitively, $\tau$ is the equivalent of the classical tem reversal operation of flipping the momenta of all particles—it transforms the state so that it evolves “backwards” by transforming “positive energy” states into “negative energy” states and vice versa.

The QCAs defined in this paper all have this time-reversal symmetry, and therefore will automatically have negative-energy solutions. This creates a dilemma in how to go from free theories to interacting theories without allowing negative-energy solutions to be excited.

7. Including Interaction Terms in 1D

QCAs, as described in the earlier sections, can produce free quantum field theories in their long-wavelength (or low-energy) limit. Due to their time-reversal symmetry, these QCAs have matching negative-energy and positive-energy solutions. For free theories, the existence of negative-energy solutions has no important consequences. One can restrict the initial state to include only positive-energy states. For fermionic theories there is the additional strategy of redefining the vacuum state to be the Dirac Sea, with all negative-energy states occupied and all positive-energy states vacant. Excitation of a positive-energy state would leave a negative-energy hole: a particle-antiparticle pair. This approach has the added advantage that the dynamics can be chosen to conserve particle number.

Having two interacting QCAs, however, will raise an immediate issue. The eigenstates of $\widehat{U}$—the negative- and positive-energy states—are delocalized states of definite momentum. If the interaction is local—as would be expected in a QCA—then (as we shall see) those local operators must include a non-zero component of negative-energy operators. For fermionic QCAs with a Dirac Sea vacuum this is not necessarily a problem; it merely means that interactions can give rise to particle-antiparticle pairs. But bosonic QCAs have no equivalent of the Dirac Sea; the negative energy states can have any number of excitations. So it is difficult to avoid an interaction that produces a cascade of particle-antiparticle pairs, paid for by emission of negative-energy bosons.

To simplify the analysis we will consider a 1D fermionic QCA coupled to a 1D bosonic QCA. One advantage of 1D is that it is possible to construct a fermionic QCA using only low-dimensional (qubit) local subsystems [18]. The tension between locality and positive energy already exists in 1D, and the solutions are easier to find and numerically investigate.

7.1. The 1D QCAs

For our 1D QCAs, each site contains two local subsystems labeled $x,+$ and $x,-$, at positions $x=j\Delta x$ for integers $0\le j<N$ with periodic boundaries. The evolution unitary alternately couples two qubits at neighboring sites and the two qubits at a single site, as shown in Figure 1.

In the fermionic case [18], the local subsystems are qubits, and the evolution unitary alternates two unitaries that act on pairs of local subsystems:

$${\widehat{U}}_F={\widehat{C}}_F{\widehat{\Sigma }}_F.$$

(64)

In this equation, the coin operator is

$${\widehat{C}}_F=\cdots \otimes C_{x-\Delta x}\otimes C_x\otimes C_{x+\Delta x}\otimes \cdots ,$$

(65)

where each of these local unitaries $C_x$ acts on the qubits at $x,+$ and $x,-$ as follows:

$$\begin{array}{ccc}\hfill C|00\rangle & =& |00\rangle ,\hfill \\ \hfill C|01\rangle & =& cos\left(\theta \right)|10\rangle +sin\left(\theta \right)|01\rangle ,\hfill \\ \hfill C|10\rangle & =& cos\left(\theta \right)|01\rangle -sin\left(\theta \right)|10\rangle ,\hfill \\ \hfill C|11\rangle & =& -|11\rangle ,\hfill \end{array}$$

(66)

where 1 indicates the presence (and 0 the absence) of a particle in the corresponding state, and $\theta$ is the small dimensionless parameter corresponding to the particle mass. The shift operator is

$$\widehat{\Sigma }=\cdots \otimes S_{x-\Delta x}\otimes S_x\otimes S_{x+\Delta x}\otimes \cdots ,$$

(67)

where each of these local unitaries $S_x$ acts on the qubits at $x,+$ and $x+\Delta x,-$ as follows:

$$\begin{array}{ccc}\hfill S|00\rangle & =& |00\rangle ,\hfill \\ \hfill S|01\rangle & =& |10\rangle ,\hfill \\ \hfill S|10\rangle & =& |01\rangle ,\hfill \\ \hfill S|11\rangle & =& -|11\rangle .\hfill \end{array}$$

(68)

In both cases, the − sign on the $|11\rangle$ state produces a phase of $-1$ when two particles cross each other, producing Fermi statistics. As shown in [18] one can introduce creation and annihilation operators $b_{x,\pm }^{†}$ and $b_{x,\pm }$ that obey the usual anticommutation properties, and rewrite the evolution unitaries in terms of them:

$$\begin{array}{ccc}\hfill C_x& =& exp\left\{-\left({\displaystyle \frac{\pi }{2}}+\theta \right)\left(b_{x,+}^{†}{b}_{x,-}+b_{x,-}^{†}{b}_{x,+}\right)\right\},\hfill \\ \hfill S_x& =& exp\left\{-{\displaystyle \frac{i\pi }{2}}\left(b_{x,+}^{†}{b}_{x+\Delta x,-}+b_{x+\Delta x,-}^{†}{b}_{x,+}\right)\right\}.\hfill \end{array}$$

(69)

In the bosonic case, each of the local subsystems is a harmonic oscillator mode with creation and annihilation operators $a_{x,\pm }^{†}$ and $a_{x,\pm }$. The evolution unitary ${\widehat{U}}_B$ has the same local structure as in the fermionic QCA, but the local unitaries become

$$\begin{array}{ccc}\hfill C_x& =& exp\left\{-{\displaystyle \frac{i\pi }{2}}\left(a_{x,+}^{†}{a}_{x,-}+a_{x,-}^{†}{a}_{x,+}\right)\right\},\hfill \\ \hfill S_x& =& exp\left\{-{\displaystyle \frac{i\pi }{2}}\left(a_{x,+}^{†}{a}_{x+\Delta x,-}+a_{x+\Delta x,-}^{†}{a}_{x,+}\right)\right\}.\hfill \end{array}$$

(70)

We do not include a parameter analogous to $\theta$ in this case; these are massless bosons, which propagate in the positive or negative direction at the speed of light $c=\Delta x/\Delta t$.

We can prove that these QCAs have time-reversal symmetry by explicitly constructing the unitary operator $\tau$ in Equation (63). Let us first make the construction in the bosonic case. Examining the unitaries defined in Equation (70), we see that every term inside the exponential includes exactly one + mode and one − mode. If one could, for example, flip the sign of all the operators $a_{x,+}^{†},a_{x,+}$, this would have the effect of flipping the sign of all the exponents, which would in turn transform ${\widehat{C}}_B\to {\widehat{C}}_B^{†}$ and ${\widehat{\Sigma }}_B\to {\widehat{\Sigma }}_B^{†}$. Define the unitary operator

$$V\equiv exp\left\{i\pi \sum _x{a}_{x,+}^{†}{a}_{x,+}\right\}.$$

(71)

It is easy to check that this operator satisfies

$$V^{†}{a}_{x,\pm }V=\mp a_{x,\pm },V^{†}{a}_{x,\pm }^{†}V=\mp a_{x,\pm }^{†}.$$

(72)

This equation in turn implies that

$$V^{†}{\widehat{U}}_{B}V=V^{†}{\widehat{C}}_{B}V{V}^{†}{\widehat{\Sigma }}_{B}V={\widehat{C}}_B^{†}{\widehat{\Sigma }}_B^{†}.$$

(73)

To complete the construction, define

$$\tau \equiv {\widehat{C}}_B{V}^{†},$$

(74)

which satisfies

$$\tau {\widehat{U}}_B{\tau }^{†}={\widehat{C}}_B{\widehat{C}}_B^{†}{\widehat{\Sigma }}_B^{†}{\widehat{C}}_B^{†}={\widehat{\Sigma }}_B^{†}{\widehat{C}}_B^{†}={\widehat{U}}_B^{†}.$$

(75)

An exactly analogous construction, using operators $b_{x,+}^{†}{b}_{x,+}$, works for the fermionic QCA.

7.2. Form of Interactions

The evolution operators ${\widehat{U}}_F$ and ${\widehat{U}}_B$ represent free evolution of the two QCAs. For interacting QCAs, we could include an additional interaction unitary, so that the full evolution would be

$$\widehat{U}={\widehat{U}}_I{\widehat{U}}_B{\widehat{U}}_F.$$

(76)

For a purely local interaction,

$${\widehat{U}}_I=\cdots \otimes U_{x-\Delta x}\otimes U_x\otimes U_{x+\Delta x}\otimes \cdots$$

(77)

These local unitaries $U_x$ act on the subsystems $x,\pm$ of both the bosonic and fermionic QCAs.

For this paper, we will assume a particular form for this interaction; but we believe that the argument regarding locality and negative energy is very general, irrespective of the particular form of the interaction. Based on the Dirac Sea construction referred to above, we choose $U_x$ to conserve the number of particles in the fermionic QCA. However, based on the examples from quantum field theory, it should not conserve bosonic particle number. Using QFT as a model, we consider an interaction unitary of the form:

$${\widehat{U}}_I=exp\left\{-i\sum _x\sum _{\ell ,m,n=\pm }\left({\alpha }_{\ell ,m,n}{a}_{x,\ell }{b}_{x,m}^{†}{b}_{x,n}+\mathrm{h}.\mathrm{c}.\right)\right\}.$$

(78)

Here, the coefficients $\left\{{\alpha }_{\ell ,m,n}\right\}$ are parameters that represent coupling between particles in different internal states. (Generally, we will assume that these parameters are small in magnitude, $|{\alpha }_{\ell ,m,n}|\ll 1$.) A term like the one shown in the exponential can annihilate a boson and cause a fermion to make a transition from one internal state to another; the hermitian conjugate can create a boson while causing a transition. If one thinks of states close to the Dirac Sea, where most negative-energy states are occupied and most positive-energy states are vacant, terms like this would represent a boson giving rise to a fermion-antifermion pair, or a fermion-antifermion pair annihilating to produce a boson.

However, in this local picture there is no way of knowing whether a boson produced would be positive-energy or negative-energy, and we will now see that in general the amplitude for both is nonzero.

7.3. Coupling to Negative-Energy States

For the 1D massless bosonic QCA, the “energy” eigenstates are just the momentum states $k,\pm$:

$$a_{k,\pm }={\displaystyle \frac{1}{\sqrt{N}}}\sum _x{e}^{ixk}{a}_{x,\pm }\iff a_{x,\pm }={\displaystyle \frac{1}{\sqrt{N}}}\sum _k{e}^{-ixk}{a}_{k,\pm }.$$

(79)

The “energies” of ${\widehat{U}}_B$ are $\pm {\varphi }_k$ for ${\varphi }_k=\left|k\right|\Delta x$. If we make the convenient choice of range for k to be $-\pi /2\Delta x<k\le \pi /2\Delta x$, with $k=j\pi /N\Delta x$ for some integer $-N/2<j\le N/2$, then the result is very intuitive: the positive-energy states are $k,+$ for $k>0$ and $k,-$ for $k\le 0$; the negative energy states are $k,-$ for $k>0$ and $k,+$ for $k\le 0$. For each k there is a positive- and a negative-energy state; we can label their annihilation operators $a_{k,\mathrm{pos}}$ and $a_{k,\mathrm{neg}}$ when we want to clarify which one we mean for a general momentum k.

For the fermionic QCA the “energy” eigenstates are also momentum eigenstates, but when $\theta \ne 0$ the internal basis depends on the value of k. However, it will suffice just to work with momentum eigenstates for the fermions.

We can use the inverse transform from Equation (79) to rewrite the form of the interaction unitary in Equation (78) in terms of the momentum eigenstates:

$$\begin{array}{ccc}\hfill {\widehat{U}}_I& =& exp\left\{-{\displaystyle \frac{i}{N^{3/2}}}\sum _x\sum _{k_1,k_2,k_3}\sum _{\ell ,m,n=\pm }\left({\alpha }_{\ell ,m,n}{e}^{-i(k_1-k_2+k_3)x}{a}_{k_1,\ell }{b}_{k_2,m}^{†}{b}_{k_3,n}+\mathrm{h}.\mathrm{c}.\right)\right\}\hfill \\ & =& exp\left\{-{\displaystyle \frac{i}{N^{1/2}}}\sum _{k_1,k_3}\sum _{\ell ,m,n=\pm }\left({\alpha }_{\ell ,m,n}{a}_{k_1,\ell }{b}_{k_1+k_3,m}^{†}{b}_{k_3,n}+\mathrm{h}.\mathrm{c}.\right)\right\},\hfill \end{array}$$

(80)

where we derived this form by summing over x, which produces a Kronecker delta $N{\delta }_{k_2,k_1+k_3}$. From the form of this interaction we can immediately conclude two things:

• Momentum is conserved. The term shown in the exponent annihilates a boson with momentum $k_1$ and causes a fermion transition from a state with momentum $k_3$ to a state with momentum $k_2=k_1+k_3$. The hermitian conjugate term would cause a fermion transition from momentum $k_2$ to $k_3$ and produce a boson with momentum $k_1=k_2-k_3$.
• The coefficients ${\alpha }_{\ell ,m,n}$ are independent of the momentum k; so if it is possible to create a positive-energy boson (say with $k>0$ and $\ell =+$ then it must also be possible to create a negative-energy boson (with $k<0$ and $\ell =+$).

These results are not surprising: the energy eigenstates of the boson are completely delocalized, and cannot be distinguished at a single site x.

7.4. Interaction Range vs. Negative-Energy Coupling

7.4.1. Finite-Range Interactions

While the above argument shows that interactions localized at a single site must allow the creation of negative-energy bosons, it does not (quite) show that long-range interactions are required to prevent the production of negative-energy states. One can generalize local interactions to allow coupling between fermions and bosons over a longer, but still finite, range. A system with such finite-range interactions would still be a QCA.

To make this more concrete, consider an interaction unitary that is a product of unitaries of the following form:

$${\widehat{U}}_{I,x}=exp\left\{i\sum _{y=-{\displaystyle \frac{M\Delta x}{2}}}^{{\displaystyle \frac{M\Delta x}{2}}}\sum _{\ell ,m,n=\pm }\left({\alpha }_{y,\ell ,m,n}{a}_{x+y,\ell }{b}_{x,m}^{†}{b}_{x,n}+\mathrm{h}.\mathrm{c}.\right)\right\},$$

(81)

where M is an even integer and $M/2$ is the range of the interaction. Since the coefficients ${\alpha }_{y,\ell ,m,n}$ now depend on the relative position y, when we transform to the momentum picture the transformed coefficients can depend on momentum. Is it possible to eliminate coupling to negative-energy states?

Let’s pull out the bosonic piece of the interaction. A finite-range interaction as in Equation (81) involves bosonic operators of the form

$$O_x=\sum _{y=-{\displaystyle \frac{M\Delta x}{2}}}^{{\displaystyle \frac{M\Delta x}{2}}}\left({\alpha }_{y,+}{a}_{x+y,+}+{\alpha }_{y,-}{a}_{x+y,-}\right),$$

(82)

(where we have fixed the internal states $m,n$ of the fermions and suppressed those indices). If we are allowed to choose the coefficients $\left\{{\alpha }_{y,\pm }\right\}$ arbitrarily, how small can we make the component of negative-energy operators in $O_y$ relative to the positive-energy operators? Writing the operators $a_{x+y,\pm }$ in terms of momentum operators, we get

$$O_x\equiv \sum _k\left({\tilde{\alpha }}_{k,\mathrm{pos}}{a}_{k,\mathrm{pos}}+{\tilde{\alpha }}_{k,\mathrm{neg}}{a}_{k,\mathrm{neg}}\right).$$

(83)

We can calculate the coefficients $\left\{{\tilde{\alpha }}_{k,\mathrm{neg}}\right\}$ by rewriting the original expression for $O_x$ in terms of momentum operators:

$${\tilde{\alpha }}_{k,\mathrm{neg}}={\displaystyle \frac{1}{\sqrt{N}}}\sum _y\left\{\begin{array}{cc}{e}^{-i(x+y)k}{\alpha }_{y,-},& k>0,\\ e^{-i(x+y)k}{\alpha }_{y,+},& k\le 0.\end{array}\right.$$

(84)

We want to choose the coefficients $\left\{{\alpha }_{y,\pm }\right\}$ to minimize the $\left\{{\tilde{\alpha }}_{k,\mathrm{neg}}\right\}$. We can quantify this by the sum

$$\begin{array}{ccc}\hfill \sum _k{\left|{\tilde{\alpha }}_{k,\mathrm{neg}}\right|}^2& =& {\displaystyle \frac{1}{N}}\sum _{y,y^{\prime }}\left[{\alpha }_{y,-}^{*}{\alpha }_{y^{\prime },-}\sum _{k>0}{e}^{-i(y-y^{\prime })k}+{\alpha }_{y,+}^{*}{\alpha }_{y^{\prime },+}\sum _{k\le 0}{e}^{-i(y-y^{\prime })k}\right]\hfill \\ & \equiv & \sum _{j,j^{\prime }=1}^{M+1}\left({\alpha }_{(j-M/2)\Delta x,-}^{*}{\alpha }_{(j^{\prime }-M/2)\Delta x,-}{d}_{j{j}^{\prime }}+{\alpha }_{(j-M/2)\Delta x,+}^{*}{\alpha }_{(j^{\prime }-M/2)\Delta x,+}{d}_{j{j}^{\prime }}^{*}\right),\hfill \end{array}$$

(85)

where the $\left\{d_{j{j}^{\prime }}\right\}$ are elements of an $(M+1)\times (M+1)$ matrix D, and we have changed from summing over y to summing over an integer index $j=1,\dots ,M+1$, where $y=(j-M/2)\Delta x$. Note that this sum is independent of x (which reflects the translation symmetry of the interaction) and takes the form

$$\sum _k{\left|{\tilde{\alpha }}_{k,\mathrm{neg}}\right|}^2=v_{-}^{†}D{v}_{-}+v_{+}^{†}{D}^{*}{v}_{+},$$

(86)

where D is the matrix with elements $d_{j{j}^{\prime }}$ that we’ve just defined, and $v_{\pm }$ are vectors whose $M+1$ elements are the coefficients ${\alpha }_{y,\pm }$. This sum is minimized by choosing $v_{-}$ and $v_{+}$ to be the eigenvectors of D and $D^{*}$, respectively, with the smallest eigenvalue. The eigenvalues of D and $D^{*}$ are equal, and they must all be $\ge 0$, so we minimize the sum by choosing the eigenvalue closest to zero.

These matrix elements $d_{j{j}^{\prime }}$ are

$$\begin{array}{ccc}\hfill d_{j{j}^{\prime }}& =& {\displaystyle \frac{1}{N}}\sum _{\ell =1}^{N/2}{e}^{-2i\pi (j-j^{\prime })\ell /N}={\displaystyle \frac{e^{-2i\pi (j-j^{\prime })/N}}{N}}{\displaystyle \frac{1-e^{-i\pi (j-j^{\prime })}}{1-e^{-2i\pi (j-j^{\prime })/N}}}\hfill \\ & =& e^{-i\pi (j-j^{\prime })(N+1)/2N}{\displaystyle \frac{sin(\pi (j-j^{\prime })/2)}{Nsin(\pi (j-j^{\prime })/N)}}\hfill \\ & \approx & e^{-i\pi (j-j^{\prime })/2}{\displaystyle \frac{sin(\pi (j-j^{\prime })/2)}{\pi (j-j^{\prime })}},N\gg M.\hfill \end{array}$$

(87)

We can make a few observations about this matrix D:

• D is a Toeplitz matrix: its elements $d_{j{j}^{\prime }}$ depend only on the difference $j-j^{\prime }$, so they are constant along the diagonals.
• D must be a positive matrix, and we will see that its eigenvalues lie between 0 and 1.
• The diagonal elements $j=j^{\prime }$ are $d_{jj}=1/2$.
• The phase factor $e^{-i\pi (j-j^{\prime })/4}$ has no effect on the eigenvalues of D, since that phase can be absorbed into the eigenvectors $v_{-}$. We can instead consider the real matrix $\overline{D}$.

$${\overline{d}}_{j{j}^{\prime }}={\displaystyle \frac{sin(\pi (j-j^{\prime })/2)}{Nsin(\pi (j-j^{\prime })/N)}}\approx {\displaystyle \frac{sin(\pi (j-j^{\prime })/2)}{\pi (j-j^{\prime })}}.$$

(88)

The numerator takes repeating values $0,1,0,-1,0,1,\dots$, so except for $j=j^{\prime }$ the even diagonals all vanish, and the odd diagonals take values $1/\pi$, $-1/3\pi$, $1/5\pi$, $-1/7\pi$, and so forth. Interestingly, these are the coefficients of the Taylor expansion of $(1/\pi )arctan\left(1\right)$.

For example, for $M=4$ the matrix $\overline{D}$ is

$$\overline{D}={\displaystyle \frac{1}{\pi }}\left(\begin{array}{ccccc}\pi /2& 1& 0& -1/3& 0\\ 1& \pi /2& 1& 0& -1/3\\ 0& 1& \pi /2& 1& 0\\ -1/3& 0& 1& \pi /2& 1\\ 0& -1/3& 0& 1& \pi /2\end{array}\right).$$

(89)

7.4.2. Toeplitz Matrices and Their Eigenvalues

We can define a family of these Toeplitz matrices $\overline{D}$ as a function of the interaction range. Define a function $f\left(n\right)$ that takes an integer argument:

$$f\left(n\right)=\left\{\begin{array}{cc}1/2,\hfill & n=0,\hfill \\ 0\hfill & \mathrm{even}n\ne 0,\hfill \\ {(-1)}^{(n-1)/2}/n\pi ,\hfill & \mathrm{odd}n.\hfill \end{array}\right.$$

(90)

The matrix elements of $\overline{D}$ are ${\overline{d}}_{j{j}^{\prime }}=f(j-j^{\prime })$. We have not (yet) found a closed-form expression for the eigenvalues of this family of Toeplitz matrices, but we have found an analytical bound on the smallest eigenvalue of this matrix [32,33]. First, we introduce the discrete Fourier transform of $f\left(n\right)$:

$$\tilde{f}\left(k\right)=\sum _{n=-\infty }^{\infty }{e}^{ink}f\left(n\right),k\in [-\pi ,\pi ].$$

(91)

We can check that $\tilde{f}$ is convergent and real, and has the properties $\tilde{f}\left(k\right)=\tilde{f}(-k)$, $\tilde{f}(k+2\pi )=\tilde{f}\left(k\right)$. Using the above observation connecting the matrix elements of $\overline{D}$ to the Taylor expansion of the arctangent, we can explicitly calculate

$$\tilde{f}\left(k\right)={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{\pi }}\left(arctan\left(e^{ik}\right)+arctan\left(e^{-ik}\right)\right),$$

(92)

which produces a striking square wave, as shown in Figure 2.

The bound on the smallest eigenvalue follows from the following theorems.

Theorem 1

([32,33]). . If $\left\{A_M\right\}$ is a sequence of $M\times M$ self-adjoint Toeplitz matrices (defined by a function $f\left(n\right)$ with well-defined transform $\tilde{f}\left(k\right)$), and if ${\sigma }_M$ is the minimum singular value of $A_M$, then $|{\sigma }_M-{\sigma }_{\infty }|<e^{-\alpha M}$ for some $\alpha >0$. (The smallest singular value of a self-adjoint matrix is the smallest absolute value of an eigenvalue, $\left|\lambda \right|$. Since the matrices $\left\{A_M\right\}$ are positive, this is the same as the smallest eigenvalue.)

Theorem 2

([32,33]). . The limit ${\sigma }_{\infty }$ of the smallest singular value of the matrix family $\left\{A_M\right\}$ is ${min}_k\tilde{f}\left(k\right)$.

As we can see in Figure 2, the minimum value of $\tilde{f}\left(k\right)$ is 0. So the minimum eigenvalues of the matrices $\overline{D}$ approach 0 exponentially quickly as M becomes large. (Note that the form of $\tilde{f}\left(k\right)$ also implies that the maximum eigenvalue approaches 1 exponentially quickly from below.) While we do not have an analytical expression for α, we can estimate it numerically.

7.4.3. Numerical Results

We have numerically calculated the minimum eigenvalue for the $(M+1)\times (M+1)$ matrix $\overline{D}$ for an interaction range up to $M=100$. As shown in Figure 3, these eigenvalue drop exponentially with the range M. (Because these eigenvalues are so tiny, high-precision numerical calculations are required, with up to 110 digits of precision.) Numerical fitting suggests that these eigenvalues fall off as $e^{-\alpha M}$ for $\alpha \approx 1.755$. (We note that this is close to $\sqrt{\pi }$, but in the absence of an analytical formula this is likely to be a coincidence.)

The behavior of this eigenvalue tells us two things about locality and negative energies. First, the coupling to negative energy modes is nonzero for every finite range. Any pair of coupled QCAs of this sort will have an amplitude to produce negative-energy bosons; to eliminate this amplitude requires non-vanishing interactions at arbitrary distances. But second, since the minimum eigenvalue falls off so rapidly with range, it is possible to suppress these amplitudes to the point where they are essentially negligible over any desired length of time. If we let the lattice size N become very large, with a very small lattice spacing $\Delta x$ we can have the interaction range $M\Delta x$ be imperceptibly small; effectively indistinguishable from a point interaction at longer length scales. If we take the limit $\Delta x\to 0$ and $N\to \infty$ to recover QFT, these become zero-range interactions.

In fact, this near-locality is even stronger than that. In Figure 4 we plot the absolute value of the eigenvector component for the minimum value of $\overline{D}$ when $M=100$. To the eye it appears to be a near-perfect Gaussian. Numerical fitting shows that it is, in fact, extremely close to a perfect Gaussian. This means that the amplitude to produce negative-energy bosons can be exponentially suppressed by an interaction over a finite range that also falls off rapidly away from the center.

One final note about the accuracy of these numerical calculations. The form of the matrix used to calculate the results shown in Figure 3 and Figure 4 is given by the large N approximation in Equation (88), and is exact only in the limit $N\to \infty$. However, any corrections should scale like $O\left({(M/N)}^2\right)$, and so also be extremely small.

8. Conclusions

We have shown in this paper that it is possible to construct a QCA that recovers free QED in the long-wavelength limit. In the case of bosons, this is accomplished by making the degrees of freedom of the local sites harmonic oscillator modes, and coupling them by alternating quadratic interactions. This construction produces the Maxwell QFT by using a six-dimensional internal space. The result is a two-dimensional positive-energy sector (representing the two helicity states), a two-dimensional negative-energy counterpart, and a two-dimensional zero-energy sector that is static and decoupled from the rest.

Using oscillator modes makes the HIlbert space of the local sites infinite dimensional. This is an interesting comparison with our earlier derivation of QCAs with fermionic statistics in 2 or more spatial dimensions, which required us to confine the state to the totally antisymmetric subspace of a QCA with high local dimension. It seems that to recover fermionic or bosonic statistics requires high-dimensional subsystems, albeit for different reasons.

The existence of negative-energy solutions for both the fermionic and bosonic cases follows from the time-reversal symmetry of the QCAs. This creates a dilemma, since a purely local interaction between two QCAs will then always be able to excite negative-energy bosons, which could produce a cascade of particle-antiparticle creation. We examined this problem in a 1D model and showed that amplitudes for negative-energy states cannot be eliminated by interactions with finite range, but they can be suppressed without limit, with interactions that appear effectively local at longer length scales.

This demonstration that the amplitude of negative-energy solutions can be suppressed exponentially in the interaction range was shown for a particular form of interaction in a 1D model. Numerical studies are more challenging in higher spatial dimensions, since the size of the analogous matrix will grow like a higher power of the range. It also will no longer be a Toeplitz matrix, which may complicate the analysis. While we do not know, we expect that the suppression will most likely be even greater in higher dimensions, given the larger number of degrees of freedom available. It also raises the question of how zero-energy degrees of freedom, which arise in the bosonic QCA models derived above in higher dimensions, play a role in interacting theories. This is a subject for future research.