\section{Introduction}

Surface acoustic waves (SAWs) \cite{Landau1986} have been used for decades as a valuable scientific and technological tool.
In the context of electronics,
they are often excited in piezoelectric materials
\cite{Hutson1962,Morgan2010,Auld1990,Royer2000,Royer2000a,Newnham2004},
where the mechanical and electrical fields are coupled. In particular, they have been applied as experimental probes of the quantum Hall effects in
two-dimensional electron gases (2DEG) \cite{Wixforth1986}.

On the other hand, there is an increased interest in 2DEG since the isolation of graphene in 2004 and the production of other two-dimensional (2D) materials which followed it.
Due to their unusual character, the properties of graphene electrons have been intensively studied during the last decade
\cite{CastroNeto2009}. Although graphene on a substrate has received considerable attention, relatively few studies have been devoted to the case of graphene in contact with a piezoelectric material.
These include the propagation of surface acoustic waves on graphene \cite{Thalmeier2010}, some acoustoelectric effects \cite{Miseikis2012} \cite{Bandhu2013}, and the relaxation induced by the
surface acoustic wave quanta on graphene electrons \cite{Zhang2013}, among others.
Recently, it has been proposed that surface acoustic waves (SAW) can provide a
diffraction grating for the conversion of light into graphene plasmons \cite{Schiefele2013}.

The coupling of the piezoelectric SAW to the electrons in
a 2DEG or in graphene has been computed, within certain simplifying assumptions and for definite substrate crystal structures, in Ref. \cite{Simon1996,Knabchen1996}. The derivation of the electron-surface phonon interaction for a piezoelectric material has been performed only within a purely elastic Rayleigh wave approximation \cite{Knabchen1996} or for definite propagating directions \cite{Simon1996}. But these methods fail in stronger piezoelectric materials and for other crystal symmetries. For instance, the isotropic Rayleigh wave approximation in the case of lithium niobate leads to surface acoustic wave velocities about 15 \% too low and lacking the correct angular dependence \cite{Morgan2010,Farnell1970}, and this material is not the one with the largest electromechanical couplings at all. Moreover, the obtained vertex are expressed in terms of a matching constant whose physical interpretation is rather obscure, allowing for just order of magnitude estimates.

In the present work we calculate a general electron-phonon interaction [see Eq. (\ref{eqn:PAVertex})], which is expressed solely in terms of physical quantities characterizing the response of the substrate surface. We emphasize that all the quantities appearing in the vertex are both computable from linear piezo-elasticity theory and experimentally measurable. One of them, the {\it electromechanical coupling coefficient}, $K_R$, will turn out to be central to all computations, serving as a natural dimensionless parameter which provides the scale for the effect of the substrate piezoelectricity on the 2D electron system. Moreover, from very general considerations explained in Appendix \ref{app:PiezoElectricVertex}, we are able to provide bounds on its size: $0\le K_R<1$ [see Eq. (\ref{eqn:interactionsRatio})]. It is important to note that the vertex written here is derived within the framework of linear piezo-elastic theory, which means that its validity should be restricted to low amplitude, low frequency and long wavelength phenomena. Bulk modes are also left aside in this work. On the other hand, our study is not restricted to any approximation based on the symmetry or the piezoelectric softness of the substrate.

Equipped with the effective electron-electron interaction which results from taking into account the exchange of these acoustic phonons between the electrons in graphene (or other 2D materials), the question can be raised of whether these interactions might be attractive and, depending on some material parameters and the tunable electronic density of graphene, perhaps strong enough to generate electron Cooper pairing and superconductivity \cite{Mahan2013}. We can further ask whether such a superconductivity could be observed at temperatures attainable in a laboratory without the recurring to huge non gate-achievable doping levels, as predicted for intrinsic graphene phonons \cite{Einenkel2011} or for Kohn-Luttinger or electronic superconductivity in other graphene heterostructures with repulsive interactions \cite{Guinea2012}. From the interaction vertex derived in the present work, it can be shown that the relative size of the static phonon-mediated electron-electron interaction with respect to the original Coulomb repulsion turns out to be exactly $K_R^2<1$, because of the aforementioned general inequality. However, by applying the Eliashberg formalism to graphene \cite{Einenkel2011}, we are able to assess, in terms of $K_R$, the influence that these low-frequency and long-wavelength phonons have on possible BCS type instabilities. The conclusion is that  present piezoelectric materials are not able to either induce $s-$wave pairing by themselves or affect in a significant way any pairing instability which could be already present in graphene. We note, however, that this conclusion could be substantially changed in case new hard piezoelectric materials would be found.

\section{Electron-phonon interaction}
\begin{figure}
	\centering
	\includegraphics[width=0.9\linewidth]{./PiezoSAWSchema}
	\caption{A flat piezoelectric substrate with an interface to air or vacuum, with $\mathbf{b}=\left[\cos(\theta),\sin(\theta)\right]$ the propagation direction of the piezoelectric SAW.}
	\label{fig:PiezoSAWSchema}
\end{figure}

The system to be considered is depicted in Fig.~\ref{fig:PiezoSAWSchema}. The $z=0$ interfacial surface is supposed to be free of tension and, when acting as a substrate to a deposited graphene sheet (or any other 2D charged system), free of electrodes as well. However, by introducing a surface charge, one can express the 2D response of the piezoelectric substrate, the {\it piezoelectric surface permittivity}, as the ratio of an electric displacement to an electric field. To be precise, let us allow for a surface charge with a harmonic dependence along $\textbf{q}=q[\cos \theta,\sin\theta]$ (here $\textbf{R}=(x,y)$)
\begin{align}
	\sigma(\mathbf{R},t)=\sigma(\mathbf{q},\omega) e^{i (\textbf{q}\cdot\textbf{R}-\omega t)},
\label{eqn:surfaceChargeDef}
\end{align}
and, from linearity, all quantities evaluated at the surface have the same 2D space-time dependence. We use SI units throughout this work. Because the medium at $z>0$ has a dielectric constant $\varepsilon_\text{vac}$, then, from Poisson's equation it follows that
\begin{align}
	\varepsilon_{\text{vac}}= \frac{D_3^+(\mathbf{q},\omega)}{q \varphi(\mathbf{q},\omega)},
\label{eqn:vaccumPermittivity}
\end{align}
where $D_3^+$ is the normal to the surface Fourier component of the electric displacement field over the surface (on the $z \rightarrow 0^+$ side) and $\varphi$ is the electric potential, which is continuous because we do not allow for anything more singular than surface charges. The same ratio taken below the surface ($z \rightarrow 0^-$ ) allows us to introduce the (relative) {\it piezoelectric surface permittivity},
\begin{align}
	\tilde{\varepsilon}(\mathbf{q},\omega) := -\frac{1}{\varepsilon_{\text{vac}}}\frac{D_3^-(\mathbf{q},\omega)}{q \varphi(\mathbf{q},\omega)} \, ,
\label{eqn:piezoSurfacePermittivityDef}
\end{align}
which can also be straightforwardly expressed in terms of the surface impedance tensor \cite{Ingebrigsten1969}.
From Poisson's equation and Eqs. (\ref{eqn:surfaceChargeDef}-\ref{eqn:piezoSurfacePermittivityDef}) it follows
\begin{equation}\label{eqn:surfaceChargePotentialRelation}
	\sigma(\mathbf{q},\omega)=D_3^+(\mathbf{q},\omega)-D_3^-(\mathbf{q},\omega)=q \varphi(\mathbf{q},\omega) [1+\tilde{\varepsilon}(\textbf{q},\omega)]\varepsilon_\text{vac}.
\end{equation}

Further analysis summarized in Appendix \ref{app:PiezoElectricVertex} shows that $\tilde{\varepsilon}(\mathbf{q},\omega)$ has a dependence of the form
$\tilde{\varepsilon}(\mathbf{q},\omega)=\tilde{\varepsilon}(\mathbf{q}/\omega)$. An immediate conclusion from Eq. (\ref{eqn:surfaceChargePotentialRelation}) is that a purely piezoelectric wave (i.e., without sources, $\sigma(\mathbf{q},\omega)=0$), can propagate without damping if and only if $\tilde{\varepsilon}(\mathbf{q}/\omega)+1=0$. Thus, if the phase velocity is $v_s(\theta)$, then $\tilde{\varepsilon}(\mathbf{b}/v_s(\theta))+1=0$ and the dispersion relation of the obtained, called piezoelectric Rayleigh waves (here referred to as SAW) is $\omega=v_s(\theta) q$.

For the high-frequency limit we introduce $\tilde{\varepsilon}_{\text{HF}}(\theta):=\tilde{\varepsilon}(\mathbf{q}/\omega), \omega\rightarrow\infty$. This should be the anisotropic dielectric function valid all the way up to the optical region, and should take into account all screening processes in the substrate except for the slow piezoelectric ones, which are estimated below for the substrate of the 2D electronic material.

A central quantity in the evaluation of devices which use piezoelectric Rayleigh waves is the {\it SAW electromechanical coupling coefficient} $K_R(\theta)$, introduced through the relation at $1+\tilde{\varepsilon}(\textbf{b}/v_s(\theta))=0$:
\begin{equation}
	\frac{K_R^2(\theta)/2}{\tilde{\varepsilon}_{\text{HF}}(\theta)+1} = \left[v_s(\theta)\frac{\partial \tilde{\varepsilon}(\textbf{b}/v)}{\partial v}\Bigr|_{v=v_s(\theta)} \right]^{-1} ~.
\label{eqn:emCouplingCoefficient}
\end{equation}

In Appendix \ref{app:PiezoElectricVertex} we show that very general considerations require
\begin{equation} \label{eqn:central-inequality}
0 \le K_R(\theta)<1\, ,
\end{equation}
which is one of the central results of this work. The inequality $(\ref{eqn:central-inequality})$  is crucial because, as we will show, it implies that at small frequencies piezoelectric phonons cannot provide the sufficient screening to overcome the bare Coulomb repulsion.

In Ref. \cite{Royer2000} it is shown that there is a relation between the amplitude of electric potential at the surface, $\varphi_0=\varphi(\mathbf{q},\omega)$ and the total energy $H_{\text{harm}}$, see Eq. (\ref{eqn:HamiltonianHarmonicPSAW}). Hence, standard quantization procedure (see the Appendix, subsection \ref{app:HamiltonianInterVertex}) shows that the interaction between the 2D electronic material sheet and the spontaneous piezoelectric Rayleigh waves can be written as
\begin{align}
	H_{\text{e-ph}}^{\text{PA}} &= \frac{1}{\sqrt{A}} {\sum_{\mathbf{k,q},\sigma}}\gamma_{\mathbf{q}}^{\text{PA}} \, a_{\mathbf{k+q},\sigma}^{\dagger}a_{\mathbf{k},\sigma}b_{\mathbf{q}}+\text{H.c.}~,
	\nonumber \\
	\gamma^{\text{PA}}_\mathbf{q} & = \frac{K_R(\theta)}{2} \sqrt{\frac{ \hbar e^2 v_s(\theta)}{\overline{\varepsilon}_{\text{HF}}(\theta)\varepsilon_\text{vac}}}=K_R(\theta)\sqrt{\frac{\pi\alpha_\text{fs}\hbar^2 v_F v_s(\theta)}{\overline{\varepsilon}_\text{HF}(\theta)}},
\label{eqn:PAVertex}
\end{align}
where $A$ is the area of the sample, $a_{\mathbf{k},\sigma},a^\dagger_{\mathbf{k},\sigma}$ the electron operators with $\sigma=\pm$ the electron spin, $b_{\mathbf{k}},b^\dagger_{\mathbf{k}}$ are the piezoelectric phonon operators, $\alpha_\text{fs}:=e^2/4\pi\varepsilon_\text{vac}\hbar v_F$ and  $\overline{\varepsilon}_{\text{HF}}(\theta) := (\tilde{\varepsilon}_{\text{HF}}(\theta)+1)/2$. The validity of the shown interaction hamiltonian Eq. (\ref{eqn:PAVertex}) requires two further assumptions: first, the 2DEG or multilayered graphene sample should be thin enough so that in Eq. (\ref{eqn:PSAWWaveVacuum}), $kd\ll 1$, where $z\equiv d$ is the width of the sample and $k$ is the maximum allowed phonon momentum. And second, this maximum allowed momentum should be sufficiently small for the classical piezo-elasticity theory, as shown in Eqs. (\ref{eqn:elasticNewton}-\ref{eqn:piezoConstitutive}), to remain valid. We assume that a maximum momentum on the order of $k_F \sim 10^6-10^7 \, \text{cm}^{-1}$ does not violate this last restriction.

The resulting total Hamiltonian for the combined system of 2D electron gas and piezoelectric Rayleigh phonon is
\begin{eqnarray}
	H= \sum_{\mathbf{k},\sigma}E_{\mathbf{k}} a_{\mathbf{k},\sigma}^{\dagger} a_{\mathbf{k},\sigma} + \hbar \sum_{\mathbf{q}}\omega_{\mathbf{q}}\, b_{\mathbf{q}}^{\dagger} b_{\mathbf{q}} \nonumber \\
	+  H_{\text{e-ph}}^{\text{PA}} + \frac{1}{2A} \sum_{\bf{q}} v_{\mathbf{q}}^{(0)} \rho(\mathbf{q}) \rho(\mathbf{-q} )~,
\label{eqn:PAHamiltonian}
\end{eqnarray}
where $E_\mathbf{k}$ is the electron energy for a 2D wave-vector $\mathbf{k}$,
$\omega_\mathbf{q}=v_s(\theta) q$ is the dispersion relation for the acoustic piezoelectric SAW phonon of 2D wave-vector $\mathbf{q}$ and $v_s(\theta)$ is the SAW propagation velocity, and
$\rho(\mathbf{q}) = \sum_{\mathbf{k},\sigma} a_{\mathbf{k+q},\sigma}^{\dagger} a_{\mathbf{k},\sigma}$ is the Fourier transform of the electron density.

We use the bare Coulomb electron-electron interaction as
\begin{align}
	v_{\mathbf{q}}^{(0)} &= \frac{ e^{2}}{2\overline{\varepsilon}_{\text{HF}}(\theta) \varepsilon_\text{vac} q},
\end{align}
which contains all high-frequency screening processes except for piezoelectric ones.

The bare electron-electron interaction mediated by phonons is \cite{Hwang2010,Mahan2013}
\begin{equation}
	V_{\text{ph}}^{\text{PA}}(\mathbf{q},\omega) = |\gamma_{\mathbf{q}}^{\text{PA}}|^{2} G_{0}^{\text{PA}}(\mathbf{q},\omega) \, ,
	\label{eqn:bareElectronPhononElectronInteraction}
\end{equation}
where
\begin{equation}
G_{0}^{\text{PA}}(\mathbf{q},\omega) = \frac{2\omega_\mathbf{q}}{\hbar\left(\omega^2-\omega_\mathbf{q}^2+i\eta\right)}.
\end{equation}
is the bare piezoelectric acoustic phonon propagator ($\eta \rightarrow 0^+$).  The resulting RPA-type approximation to the dielectric function and effective interaction are:
\begin{align}
	V_{\text{eff}}(\mathbf{q},\omega) &=\frac{ e^2}{2\varepsilon(\mathbf{q},\omega)\varepsilon_\text{vac}q} \nonumber \\
	 &=\frac{v_{\mathbf{q}}^{(0)}+V_{\text{ph}}^{\text{PA}}(\mathbf{q},\omega)}{1-\left[v_{\mathbf{q}}^{(0)}+V_{\text{ph}}^{\text{PA}}(\mathbf{q},\omega)\right]\,\Pi_{0}(\mathbf{q},\omega)}
	\label{eqn:V_eff}
\end{align}
which can also be written as
\begin{equation} \label{eqn:V-eff-RPA}
V_{\text{eff}}(\mathbf{q},\omega)
= \frac{v_{\mathbf{q}}^{(0)}}{\varepsilon_{\text{RPA}}(\mathbf{q},\omega)} + \left|\frac{\gamma_{\mathbf{q}}^{\text{PA}}}{\varepsilon_{\text{RPA}}(\mathbf{q},\omega)} \right|^{2} \tilde{G}^{\text{PA}}(\mathbf{q},\omega) ~.
\end{equation}
where $\varepsilon_{\text{RPA}}(\mathbf{q},\omega) \equiv 1-v_{\mathbf{q}}^{(0)}\Pi_{0}$, with  $\Pi_{0}$ the irreducible polarization function,
\begin{equation}
\tilde{G}^{\text{PA}}(\mathbf{q},\omega) =
\frac{G_{0}^{\text{PA}}(\mathbf{q},\omega)}
{1-\frac{|\gamma_{\mathbf{q}}^{\text{PA}}|^{2}G_{0}^{\text{PA}}(\mathbf{q},\omega)\Pi_{0}(\mathbf{q},\omega)}{\varepsilon_{\text{RPA}}(\mathbf{q},\omega)}}\, .
\end{equation}
In the low frequency limit, $\Pi_{0}(\mathbf{q},\omega\simeq 0)\simeq -D(E_F)=-2k_F/\pi \hbar v_F$, for monolayer graphene \cite{Wunsch2006}, or $\Pi_{0}(\mathbf{q},\omega\simeq 0)\simeq -D(E_F)=-m/2\pi \hbar^2$ for a 2DEG with effective mass $m$ \cite{Giuliani2005}. These two last static limits are exact for $q<2k_F$. In Eq. (\ref{eqn:V-eff-RPA}), the total interaction has been rewritten as the sum of a purely electronically screened Coulomb repulsion and a phonon-induced effective part in which the vertex and phonon-propagator are also screened by just the conducting electrons of the 2DEG or graphene \cite{Mahan2013,Mattuck2012}.
For frequencies small in the scale of the acoustic phonons (or the Bloch-Gr\"uneisen temperature $k_B T_{\rm BG}:=2\hbar v_s k_F $), the bare electron-phonon-electron interaction contributes to the long-range part of the total interaction with a ${\bf q}$-dependence similar to that of the Coulomb repulsion:
\begin{equation}
	V_{\text{ph}}^{\text{PA}}(\mathbf{q},\omega \simeq 0)=|\gamma_{\mathbf{q}}^{\text{PA}}|^{2}G_{0}^{\text{PA}}(\mathbf{q},\omega \simeq 0) = - \frac{2|\gamma_{\mathbf{q}}^{\text{PA}}|^{2}}{\hbar v_{s}q}~,
\label{eqn:Vstatic}
\end{equation}
 Note that it is the acoustic phonon propagator $G_{0}^{\text{PA}}(\mathbf{q},\omega \simeq 0)$ that introduces the coulombic long-range dependence in $q$ via the dispersion of the modes. In the next subsection we shall see that a similar final $q$ dependence has a different origin.

In the limit of low frequencies, $\omega \simeq 0$,  there can be no effective attraction for electrons close to the Fermi-surface because, as shown in Eq. (\ref{eqn:kSquareBound}), the following inequality is satisfied
\begin{equation}
	\frac{-V_{\text{ph}}^{\text{PA}}(\mathbf{q},\omega \simeq 0)}{v_{\mathbf{q}}^{(0)}}=K_R^2(\theta)<1 \label{eqn:interactionsRatio} \,
\end{equation}
which is, in conjunction with the interaction vertex given by Eq. (\ref{eqn:PAVertex}), a central result of this paper.

\subsection{Comparison with optical phonons}

For simplicity, we focus a single branch of the longitudinal optical (LO) for which we assume a constant frequency $\omega_0$. The total Hamiltonian reads as in Eq. (\ref{eqn:PAHamiltonian}) except for the replacements:
\begin{align}
\omega_{\mathbf{q}} & \rightarrow \omega_0 \\
v^{(0)}_{\mathbf{q}} & \rightarrow v^{(\infty)}_{\mathbf{q}} := \frac{e^2}{2\overline{\varepsilon}_\infty \varepsilon_\text{vac} q} \\
\overline{\varepsilon}_\infty & := \frac{\varepsilon_\infty+1}{2} \\
\gamma^{\text{PA}}_\mathbf{q} & \rightarrow \gamma^{\text{OP}}_\mathbf{q} := \sqrt{ g \frac{e^ 2 \hbar \omega_0}{2\varepsilon_\text{vac}q}} \\
g & := \left(\frac{1}{\varepsilon_\infty+1} - \frac{1}{\varepsilon_0+1}\right) >0\, ,
\end{align}
where standard notation for dielectrics is used: $\varepsilon_\infty$ the dielectric constant coming from very high frequency interband electronic transitions and $\varepsilon_0$ would be static dielectric constant in the absence of the piezoelectric phonons at frequencies much smaller than $\omega_0$.

Again, as shown in the discussion around Eq. (\ref{eqn:Vstatic}), for small frequencies ($\omega \ll \omega_0$), the bare phonon-mediated electron-electron interaction contributes to the long-range part of the total interaction like the Coulomb repulsion:
\begin{equation}
V_{\text{ph}}^{\text{OP}}(\mathbf{q},\omega \simeq 0)=|\gamma_{\mathbf{q}}^{\text{OP}}|^{2}G_{0}^{\text{OP}}(\mathbf{q},\omega \simeq 0) = - g \frac{e^2}{\varepsilon_\text{vac} q}~. \label{eqn:VOPstatic}
\end{equation}
However, in contrast to the piezoelectric case, here it is the vertex that introduces the coulombic dependence in $q$.

At small frequencies, $\omega \ll \omega_0$, a single optical phonon is not enough to provide over-screening, because
\begin{equation}
	\frac{-V_{\text{ph}}^{\text{OP}}(\mathbf{q},\omega \simeq 0)}{v_{\mathbf{q}}^{(\infty)}}= \frac{\varepsilon_0-\varepsilon_\infty}{\varepsilon_0+1}<1~.
	\label{eqn:optInteractionsRatio}
\end{equation}

\section{Effect of piezoelectric phonons on superconducting instabilities}

From (\ref{eqn:interactionsRatio}) and (\ref{eqn:V_eff}), we see that, in the static limit ($\omega \simeq 0$), and for $q<2k_F$,
$V_{\text{eff}}$ can be written in the form
\begin{equation}
V_{\text{eff}} (\mathbf{q},0) = \frac{[1-K_R^2(\theta)]v_{\mathbf{q}}^{(0)}}{1+\left[1-K_R^2(\theta)\right]v_{\mathbf{q}}^{(0)}D(E_F)}~,
\label{eqn:Pairing2dpiezo}
\end{equation}
where we note that we have not assumed $q\ll k_F$, as discussed in the paragraph following (\ref{eqn:PAVertex}). From the inequality in (\ref{eqn:interactionsRatio}), we are led to conclude
that over-screening of the Coulomb repulsion by the phonon-mediated attraction is not possible.
Moreover, and following standard textbook reasoning (see for example \cite{Ashcroft2011}), we conclude that BCS-type instabilities must also be ruled out.
A similar result holds for a single branch of optical phonons, as can be seen from Eq. (\ref{eqn:optInteractionsRatio}) (see however \cite{Gorkov2015} for the effect of multiple optical phonon branches from the substrate on superconducting instabilities).

Moreover, in case such over-screening occurred, the static dielectric constant from Eq. (\ref{eqn:V_eff}) would predict unphysical features such as unstable phononic modes with $\tilde{\omega}(q_{c}) = 0$ for some $q_{c} \neq 0$ and even imaginary frequencies for $q<q_{c}$. No matter how small the absolute difference $|1-K_R^2(\theta)|$ happened to be, there would always exist a pole for the static (\ref{eqn:V_eff}) at small enough $q$ (what cannot occur in standard BCS metals), signaling a different type of instability, possibly a charge density wave.

On the other hand, the result (\ref{eqn:interactionsRatio}) for the vertex could still lead to higher angular momentum pairing instabilities (as in the Kohn-Luttinger mechanism \cite{kohn1965new}) provided that $K_R^2(\theta)$ is sufficiently large and anisotropic, a case not considered by us.

\subsection{Eliashberg formalism \cite{Einenkel2011}}

The previous reasoning about the absence of superconducting instabilities, is incomplete and somewhat oversimplified. Three reasons support this claim: (i) Long-wave piezoelectric phonon excitations (as considered in the present work) can never be the only source of effective electron-electron interactions; in particular, we have not taken into account the short range electric fluctuations of the substrate. (ii) There is definitely some dynamic over-screening at high frequencies [see Eq. (\ref{eqn:bareElectronPhononElectronInteraction})]. And (iii) Coulomb interaction has to be properly renormalized by taking into account collisions with high momentum transfer, which diminishes the Coulomb repulsion and thus comparatively strengthens the other attraction mechanisms.

Leaving aside the first objection momentarily, we can use the Eliashberg formalism, as applied to graphene in Ref. \cite{Einenkel2011}, to deal with the other two objections. The effective interaction could cause superconducting instabilities if a dimensionless electron-phonon coupling $\lambda^\text{PA}$ happened to be greater than an also dimensionless Coulomb pseudopotential $\mu^*$ coming from high-energy renormalizations \cite{Morel1962,Einenkel2011}. The coupling constant $\lambda^\text{PA}$ in the Eliashberg formalism is the same appearing in (the real part of the) self-energy calculations to renormalize the Fermi velocity \cite{Gonzalez2015} and is given by:
\begin{align}
&\alpha_\text{PA}^2F(\omega)=\frac{|\gamma^\text{PA}|^{2}}{2 \pi^2 \hbar^2 v_s v_F}\frac{\sqrt{1-(\frac{\omega/v_s}{2k_F})^2}}{(1+\frac{k_\text{TF}}{\omega/v_s})^2}~, \nonumber \\
&\lambda^\text{PA}=2\int\limits_{0}^{\infty}\frac{\alpha_\text{PA}^2F(\omega)}{\omega}d\omega = \frac{r_s}{\pi} K_R^2\, F(2r_s)~, \label{eqn:lambda}\\
& F(x) = \int\limits_{0}^{1}\frac{t\sqrt{1-t^2}\,dt}{(t+x)^2}=-2+x\pi + \frac{(1-2x^2)\, \text{acosh}(x^{-1})}{\sqrt{1-x^2}}~, \nonumber
\end{align}
where $r_s(\theta):=\alpha_{\text{fs}}/\varepsilon(\theta)$, and the symbols $r_s, v_s, \gamma^{\text{PA}}, K_R^2$ stand for the Fermi surface angle-averaged quantities of the same name. The constant $\mu^*$ equals
\begin{equation}
\mu^* = \frac{\frac{1}{4}D(E_F)V}{1+\frac{1}{4}D(E_F)V \log(\frac{E_F}{\hbar \omega_c})}~,
\label{eqn:PseudoPotential}
\end{equation}
where $\omega_c$ is some energy cutoff which should satisfy $\omega_\text{Debye} \ll \omega_c \ll E_F/\hbar $ \cite{Einenkel2011} and $V$ comes from the Fermi surface average of the Thomas-Fermi renormalized Coulomb repulsion $v_\mathbf{q}^{(0)}/(1+\frac{k_\text{TF}}{q})$. We have
\begin{align}
&\frac{1}{4}D(E_F)V=\frac{r_s}{\pi}G(2r_s), \\
& G(x)=\int\limits_{0}^{1}\frac{\sqrt{1-t^2}\,dt}{t+x}=-1+\frac{\pi x}{2}+\sqrt{1-x^2}\,\text{acosh}(x^{-1})~, \nonumber
\end{align}
and therefore, provided that one takes $\hbar\omega_c \simeq k_B T_\text{BG}$, so that $\log\left(\frac{E_F}{\hbar \omega_c}\right) \simeq \log(\frac{v_F}{2v_s}) \simeq 5$. Thus, an estimate of the effective pseudo-potential is
\begin{equation}
\mu^*\simeq\frac{\frac{r_s}{\pi}G(2r_s)}{1+\frac{5r_s}{\pi}G(2r_s)} ~.
\end{equation}
There could be intravalley \footnote{In this analysis, we are taking into account only long-wavelength piezoelectric phonons, hence no intervalley pairing instability could occur. We address this question in the next paragraph.} superconducting instabilities provided that
\begin{equation}
1<\frac{\lambda^\text{PA}}{\mu^*}=K_R^2 \frac{F(2r_s)}{G(2r_s)}\left[1+\frac{5r_s}{\pi}G(2r_s)\right]~,
\end{equation}
which imposes a constraint on the value of $K_R^2$ from the piezoelectric substrate with respect to quantities depending on $r_s$. The coupling $K_R$ should be very large and actually greater than 1 for this choice of $\omega_c$, although there could exist superconductivity in this idealized case of a system consisting just of the graphene electrons and long wavelength piezoelectric phonons, provided that $E_F/\omega_c$ is larger and $K_R^2$ close to 1.

In order to amend the first objection, we have to consider proper phonons of the electronic system (here we go on considering graphene), in conjunction with the short range of the piezoelectric ones. Then, pairing instabilities due to intervalley scattering have to be considered as well, because intravalley scattering terms contribute also to the intervalley pairing gap.
With the notation in Ref. \cite{Einenkel2011}, an estimate on the critical temperature for the intravalley pairing is \cite{Einenkel2011} $T_c^\text{intra} = 1.13\omega_\text{Debye}\,\exp\left(-\frac{1+\lambda}{\lambda_{11}-\mu_{11}^*}\right)$, and a very similar is obtained for the intervalley transition
$T_c^{\text{inter}} = 1.13\omega_\text{Debye}\,\exp\left(-\frac{1+\lambda}{\lambda-\mu_{12}^*}\right)$, with $\lambda = \lambda_{11}+\lambda_{12}$ and the previously computed $\lambda^{\text{PA}}$ included into the intravalley term $\lambda_{11}$ ($\lambda_{12}$ denotes the contribution from all intervalley terms). Here the pseudo-potential $\mu_{12}^*$ is only slightly larger than $\mu_{11}$, and both are given by similar formulas as in Eq. (\ref{eqn:PseudoPotential}), but with $\omega_c\rightarrow \omega_\text{Debye}$.

The upshot of this discussion is that the long-wavelength piezoelectric phonons work in favor of pairing instabilities, as shown in Fig. \ref{fig:tcPlot}. We emphasize, however, that we are not claiming that a piezoelectric substrate {\it per se} necessarily increases the critical temperature, since it could be the case that other piezoelectric fluctuations not considered in the present study (e.g. shorter wavelength modes) could work against pairing instabilities.

\begin{figure}
\raggedleft
\includegraphics[width=0.92\linewidth]{./tcPlot}
\includegraphics[width=0.9\linewidth]{./deltatcPlot}
\caption{(Color online) Critical temperature and variation from the \textquotedblleft bare\textquotedblright~one in \cite{Einenkel2011} for graphene on a piezoelectric substrate, as a function of the conduction band density. Three pairs of plots are given for three different values of $K_R$ and the two values of the constant $C = 3.5$ and 5 in Eq. (5.1) for $\lambda$ in the previous reference \cite{Einenkel2011}. }
\label{fig:tcPlot}
\end{figure}

\section{Conclusions}
In conclusion, we have derived a general expression for the two-dimensional electron-phonon piezoelectric interaction valid for any piezoelectric substrate covered by a two-dimensional electron system, as in the classical 2D Fr\"ohlich hamiltonian for the optical phonons, and characterized the magnitude of the interaction. Our results show that electron overscreening cannot be achieved just with the strongest piezoelectric phonons within our assumptions because $K_R^2<1$ is always satisfied. Nevertheless, these phonons could help further in other contexts where the 2D superconductivity is known to exist, for example in bulk few-layer $\text{MoS}_{2}$ with most of the carriers confined to the first layer \cite{Ye2012,Roldan2013}; or postulated to exist but not yet observed due to experimental difficulties (e.g. very heavily doped graphene \cite{Einenkel2011}). Other example is the recent high-temperature superconductor system of 2D FeSe on top of the ferroelectric $\text{SrTiO}_{3}$, whose optical phonons have been analyzed arriving at conclusions similar to ours \cite{Gorkov2015}, and where the strong piezoelectric phonons could play a role as well.
\begin{acknowledgments}
We want to thank Fernando Calle and J\"urgen Schiefele for valuable discussions. DGG acknowledges financial support from Campus de Excelencia Internacional (Campus Moncloa UCM-UPM). This work has been supported by the Spanish Ministry of Economy (MINECO) through Grants No. FIS2011-23713, FIS2013-41716-P; the European Research Council Advanced Grant (contract 290846), and the European Commission under the Graphene Flagship, contract CNECTICT- 604391.

\end{acknowledgments}


\section{Introduction}

The notion of \emph{cellular automaton} was introduced by J. von Neumann in his seminal paper \cite{neumann1966theory} where he aimed at modeling a self-reproducing entity. The idea behind the concept of a cellular automaton is that the richness of states exhibited by the evolution of a macroscopic system could emerge from a simple local interaction rule among its elementary constituents. More precisely, a cellular automaton is a lattice of cells that can be in a finite number of states, together with a rule for the update of cell states from time $t$ to time $t+1$. The principal requirement for such rule is \emph{locality}: The state of the cell $\bvec{x}$ at step $t + 1$ depends on the states of a finite number of neighboring cells at step $t$. The use of classical cellular automata for simulation of quantum mechanics was proposed by 'tHooft \cite{hooft2014cellular}, followed by other authors \cite{elze2014action}.

The first author to suggest the introduction of the quantum version of cellular automata was R. Feynman in the celebrated paper of Ref.~\cite{feynman1982simulating}. Since then, the interest in \emph{quantum cellular automata}
(QCAs), has been rapidly growing, especially in the Quantum Information community, leading to many results about their general structure (see e.g Refs.\cite{schumacher2004reversible,arrighi2011unitarity,gross2012index} and references therein).
Special attention is devoted in the literature to QCAs with linear evolution, known as \emph{Quantum Walks} (QWs)
\cite{grossing1988quantum,aharonov1993quantum,ambainis2001one,reitzner2011quantum}, which were especially applied in the design of quantum algorithms
\cite{childs2003exponential,ambainis2007quantum,magniez2007quantum,farhi2007quantum}, providing a speedup for relevant computational problems.

More recently QCAs have been considered as a new mathematical framework for Quantum Field Theory
\cite{darianovaxjo2010,darianopla,Bisio2015244,PhysRevA.90.062106,bisio2014quantum,arrighi2014dirac,arrighi2013decoupled,farrelly2014causal,farrelly2013discrete}.
Within this approach, each cell of the lattice corresponds to the
evaluation $\psi(\bvec{x})$ of a quantum field at the site $\bvec{x}$ of a lattice, with the dynamics updated in discrete time steps by a local unitary evolution.  Assuming that the lattice spacing corresponds to an hypothetical discrete Planck scale\footnote{Other approaches to discrete space-time based on p-adic numbers were studied in Refs.~\cite{albeverio2009p}.}, the usual quantum field evolution should emerge as a large scale approximation of the automaton dynamics. On the other hand, the QCA dynamics will exhibit a different behaviour at a very small scale, corresponding to ultra-relativistic wave-vectors.

The analysis of this new phenomenology is of crucial importance in providing the first step towards an experimental test of the theory as well as a valuable insight on the distinctive features of the QCA theory. Until now the research was mainly focused on linear QCAs which describe the dynamics of free field. By means of a Fourier transform the linear dynamics can be easily integrated and then, as we will show in section \ref{sec:interp-hamilt-diff},
an approximated model for the evolution of particle states (i.e. state of the dynamics narrow-band in wave-vector) can be obtained.
Moreover, it is also possible to derive an analytical solution of the evolution in terms  of a path sum in the position space, thus giving the QCA analog of the Feynman propagator \cite{d2014path,d2014discrete}.
In section \ref{sec:phenomenology} we will exploit these tools
to explore many dynamical features of the QCA models for the Weyl and Dirac fields and to compare them with the corresponding counterpart emerging from the Weyl and Dirac Equation. We will see that, when considering massive Fermionic fields (e.g electrons) the deviations from the usual field dynamics cannot be reached by present day experiments, contrarily to the case of the QCA theory of the free electromagnetic field.
In section \ref{sec:phen-qca-theory} we will review the main
phenomenological aspects of the QCA model for free photons
(that in this framework become composite particles) with special emphasis of the emergence of a frequency-dependent speed of light, a Planck-scale effect already considered by other authors in the Quantum Gravity community
\cite{ellis1992string,lukierski1995classical,Quantidischooft1996,amelino2001testable,PhysRevLett.88.190403}.
In the final Section of this paper we address two issues of the
QCA theory that are still under investigation. The first one concerns
the notion of Lorentz covariance: Because of its intrinsic discreteness, a QCA model cannot enjoy a notion of
Lorentzian space-time and the usual Lorentz covariance must break down
at very small distances. One way of addressing the problem of changing
the reference frame
is to assume that every inertial observer must observe the same dynamics.
Then one can look for a set of modified Lorentz transformtion which
keep the QCA dispersion relation invariant. The first step of this
analysis are reported in Section  \ref{sec:lorentz-covar-deform}.
The second issue we will briefly address in Section
\ref{sec:therm-free-ultr}
 are thermodynamical effects that could emerge from modified QCA dynamics.

\section{Weyl, Dirac and Maxwell automata}
\label{sec:weyl-dirac-maxwell}

A Quantum Cellular Automaton (QCA) describes the discrete
time evolution of a set of cells, each one containing an array of
quantum modes. In this section we review the QCA models for the free fermions and for
the free electromagnetic field. For a complete presentation of these
results we refer to Refs. \cite{Bisio2015244,PhysRevA.90.062106,bisio2014quantum}.
  Within our
framework we will consider \emph{Fermionic} fields, our choice being
motivated by the requirement that amount of information in finite
number of cells must be finite. Then, each cell $\bvec{x}$ of the
lattice is associated with the Fermionic algebra generated by the
field operators $\{ \psi(\bvec{x}), \psi^{\dag}(\bvec{x})\}$ which
obey the canonical anticommutation relation $[\psi(\bvec{x}),
\psi^{\dag}(\bvec{x}')]_+ = \delta_{\bvec{x},\bvec{x}'}$ and
$[\psi(\bvec{x}), \psi(\bvec{x}')]_+ =0$ \footnote{We denote as
  $[A,B]_+$ the anticommutator $AB+BA$.  The commutator $AB-BA$ will
  be denoted as $[A,B]_-$.}.  With a slight genealization, we consider
the case in which each cell correspond to more than one Fermionic
mode. Different Fermionic modes will be denoted by an additional
label, e.g. $\psi_i(\bvec{x})$.  The automaton evolution will be
specified by providing the unit-step update of the Fermionic field
operators.  This rule defines the primitive physical law, and must
then be as simple and universal as possible. This principle translates
into a minimization of the amount of mathematical parameters
specifying the evolution. In particular we constrain the automaton to
describe a \emph{unitary} evolution which is \emph{linear} in the
field. We
notice that the linearity of the QCA restrict the scenario to
non-interacting field dynamics. Then we require the evolution to be
\emph{local}, which means that at each step every cell interacts with
a finite number of neighboring cells, and \emph{homogeneous}, meaning
that all the steps are the same, all the cells are identical systems
and the interactions with neigbours is the same for each cell (hence
also the number of neigbours, and the number of Fermionic modes in
each cell).  The neighboring notion also naturally defines a graph $\Gamma$
with $\bvec{x}$ as vertices and the neighboring
couples as edges.
We also assume \emph{transitivity}, i.e. that every two cells are
connected by a path of neighbors and \emph{isotropy} which means that
the neighboring relation os symmetric and there exist a group of automorphisms for the graph under which the automaton is covariant.  From these assumptions one can show\footnote{This step
 would requires a more precise mathematical characterization (which we omit) of the
 presented assumptions. See Ref. \cite{PhysRevA.90.062106} for the details.} that graph $\Gamma$ is a \emph{Cayley graph} of a group $G$. In the following, we consider the Abelian case $G = \mathbb{Z}^3$.

Let $\mathsf{S}_+$ denote the set of generators of
$\mathbb{Z}^3$ corresponding to the Cayley graph $\Gamma$
and let $\mathsf{S}_-$ be the set of inverse generators.
 For a given cell $\bvec{x}$ the set of neighboring cells
is given by the set $\mathcal{N}_\bvec{x} := \{ \bvec{x} + \bvec{z} \,|\,
\bvec{z} \in  \mathsf{S} := \mathsf{S}_+ \cup \mathsf{S}_-\}$,
where we used the additive notation for the group composition.
If $s$ is the number of Fermionic modes in each cell,
the single step evolution can then be represented in terms of
$s \times s $ transition matrices $A_{\bvec{z}}$ as follows
\begin{align}
  \boldsymbol{\psi}(\bvec x, t+1 ) = \sum_{\bvec{z} \in \mathsf{S}} A_{\bvec{z}}
  \boldsymbol{\psi}(\bvec x + \bvec{z}, t).
\end{align}
where $\boldsymbol{\psi}(\bvec x , t)$ is the array of field operators at $\bvec x$
at step $t$.
Upon introducing the Hilbert space $\ell^2( \mathbb{Z}^3 )$,
the automaton evolution can be described by the unitary matrix
$A$ on   $\ell^2( \mathbb{Z}^3 ) \otimes \mathbb{C}^s$ given by
 \begin{align}
A :=  \sum_{\bvec{z} \in \mathsf{S}} T_{\bvec{z}} \otimes A_{\bvec{z}},
\end{align}
where $T_\bvec{x}$ denotes the unitary representation of
$ \mathbb{Z}^3$ on $\ell^2( \mathbb{Z}^3 )$,
$T_\bvec{y} \ket{\bvec x} = \ket{\bvec x + \bvec{y}} $.
If $s=1$ , i.e. there only one Fermionic mode in each cell,
one can prove that the only evolution which obeys our set of
assumptions is the trivial one ($A$ is the identity matrix).
Then we are led to consider the $s=2$ case and we denote the two
Fermionic modes as $\psi_L(\bvec x,t)$ and $\psi_R(\bvec x,t)$.
Moreover in the $s=2$ case one can show that our assumptions\footnote{In order to
  prove this step one need a stronger isotropy condition than the one
  presented in the text. See Ref. \cite{PhysRevA.90.062106} for the
 details.} imply that the only lattice which admits a nontrivial
evolution is the body centered cubic (BCC) one.
Being $\mathbb{Z}^3 $ an abelian group, the Fourier transform is well
defined and the operator $A$ can be block-diagonalized as follows
\begin{align}
\label{eq:autofourier}
  A = \int_B \d^3 \!\!\bvec k \,\, \ketbra{\bvec k}{\bvec k} \otimes A_\bvec k,
\end{align}
where
$|\bvec k\>:=(2\pi)^{-\tfrac{3}{2}}\sum_{\bvec x\in
  \mathbb{Z}^3}e^{i\bvec k\cdot\bvec x}|\bvec x\>$,
$B$ is the first Brillouin zone of the BCC lattice and
$A_\bvec k:=\sum_{\bvec{z}\in S}\bvec k\,e^{-i\bvec
  k\cdot\bvec{z}}A_\bvec{z}$ is a $2 \times 2$ unitary for every
$\bvec k$.
We have only two (up to a local change of basis) non trivial QCAs corresponding to the unitary
matrices
\begin{equation}\label{eq:weyl3D}
A^{\pm}_{\bvec k} := d^{\pm}_{\bvec k} I+\tilde{\bn}^{\pm}_{\bvec k}\cdot\boldsymbol{\sigma}
=\exp[-i\bvec{n}^{\pm}_{\bvec k} \cdot \boldsymbol{\sigma}],
\end{equation}
where $ \boldsymbol{\sigma}$ is the array
$(\sigma_x,\sigma_y,\sigma_z)$
of Pauli matrices and we defined
\begin{align}
&\tilde{\bn}^{\pm}_{\bvec k} :=
\begin{pmatrix}
s_x c_y c_z \mp c_x s_y s_z\\
\mp c_x s_y c_z - s_x c_y s_z\\
c_x c_y s_z \mp s_x s_y c_z
\end{pmatrix}\!\!,\,
{\bn}^{\pm}_{\bvec k}:=\frac{\lambda^{\pm}_{\bvec k}\tilde{\bn}^{\pm}_{\bvec k}}{\sin\lambda^{\pm}_{\bvec k}},\nonumber\\
&\d^{\pm}_{\bvec k} :=  (c_x c_y c_z \pm s_x s_y s_z ),\;
\lambda^{\pm}_{\bvec k}:=\arccos(d^{\pm}_{\bvec k}),\nonumber\\
&c_\alpha := \cos({k}_\alpha/\sqrt{3}),\;s_\alpha:= \sin({k}_\alpha/\sqrt{3}),\;\alpha = x,y,z.\nonumber
\end{align}
The matrices $A_{\bvec k}^\pm$ in Eq. \eqref{eq:weyl3D} describe the evolution of a two-component
Fermionic field,
\begin{align}
  {\psi} ({\bvec k},t+1) =
A_{\bvec k}^\pm {\psi} ({\bvec k},t),
\quad
  {\psi} ({\bvec k},t) : =
\begin{pmatrix}
 {\psi}_R ({\bvec k},t)\\
 {\psi}_L ({\bvec k},t)
  \end{pmatrix}.
  \end{align}
  The adimensional framework of the automaton corresponds to measure
  everything in Planck units. In such a case the limit $|{\bvec k}|\ll 1$
  corresponds to the relativistic limit, where on has
\begin{equation}
\bn^{\pm}({\bvec k})\sim\tfrac{{\bvec k}}{\sqrt{3}},\quad A^{\pm}_{\bvec k}\sim\exp[-i\tfrac{{\bvec k}}{\sqrt{3}} \cdot\boldsymbol{\sigma}],
\end{equation}
corresponding to the Weyl's evolution, with the
rescaling  $\tfrac{{\bvec k}}{\sqrt{3}} \to \bvec k$.
Since the QCAs $A^{+}$ and $A^{-}$ reproduce the dynamics of the
Weyl equation in the  limit $|{\bvec k}|\ll 1$, we refer to them as
\emph{Weyl automata}. For sake of simplicity,
in the following we will consider only one Weyl automaton, i.e. we
define $A_{\bvec k} := A_{\bvec k}^-$ and we similarly drop all the others
$\pm$ superscripts.
 This choice is completely painless since
all the methods that we will use can be easily adapted to the choice
$A_{\bvec k} = A_{\bvec k}^-$. However the two automata, beside giving the
Weyl equation for small $\bvec k$, exhibit a different behaviour at high
$\bvec k$ and we will point out those differences whenever it will be relevant.

The derivation that we sketch previously can be carried on also in the
two dimensional case (considering QCA on Cayley graphs  of
$\mathbb{Z}^2$)
 and in the one dimensional case (considering QCA on Cayley graphs  of
$\mathbb{Z}$).
In the 2-dimensional case we obtain a unique (up to a local change of basis) the QCA on the square latticeand it leads to
\begin{align}
\label{eq:weyl2D}
   A^{(2D)}_{\bvec k}=I d^A_\bvec k-i\boldsymbol\sigma\cdot\bvec a^A_\bvec k,
\end{align}
where the functions $\bvec a_\bvec k$ and $d_\bvec k$ are expressed in terms of $k_x:=\frac{k_1+k_2}{\sqrt2}$ and $k_y:=\frac{k_1-k_2}{\sqrt2}$ as
$
 (a^A_\bvec k)_x:= s_xc_y
$,
$
  (a^A_\bvec k)_y:=c_xs_y
$,
$
  (a^A_\bvec k)_z:=s_xs_y
$,
$
  d^A_\bvec k:=c_xc_y.
$
where $c_i =\cos\tfrac{k_i}{\sqrt2}$ and $s_i=\sin\tfrac{k_i}{\sqrt2}$.
In the one dimensional case we find
\begin{align}
\label{eq:weyl1D}
   A^{(1D)}_{k}=
   \begin{pmatrix}
     e^{-i k} & 0 \\
      0&  e^{i k}  \\
   \end{pmatrix}
\end{align}
Both in the 2-dimensional and 1-dimensional cases the limit
$|{\bvec k}|\ll 1$ gives the 2-dimensional and 1-dimensional Weyl equation
respectively (in the 2-dimensional we need the rescaling  with the
rescaling  $\tfrac{{\bvec k}}{\sqrt{2}} \to \bvec k$.).
The QCA in Eqs. (\ref{eq:weyl3D},~(\ref{eq:weyl2D}) and~(\ref{eq:weyl1D}))
describe the dynamics of free massless Fermionic fields. If we couple
two Weyl automata $A_\bvec k$ with a mass term we obtain a new QCA $U_\bvec k$ given by
\begin{align} \label{eq:diracauto}
  U_\bvec k = \begin{pmatrix}
     n A_\bvec k & im I \\
      im I & n A_\bvec k^\dagger  \\
   \end{pmatrix}
\quad n^2 + m^2 =1.
\end{align}
Clearly this construction can be done in the $1$,$2$ and
$3$-dimensional cases and the resulting QCA is alway unitary and local.
One can easily see that in the limit $|{\bvec k}|\ll 1$  and $m \ll 1$
Eq. \ref{eq:diracauto}  (with the appropriate rescaling of $\bvec k$ in
$2$ and $3$ dimensions) gives the same evolution as Dirac equation
and then we denote the automata of Eq. (\ref{eq:diracauto})
\emph{Dirac automata}.

The $3$-dimensional Weyl QCA can also be use as a building block for
QCA model of free electrodynamics. The basic idea is to interpret the
photon as a pair of Weyl fermions that are suitably correlated in
wave-vector. Then one can show that, in an appropriate regime, this field
obeys the dynamics dictated by the Maxwell equations and the bosonic
commutation relation are recovered. This approach recall the so-called
beutrino theory of light of De Broglie
\cite{de1934nouvelle,jordan1935neutrinotheorie,kronig1936relativistically,perkins1972statistics,perkins2002quasibosons} which suggested that
the photon could be a composite particle made of of a
neutrino-antineutrino pair.  Within our framework (we omit the details
of this construction that can be found in
Ref. \cite{bisio2014quantum}) the electric and magnetic field are
given by
\begin{align}
  \label{eq:electric and magnetic field}
  &\bvec{E}:=|{\bn}_{\tfrac\bk2}|(\bvec{F}_T+\bvec{F}_T^\dag),\quad\bvec{B}:=i|{\bn}_{\tfrac\bk2}|(\bvec{F}_T^\dag-\bvec{F}_T),\\
  &2|{\bn}_{\tfrac\bk2}|\bvec{F}_T=\bvec{E} + i \bvec{B}
 \nonumber\\
&\bvec{F}_T(\bvec k) := \bvec{F}(\bvec k)  -
\left(\frac{\bvec{n}_{\frac{\bvec k}{2}}}{|\bvec{n}_{\frac{\bvec k}{2}}|} \cdot
{\bvec{F}}(\bvec k) \right)
\frac{\bvec{n}_{\frac{\bvec k}{2}}}{|\bvec{n}_{\frac{\bvec k}{2}}|}
\nonumber\\
&\bvec{F}(\bvec k) :=
(
F^{1}(\bvec k), F^{2}(\bvec k), F^{3}(\bvec k)
) ^T, \nonumber\\
&F^{j}(\bvec k) :=
 \int \frac{ d \bvec{q}}{(2 \pi)^3}
f_{\bvec k}(\bvec{q})
\phi
\left(\tfrac{\bvec k}{2}-\bvec{q}\right)
\sigma^{j}
\psi
 \left(\tfrac{\bvec k}{2}+\bvec{q}\right) \nonumber
\end{align}
where
$\int\frac{d\bvec{q}}{(2\pi)^3} |f_{\bvec k}(\bvec{q})|^2 =1, \forall \bvec k$
and $\phi
\left(\bvec k  \right)$, $\psi
\left(\bvec k  \right)$
are two massless Fermionic fields whose evolution is dictated by the
automaton
$A_\bvec k^*$ and $A_\bvec k$ respectively\footnote{We denote as $A^*$ the
  complex conjugate of $A$}, i.e.
\begin{align}
  \phi
\left(\bvec k , t \right) &=
A_\bvec k^{* t}  \phi
\left(\bvec k   \right)
\qquad
  \psi
\left(\bvec k , t \right) =
A_\bvec k^{t}  \phi
\left(\bvec k   \right)  \nonumber \\
  \phi
\left(\bvec k \right) &=
\begin{pmatrix}
   \phi_R
\left(\bvec k \right)  \\
\phi_L
\left(\bvec k \right)
\end{pmatrix}
\qquad
  \psi
\left(\bvec k \right) =
\begin{pmatrix}
   \psi_R
\left(\bvec k \right)  \\
\psi_L
\left(\bvec k \right)
\end{pmatrix}
\nonumber .
\end{align}
For an appropriate choice of the functions $f_{\bvec k}(\bvec{q})$ (see
Ref. \cite{bisio2014quantum}), one can prove that the evolution of the
Electric and Magnetic fields which are defined in
Eq.\eqref{eq:electric and magnetic field} is given by the following
equations\footnote{Since a QCA describe an evolution which discrete in
  time, the derivative with respect time is not defined in this
  context. However we can imagine
$A_\bvec k^{t}$ to be defined for any
  real value of $t$ and then derive with respect the continuous
  variable $t$. This is the construction which underlies Eq.~(\ref{eq:maxwell2})}
\begin{align}
  \begin{split} \label{eq:maxwell2}
    &\partial_t\bvec{F}_T(\bvec k,t) =
2\bvec{n}_{\frac{\bvec k}{2}} \times \bvec{F}_T(\bvec k,t)\\
&2\bvec{n}_{\frac{\bvec k}{2}} \cdot \bvec{F}_T(\bvec k,t) = 0,
  \end{split}
\end{align}
which in the limits $|{\bvec k}|\ll 1$ and with the rescaling
$\tfrac{\bvec k}{\sqrt{3}} \to \bvec k$
 become the Fourier trasform of the usual vacuum Maxwell equations
in position space
\begin{align}
      \label{eq:maxwellstandard}
  \begin{array}{lcl}
    \nabla \cdot \bvec{E}=0    & &\nabla \cdot \bvec{B} =0\\
    \partial_t \bvec{E} = c\nabla \times \bvec{B}  && \partial_t \bvec{B} = -c\nabla \times \bvec{E} \;\;.
  \end{array}
\end{align}
Because of this result, we refer to the construcion of Eq. (\ref{eq:electric
  and magnetic field}) as the \emph{Maxwell automaton}. This is a
slight abuse of notation since Eq. (\ref{eq:electric and magnetic
  field}) does not introduce any new QCA model but it defines a field
of bilinear operators, each one of them evolving with the Weyl
automaton, that  can be interpreted as the electromagnetic field in
vacuum. In this sense the expression ``Maxwell automaton'' actually
means QCA model for the Maxwell equations (in vacuum).

\section{Analysis of the dynamics}

The aim of this section is to analyse the dynamics of the QCA models
presented in the previous section.  The material of this section can
be found in Refs.
\cite{Bisio2015244,PhysRevA.90.062106,bisio2013dirac,d2014path,d2014discrete}.

In this section we fosus on the single particle sector of the QCA. Since the
QCAs we are considering are linear in the fields the single particle
sector contains all the information of the dynamics (it is a free
theory). We can then write $\ket{\psi (t)} = A^t \ket{\psi(0)}$ where
$\ket{\psi (t)} := \sum_x \ket{\psi(x,t)} \ket{x}$ is generic  one particle state. This kind of framework is
better known in the literature under the name \emph{Quantum Walk}
\cite{grossing1988quantum,aharonov1993quantum,ambainis2001one,succi1993lattice,bialynicki1994weyl,meyer1996quantum}.

\subsection{Interpolating Hamiltonian and differential equation for single-particle wave-packets} \label{sec:interp-hamilt-diff}

Since a QCA (and a Quantum Walk) describes a discrete evolution on a lattice,
the notion of Hamiltonian (like any other differential operators)
is completely deprived of physical meaning.
However it is useful to introduce a Hamiltonian operator $H_\bvec k$, that we
call \emph{interpolating Hamiltonian} that obeys the following equation
\begin{align}
\label{eq:interpolhamil}
  A^t_\bvec k = e^{it H_\bvec k}.
\end{align}
where $A_\bvec k^{t}$ is defined for any real value of $t$ (the automaton
$A_\bvec k$ can be any of the QCA models of Section
\ref{sec:weyl-dirac-maxwell}).  It is clear that $ H_\bvec k$ is the
genrator of the continuous time evolution which interpolate the QCA
dynamics between the integer steps.  The eigenvalue of $H_\bvec k$ have
the same modulus and its analytical expression, denoted by
$\omega(\bvec k)$, is the \emph{dispersive relation} of the automaton and
it provides a lot of information about the dynamics.  In analogy with
what we do in Quantum Field Theory, states of the dynamics
corresponding to the positive eigenvalue $\omega(\bvec k)$ are called
\emph{particle states} while eigenstate with negative eigenvalue
$-\omega(\bvec k)$ are called \emph{antiparticle states}.
If we denote with $\ket{u}_\bvec k$ a positive frequency eigenstate of
 $H_\bvec k$ (i.e. $H_\bvec k  \ket{u}_\bvec k = \omega(\bvec k) \ket{u}_\bvec k$)
a generic particle state with positive frequency
\begin{align}\label{eq:particlestates}
\ket{\psi}_+ =  \int  \frac{d\bvec{k}}{(2\pi)^n}  g(\bvec k)
\ket{u}_\bvec k \ket{\bvec k}
\end{align}
where $n$ is the dimension of the lattice and $g(\bvec k)$ is a normalized
probability amplitude. We remind that for the $2$ and $3$-dimensional
Dirac QCA $H_\bvec k$ has dimension $4$ and then both the eigenvalues
have degeneracy (corresponding to the spin degree of freedom).
The construction of Eq. (\ref{eq:particlestates}) can be
straightforwardly applied also for the definition general antiparticle
states with negative frequency $\ket{\psi}_-$.

One can use the
interpolating Hamiltonian $H_\bvec k$ in order to rephrase the continuous
evolution in terms of a differential Equation, i.e.
 \begin{align}
\label{eq:schroding1}
  i \partial_t \ket{\psi(\bvec k,t)} =  H_\bvec k \ket{\psi(\bvec k,t)}.
 \end{align}
where $\ket{\psi(\bvec k,t)}$ is the wave-vector representation of a
one-particle state, i.e. $\ket{\psi} = \int  \tfrac{d\bvec{k}}{(2\pi)^3}   \ket{\psi(\bvec k,t)}
\ket{\bvec k}$ .
When  the initial state has positive frequency (see
 Eq. (\ref{eq:particlestates})) and its distribution  $g(\bvec k)$ is
 smoothly peaked around a given $\bvec k_0$, the evolution of Eq.~(\ref{eq:schroding1})
can be approximated by the following dispersive equation
\begin{equation}
\label{eq:schroding2}
i\partial_t \tilde{g}(\bvec x,t)=\pm[\bvec v\cdot\bvec \nabla+\tfrac{1}{2}\bvec D\cdot\bvec
\nabla\bvec \nabla]\tilde{g}(\bvec x,t),
\end{equation}
where $\tilde{g}(\bvec x,t)$ is the Fourier transform of
$\tilde{g}(\bvec k,t):=e^{-i\bvec k_0\cdot\bvec x+i\omega(\bvec k_0) t}\psi(\bvec k,t)$, and $\bvec v$ and $\bvec D$ are the drift
vector $\bvec v=\left(\bvec\nabla_{\bvec k}\omega\right)(\bvec k_0)$ and diffusion tensor
$\bvec
D=\left(\bvec\nabla_{\bvec k}\bvec\nabla_{\bvec k}\omega\right)(\bvec k_0)$,
respectively.
Intuitively the vector $\bvec v$ represent the velocity of the
wavepacket and the tensor $\bvec D$ tells us how the
wavepacket spreads during the evolution.
The accuracy of the approximation can be analytical evaluated (see
Ref. \cite{Bisio2015244}) and compared with comuter simulation as in Fig.~\ref{fig:Hermite}.

\begin{figure}[h!]
  \centering
\includegraphics[width=0.45\textwidth]{Hermite4times1.pdf}
\qquad
\includegraphics[width=0.45\textwidth]{Gaussian4times1.pdf}\\
\includegraphics[width=0.45\textwidth]{Hermite4times3.pdf}
\qquad
\includegraphics[width=0.45\textwidth]{Gaussian4times3.pdf}\\
\includegraphics[width=0.45\textwidth]{Hermite4times4.pdf}
\qquad
\includegraphics[width=0.45\textwidth]{Gaussian4times4.pdf}
\caption{(Colors online) Test of the approximated evolution of Eq.
\eqref{eq:schroding2} of the one dimensional Dirac automaton evolution. {\bf Left figure:} here the
  state is a superposition of Hermite functions
  (the polynomials $H_j(x)$ multiplied by the Gaussian) peaked around $k_0=3\pi/10$. {\bf Right
    figure:} here the initial state is Gaussian profile peaked around
  $k_0=0.1$. This figure is published in Ref. \cite{Bisio2015244}.}\label{fig:Hermite}
\end{figure}

\section{Phenomenology}
\label{sec:phenomenology}
This section is devoted to the study of the various phenomenological
effects of the QCA model presented in Section
\ref{sec:weyl-dirac-maxwell}. The aim of this analysis is to
understand the properties of the QCA dynamics and compare its features
to the known results about the dynamics of free quantum fields. The
ultimate goal to identify experimental situation in which it is
possible to falsify the validity of the QCA theory.

\subsection{Zitterbewegung}
The first feature of the QCA dynamics we are going to explore (for a
more complete presentation see Ref.\cite{bisio2013dirac}) is the
appearence of a fluctuation of the position in the particle
trajectory, the so called \emph{zitterbewegung}.

The Zitterbewegung was first recognized by Schr\"odinger in 1930
\cite{schrodinger1930kraftefreie} who noticed that in the Dirac
equation describing the free relativistic electron the velocity
operator does not commute with the Dirac Hamiltonian: the evolution of
the position operator,exhibits a
very fast periodic oscillation around the mean position with frequency $2mc^2$ and amplitude
equal to the Compton wavelength $\hbar/mc$ with $m$ the rest mass of
the relativistic particle.  Zitterbewegung oscillations cannot be
directly observed by current experimental techniques for an
electron since the amplitude is very small $\approx 10^{-12}$
m.  However, it can be seen in a number of solid-state,
atomic-physics, photonic-cristal and optical waveguide simulators
\cite{lurie1970zitterbewegung,cannata1991effects,ferrari1990nonrelativistic,cannata1990dirac,PhysRevLett.100.113903}.

Here we focus on the one-dimensional Dirac QCA
whose epression,  introduced in Section
\ref{sec:weyl-dirac-maxwell}, is easily obtained as special case
of Eq.~\eqref{eq:diracauto}\footnote{More precisely,
  Eq.~\eqref{eq:weyl1D} leads to two identical copies of Eq.~\eqref{eq:Dirac1d}}
\begin{align}
  \label{eq:Dirac1d}
  U = \int_{-\pi}^{\pi} dk \ketbra{k}{k} \otimes U_k \quad
U_k =
\begin{pmatrix}
  n e^{-ik}  & im \\
  im  &n e^{ik}
\end{pmatrix}
\end{align}
The ``position'' operator $X$ corresponding to the representation $|x\>$ (i.e. such that
$X|s\>|x\>=x|s\>|x\>$, $x \in \mathbb{Z}$) is defined as follows
\begin{align}
X=\sum_{x\in\mathbb{Z}}x(I\otimes\ketbra{x}{x}).
\end{align}
and it provides the average location of a wavepacket in terms
of $\<\psi|X|\psi\>$.
If we write the single particle
in terms of its positive frequency and negative frequency components, i.e.
$|\psi\> = c_+ |\psi\>_+  +  c_- |\psi\>_-$,
yhe time evolution of the mean value of the position operator
$\bra{\psi}X(t)\ket{\psi}$ is given by
\begin{align}\nonumber
x_\psi(t) :=  \bra{\psi} X(t) \ket{\psi} =
x_{\psi}^+(t) + x_{\psi}^-(t) + x_{\psi}^{\rm{int}}(t)\\\nonumber
x_{\psi}^{\pm}(t) : = \bra{\psi _\pm} X(0) + Vt \ket{\psi_\pm}\\\label{e:zitterbewegung}
x_{\psi}^{\rm{int}}(t) :=
2 \Re [\bra{\psi_+} X(0) - {Z}_{{X}}(0) + {Z}_{{X}}(t) \ket{\psi_-}]
\end{align}
where $V$ is a time independent operator corresponding to the group
velocity and $Z_X(t)$ is the operator that gives the oscillatory
motion (see Ref. \cite{bisio2013dirac} for the details).
We notice that the interference between positive
and negative frequency is responsible of the oscillating term
$x_{\psi}^{\rm{int}}(t)$ whose magnitude is bounded by $1/m$ which in the usual dimensional units
corresponds to the Compton wavelength $\hbar/ m c$.  These results
show that $x_{\psi}^{\rm{int}}(t)$ is the automaton analogue of the
 Zitterbewegung for a Dirac particle.  for $t \to \infty$ the term $2 \Re [\bra{\psi_+} {Z}_{{X}}(t)
\ket{\psi_-}]$, which is responsible of the oscillation, goes to $0$
as $1/\sqrt{t}$ and only the additional shift contribution given by $2 \Re
[\bra{\psi_+} X(0) - {Z}_{{X}}(0) \ket{\psi_-}]$ survives.
In Fig. \ref{fig:Zitt11-12-21-22} one can se the simulation of the evolution of states with particle
and antiparticle components smoothly peaked around some $k_0$.
\begin{figure}[t!]
\includegraphics[width=.23\textwidth]{Zitt11.pdf}\quad
\includegraphics[width=.23\textwidth]{Zitt12.pdf}
\includegraphics[width=.23\textwidth]{Zitt31.pdf}\quad
\includegraphics[width=.23\textwidth]{Zitt32.pdf}
\caption{ Zitterbewegung in the one dimensiona Dirac QCA.  {\bf
    Top:} The mass of the particle is $m=0.15$. The amplitudes of the
superposition between
  positive and negative frequency states are $c_+=1/ \sqrt{2}$
  $c_-=i/\sqrt{2}$ respectively. The wavepacket is peaked around
  $k_0=0$.
 The shift and oscillation frequency are
  respectively $\bra{\psi} X(0)+ Z_{X}(0)\ket{\psi}=3.2$ and
  $\omega(0)/\pi=0.05$.  {\bf Middle:}
  $m=0.15$, $c_+=1/\sqrt{2}$, $c_-=1/\sqrt{2}$, $k_0=0$,
  $\sigma=40^{-1}$.  The shift and oscillation frequency
  are $0$ and $0.13$, respectively.  {\bf Bottom:} $m=0.13$,
  $c_+=\sqrt{2/3}$, $c_-=1/\sqrt{3}$, $k_0=10^{-2}\pi$,
  $\sigma=40^{-1}$.  In this case the particle and antiparticle
  contribution are not balanced and the average position drift
  velocity is thus $\bra{\psi _+} V \ket{\psi_+}+\bra{\psi _-} V
  \ket{\psi_-}=(|c_+|^2-|c_-|^2)v(k_0)=0.08$, corresponding to an
  average position $x_{\psi}^+(800) + x_{\psi}^-(800)=464$.
Notice that for $t \to \infty$ the
  term $2 \Re [\bra{\psi_+} {Z}_{{X}}(t) \ket{\psi_-}$, which is
  responsible of the oscillation, goes to
  $0$. This figure is published in Ref. \cite{bisio2013dirac}.}\label{fig:Zitt11-12-21-22}
\end{figure}

\subsection{Scattering against a potential barrier}\label{s:barrier}
In this section we study the dynamics of the one dimensional Dirac
automaton in the presence of a potential.
In the position representation the one particle evolution of the
one dimensioanal Dirac QCA
reads as follows:
\begin{align}
\label{eq:dirac1dposition}
  U := \sum_x
\left(
\begin{array}{ll}
   n \ketbra{x-1}{x} &   -im \ketbra{x}{x} \\
  -im \ketbra{x}{x} &   n \ketbra{x+1}{x}
  \end{array}
\right).
\end{align}
The presence of a  potential $\phi(x)$, modifies the unitary evolution
of Eq. \eqref{eq:dirac1dposition} with a position dependent phase as
follows (see also Ref \cite{kurzynski2008relativistic,meyer1997quantum}):
\begin{align*}
  U_\phi := \sum_x e^{-i \phi(x)}
\left(
\begin{array}{ll}
   n \ketbra{x-1}{x} &   -im \ketbra{x}{x} \\
  -im \ketbra{x}{x} &   n \ketbra{x+1}{x}
  \end{array}
\right).
\end{align*}
We now review the analysis (carried on in Ref. \cite{bisio2013dirac}) of the case in which
 $\phi(x) := \phi \, \theta(x)$
($\theta(x)$ is the Heaviside step function) that is a potential step  which is $0$
for $x <0 $ and has a constant value $\phi \in
[0, 2\pi]$ for
$x \geq 0$.
Let us consider the situation in which, for $t \ll 0$, the state is
a positive frequency
wavepacket peaked around $\bvec k_0$ that moves at group velocity
$v(k_0)$ and hits the barrier form the left.
Then and one can show that for $t \gg 0$ the state is evolved into a  superposition of a reflected
and a transmitted wavepacket as follows (we use the notation of
Eq.~\eqref{eq:particlestates} adapted at the one-dimensional case):
\begin{align*}
  \begin{split}
    \ket{\psi(t)} \xrightarrow{ t \gg 0} \beta(k_0) \int \df{k}{\sqrt{2\pi}}  g_{k_0}(k) e^{-i \omega(k)t}  \ket{u}_{-k}\ket{k} +\\
 {}+ \tilde{\gamma}(k_0) e^{-i \phi t} \int \df{k}{\sqrt{2\pi}}  \tilde{g}_{k'_0}(k') e^{-i
    \omega(k')t} \ket{u}_{k'}\ket{k'}
  \end{split}
\end{align*}
where we defined
\begin{align*}
&k'_0 \mbox{ s.t. }    \omega(k'_0) = \omega(k_0) -\phi, \\
&\tilde{\gamma}(k_0) := {\gamma}(k_0)
\sqrt{\frac{v(k'_0)}{v(k_0)}}, \qquad
\tilde{g}_{k'_0}(k') =
\sqrt{\frac{v(k'_0)}{v(k_0)}}{g}_{k'_0}(k')
\end{align*}
(one can check $ \int \df{k}{\sqrt{2\pi}}
|\tilde{g}_{k'_0}(k')|^2 = 1$), whose group velocities are
$-v(k_0)$ for the reflected wave packet and $v(k^\prime_0)$ for the
transmitted wave packet.

The probability of finding the particle in the reflected wavepacket is
then $R = |\beta(k_0)|^2$
(reflection coefficient) while the probability of finding the particle in the transmitted
wavepacket is $ T= |\tilde{\gamma}(k_0)|^2$ (trasmission coefficient). The consistency of the result
can be verified by checking that $R+T = 1$.
Clearly $\phi=0$ implies
$R = 0$ and increasing $\phi$ for a fixed $k$ increases the value of
$R$ up to $R=1$.
By further increasing $\phi$ a
transmitted wave reappears and the reflection coefficient decreases. This is the so called ``Klein
paradox'' which is originated by the presence of positive and negative frequency eigenvalues of the
unitary evolution.  The width of the $R=1$ region is an increasing function of the mass equal to
$2\arccos(n)$ which is the gap between positive and negative frequency solutions.
\begin{figure}[t!]
\includegraphics[width=.4\textwidth]{R-m0_4-K2.pdf}\qquad
\includegraphics[width=.4\textwidth]{Vt-m0_4-K2.pdf}
\caption{Reflection coefficient for $m=0.4$ and wave-vector of the
  incident particle $k_0=2$ as a function of the potential barrier
  height $\phi$. This Figure is published in Ref. \cite{bisio2013dirac}.}\label{fig:Section}
\end{figure}
\begin{figure}[ht!]
\includegraphics[width=.4\textwidth]{R0_25-m0_4-V1_42-k2.pdf}\qquad
\includegraphics[width=.4\textwidth]{R0_75-m0_4-V1_55-k2.pdf}\\
\includegraphics[width=.4\textwidth]{R1-m0_4-V2-k2.pdf}\qquad
\includegraphics[width=.4\textwidth]{R0_5-m0_4-V2_4-k2.pdf}
\caption{Simulations of the one dimensional Dirac automaton evolution with a square
  potential barrier. Here the automaton mass is $m=0.2$ while the
  barrier turns on at $x=140$. In the simulation the incident state
  is a smooth state of the form $\ket{\psi(0)}=  \int
  \df{k}{\sqrt{2\pi}} g_{k_0}(k) \ket{+}_k$ peaked around the positive frequency eigenstate
  $\ket{+}_{k_0}$ with $k_0=2$ and
  with $g_{k_0}$ a Gaussian having width
  $\sigma=15^{-1}$. The incident group velocity is $v(k_0)=0.90$. The
  simulation is run for four increasing values of the potential $\phi$.
  {\bf  Top-Left:}
 Potential barrier height $\phi=1.42$, reflection coefficient
  $R=0.25$, velocity of the transmitted particle $v(k^\prime_0)=0.63$.
  {\bf Top-Right:} $\phi=1.55$, $R=0.75$, $v(k^\prime_0)=0.1$.
  {\bf Bottom-Left:} $\phi=2$, $R=0.1$, $v(k^\prime_0)=0$.
  {\bf Bottom-Right:} $\phi=2.4$, $R=0.50$, $v(k^\prime_0)=0.33$. This Figure is published in Ref. \cite{bisio2013dirac}.} \label{fig:Simulation}
\end{figure}

In Fig. \ref{fig:Section} we plot the reflection $R$ coefficient and the transmitted wave velocity
group $v(k_0')$ as a function of the potential barrier height $\phi$ with the incident wave packet
having $k_0=2$ and $m=0.4$. From the figure it is clear that after a plateau with $R=1$ the
reflection coefficient starts decreasing for higher potentials. In Fig. \ref{fig:Simulation} we show
the scattering simulation for four increasing values of the potential, say
$\phi=1.42,\,1.55,\,2,\,2.4$.

\subsection{Travel-time and Ultra-high energy cosmic rays}

The approximated evolution studied in
Section~\ref{sec:interp-hamilt-diff}
provide a useful analytic tool for evaluating the
macroscopic evolution of the automaton.
We now consider an elementary experiment,
based on
particle fly-time, which compares
the Dirac automaton evolution with the one given by the Dirac equation.

Consider a protonwith $m_p\approx 10^{-19}$ and wave-vector peaked around
$k_{CR}\approx 10^{-8}$ in Planck units\footnote{As for order of
  magnitude, we consider numerical values corresponding to ultra high
  energy cosmic rays (UHECR) \cite{takeda1998extension}} , with a spread
$\sigma$ of the wave-vector. We ask what is the minimal time $t_{CR}$ for observing a
complete spatial separation between the trajectory predicted by the
cellular automaton model and the one described by the usual Dirac
equation. Thus we require the separation between the two trajectories
to be greater than $\hat\sigma=\sigma^{-1}$ the initial proton's width
in the position space. We approximate the state evolution of the
wave-packet of the proton using the differential equation
\eqref{eq:schroding2} for an initial Gaussian state. The time required to
have a separation $\hat\sigma$ between the automaton and the Dirac
particle
is
\begin{align}\label{eq:flying-time2}
t_{CR}\approx 6\frac{\hat\sigma}{m_p^2}.
\end{align}
and for $\hat\sigma=10^2\text{fm}$ (that is reasonable for a proton wave-packet) the flying time request for complete
separation between the two trajectories is
$
t_{CR}\approx 6\times 10^{60}$ Planck times, i.e. $\approx 10^{17}s$,
a value that is comparable with the age of the universe and then
incompatible with a realistic setup.

\subsection{Phenomenology of the QCA Theory of Light}
\label{sec:phen-qca-theory}
In this section we present an overview of the new phenomenology
emerging
 from QCA theory of free electrodynamics presented in Section
 \ref{sec:weyl-dirac-maxwell}. For a more detailed presentation we
 refer to Ref. \cite{bisio2014quantum}.

\subsubsection{Frequency dependent speed of light}
From Eq.~\eqref{eq:maxwell2} one has that
the angular frequency of the electromagnetic waves
is given by the modified dispersion relation
\begin{align}
\label{eq:modifieddisprelmax}
\omega(\bvec k) = 2 | \bvec{n}_{\tfrac{\bvec k}{2}} |   .
\end{align}
and the usual relation $\omega(\bvec k) =  | \bvec k  |  $
is recovered in only the $| \bvec k  | \ll 1$ regime.
The speed of light is the group velocity of the electromagnetic
waves, i.e.~the gradient of the dispersion relation. The major consequence
of Eq. \eqref{eq:modifieddisprelmax} is that the speed of light depends on
the value of $\bvec k$, as if the vacuum were a dispersive medium.

The phenomenon of a $\bvec k$-dependent speed of light is also studied in the
quantum gravity literature where many authors considered the hypothesis that the existence of an invariant
length (the Planck scale) could manifest itself in terms of dispersion
relations that differ from the usual relativistic one
\cite{ellis1992string,lukierski1995classical,Quantidischooft1996,amelino2001testable,PhysRevLett.88.190403}.
In these models the $\bvec k$-dependent speed of light $c(\bvec k)$, at the leading order in $k :=| \bvec k |$,
is expanded as $c(\bvec k) \approx 1 \pm \xi k^{\alpha}$, where $\xi $ is a numerical factor of order
$1$, while $\alpha$ is an integer.  This is exactly what happens in our framework, where the
intrinsic discreteness of the quantum cellular automata leads to the dispersion relation of
Eq. \eqref{eq:modifieddisprelmax} from which the following $\bvec k$-dependent speed of light
\begin{align} \label{eq:freqdepsol}
  c(\bvec k) \approx 1 \pm 3\frac{k_x k_y k_z}{|\bvec k|^2} \approx
 1 \pm \tfrac{1}{\sqrt{3}}k,
\end{align}
can be obtained by computing the modulus of the group velocity and power expanding in $\bvec k$ with the
assumption $ k_x = k_y = k_z = \tfrac{1}{\sqrt{3}} k $, $k =
|\bvec k|$. The $\pm$ sign in Eq.~\eqref{eq:freqdepsol} depends on whether
we considered the $A^{+}(\bvec k)$ or the $A^{-}(\bvec k)$
Weyl QCA. This prediction can possibly be experimentally tested
in the astrophysical domain, where tiny corrections are
magnified by the huge time of flight.  For example, observations of the arrival times of pulses
originated at cosmological distances, like in some $\gamma$-ray
bursts\cite{amelino1998tests,abdo2009limit,vasileiou2013constraints,amelino2009prospects}, are now
approaching a sufficient sensitivity to detect corrections to the relativistic dispersion relation
of the same order as in Eq. \eqref{eq:freqdepsol}.

\begin{figure}[ht]
  \begin{center}
    \includegraphics[width=.55\textwidth]{graph_vectors_QCATOL-rc3.pdf} \,
    \includegraphics[width=.42\textwidth]{graph_wave_QCATOL-beta-c.pdf}
    \caption{(colors online) {\bf Left:} the graphics shows the vector $2\bvec{n}_{\tfrac{\bvec k}{2}}$ (in green),
      which is orthogonal to the polarization plane, the wavevector $\bvec k$ (in red) and the group
      velocity $\nabla \omega (\bvec k)$ (in blue) as function of $\bvec k$ for the value $|\bvec k|= 0.8$ and
      different directions.
{\bf Right: }A
      rectilinear polarized electromagnetic wave. We notice that the
      polarization plane (in green) is sligtly tilted with respect the
      plane orthogonal to $\bvec k$ (in gray). This Figure is published in
    Ref. \cite{bisio2014quantum}.}
    \label{fig:relvectors}
  \end{center}
\end{figure}

\subsubsection{Longitudinal polarization}
A second distinguishing feature of Eq.~\eqref{eq:maxwell2} is that the polarization plane is neither
orthogonal to the wavevector, nor to the group velocity, which means that the electromagnetic waves
are no longer exactly transverse (see Fig. \ref{fig:relvectors}).
The angle $\theta$ between the polarization plane and the plane orthogonal to $\bvec k$ or
$\nabla\omega(\bvec k)$ is of the order $\theta \approx 10^{-15}\mathrm{rad}$ for a
$\gamma$-ray wavelength, a precision which is not reachable by the present technology. Since for a
fixed $\bvec k$ the polarization plane is constant, exploiting greater distances and longer times does
not help in magnifying this deviation from the usual electromagnetic theory.

\subsubsection{Composite photons and modified commuation relations }

Finally, the third phenomenological consequence of the QCA theory of light is
the deviation from the exact Bosonic statistics due to the composite
nature of the photon.  As shown in Ref.\cite{bisio2014quantum},
the choice of the function ${f}_\bvec k(\bvec{q})$ in Eq.~\eqref{eq:electric and magnetic field}
determines the regime where the composite photon can be approximately treated as a Boson.  However,
independently on the details of function ${f}_\bvec k(\bvec{q})$, one can prove that a Fermionic
saturation of the Boson is not visible, e.g. for the most powerful laser \cite{dunne2007high} one
has a approximately an Avogadro number of photons in $10^{-15}$cm${}^3$, whereas in the same volume
on has around $10^{90}$ Fermionic modes.
Another test for the composite nature of photons is provided by the prediction of
deviations from the Planck's distribution in Blackbody
radiation experiments. A similar analysis was carried out in
Ref. \cite{perkins2002quasibosons}, where the author showed that the predicted
deviation from Planck's law is less than one part over $10^{-8}$,
well beyond the sensitivity of present day experiments.

\section{Future perspectives}

We conclude this paper with an overview of the future developments of
the research program on QCA for Field Theory.

\subsection{Lorentz covariance and Deformed Relativity}
\label{sec:lorentz-covar-deform}
Because of the intrinsic discreteness of the model, a dynamical evolution described in terms of a QCA cannot satisfy the usual Lorentz covariance, which must break down at the Planck scale. Moreover the very notions of spacetime and boosted reference frame break down at small scales, and need a thoughtful reconsideration. In Ref.~\cite{bisio2015lorentz} a definition of reference frame was introduced in a background-free scenario, in terms of labelling of irreducible representations of the group $G$. The Lorentz symmetry is then recovered by imposing a generalized relativity principle on possible changes of reference frame, allowing only those changes that leave the automaton invariant. A preliminary analysis of the one-dimensional case can be found in Ref.~\cite{bibeau2013doubly}, where only the necessary condition of preserving the dispersion relation was considered. Focusing on the one dimensional Dirac QCA we have
\begin{align}
  \label{eq:dispreldirac1d}
  \omega(k) = \arccos(\sqrt{1-m^2} \cos(k))
\end{align}
and one can see that in the $k \ll 1$,$m \ll 1$ limit Eq.~\eqref{eq:dispreldirac1d} reduces to the usual relativistic
dispersion relation $\omega^2=k^2+m^2$. It is also immediate to check that the automaton dispersion relation of
Eq.~\eqref{eq:dispreldirac1d} is not invariant under standard Lorentz transformation.  In order to preserve Eq.~\eqref{eq:dispreldirac1d} one needs to introduce a non-linear representation of the Lorentz transformation in the wave-vector space---as proposed in the so
called \emph{deformed special relativity} (DSR) models
\cite{amelino2002relativity,amelino2001planck,amelino2002quantum,amelino2001testable,PhysRevLett.88.190403,magueijo2003generalized}.

In Ref.~\cite{bisio2015lorentz} the boosts preserving the three-dimensional Weyl automaton were then derived in the form of the following non-linear representation of the Lorentz group
\begin{equation}\label{eq:non-linear-boosts}
  {L}^D_\beta:= \mathcal D^{-1}\circL_\beta\circ\mathcal D,
\end{equation}
where $\mathcal D:\Reals^4\to\Reals^4$ is a non-linear map.
The specific form of $\mathcal D$ gives rise to a particular frequency/wave-vector Lorentz deformation.

These ideas can also be applied to the three-dimensional Dirac QCA. In this case one can show that a change of the rest mass should be involved in the representation of boosts, in order to obey our generalized relativity principle. Interestingly, this unexpected feature gives rise to an emergent space-time with a non-linear de Sitter symmetry instead of the Lorentz one.

Another challenging line of research is to characterize the emergent spacetime of the QCA framework. The DSR models provide a complete description of Lorentz symmetry in frequency/wave-vector space but there are heuristic
ways to extend this framework to the position-time space. Relative locality \cite{amelino2011relative,amelino2013relative},
non-commutative spacetime \cite{connes1991particle} and Hopf algebra symmetries \cite{lukierski1991q,majid1994bicrossproduct}
have been considered in order to give a real space formulation of deformed relativity.

Finally, we would like to stress that space-time emerges from: i) the structure of the group $G$, ii) the specific expression of the automaton and iii) the generalized relativity principle, while all these concepts do not require any space-time background. Thus, outside the limits in which the relativistic approximations hold, the very structure of our usual space-time break down, substituted by other counterintuitive effects. In particular this is true in all physical situations where the discrete structure of the lattice $G$ becomes relevant.

\subsection{Thermodynamics of free ultra-relativistic particles and QCA}

\label{sec:therm-free-ultr}
Most of the analysis that we presented in this paper was focused on
the dynamics of one particle state and the deviations of this
kinematics from the usual relativistic one. On the other hand, it would be
interesting to explore the QCA phenomenology when the number of
particle goes to infinity, namely a thermodynamic limit. Since the QCA
we are considering descrive a non-interacting dynamics, the
thermodynamic that will emerge will describe a gas of free particles.
However, since the dispersion relation of the QCA differs from the
relativistic one, the density of states will be different.

In the case
of free fermions this will result in a shift of the Fermi energy
that could become relevant when the number of fermions becomes very
large. One could for example analyze how the
Chandrasekhar limit of white dwarfs is modified in this context (see
Ref. \cite{amelino2012uv,camacho2006white} for a similar analysis in a
different context).

\subsection{Interacting extensions}

The theory of linear QCAs naturally leads to free quantum field theories. In order to introduce interactions, one needs to relax the linearity assumption. This can be done by splitting the computational step in two stages, the first one acting linearly, and the second one representing a nonlinear and completely local evolution. This can be motivated in terms of a time-local gauge symmetry that must preserve some local degree of freedom, in particular the local number of excitations. This simple modification of the linear automaton introduces a non-trivial interaction, making the automaton non-trivially reducible to a quantum walk. Preliminary analysis shows that this minimal relaxation of linearity is sufficient to give rise to couplings that might reproduce the phenomenology of quantum electrodynamics.

\section*{Acknowledgements}
This work has been supported in part by the Templeton Foundation under the project ID\# 43796 {\em A
  Quantum-Digital Universe}.

\bibliographystyle{apsrev4-1}