# Clifford algebra from quantum automata and unitary Wilson fermions

Pablo Arnault

► To cite this version:

Pablo Arnault. Clifford algebra from quantum automata and unitary Wilson fermions. Physical Review A, 2022, 106, pp.012201. 10.1103/PhysRevA.106.012201. hal-03577279

HAL Id: hal-03577279

https://hal.science/hal-03577279

Submitted on 16 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Clifford algebra from quantum automata and unitary Wilson fermions

Pablo Arnault*

Laboratoire Méthodes Formelles, Université Paris-Saclay, CNRS,
ENS Paris-Saclay, CentraleSupélec, INRIA, F-91190 Essonne, France

We introduce a spacetime discretization of the Dirac equation that has the form of a quantum automaton and that is invariant upon changing the representation of the Clifford algebra, as the Dirac equation itself. Our derivation follows Dirac’s original one: We required that the square of the discrete Dirac scheme be what we define as an acceptable discretization of the Klein-Gordon equation. Contrary to standard lattice gauge theory in discrete time, in which unitarity needs to be proven, we show that the quantum automaton delivers naturally unitary Wilson fermions for any choice of Wilson’s parameter.

**Introduction.** Symmetries are at the grounding of physics. Indeed, when trying to establish a law empirically, one typically starts by discarding irrelevant parameters by identifying the symmetries of the problem under study, e.g., spatial or temporal symmetries. To actually establish a symmetry empirically, one has to reproduce an experiment several times, making physically vary the parameter suspected not to be relevant in the law, e.g., changing the location of the experiment or the time at which it is carried out, and see if the result of the experiment is still the same. Thus, a given physical law has the mathematical invariances, also called symmetries, that reflect the physical symmetries of the phenomena it describes.

The Dirac equation (DE) is the law of motion of relativistic quantum particles of matter (more precisely, of so-called fermionic particles), and is one of the main equations of the standard model of particle physics. One fundamental property of the DE is that its square delivers the Klein-Gordon equation. Now, the constraints imposed by this property on the equation take the form of algebraic relationships, the so-called Clifford algebra, between objects that enter the definition of the DE, the so-called gamma matrices. Any particular family of gamma matrices that satisfy the Clifford algebra is thus a good candidate to write down the DE and describe with it the appropriate physics; such a family is called a representation of the Clifford algebra. Hence, the DE has the following symmetry: it is invariant upon changing the Clifford-algebra representation. This symmetry is fundamental: the physics is captured by the Clifford algebra, and not by a single, particular family of gamma matrices chosen. Moreover, the Clifford algebra is at the grounding of modern geometric algebra [1], which shows its fundamental character beyond the DE, for understanding the geometry of our world.

It is well-known that numerically simulating real-time dynamics of quantum multiparticle systems is exponentially hard, and that quantum computers could overcome this difficulty. The DE is a central ingredient of these quantum multiparticle systems in the relativistic regime. In order to perform simulations involving the DE, one starts by discretizing it on a spacetime grid [2, 3]. But, if one does so without care, one breaks certain symmetries satisfied by the DE. Now, discrete-time quantum walks (DQWs) are quantum transport schemes in discrete spacetime that have been the subject of much attention in the last 30 years, in particular because of their suc-

cess as discretizations of the DE that preserve various of its symmetries [4–6], also when coupled to the fundamental force fields of nature [7–14]. There are two main such symmetries that DQWs preserve: the first is unitarity; the second, which is actually rather a property, is relativistic locality, i.e., the group velocity in vacuum cannot overcome the speed of light. We call Dirac DQW a DQW corresponding to a discretization of the DE that exhibits the two aforementioned properties.

Up to our knowledge, no existing Dirac DQW has yet managed to preserve the invariance of the DE upon changing the Clifford-algebra representation. More precisely, the representation chosen when discretizing the DE with a Dirac DQW is actually always the same, e.g., in 1 + 1 dimensions, the first alpha matrix always equals the third Pauli matrix. If one chooses initially a different representation, one will need a different discretization method. So, not only the structure of the resulting discretization, but also the method used to discretize the equation, depend on the Clifford-algebra representation, which is strongly unsatisfying.

In the present work, we give a solution to this issue, by introducing a Dirac-DQW-based discretization method of the DE that can be carried out whatever the Clifford-algebra representation. More precisely, we make it possible the emergence of the Clifford algebra out of the operators defining the Dirac DQW. This is achieved by choosing appropriate such operators, but still independent of any choice of basis; they only need to satisfy the above-mentioned Clifford algebra. This Clifford algebra is obtained, following Dirac’s own procedure in continuous spacetime, by requiring that the Dirac DQW squares a spacetime-discretized version of the Klein-Gordon equation that we consider valid, i.e., it must not only deliver the correct Klein-Gordon equation in the continuum limit, but also (i) be applicable to scalar state sequences *at the discrete level*, as well as satisfying an extra condition, namely, (ii) the vanishing of the crossed term *at the discrete level*. We then show that DQWs also contain, “naturally” (provided we make certain appropriate choices), a Wilson term that makes fermion doubling be avoided, which is well-known in lattice gauge theory (LGT). Unitarity is maintained at each step of the derivation, and this *for any choice* $r \in \mathbb{R}$ of Wilson’s parameter.

**Discrete-time quantum walks.** A DQW is a unitary automaton with an ultralocal evolution operator. The system is called a *walker*, and its state at time $j \in \mathbb{N}$ is given by a sequence $\Psi_j : (j, p) \mapsto \Psi_{j,p}$ defined over a lattice with

<footnote>* pablo.arnault@inria.fr</footnote>
sites labeled by $p \in \mathbb{Z}$; the lattice is considered 1D for simplicity. The dynamics of this system is called a *walk*, written $\Psi_{j+1} = \mathcal{U}\Psi_j$, where $\mathcal{U}$ is the unitary one-step evolution *walk operator*. That $\mathcal{U}$ is *ultralocal* means, by definition, that the internal state $\Psi_{j+1,p}$ solely depends on internal states $\Psi_{j,p'}$ that belong to a bounded spatial-lattice neighborhood around $p$. In the context of automata, ultralocality of $\mathcal{U}$ is often implicit. In quantum computation it is frequent that $\mathcal{U}$ is not ultralocal: this can be either because we purposely apply, in discrete time, a gate which is non-local [15], or because we are in continuous time and evolve a system via a nearest-neighbors Hamiltonian, which yields approximate but not exact ultralocality via Lieb-Robinson bounds [16].

To specify the nature of $\Psi_j$, it is useful to invoke Meyer's no-go result [5, 17]: No non-trivial ultralocal unitary automaton with a one-step evolution operator that is homogeneous in space (i.e., translationally invariant) can have a scalar walker; the minimum number of internal components for $\Psi_{j,p}$ is thus 2. Hence, we consider $\Psi_{j,p} \in \mathcal{H}_2$, a complex Hilbert space of dimension 2. In formal parallel with classical random walks [18], in which a coin is tossed to determine the direction taken by the walker, the internal state $\Psi_{j,p}$ is called *coin state*, and $\mathcal{H}_2$ is called *coin space*.

**Transport coin operators of a DQW.** Consider the DQW defined initially, $\Psi_{j+1} = \mathcal{U}\Psi_j$. The one-step evolution operator $\mathcal{U}$ can be given under a multiplicative [19, 20] or an additive form [5, 6]. We here give $\mathcal{U}$ under the generic additive form $\mathcal{U} := W_{-1}\mathcal{T}_1^{-1} + W_1\mathcal{T}_1 + W_0$. In this equation, (i) $\mathcal{T}_1$ is the translation operator by one lattice site in the direction of growing ps (that is, $(\mathcal{T}_1\Psi)_{j,p} := \Psi_{j,p-1}$), and (ii) the $W_i$s, $i = -1, 0, 1$, are operators acting solely on the coin space, that we call *jump coin operators*. While one may view the $W_i$s as $2 \times 2$ complex matrices, viewing them abstractly is actually the purpose of the present work, and we *will not introduce any basis of* $\mathcal{H}_2$. The translation operator is by construction related to the momentum operator, $\mathcal{K}$, via $\mathcal{T}_1 = e^{-i\mathcal{K}}$, where $\mathcal{K} \in \mathbb{R}^{+,*}$ will be identified further on with the spatial lattice spacing. In the Supplemental Material, we derive the constraints that the unitarity of $\mathcal{U}$ imposes on the $W_i$s [5, 6].

We define the following coin operators, that we call *transport coin operators*,

$$B := W_1 - W_{-1} \qquad (1a)$$

$$V := W_1 + W_{-1} \qquad (1b)$$

$$M := \sum_{i=-1,0,1} W_i = V + W_0. \qquad (1c)$$

In terms of these operators, $\mathcal{U}$ reads $\mathcal{U} = \frac{1}{2}(V - B)\mathcal{T}_1^{-1} + \frac{1}{2}(V + B)\mathcal{T}_1 + M - V$. In the Supplemental Material, we translate on the transport coin operators the constraints imposed on the *jump coin operators* by the unitarity of $\mathcal{U}$.

**Local Hamiltonian of a DQW.** From the one-step dynamics $\Psi_{j+1} = \mathcal{U}\Psi_j$, one can conceive a dynamics $i(\Psi_{j+1} - \Psi_{j-1})/2 = \mathcal{H}\Psi_j$ determined by the Hermitean operator $\mathcal{H} := \frac{1}{2}(\mathcal{U} - \mathcal{U}^\dagger)$, which is (ultra)local since $\mathcal{U}$ is ultralocal, and that we call *local Hamiltonian of the DQW*. This dynamics is “two-step”, meaning that while the one-step dynamics takes as initial condition $\Psi_{j=0}$, the two-step one takes as initial condition both $\Psi_{j=0}$ and $\Psi_{j=1}$. The two-step dynamics is

equivalent to the one-step one provided that $\Psi_1 = \mathcal{U}\Psi_0$. That the Hamiltonian $\mathcal{H}$ is (ultra)local is in contrast with the case of the well-known effective Hamiltonian of the DQW. Note that both Hamiltonians are related by a proportionality constant in Fourier space [21].

The operator $\mathcal{H}$ can be written in terms of the transport coin operators as

$$\mathcal{H} = \mathcal{H}'_Q := A^1(-i\mathcal{D}_1) + \frac{r}{2}Q(-\mathcal{L}) + \epsilon m A^0, \qquad (2)$$

where $m \in \mathbb{R}^+$ and $r \in \mathbb{R}$ are two parameters that we force to appear (we also force $\epsilon$ to appear in $\epsilon m$). In Eq. (2), we have introduced (i) the following operators acting on the position space solely, $\mathcal{D}_1 := \frac{1}{2}(\mathcal{T}_1^{-1} - \mathcal{T}_1)$ and $\mathcal{L} := \mathcal{T}_1^{-1} - \mathcal{T}_1 - 2$, and the following coin operators made out of the transport coin operators, $A^1 := (B+B^\dagger)/2$, $Q := -\frac{1}{r}(V-V^\dagger)/2$, and $A^0 := \frac{i}{\epsilon m}(M-M^\dagger)/2$. We have used the notation $\mathcal{H} = \mathcal{H}'_Q$ because we are going to consider both the case $Q = 0$ and the case $Q \neq 0$.

**Dirac-continuum-limit requirement.** The two-step dynamics can be written

$$i(\mathcal{D}_0\Psi)_j = \mathcal{H}\Psi_j, \qquad (3)$$

where $\mathcal{D}_0 := (\mathcal{T}_0^{-1} - \mathcal{T}_0)/2$, with $\mathcal{T}_0$ the shift by one lattice site forward in time, i.e., $(\mathcal{T}_0\Psi)_j := \Psi_{j-1}$. We wish that the two-step dynamics deliver the $(1+1)$D DE in the continuum limit. In order to take the continuum limit, we introduce continuous time and space coordinates $t$ and $x$, as well as a function of these continuous coordinates, $\Psi(\cdot, \cdot) : (t, x) \mapsto \Psi(t, x)$, that is as smooth as wished, and that coincides at coordinates $(t_j := j\epsilon, x_p := pa)$ with the “value” taken by the coin state at point $(j, p)$, that is, $\Psi(t_j, x_p) := \Psi_{j,p}$. We then consider the ballistic scaling $\epsilon = a$ [22], and Taylor expand the two-step dynamics in time and space around the point $(t_j, x_p)$, at second order in $\epsilon$.

For the continuum limit of the two-step dynamics, Eq. (3) (divide Eq. (3) by $\epsilon$ and let $\epsilon \to 0$ after the Taylor expansion at second order), to coincide with the $(1+1)$D DE, it is sufficient that the two following *Dirac-continuum-limit constraints* be satisfied, $A^1 := \frac{B+B^\dagger}{2} \sim \alpha^1$ and $\epsilon m A^0 := i \frac{M-M^\dagger}{2} \sim \epsilon m \alpha^0$, where $\alpha^1$ and $\alpha^0$ must satisfy $(\alpha^0)^2 = (\alpha^1)^2 = 1$ and $\alpha^0\alpha^1 + \alpha^1\alpha^0 = 0$, in order for them to correspond to the well-known $\alpha$ operators of the DE.

Note that no continuum-limit constraint is imposed on $Q$ and thus neither on $V$ (apart from $Q$ not scaling as $\epsilon^\delta$ with $\delta \le -1$, i.e., we must have $\delta > -1$), because as $\epsilon \to 0$, we have that $\mathcal{L} \sim \epsilon^2 \partial_x^2$ while $\mathcal{D}_1 \sim \epsilon \partial_x$, so that the term $\frac{r}{2}Q(-\mathcal{L})$ in Eq. (2) vanishes in the continuum (again, provided $\delta > -1$); we call this term *Wilson $Q$ term*. The Wilson $Q$ term can be chosen non-vanishing and the scheme still deliver the DE. This observation will be useful further on, but for now let us consider that $Q = 0$. At this point one may be tempted to make the trivial choice $V = 0$ to obtain $Q = 0$, but this is unlikely to be possible due to the unitarity constraints (see the Supplemental Material). The constraint we have in order to obtain $Q = 0$ is $V = V^\dagger$. We will see below what choice we finally make for $V$.
**Klein-Gordon-square requirement and Clifford algebra from quantum automata.** Squaring the equation $i\mathcal{D}_0\Psi|_j = \mathcal{H}'_{Q=0}\Psi_j$ delivers $\mathcal{D}_0^2\Psi|_j = -(\mathcal{H}'_{Q=0})^2\Psi_j$, with $(\mathcal{H}'_{Q=0})^2 = -(A^1)^2\mathcal{D}_1^2 + (\epsilon m)^2(A^0)^2 + \epsilon m[A^0A^1 + A^1A^0](-i\mathcal{D}_1)$. Thanks to the Dirac-continuum-limit constraints, we obtain the correct Klein-Gordon equation in the continuum limit (divide by $\epsilon^2$ and let $\epsilon \to 0$ after a Taylor expansion at second order).

Now, we wish that the discrete scheme obtained by squaring Eq. (3) be a valid discretization of the Klein-Gordon equation, and by valid we mean, not only that it delivers the correct continuum limit, but that it be applicable, at the discrete level, to scalar state sequences (as it is the case in the continuum setting), and for this we need to impose $(A^0) \propto 1$ and $(A^1)^2 \propto 1$, i.e., not operator-valued<fnref n="1" />, as well as $A^0A^1 + A^1A^0 \propto 1$ or equal to 0. We choose to impose that the A operators $A^0$ and $A^1$ satisfy $A^0A^1 + A^1A^0 = 0$, i.e., we impose the crossed term to vanish (which may be taken as part of the definition of what we consider a valid discrete Klein-Gordon scheme). To sum up, what we impose is that the A operators satisfy the same algebra as that of the $\alpha$ operators of the DE up to multiplicative factors, that is, we impose

$$(A^0)^2 \propto 1 \qquad (4a)$$

$$(A^1)^2 \propto 1 \qquad (4b)$$

$$A^0A^1 + A^1A^0 = 0. \qquad (4c)$$

We now have to find A operators which satisfy (i) this algebra, (ii) the Dirac-continuum-limit constraints, and (iii) the unitarity constraints. Can we find such operators? We are going to see that “yes”, a result which is a priori not trivial at all. The trivial choice $A^1 \propto \alpha^1$, more concretely, $B \propto \alpha^1$, is going to make us reach our purpose, so that we make this choice. What is untrivial is the proportionality constant, as well as the choice of M (which determines $A^0$). Indeed, a standard choice for the mass term $\epsilon m A^0$ in the literature is $e^{-i\epsilon m \alpha^0}$, but this term together with $A^1 \propto \alpha^1$ makes Constraint (4c) unsatisfied. The term that is going to make us reach our purpose is a suggestion from both Feynman’s original scheme [23], and Succi’s quantum lattice Boltzmann schemes [24, 25] (while $e^{-i\epsilon m \alpha^0}$ is a suggestion from multiplicative constructions, while here our construction is additive): we choose $M = \mu_\epsilon(1 - i\epsilon m \alpha^0)$, where $\mu_\epsilon$ is a factor that is imposed to us by one of the unitarity constraints (see the Supplemental Material), namely, $M^\dagger M = 1$, so that $\mu_\epsilon := \frac{1}{\sqrt{1+\epsilon^2m^2}}$. This $\mu_\epsilon$ is also the proportionality factor evoked above: we choose $B := \mu_\epsilon \alpha^1$. The algebra (4) is thus now satisfied, since $A^0 = \mu_\epsilon \alpha^0$ and $A^1 = \mu_\epsilon \alpha^1$, and this algebra is equivalent to the following Clifford algebra, $\{\Gamma^\mu, \Gamma^\nu\} = 2\tilde{\eta}^{\mu\nu}$, where $\{\cdot, \cdot\}$ is the anticommutator, and where we have introduced the following modified Minkowski metric, $\tilde{\eta} := \text{diag}(1/\mu_\epsilon^2, -1)$, as well as the $\Gamma$ operators, $\Gamma^0 := (A^0)^{-1} = \alpha^0/\mu_\epsilon$, and $\Gamma^1 := (A^0)^{-1}A^1 := \alpha^0\alpha^1 := \gamma^1$. One can check that the other and last unitarity constraint involving M but not V, namely, $B^\dagger M = M^\dagger B$, is satisfied.

Let us recap how M is built: the “$\alpha^0$” in M is both for the unitarity constraint $B^\dagger M = M^\dagger B$ to be satisfied and so that $A^0A^1 + A^1A^0 = 0$; the “$\epsilon m$” in front of the $\alpha^0$ is for the continuum limit to be the good one; the “i” is because $A^0 \propto i(M - M^\dagger)$; and finally the additional term “1” is for the unitarity constraint $M^\dagger M = 1$ to be satisfied. Now, finally, for the unitarity constraints on V to be satisfied, we can choose $V := \mu_\epsilon$. With these choices for B, M, and V, our DQW with $Q=0$ is invariant under unitary transformations of the coin state, in such a way that the algebra (4) is preserved, in exact parallel with the continuum situation.

Note that the operators $A^0$ and $A^1$ depend on the mass, while in the continuum limit none of the $\alpha$ operators do. Note also that $\mathcal{H}'_{Q=0} = \mu_\epsilon[\alpha^1(-i\mathcal{D}_1) + \epsilon m \alpha^0]$ (and $(\mathcal{H}'_{Q=0})^2 = \mu_\epsilon^2[-\mathcal{D}_1^2 + (\epsilon m)^2]$, with $\mu_\epsilon \neq 1$ if $\epsilon m \neq 0$ (although $\mu_\epsilon \to 1$ as $\epsilon m \to 0$)), which is the price to pay to obtain a unitary discretization while having nevertheless discretizing the transport term naively, by a symmetric finite difference $\mathcal{D}_1$. Finally, let us mention the following: in the Supplemental Material, we have only written the unitarity constraints ensuing from $U^\dagger U = 1$, but one can check that the unitarity constraints ensuing from $UU^\dagger = 1$ are also satisfied with our choices for B, M and V.

**Towards avoiding fermion doubling with the Wilson Q term.** Consider the two-step Hamiltonian of Eq. (2) for a non-vanishing Wilson Q term; it is convenient in order to further take a continuum limit to rather consider the Hamiltonian $h := \mathcal{H}/\epsilon$ (from the continuum point of view, this is actually the correct Hamiltonian dimensionally). The reason we have called the Wilson Q term that way is because its spatial-operator part, $-\mathcal{L}$, is that of the well-known Wilson term of LGT, namely, $(r\epsilon/2)\alpha^0(-\mathcal{L}/\epsilon^2)$ [26], i.e., the same as ours but with $\alpha^0$ instead of Q. In LGT, this term enables to avoid the so-called fermion doubling problem, which appears when discretizing naively the DE; these facts are recalled in the Supplemental Material. It is remarkable that such a Wilson term is naturally contained in the decomposition of a generic DQW. What Q we can take in order for the fermion doubling to be avoided, while satisfying the unitarity constraints? Is this even possible? We are going to see that the answer to the last question is “yes”, a result that is a priori not trivial at all.

**Klein-Gordon-square requirement with the Wilson Q term.** Assuming $A^0A^1 + A^1A^0 = 0$, we have that $h^2 = -(A^1)^2\frac{\mathcal{D}_1^2}{\epsilon^2} + m^2(A^0)^2 + \frac{r^2}{4\epsilon^2}Q^2\mathcal{L}^2 + \frac{r}{2\epsilon^2}[QA^0 + A^0Q](-\mathcal{L}) + \frac{r}{2\epsilon^2}[QA^1 + A^1Q](-\mathcal{L})(-i\mathcal{D}_1)$. Notice that the three last terms vanish in the continuum limit, so that we recover the Klein-Gordon equation. Now, the validity of the discrete Klein-Gordon scheme, which we define, as before, as the possibility of applying the scheme to scalar state sequences, can be obtained if we require, in addition to the previous Clifford-algebra requirement, Eqs. (4), that as $\epsilon \to 0$, (i) $Q^2 \sim 1$, and (ii) each of the two anticommutators involving Q either vanishes or is proportional to the identity; these requirements suggest to use for the choice of V the same trick as that used for the choice of M.

Hence, we choose the ansatz $V = \nu_\epsilon(1 + i\epsilon^\rho r\alpha^\lambda)$, where

<footnote n="1">¹ We will see further down that requiring mere proportionality and not equality turns out to be indeed necessary, for unitarity to be preserved.</footnote>
(i) $\alpha^\lambda$ is an $\alpha$ operator satisfying the usual algebra of the $\alpha$ operators of the DE, (ii) $\nu_\epsilon$ is a normalization factor to be determined, and (iii) $\rho$ is an exponent to be determined. Since we are in one spatial dimension, $\lambda$ should only take two possible values, $\lambda = 0$ and $\lambda = 1$ (since $\alpha^\lambda$ has been defined as an $\alpha$ operator of the DE and we are in one spatial dimension). That being said, and we mention this since it is going to be used further down, in general terms: let us say we are in $n$ spatial dimensions (so that $\lambda$ can take values from 0 to $n$) with a given representation of the Clifford algebra; then $\lambda$ can take an additional value $n+1$ provided that in $n+1$ spatial dimensions one can still find a representation of the Clifford-algebra which has the same dimension as the representation found in dimension $n^2$. We know that this is precisely the case when going from one spatial dimension to two: e.g., if $\alpha^0$ and $\alpha^1$ are represented by two Pauli matrices, then we can choose an $\alpha^2$ that is represented by the last Pauli matrix. Hence, we consider that $\lambda = 0, 1$, or 2.

The unitarity constraint $B^\dagger V = V^\dagger B$ requires $\lambda \neq 1$, which we assume, so that the anticommutator $QA^\dagger + A^\dagger Q$ vanishes. We do not change our choice of $M$, neither its normalization factor (the latter enables to satisfy $M^\dagger M = 1$). Now, both $\lambda = 0$ and $\lambda = 2$ are compatible with the unitarity constraint $2V^\dagger V = V^\dagger M + M^\dagger V$, which determines the following associated normalization factors: $\nu_\epsilon^{\lambda=0} := \mu_\epsilon \frac{1-\epsilon^{1+\rho}mr}{1+(\epsilon^\rho r)^2}$, and $\nu_\epsilon^{\lambda=2} := \frac{\mu_\epsilon}{1+(\epsilon^\rho r)^2}$. Finally, for the unitarity constraint $V^\dagger V = B^\dagger B$ to be satisfied, we need to change the normalization factor for $B$: we choose $B := \eta_\epsilon^\lambda \alpha^1$, and $\eta_\epsilon^\lambda$ is determined by the constraint just mentioned, yielding $\eta_\epsilon^\lambda := \nu_\epsilon^\lambda \sqrt{1+(\epsilon^\rho r)^2}$. All unitarity constraints are now satisfied, including those coming from $UU^\dagger = 1$. Now, since $\mu_\epsilon$ goes to 1 in the continuum limit, we see that for both $\nu_\epsilon^\lambda$ and $\eta_\epsilon^\lambda$ to go to 1 in the continuum limit, we must choose $\rho > 0$: this finally yields $Q = \nu_\epsilon^\lambda \epsilon^\rho \alpha^\lambda$, which in the continuum limit goes as $\epsilon^\rho \alpha^\lambda$ – that is, in the case $\lambda = 0$, as $\epsilon^\rho \alpha^0$ and not as $\alpha^0$ as in standard LGT. We finally have

$$h^\lambda = \eta_\epsilon^\lambda \alpha^1 \left(\frac{-i\mathcal{D}_1}{\epsilon}\right) + \mu_\epsilon m \alpha^0 + \nu_\epsilon^\lambda \epsilon^\rho \frac{r}{2\epsilon} (-\mathcal{L}) \alpha^\lambda. \quad (5)$$

We will keep working with both models $\lambda = 0$ and $\lambda = 2$ and see if one performs better than the other at avoiding fermion doubling (avoidance which for now is not a given, we are precisely going to explain further down at which condition fermion doubling is avoided).

**Avoiding fermion doubling.** To find solutions of our two-step scheme, Eq. (3), we consider a superposition-of-plane-waves ansatz (since the solution of the DE, that we seek to simulate, has this form): $\Psi(t, x) = \frac{1}{2\pi N} \sum_{i=+,-} \int_{-\infty}^{+\infty} dk \tilde{\Psi}_i(0, k) e^{-i(\omega_i(k)t - kx)}$, where $(t, x) = (t_j, x_p)$. Inserting the plane-wave ansatz into Eq. (3) with $h = \mathcal{H}/\epsilon$ given by Eq. (5), leads, after squaring, to a dis-

persion relation

$$\frac{\sin^2(\omega_i(k)\epsilon)}{\epsilon^2} = F^{\text{DQW},\lambda}(k), \quad (6)$$

where

$$F^{\text{DQW},\lambda=0}(k) := (\eta_{\epsilon}^{0})^{2} \frac{\sin^{2}(k\epsilon)}{\epsilon^{2}} + \left[ \mu_{\epsilon}m + \nu_{\epsilon}^{0} \epsilon^{\rho} \frac{r}{\epsilon} (1 - \cos(k\epsilon)) \right]^{2} \quad (7a)$$

$$F^{\text{DQW},\lambda=2}(k) := (\eta_{\epsilon}^{2})^{2} \frac{\sin^{2}(k\epsilon)}{\epsilon^{2}} + (\mu_{\epsilon}m)^{2} + \left[ \nu_{\epsilon}^{2} \epsilon^{\rho} \frac{r}{\epsilon} (1 - \cos(k\epsilon)) \right]^{2}. \quad (7b)$$

Now, discrete-time LGT is usually formulated in a Lagrangian way [27, 28], whereas we formulated our discrete-time scheme in a Hamiltonian one. In Lagrangian LGT, a term is naturally added to the action to remove the temporal doublers along with the term added to remove the spatial doublers. Here, this is not the case, and we only treat spatial doublers. Considering a low-frequency limit, $\omega_i(k) \ll 1$, of the dispersion relation, finally leads to the solutions

$$\omega_{\pm}^{\text{DQW},\lambda}(k) = \pm \sqrt{F^{\text{DQW},\lambda}(k)}. \quad (8)$$

Expressions (7) are to be compared with the dispersion relation of standard LGT,

$$F^{\text{LGT}}(k) := \frac{\sin^2(k\epsilon)}{\epsilon^2} + \left[ m + \frac{r}{\epsilon} (1 - \cos(k\epsilon)) \right]^2. \quad (9)$$

The difference between $F^{\text{DQW},\lambda=0}(k)$ and $F^{\text{DQW},\lambda=2}(k)$ is, apart from the normalization factors $\nu_\epsilon^\lambda$ and $\eta_\epsilon^\lambda$, the *crossed term* of the square of the sum in $F^{\text{DQW},\lambda=0}(k)$. In standard LGT, this crossed term is also present. Since the upcoming discussion is the same whether there is a crossed term or not, let us forget about the latter for now, we will come back to it later. Here comes the discussion. In standard LGT, what makes fermion doubling be avoided is that the function $1 - \cos(k\epsilon)$ raises the value of the dispersion relation at the edges of the Brillouin zone, i.e., at $\pm\pi/\epsilon$. More precisely, if the amplitude by which this dispersion relation is raised, which we call *raising amplitude*, does not go to zero with $\epsilon$, then for sure fermion doubling is avoided, and this is indeed the case in standard LGT since the associated raising amplitude is $(r/\epsilon)^2$, see Eq. (9). In our model, the raising amplitude is $(\nu_\epsilon^\lambda r \epsilon^{\rho-1})^2$, see Eqs. (7). Since $\nu_\epsilon^\lambda$ goes to 1 in the continuum limit, we are fine as long as

$$\rho < 1, \quad (10)$$

i.e., this condition ensures that fermion doubling is avoided, both for $\lambda = 0$ and $\lambda = 2$.

**Initial slopes.** The limit to the continuum is obtained by letting $\epsilon$ go to zero. Let us consider the gapless frequencies

$$f^{\text{M}}(k) := \sqrt{(\omega_{\pm}^{\text{M}}(k))^2 - M^{\text{M}}} = \sqrt{F^{\text{M}}(k) - M^{\text{M}}}, \quad (11)$$

where “M” is the considered model, and takes for now the “values” $M = (\text{DQW}, \lambda = 0)$ or $M = (\text{DQW}, \lambda = 2)$, and where $M^M$ is the central $(k = 0)$ gap of the model.

<footnote n="2">² Actually, this condition may be necessary only for irreducible representations: if by going from dimension $n$ to $n+1$ the dimension of irreducible representations is augmented, one may be able to simply consider in dimension $n$ a reducible representation having the same dimension as the representation of dimension $n+1$.</footnote>
Let us consider the Taylor expansion of the dispersion relations $(f^M(k))^2$ at next-to-lowest order in $\epsilon$. A priori, the two Taylor-expanded dispersion relations (one for $\lambda = 0$ and one for $\lambda = 2$) are different: more precisely, there are a priori more terms in the case $\lambda = 0$, namely, the crossed terms. Now, it turns out that these additional terms actually cancel each other, at least at next-to-lowest order in $\epsilon$, so that both the model $\lambda = 0$ and $\lambda = 2$ have, at next-to-lowest order in $\epsilon$, the same Taylor expansion in $\epsilon$, namely,

$$ (f^{\text{DQW},\lambda=0}(k))^2 \approx (f^{\text{DQW},\lambda=2}(k))^2 \approx (1 - \frac{1}{2}r^2\epsilon^{2\rho})k^2. \qquad (12) $$

Actually, even if these above-mentioned additional terms had not cancelled each other, they would be of higher order because $\rho < 1$, so both models $\lambda = 0$ and $\lambda = 2$ would still have the same small-$\epsilon$ expansion, we leave the verification of this observation to the reader. Moreover, the criterion $\rho < 1$ for the avoidance of fermion doubling in our model, is the same whether there is a crossed term or not. Hence, we consider both models $\lambda = 0$ and $\lambda = 2$ equivalently good for our task of avoiding fermion doubling and converging as fast as possible to the continuum limit. We choose $\lambda = 0$, in order for comparisons to LGT to be simpler, and from now on we call this model M = DQW. At this point the reader may wonder why we have kept the discussion about $\lambda = 0$ and $\lambda = 2$ if in the end we consider both models equivalently good and choose e.g. $\lambda = 0$: apart from a mere informative reason, this discussion is going to be useful just below in the next paragraph, to understand the role of the crossed term in the LGT model at small $\epsilon$, and the subsequent interest of getting rid of this crossed term by choosing $\lambda = 2$ in this LGT model.

We call the factor in front of $k^2$ in Eq. (12) the (squared) initial slope of the model. In standard LGT, the lowest-order modification of the initial slope comes, this time, from the crossed term (this is detailed in the Supplemental Material): at next-to-lowest order in $\epsilon$ we have that

$$ (f^{\text{LGT}}(k))^2 \approx (1 + \epsilon mr)k^2. \qquad (13) $$

Hence, for the initial slope of our model to converge faster to the continuum limit than that of the LGT model, we need to choose

$$ \rho > 0.5. \qquad (14) $$

In the end, we have, in addition to in-built unitarity, a model that performs better at reaching the continuum limit than that of LGT if we choose $\rho \in ]0.5, 1[$. That being said, if we choose in the LGT model $\lambda = 2$, then the initial slope is exactly that of the continuum, that is, 1 (we leave the verification of this to the reader), and in that case the LGT model performs better than ours at reaching the continuum limit, whatever value we choose for $\rho$.

**Figure.** In Fig. 1, we plot the gapless frequencies defined above in Eq. (11), where this time the model M takes the “values” M = Dirac, naive, LGT, DQW. We see that the naive discretization of the DE leads to two extra poles in the gapless frequency (gold plot), which causes the fermion doubling problem (see the Supplemental Material), while the

LGT (green plot) and DQW (red plot) models avoid it by creating gaps at the edges. In-built unitarity for any $r \in \mathbb{R}$ is an important advantage of our model with respect to LGT.

Take a closer look at the initial slopes (i.e., around $k = 0$). One can appreciate, though weakly, the fact that the initial slope of the DQW (resp. LGT) model is a bit smaller, see Eq. (12) (resp. bigger, see Eq. (13)) than that of Dirac fermions (although the two former slopes obviously converge to the latter in the continuum limit i.e. for $\epsilon$ going to 0). We could have made this fact more visible on the figure, but then the plots do not give the impression (which is anyways of course a truth, as shown above) that in the continuum limit one does get the correct initial slope, while in Fig. 1 this impression is given.

Let us make a final remark. In the figure, we have chosen $r = 1$ because it is a standard choice. That being said, let us mention that for $r = 1$ the discrete-time LGT model of Wilson fermions (which is a Lagrangian model) can actually be proven to be unitary. This proof does however not hold for $r \neq 1$ [27], while our model has in-built unitarity for any $r$. Pay attention that if we choose, for $r = 1$, $\lambda = 2$ in the discrete-time LGT model, unitarity would again have to be proven, and it may actually not hold.

<figure><img src="image_2.png" /><figcaption>Figure 1. Gapless frequency $f^M(k)$ of Model “M”, for “M” = “Dirac” (blue), “naive” (gold), “LGT” (green), and “DQW” (red), where $f^{\text{Dirac}}(k) = |k|$, $f^{\text{naive}}(k) = \frac{|\sin(k\epsilon)|}{\epsilon}$, $f^{\text{LGT}}(k) = \sqrt{\frac{\sin^2(k\epsilon)}{\epsilon^2} + [m + \frac{r}{\epsilon}(1 - \cos(k\epsilon))]^2 - m^2}$, and $f^{\text{DQW}}(k) = \sqrt{(\eta_{\epsilon}^0)^2 \frac{\sin^2(k\epsilon)}{\epsilon^2} + [\mu_{\epsilon}m + \nu_{\epsilon}^0 \epsilon^{\rho} \frac{r}{\epsilon}(1 - \cos(k\epsilon))]^2 - (\mu_{\epsilon}m)^2}$, for $\epsilon = 0.1$, $m = 1$, $r = 1$ (Wilson’s choice), and $\rho = 0.6$.</figcaption></figure>

**Conclusion.** We showed that a Clifford algebra emerges out of the internal-state operators defining a quantum-automaton discretization of the Dirac equation. This discretization, which is a DQW, is unitary by construction, while discrete-time versions of lattice gauge theory are usually Lagrangian so that unitarity is not in-built. Our DQW is invariant under representation changes of the above-mentioned Clifford algebra, which parallels exactly the continuum situation. Moreover, we show that DQWs naturally contain a Wilson term that makes spatial fermion doubling be avoided, and this without breaking unitarity and for any choice $r \in \mathbb{R}$ of Wilson’s parameter.

**Acknowledgements.** The author thanks C. Cedzich for making him notice that the two-step Hamiltonian is (ultra)local, while the effective Hamiltonian, as a logarithm of the one-step evolution operator, is not (infinite series of grow-
ing powers of the one-step evolution operator) This is an in-
teresting feature which is in the prolongation of the ultralocal-
ity of the one-step evolution operator that is seeked in the de-

velopment of quantum cellular automata and their one-particle
versions, DQWs. The author also thanks A. Grinbaum for his
stylistic and grammatical energetic advice and corrections.

[1] D. Lundholm and L. Svensson, arXiv:0907.5356 (2009).

[2] K. G. Wilson, Phys. Rev. D 10, 2445 (1974).

[3] J. B. Kogut and L. Susskind, Phys. Rev. D 11, 395 (1975).

[4] I. Bialynicki-Birula, Phys. Rev. D 49, 6920 (1994).

[5] D. A. Meyer, J. Stat. Phys. 85, 551 (1996).

[6] A. Bisio, G. M. D'Ariano, and A. Tosini, Ann. Phys. (N. Y.) 354, 244 (2015).

[7] P. Arnault and F. Debbasch, Physica A 443, 179 (2016).

[8] P. Arnault and F. Debbasch, Phys. Rev. A 93, 052301 (2016).

[9] C. Cedzich, T. Rybár, A. H. Werner, A. Alberti, M. Genske, and R. F. Werner, Phys. Rev. Lett. 111, 160601 (2013).

[10] G. Di Molfetta, M. Brachet, and F. Debbasch, Phys. Rev. A 88, 042301 (2013).

[11] P. Arnault and F. Debbasch, Ann. Phys. (N. Y.) 383, 645 (2017).

[12] P. Arrighi, S. Facchini, and M. Forets, Quantum Inf. Process. 15, 3467 (2016).

[13] P. Arrighi and S. Facchini, *Quantum Info. Comput.* **17**, 810 (2017).

[14] P. Arnault, G. Di Molfetta, M. Brachet, and F. Debbasch, Phys. Rev. A 94, 012335 (2016).

[15] N. M. P. Neumann, R. van Houte, and T. Attema, in *Computational Science – ICCS 2020*, edited by V. V. Krzhizhanovskaya, G. Závodszky, M. H. Lees, J. J. Dongarra, P. M. A. Sloot, S. Brissos, and J. Teixeira (Springer International Publishing, Cham, 2020) pp. 605–615.

[16] E. H. Lieb and D. W. Robinson, Commun. Math. Phys. 28, 251 (1972).

[17] D. A. Meyer, Phys. Lett. A 223, 337 (1996).

[18] J. Kempe, Contemporary Physics 44, 307 (2003).

[19] H. Vogts, Ph.D. thesis, Technischen Universität zu Braunschweig (2013).

[20] C. Cedzich, T. Geib, and R. F. Werner, arXiv:2102.12951 (2021).

[21] P. Arnault, B. Pepper, and A. Pérez, Phys. Rev. A 101, 062324 (2020).

[22] F. Debbasch, G. Di Molfetta, D. Espaze, and V. Foulonneau, Phys. Scripta 151, 014044 (2012).

[23] S. S. Schweber, Rev. Mod. Phys. 58, 449 (1986).

[24] S. Succi and R. Benzi, Physica D 69, 327 (1993).

[25] S. Succi, F. Fillion-Gourdeau, and S. Palpacelli, Eur. Phys. J. Quantum Technol. 2, 12 (2015).

[26] L. Susskind, Phys. Rev. D 16, 3031 (1977).

[27] M. P. Hernández, in *Modern Perspectives in Lattice QCD: Quantum Field Theory and High Performance Computing, Lecture Notes of the Les Houches Summer School: Volume 93*, August 2009 (Oxford University Press, 2011) pp. 1–91.

[28] H. J. Rothe, *Lattice Gauge Theories: An Introduction*, 4th Edition (World Scientific, 2012).
# Supplemental Material

In the first section, we translate, on the transport coin operators $B, V$ and $M$, the unitarity constraints imposed on the jump coin operators $W_i$, $i = -1, 0, 1$, by the unitarity of the one-step evolution operator $U$. In the second section, we explain the fermion-doubling problem of lattice gauge theory, and present Wilson's method to solve it, with so-called Wilson fermions.

## CONTENTS

I. Unitarity constraints 7

II. Fermion doubling and Wilson fermions in continuous and discrete time 8

A. The Schrödinger equation for translationally invariant systems 8

1. The generic Schrödinger equation and its spectral solution 8

2. Fourier analysis 9

3. Case of translationally invariant systems 9

B. The solution of the Schrödinger equation for translationally invariant systems: A superposition of plane waves 9

1. No internal structure for $\Psi(t, x)$ 9

2. Internal structure for $\Psi(t, x)$ 10

C. Dirac fermions (continuous spacetime) 11

D. The doubling problem when discretizing space but keeping time continuous: Spatial doublers 11

E. The doubling problem in discrete spacetime: Temporal doublers in addition to the spatial doublers 12

F. Removing the doublers with Wilson fermions 13

1. We limit ourselves to spatial doublers 13

2. Treatment of spatial doublers in the naive continuous-time scheme, via Wilson fermions 14

## I. UNITARITY CONSTRAINTS

We impose the unitarity of $U$, that is, $U^\dagger U = 1$. A couple of computation lines lead to the following conditions on the jump coin operators $W_i$, $i = -1, 0, 1$, that we call unitarity constraints on the jump coin operators:

$$W_{-1}^{\dagger} W_{-1} + W_{1}^{\dagger} W_{1} + W_{0}^{\dagger} W_{0} = 1 \qquad (15a)$$

$$W_{-1}^{\dagger} W_{0} + W_{0}^{\dagger} W_{1} = 0 \qquad (15b)$$

$$W_{-1}^{\dagger} W_{1} = 0. \qquad (15c)$$

Let us translate these constraints on the transport coin operators $B, V$ and $M$. We start by Constraint (15c), which is the simplest to translate in the sense that it does not involve $M$ while both Constraints (15a) and (15b) do. Constraint (15c) yields

$$0 = (V^{\dagger} - B^{\dagger})(V + B) \qquad (16a)$$

$$= V^{\dagger}V - B^{\dagger}B - B^{\dagger}V + V^{\dagger}B. \qquad (16b)$$

Computing the sum and the difference, $(16b)^{\dagger} + (16b)$ and $(16b)^{\dagger} - (16b)$ yields, respectively,

$$V^{\dagger}V = B^{\dagger}B \qquad (17a)$$

$$B^{\dagger}V = V^{\dagger}B. \qquad (17b)$$

Conversely, it is trivial to check that these two constraints imply Constraint (15c).

We now proceed with translating Constraint (15b), because we will use one of the two resulting constraints to translate Constraint (15a). Constraint (15b) yields

$$0 = (V^{\dagger} - B^{\dagger})(M - V) + (M^{\dagger} - V^{\dagger})(V + B) \qquad (18a)$$

$$= V^{\dagger}M - B^{\dagger}M - V^{\dagger}V + B^{\dagger}V + M^{\dagger}V - V^{\dagger}V + M^{\dagger}B - V^{\dagger}B. \qquad (18b)$$
Computing the sum and the difference (18b)$^\dagger$ + (18b) and (18b)$^\dagger$ - (18b), and using Eq. (17b) in the difference, yields, respectively,

$$2V^\dagger V = V^\dagger M + M^\dagger V \qquad (19a)$$

$$B^\dagger M = M^\dagger B. \qquad (19b)$$

Conversely, it is trivial to check that these two constraints together with Constraint (17b) imply Constraint (15b). We finally proceed with translating Constraint (15a), which can be written

$$4 = (V^\dagger - B^\dagger)(V - B) + (V^\dagger + B^\dagger)(V + B) + 4W_0^\dagger W_0 \qquad (20a)$$

$$= V^\dagger V + B^\dagger B - B^\dagger V - V^\dagger B + V^\dagger V + B^\dagger B + B^\dagger V + V^\dagger B + 4W_0^\dagger W_0. \qquad (20b)$$

Inserting both Constraint (17a) and Constraint (17b) in the preceding one, Constraint (20b), delivers

$$1 = V^\dagger V + W_0^\dagger W_0 \qquad (21a)$$

$$= V^\dagger V + (M^\dagger - V^\dagger)(M - V) \qquad (21b)$$

$$= 2V^\dagger V + M^\dagger M - M^\dagger V - V^\dagger M. \qquad (21c)$$

Finally, inserting Constraint (19a) into the preceding one, Constraint (21c), results in

$$M^\dagger M = 1. \qquad (22)$$

## II. FERMION DOUBLING AND WILSON FERMIONS IN CONTINUOUS AND DISCRETE TIME

### A. The Schrödinger equation for translationally invariant systems

#### 1. The generic Schrödinger equation and its spectral solution

The generic Schrödinger equation is a partial differential equation (PDE) of the form

$$i\partial_0 \Psi|_t = h\Psi(t), \qquad (23)$$

where $h$ is a Hermitean linear operator acting on the function $\Psi(t) : x \mapsto \Psi(t, x)$, where here we choose $x \in \mathbb{R}$. Since this equation is linear, we solve it spectrally, i.e., by finding the eigen-elements $(\omega_\sigma, \Phi_\sigma)_{\sigma \in \Sigma}$ of $h$, where $\Sigma$ is a certain indexing space; By definition, these eigen-elements satisfy

$$h\Phi_\sigma = \omega_\sigma \Phi_\sigma, \qquad (24)$$

where the eigen-values $\omega_\sigma$ are real since $h$ is Hermitean.

The method is the following. Assume we have determined the eigen-elements of $h$. Since $\Psi(t)$ belongs to a Hilbert space, we can decompose it on the eigen-basis $(\Phi_\sigma)_{\sigma \in \Sigma}$ at an arbitrary time $t$:

$$\Psi(t) = \sum_{\sigma \in \Sigma} C_{\sigma}(t) \Phi_{\sigma}, \qquad (25)$$

where the $C_\sigma(t)$s are the coefficients of $\Psi(t)$ on the eigen-basis. Now, using Eq. (24), the generic Schrödinger equation on $\Psi$, Eq. (23), which is a PDE, can be translated into a family of ordinary differential equations (ODEs) – indexed by $\sigma$ – for the coefficients $C_\sigma$, that is,

$$i\partial_0 C_\sigma|_t = \omega_\sigma C_\sigma(t), \qquad (26)$$

whose solution is well-known:

$$C_{\sigma}(t) = C_{\sigma}(0)e^{-i\omega_{\sigma}t}. \qquad (27)$$

Hence, the solution seeked is

$$\Psi(t) = \sum_{\sigma \in \Sigma} C_{\sigma}(0)e^{-i\omega_{\sigma}t}\Phi_{\sigma}. \qquad (28)$$

Because they intervene in the periodic functions $t \mapsto e^{-i\omega_\sigma t}$, the $\omega_\sigma$ are called *frequencies*; More precisely, they are the *eigen-frequencies of h*. To be more definite, one should actually use the denomination “angular frequency” rather than “frequency”.
## 2. Fourier analysis

We will use the more definite notation $\Psi(t, \cdot)$ for $\Psi(t)$ when needed. Let us take the Fourier transform of $\Psi(t, \cdot)$ at a given time $t$:

$$ \tilde{\Psi}(t, k) := \frac{1}{\sqrt{2\pi}} \int_{\mathbb{R}} dx \, \Psi(t, x) e^{-ikx}. \qquad (29) $$

Inverting this equation, we obtain the decomposition of $\Psi(t)$ into its Fourier components:

$$ \Psi(t, x) = \frac{1}{\sqrt{2\pi}} \int_{\mathbb{R}} dk \, \tilde{\Psi}(t, k) e^{ikx}. \qquad (30) $$

To be precise, the function $x \mapsto \tilde{\Psi}(t, k)e^{ikx}$ is the Fourier component of $\Psi(t, \cdot)$ associated to the value $k$ of the Fourier variable, and $\tilde{\Psi}(t, k)$ is the Fourier coefficient, or Fourier amplitude of $\Psi(t, \cdot)$ associated to the value $k$. Because $x$ is a spatial position, the Fourier variable $k$ is a spatial frequency; Again, as in the case of $\omega_\sigma$ above, to be more definite one should use the denomination “angular pulsation” rather than “frequency”.

## 3. Case of translationally invariant systems

If $h$ does not depend on the point $x$, i.e., if $h$ is translationally invariant, then one can check by considering Eq. (23) for Expression (30) that each Fourier coefficient satisfies the equation

$$ i\partial_0 \tilde{\Psi}(\cdot, k)|_t = \tilde{h}(k) \tilde{\Psi}(t, k), \qquad (31) $$

where $\tilde{h}(k)$ is the expression obtained when replacing, in $h$, the operator $-i\partial_1$ by the real number $k$.

Hence, if $h$ does not depend on the point $x$, each Fourier coefficient evolves independently of the others, while this is not the case if $h$ does depend on $x$. Moreover, the “Schrödinger equation in Fourier space”, Eq. (31), is simpler to solve than the original Schrödinger equation, (23), because the operator $-i\partial_1$ has been replaced by a real number $k$, so that Eq. (31) is not a PDE as the original Schrödinger equation, but a family of ODEs indexed by $k$. Let us now make the link with Sec. II A: One can actually mathematically show (via, e.g., “mere” constructive proofs) for a large class of operators $h$, that there exist an indexing space $\Sigma$ such that $k$ is one of the indices, i.e., $k \in \sigma^3$. In the language of quantum mechanics, and index $i \in \sigma$ is refered to as a good quantum number; It is an eigenvalue of an operator whose diagonalization serves as a partial diagonalization of $h$, i.e., $h$ is codiagonalizable with that operator.

We often speak of Fourier modes for the Fourier components; The term “mode” refers, in its most general acception, to one of the terms of a particularly relevant decomposition of an object, take, e.g., the Fourier decomposition of a function. In the case of an $x$-independent $h$, the Fourier modes are actually proper modes, because by considering the Fourier version $\tilde{h}(k)$ of $h$ we have (at least partially) “diagonalized” $h$ (which is summed up by writing $k \in \sigma$).

## B. The solution of the Schrödinger equation for translationally invariant systems: A superposition of plane waves

### 1. No internal structure for $\Psi(t, x)$

**a. Final solution.** If $\Psi(t, x)$ has no internal structure, i.e., if $\Psi(t, x) \in \mathbb{C}$, then we simply have that $k = \sigma$, and the eigenvalues are $\omega_\sigma = \omega(k) := \tilde{h}(k) \in \mathbb{R}$; The sum over $\sigma$ in Eq. (28) is an integral over $k$, and the eigen-basis is the following family of functions of $x$, $(\Phi(\cdot, k))_{k \in \mathbb{R}}$, with

$$ \Phi(x, k) := \frac{e^{ikx}}{\sqrt{2\pi}}, \qquad (32) $$

so that the solution given in Eq. (28) here reads

$$ \Psi(t, x) = \frac{1}{\sqrt{2\pi}} \int_{\mathbb{R}} dk \, C(0, k) e^{-i(\omega(k)t - kx)}. \qquad (33) $$

<footnote n="3">³ We take the liberty to consider $\sigma$ as convenient; Here, it is a set of indices, i.e., an unordered family; In practice, we will often order its elements, i.e.,</footnote>

<footnote>we will consider a list rather than a set.</footnote>
By taking $t = 0$ in this equation we realize by identification (they are unique) that the $C(0, k)$ are the Fourier coefficients of $x \mapsto \Psi(0, x)$, i.e.,

$$C(0, k) \equiv \tilde{\Psi}(0, k). \qquad (34)$$

The value $\omega(k)$ is the eigen-frequency associated to the spatial frequency $k$. Now, an important remark is that $\Psi(t, \cdot)$ is actually a supersposition of plane waves, with weights the Fourier coefficients of the initial condition. Because of this (plane-)wave structure, the spatial frequency $k$ is called *wavevector*<fnref n="4" />. The fact that the frequency of the wave $\omega(k)$ depends on the wavevector $k$ is called *dispersion*, and the expression $\omega(k)$ is called the *dispersion relation*. Notice that the fact that we have a dispersion phenomenon while we are in vacuum is specific to quantum mechanics, more precisely, to the quantum mechanics of massive bodies.

In quantum mechanics, the use of the word “wavevector” is actually extended to non translationally invariant systems, because this wavevector is, due to the wave-particle duality, in one-to-one correspondence with the momentum of the free particle associated to the wave in question, solution of a free (i.e., with $x$-independent $\hbar$) Schrödinger equation: More precisely, $k$ is the De Broglie wavevector of a particle of momentum $p = k$ (with $\hbar = 1$). There is also an analog relation for the eigen-frequencies, refered to as Einstein’s relation: The frequency $\omega_\sigma$ is the frequency associated to an energy $E_\sigma = \omega_\sigma$ (with $\hbar = 1$) for the particle in question.

b. *Final solution rederived by focusing on the linear algebraic structure.* It is the opportunity to rederive the above solution, Eq. (33), by focusing on the linear algebraic structure of the computations. This will make gentler the above replacement of $\sigma$ by $k$. We first rewrite Eq. (28) with a braket notation, which emphasizes the linear algebraic structure,

$$|\Psi(t)\rangle = \sum_{\sigma \in \Sigma} \underbrace{\langle \Phi_{\sigma} | \Psi(0) \rangle}_{\equiv C_{\sigma}(0)} e^{-i\omega_{\sigma}t} |\Phi_{\sigma}\rangle. \qquad (35)$$

By taking $t = 0$ in this equation, we see that we have simply applied the closure relation $\sum_{\sigma \in \Sigma} |\Phi_{\sigma}\rangle\langle\Phi_{\sigma}|$ to $|\Psi(0)\rangle$, and then evolved the eigen-states up to time $t$. We then consider the present case $\sigma = k$, so that $\Phi_{\sigma} = \Phi(\cdot, k)$, which means that the eigen-basis is the basis made up of the $|k\rangle := |\Phi(\cdot, k)\rangle$ with $k \in \mathbb{R}$,

$$|\Psi(t)\rangle = \int_{\mathbb{R}} dk \underbrace{\langle k | \Psi(0) \rangle}_{\equiv \tilde{\Psi}(0,k)} e^{-i\omega_{\sigma}t} |k\rangle. \qquad (36)$$

We then apply the bra $\langle x|$ to have the value at point $x$,

$$\langle x | \Psi(t) \rangle = \int_{\mathbb{R}} dk \tilde{\Psi}(0, k) e^{-i\omega_{\sigma}t} \langle x | \Phi(\cdot, k) \rangle, \qquad (37)$$

that is, using Eq. (32), exactly Eq. (33),

$$\Psi(t, x) = \frac{1}{\sqrt{2\pi}} \int_{\mathbb{R}} dk \tilde{\Psi}(0, k) e^{-i(\omega(k)t - kx)}. \qquad (38)$$

## 2. Internal structure for $\Psi(t, x)$

Now, if $\Psi(t, x)$ has an internal structure, then $k \subsetneq \sigma$<fnref n="5" />, and $\tilde{h}(k)$ can be seen as a matrix indexed by $k$, that one has to diagonalize to finally find the eigen-values of $\tilde{h}$. The eigen-values of $\tilde{h}(k)$ can be denoted $\omega_\sigma = \omega_i(k)$, where $i = 1, ..., d$, with $d$ the dimension of the matrix $\tilde{h}(k)$ (some eigen-values may be equal, e.g., $\omega_{i_a}(k) = \omega_{i_b}(k)$). The final solution is then

$$\Psi(t, x) = \sum_{i=1}^{d} \frac{1}{\sqrt{2\pi}N} \int_{\mathbb{R}} dk \tilde{\Psi}_{i}(0, k) e^{-i(\omega_{i}(k)t - kx)}, \qquad (39)$$

where $N$ is a normalization factor, needed if we want the eigenvectors $\tilde{\Psi}_i(0, \cdot)$ of $\tilde{h}(k)$ to be normalized. We see that the solution is still a superposition of plane waves, and there are $d$ dispersion relations $\omega_i(k)$, $i = 1, ..., d$.

<footnote n="4">⁴ The question that one may ask oneself is whether there exist a “physical” $\tilde{h}$ that is both linear and Hermitean but that does not lead to an ondulatory</footnote>

<footnote>phenomenology.</footnote>

<footnote n="5">⁵ We have identified $k$ with the singlet $\{k\}$.</footnote>
### C. Dirac fermions (continuous spacetime)

The Dirac Hamiltonian in 1 + 1 dimensions is

$$ \hat{h}^{\text{Dirac}} := \alpha^1 \hat{k} + m\alpha^0, \qquad (40) $$

where $\hat{k}$ is the momentum operator, i.e., the abstract version of the operator $-i\partial_1$, and the alpha matrices satisfy $(\alpha^0)^2 = (\alpha^1)^2 = 1$ and $\alpha^0\alpha^1 + \alpha^1\alpha^0 = 0$. We see that $\hat{h}^{\text{Dirac}}$ is Hermitean, and so it is a valid Hamiltonian for the generic Schrödinger equation considered in Sec. II A and II B above. Moreover it is translationally invariant.

We consider the “Hamiltonian in momentum space”,

$$ \tilde{h}^{\text{Dirac}}(k) := \alpha^1 k + m\alpha^0 \qquad (41) $$

The eigen-value equation with unknowns the eigen-elements $(\omega_i(k), V_i(k))_{i=1,...,d; k \in \mathbb{R}}$ of $\tilde{h}^{\text{Dirac}}(k)$ is, in matrix notation,

$$ \tilde{h}^{\text{Dirac}}(k)V_i(k) = \omega_i(k)V_i(k). \qquad (42) $$

Now, to find the eigen-values of $\tilde{h}^{\text{Dirac}}(k)$, there is actually a “trick”, related to the fact that the square of the Dirac equation is, by (historical) construction of the Dirac equation, the Klein-Gordon equation, and hence applicable to scalar state functions: The square $(\tilde{h}^{\text{Dirac}}(k))^2$ is proportional to the identity matrix. By squaring Eq. (42), we arrive to

$$ \omega(k)^2 = k^2 + m^2, \qquad (43) $$

so that the eigen-values are

$$ \omega_{\pm}^{\text{Dirac}}(k) := \pm \sqrt{k^2 + m^2}. \qquad (44) $$

### D. The doubling problem when discretizing space but keeping time continuous: Spatial doublers

To discretize space, we can simply perform the naive replacement of the partial derivative $\partial_1$ by a finite difference on a 1D lattice that we introduce, with sites labeled by $p \in \mathbb{Z}$ and lattice spacing $a$. The finite difference has to be symmetric if we want the resulting Hamiltonian to be Hermitean. Since the translation operator in the direction of growing $ps$ is, in abstract space, $\hat{T} = e^{-ika} {}^6$, the substitution of $\partial_1 = i\mathcal{K}$, where $\mathcal{K} := -i\partial_x$, by the announced symmetric finite difference,

$$ \mathcal{D}_1 := \frac{1}{2}(\mathcal{T}_1^{-1} - \mathcal{T}_1), \qquad (45) $$

corresponds to the following substitution in the “Dirac Hamiltonian in momentum space” $\tilde{h}^{\text{Dirac}}$, as well as in the dispersion relation, Eq. (44):

$$ k \to (-i) \frac{e^{ika} - e^{-ika}}{2a} = \frac{\sin(ka)}{a}. \qquad (46) $$

To explain the doubling problem, it is customary to consider the following function of $k$,

$$ f^{\text{naive}}(k) := \sqrt{(\omega_{\pm}^{\text{naive}}(k))^2 - m^2} \equiv \left| \frac{\sin(ka)}{a} \right|, \qquad (47) $$

to be compared with

$$ f^{\text{Dirac}}(k) := \sqrt{(\omega_{\pm}^{\text{Dirac}}(k))^2 - m^2} \equiv |k|. \qquad (48) $$

Notice first that the spatial discretization implies that now $k \in [-\pi/a, \pi/a)$. Second, notice that we of course recover the continuum situation for $ka \ll 1$, because $\sin(ka) = ka + O((ka)^3)$. We call the function $f^M(k)$ the gapless frequency of Model M.

<footnote n="6">⁶ Pay attention to the fact that, although we use the notation $\hat{k}$, on the lattice $\hat{k}$ and the position operator $\hat{x}$ do not satisfy the canonical commutation relations.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2. This figure illustrates the fermion doubling problem in discrete space but continuous time. The lattice dispersion relation under the form $f_{\text{naive}}(k)$ (in gold color) has 3 poles, $-\pi/a$, $0$, $\pi/a$, instead of a single one, $0$, for the continuum dispersion relation (in blue color), so that for each target value $\phi_0$ of the gapless frequency, we have two extra momentum solutions on the lattice with respect to the continuum situation.</figcaption></figure>

In Fig. 2, we plot both $f^{\text{Dirac}}(k)$ and $f^{\text{naive}}(k)$ over the Brillouin zone $[-\pi/a, \pi/a)$. The doubling problem is the following. For a given target value $\phi_0 := \sqrt{\omega_0^2 + m^2}$ of the gapless frequency, there are, in the naive discretization, not 2 possibilities for the momentum as in the continuum situation, $k_0$ and $-k_0$ such that $f^{\text{Dirac}}(k_0) = \omega_0$, i.e., $k_0 = \omega_0$, but 4 solutions, 2 corresponding to the low-momentum modes that we seek to simulate with the discretization, which have $k'_0 \simeq k_0$ and $-k'_0$ such that $f^{\text{naive}}(k'_0) = \phi_0$, and 2 additional, high-momentum modes, namely, $\pi/a - k'_0$ and $-(\pi/a - k'_0)$, so that frequencies and momenta are not in one-to-one correspondence anymore. In a non-interacting model, i.e., if $h$ does not depend on $x$, this is actually not a problem because the Fourier modes are independent from each other, and so the one-to-one correspondence between frequency and momentum can be tracked, e.g., fundamentally, as time evolves; More precisely and concretely: In a non-interacting model, the momentum distribution is unchanged by the dynamics, i.e., in other words, the states of fixed momentum are stationary states. In an interacting model, the Fourier modes will not evolve independently from each other, and the interaction term will cause the production of high-momentum modes from low-momentum ones – because of the 2 extra poles of $f^{\text{naive}}(k)$ with respect to $f^{\text{Dirac}}(k)$ – which can be proved rigorously in, e.g., perturbative studies of interacting models having as zeroth order $f^{\text{naive}}(k)$ [28?].

As a conclusion: In discrete space (but keeping time continuous), we will have extra, spurious modes when looking for superpositions of plane waves as Eq. (39) for the solutions. These spurious modes are called *spatial doublers*, where the specification spatial is due to the fact that what is spurious in these modes is the spatial part (high momenta, not compatible with a continuum description, even if the temporal part is – i.e., low frequencies).

### E. The doubling problem in discrete spacetime: Temporal doublers in addition to the spatial doublers

We start from the discrete-space but continuous-time situation described just above in Sec. II F, that is

$$i\partial_0\Psi(\cdot, x)|_t = h^{\text{naive}}\Psi(t, \cdot)|_x. \qquad (49)$$

As we did for space in Sec. II F, we discretize time naively, with a symmetric finite difference, in order to treat time on the same footing as space (for which we have indeed used a symmetric finite difference), which yields (replacing $h^{\text{naive}}$ by its expression):

$$\frac{i}{2\epsilon}(\Psi(t+\epsilon) - \Psi(t-\epsilon)) = \left[ \frac{-i}{2a}(e^{i\kappa a} - e^{-i\kappa a})\alpha^1 + m\alpha^0 \right] \Psi(t, \cdot)|_x. \qquad (50)$$

Notice right away that this scheme takes two initial conditions, exactly as the two-step scheme we present in our paper; The only difference with the two-step scheme we present is that in the latter there is a factor $\mu(\epsilon)$ in front of $m\alpha^0$, but this is enough to
make that scheme unitary, while the present, naive one, is not$^7$. We consider an ansatz which is a superposition of plane waves with internal components, that is, we look for solutions of the form of Eq. (39). If we insert Eq. (39) into Eq. (50), we obtain after a few computation lines the following equation in momentum space:

$$ \frac{\sin(\omega_i(k)\epsilon)}{\epsilon} \tilde{\Psi}_i(0, k) = \left[ \frac{\sin(\omega_i(k)a)}{a} \alpha^+ m \alpha^0 \right] \tilde{\Psi}_i(0, k), \qquad (51) $$

an equation which, squared, and choosing the ballistic scaling $\epsilon = a$, finally yields the following dispersion relation,

$$ \sin^2(\omega_i(k)\epsilon) = \sin^2(k\epsilon) + \epsilon^2 m^2. \qquad (52) $$

First of all, notice that for low frequencies $\omega_i(k)\epsilon \ll 1$ in Eq. (52), we recover the discrete-space but continuous-time situation, with gapless frequency given by Eq. (47) in the preceding section, Sec. II F. Now, there are solutions $\omega_i(k)$ to Eq. (52) if and only if $|\sin^2(k\epsilon) + \epsilon^2 m^2| \le 1$, which leads to

$$ \epsilon^2 m^2 \le \cos^2(k\epsilon). \qquad (53) $$

Replacing in the dispersion relation, Eq. (52), $\sin^2(A)$ by $(1 - \cos(2A))/2$, we obtain the two following solutions,

$$ \omega_{\pm}^{\text{temporal}}(k) := \pm \frac{2}{\epsilon} \arccos \left( 1 - 2 \sin^2(k\epsilon) - 2\epsilon^2 m^2 \right). \qquad (54) $$

For small enough $k$ and $m$, i.e., $k\epsilon \ll 1$ and $\epsilon m \ll 1$, we have that $|\omega_{\pm}^{\text{temporal}}(k)\epsilon| < \pi/2$, and actually that $|\omega_{\pm}^{\text{temporal}}(k)\epsilon| \ll 1$, and these two solutions approach the low frequencies $\omega_{\pm}^{\text{Dirac}}(k)$ of the continuum model.

Now, in addition to these two solutions, Eq. (54), we also have the solutions,

$$ \Omega_{\pm}^{\text{temporal}}(k) := \pm \left[ \frac{\pi}{\epsilon} - \omega_{\pm}^{\text{temporal}}(k) \right], \qquad (55) $$

which are high-frequency solutions when $|\omega_{\pm}^{\text{temporal}}(k)\epsilon| < \pi/2$. These two extra solutions are spurious because not compatible with a continuum description, but they will intervene in the dynamics in interacting models, and the modes associated to these solutions are called temporal doublers. As we have seen, these temporal doublers arise even for low momenta$^8$ (the only ones compatible with a continuum description), and this is best seen as follows. Consider low momenta in Eq. (52); This yields, replacing $k$ by the notation $\kappa(w)$ and $\omega_i(k)$ by the variable $w$,

$$ g^{\text{naive}}(w) := \sqrt{\kappa(w)^2 + m^2} = \left| \frac{\sin^2(w\epsilon)}{\epsilon} \right|, \qquad (56) $$

an expression which, apart from the fact that there is a $+m^2$ instead of a $-m^2$, corresponds exactly to the expression of the gapless frequency, Eq. (47), but having exchanged in it the roles of $\omega(k)$ and $k$, i.e., replaced the latter by $\kappa(w)$ and $w$, respectively, so that one can derive the same explanations for the temporal doublers than for the spatial doublers. Notice that, in this low-momentum context, the condition (53) for a solution $\omega_i(k)$ to exist is

$$ \epsilon^2 m^2 \le 1 - \frac{1}{2}(k\epsilon)^2. \qquad (57) $$

## F. Removing the doublers with Wilson fermions

### 1. We limit ourselves to spatial doublers

In the preceding section, Sec. II E, we have illustrated the problem of temporal doublers starting from a continuous-time description. A framework which is more appropriate to further remove temporal doublers is that of Lagrangian LGT, which starts from a Lagrangian continuum description rather than a Hamiltonian one. This leads to a modification of the action rather than the Hamiltonian to solve the problem of fermion doubling, with a Wilson term, and this procedure removes the spatial doublers, as well as the temporal ones for Wilson's choice $r = 1$ [27]. Here, however, we will stick to a Hamiltonian formulation and treat only spatial doublers. This is because in our paper the scheme is formulated in a Hamiltonian way. To treat temporal doublers in the two-step scheme of our paper one would have to modify the original, DQW scheme.

<footnote>$^7$ The present scheme is indeed not unitary in the following sense. Recall first that for the unitarity of the two-step scheme of our paper has the following meaning: That two-step scheme is equivalent to the one-step scheme,</footnote>

<footnote>for which the unitarity meaning is clear, provided that $\Psi_{j=1} = U\Psi_{j=0}$. In the present case, endowing the naive two-step scheme with the same unitarity meaning does not yield a unitary one-step scheme.</footnote>

<footnote>$^8$ We here speak of the momentum of the discrete model, i.e., the product $\epsilon k$.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 3. Gapless frequency $f^M(k)$ of Model “M” (left panel) and absolute difference $f^M(k) - f^{\text{Dirac}}(k)$ (right panel), for “M” = “Dirac” (blue), “naive” (gold), “LGT” (green), and “LGT non-crossed” (red), where $f^{\text{Dirac}}(k) = |k|$, $f^{\text{naive}}(k) = \frac{|\sin(ka)|}{a}$, $f^{\text{LGT}}(k) = \sqrt{\frac{\sin^2(ka)}{a^2} + [m + \frac{r}{a}(1 - \cos(ka))]^2 - m^2}$, and $f^{\text{LGT non-crossed}}(k) = \sqrt{\frac{\sin^2(ka)}{a^2} + \frac{r^2}{a^2}(1 - \cos(ka))^2}$, for $a = 1$, $m = 0.1$, and $r = 1$ (Wilson’s choice).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. In blue, $(f^{\text{Dirac}}(k))^2$, in gold, $(f^{\text{naive}}(k))^2$, in green, the non-crossed term $g^{\text{non-crossed}}(k)$, and in red, the crossed term $g^{\text{crossed}}(k)$. We see that the biggest contribution to the fixing of the fermion doubling comes from the non-crossed term (with respect to the crossed term). That being said, any of the two terms is actually sufficient on its own to fix the doubling problem.</figcaption></figure>

## 2. Treatment of spatial doublers in the naive continuous-time scheme, via Wilson fermions

In order to treat the problem of spatial doublers in continuous-time, described in Sec. IIF, we consider the following Hamiltonian [26],

$$h^{\text{LGT}} := h^{\text{naive}} + h^{\text{Schrö.}}, \qquad (58)$$

where

$$h^{\text{naive}} := \alpha^1 \left( -i \frac{D_1}{a} \right) + m\alpha^0, \qquad (59)$$

is the naive-discretization Hamiltonian considered in Sec. IIF, and

$$h^{\text{Schrö.}} := \alpha^0 \frac{r}{a} (-\mathcal{L}), \qquad (60)$$

is the Wilson term, whose spatial operator, $\mathcal{L}/a$, is a discrete Laplacian<fnref n="9" />,

$$\mathcal{L} = \mathcal{T}_1^{-1} + \mathcal{T}_1 - 2. \qquad (61)$$

<footnote n="9">⁹ This discrete Laplacian has been obtained by replacing $a\partial_1^2$ by $\mathcal{L} \equiv d_1^* d_1$, where $d_1$ and $d_1^*$ are respectively a forwards and a backwards finite difference (let us notice that they commute).</footnote>
Let us notice that the Laplacian is the spatial operator that intervenes in the non-relativistic Schrödinger equation, hence the superscript “Schrö.”.

The dispersion relation is now, plotted under the form of a gapless frequency,

$$f^{\text{LGT}}(k) := \sqrt{(\omega_{\pm}^{\text{LGT}})^2 - m^2} = \sqrt{\frac{\sin^2(ka)}{a^2} + \left[m + \frac{r}{a}(1 - \cos(ka))\right]^2 - m^2}. \qquad (62)$$

In the left panel of Fig. 3, we see that the doubling problem is fixed with this expression, since there are no more poles at $-\pi/a$ and $\pi/a$. To develop explanations, it is practical to consider the square of the previous expression,

$$(f^{\text{LGT}}(k))^2 = g^{\text{naive}}(k) + g^{\text{non-crossed}}(k) + g^{\text{crossed}}(k). \qquad (63)$$

where

$$g^{\text{naive}}(k) := \frac{\sin^2(ka)}{a^2} = (f^{\text{naive}}(k))^2 \qquad (64a)$$

$$g^{\text{non-crossed}}(k) := \frac{r^2}{a^2} (1 - \cos(ka))^2 \qquad (64b)$$

$$g^{\text{crossed}}(k) := 2m\frac{r}{a}(1 - \cos(ka)). \qquad (64c)$$

Let us develop this expression, Eq. (63), at next-to-next-to-lowest order in $a$,

$$(f^{\text{LGT}}(k))^2 = (1 + amr)k^2 + \left(\frac{1}{4}a^2r^2 - \frac{2}{3!}a^2\right)k^4. \qquad (65)$$

In Fig. 4, we see that the biggest contribution in fixing the doubling problem comes from the non-crossed term; This is actually also visible in the left panel of Fig. 3. That being said, any of the two terms is actually sufficient on its own to fix the doubling problem. Moreover, we see in Eq. (65) that the crossed term unfortunately adds a first-order correction (in the lattice spacing) to the initial slope, that is, the latter becomes $1 + amr$ instead of 1, which is visible in the right panel of Fig. 3. Now, one can actually make this crossed term disappear, by choosing in the Wilson term $h^{\text{Schrö}}$, Eq. (60), the operator $\alpha^2$ instead of $\alpha^0$. While in continuous time this has no impact on the unitarity of the scheme, it will a priori have an impact in discrete time [27], so this replacement cannot be made carelessly.