Giorgio Krstulovic

Stéphan Fauve

Fabio Sciarrino

Stéphane Attal

Carlo Di Franco

Christopher Cedzich

Fumiaki Matsuoka

Giuseppe Di Molfetta
email: giuseppe.dimolfetta@ens.fr

Lauchlan Honter

Ben B Luo

Tatsuaki Wada

Yutaka Shikano
email: yshikano@ims.ac.jp

G Di Molfetta

Massless Dirac equation from Fibonacci discrete-time quantum walk

Keywords: Discrete Time Quantum Walks: From synthetic gauge fields to spontaneous equilibration Discrete Time Quantum Walks, Quantum Simulations, Synthetic Gauge Fields Quantum Decoherence, Quantum Gauge Lattice, Thermalization Marches quantiques à temps discret, Simulation quantique, Champs de jauge de synthèse, Décohérence quantique, Théories de jauge sur réseau, Thermalisation Quantum Walk, Massless Dirac equation, Quantum Simulation

J'ai commencé ma thèse en septembre 2012 au Laboratoire d'Études du Rayonnement et de la Matière en Astrophysique et Atmosphères (LERMA) de l'Observatoire de Paris, dont je tiens à remercier l'exdirecteur Michel PÉRAULT et l'actuel directeur Darek LIS pour leur accueil au sein de leur unité. Je remercie également Maryvonne GERIN pour m'avoir accueilli dans son èquipe du LRA à l'Ecole Normale Supérieure où j'ai passé mes deux premières années de these, et Chantal STEHLE pour avoir fait de même à Jussieu pendant ma dernière année.

Cette thèse a été dirigée par Fabrice DEBBASCH, que j'ai rencontré pendant mon premier stage de M1 à l'Université Pierre Marie Curie. Avec lui, j'ai découvert le sujet de ce doctorat et j'ai échangé sur une infinité de sujets scientifique. Je le remercie pour la confiance et le soutien dont il a fait preuve pendant tout mon travail, jusqu'aux derniers jours de rédaction du manuscrit. Ce qu'il ma donné de plus précieux, c'est son grand enthousiasme pour la science, la recherche et la découverte, que je garderai dans mon futur professionnel. Je voudrais également remercier Marc-Etienne BRACHET pour m'avoir co-encadré pendant ces trois ans. Notamment, je lui suis reconnaissant pour sa patience, son expertise scientifique et la transmission d'une partie de son énorme connaissance. À lui va aussi ma gratitude sa profonde humanité.

Je voudrais également remercier Martine BEN AMAR pour m'avoir accueilli dans le Master de Physique des systèmes complexes et Jean-Bernard ZUBER pour sa disponibilité pendant mes premiers mois de permanence à Paris. J'exprime ma plus profonde reconnaissance à Pablo ARRIGHI et Dieter MESCHEDE qui m'ont fait l'honneur d'accepter d'être rapporteurs de la prèsente thèse (thank you Prof. P. ARRIGHI and Prof. D. MESCHEDE to have accepted to be the examiners of this thesis). J'associe à ces remerciements les autres membres de jury: Jean-Michel RAIMOND et Yutaka SHIKANO. En particulier je remercie J.-M. RAIMOND pour l'intérêt qu'il a montré envers les résultats que j'ai obtenu lors de mon premier stage. Je suis aussi reconnaissant envers Y. SHIKANO pour son accueil au sein de son unité à Okazaki. Mon séjour au Japon, qui a été la plus belle et intense expérience professionnelle pendant mon doctorat, aurait été très difficile sans la disponibilité et la sympathie de Naoko KONDO, Mayuko KATO et tous les autres membres du IMS. Je remercie parement Beatrice GUIBAL du LERMA pour avoir géré toutes mes questions administratives. Pendant ces trois années de travail, j'ai eu l'occasion de profiter de la grande expertise de

INTRODUCTION

Quantum simulation has recently established itself as an area of study in quantum physics that merges fundamental and applied questions. Such an interaction results in a more operational understanding of some aspects of quantum mechanics in terms of nature description. The idea to simulate the dynamics of a quantum system by a quantum device was first introduced by Richard Feynman in a seminal article in 1982 [19] and developed in different frameworks, from quantum optics to condensed matter physics. A way to describe quantum systems and their dynamics within a computational perspective is given by the large class of Quantum Cellular Automata (QCA), that is a grid of quantum autonomous systems, interacting through local rules [44,37,8,7,38]. Over the last years, many definitions and models of QCA have been proposed. For instance, Schumacher and Werner [37] proposed a class of reversible QCA, where the evolution of the QCA is given by local unitary operations that act on the neighborhood of a given cell, and it is further updated by a single cell unitary operation. This is in accordance with the expected microscopic fully quantum dynamics. However in recent years an agreement has been achieved and many definitions have been proven to be equivalent by Arrighi and Grattage [6]. This class of quantum systems displays a wide range of complex phenomena and it has been extensively used for the comprehension of Quantum Field Theory (QFT). Bialynicki-Birula [9] first underlines the connection between QCA and QFT. Then Meyer [29] showed that QCA mostly correspond to lattice gas models and are closely related to the Dirac equation; in addition to this, again Meyer and Shakeel [30], have recently proved that QCA with no particle interpretation exist and that they propagate information.

Quantum Walks (QWs) are a special case of reversible QCA , but in the sigle particle sector. Differently from QCA, they describe the unitary dynamics of one quantum particle. They have been introduced independently by Grossing and Zeilinger [22] and Y. Aharonov et al. [3] and then extended systematically on graphs by D. Aharonov et al. [2]. QWs have been, later, fully mathematically examined by Konno [27].

INTRODUCTION

in the continuous limit.

Nowadays QWs can be realized experimentally with a large spectrum of physical objects and setups. A first natural way to physically implement a QW is in optical networks. For instance, Knight et al. [26] realized experimentally the first optical implementation of a QW in an optical ring cavity, as a modified experimental setup of the optical linear quincunx described by Bouwmeester et al. [10], the Optical Galton Board.

A second way to implement QWs is the Nuclear Magnetic Resonance (NMR), proposed by Ryan et al. [33] for a DTQW, subsequent to the work of Du et al. [15]. This method consists in manipulating and coupling spin orientations, in constituting quantum logic gates and in performing unitary operations on the nuclei quantum state in the magnetic field.

Other possible implementations can be realized in Quantum Electrodynamics (QED). In effect, the physical scheme proposed by Y. Aharonov in its seminal paper [3] has inspired Agarwal and Pathak [1] to realize a DTQW in this framework. They implemented it by injecting a single Rydberg atom into an optical cavity and by driving it by a strong external field. The effective Hamiltonian introduced in their article reproduces perfectly the Aharonov's scheme. A more recent implementation of QW in cavity QED is due to Sanders et al. [34].

Other techniques range from quantum optics (Zhang et al. [45], [46]), to ion traps (Travaglione and Milburn [43]), and to neutral atom traps (Chandrashekar [12],Eckert et al. [17], Dür et al. [16]). In particular, quite useful for the implementation of QWs on cycle is the interacting quantum dots architecture presented by Solenov and Fedichkin [40] in a solid state framework.

Recent developments, however, have shown an increasing interest for two or more correlated QWs, due to the large class of emergent quantum phenomena (Peruzzo et al. [32], Sansoni et al. [35], Schreiber et al. [36] and more recently Defienne et al. [14]). However, this thesis is dedicated to the study of the single-particle dynamics and we will not take into account multi-particle systems.

As for the physical implementations, the potential applications of QWs are too numerous to be listed here. They spread from search algorithms [13] and graph isomorphism [4,21] to modeling and simulating quantum dynamics [36] and classical system [41]. We can redirect the reader to Lovett et al. [28] for an introductory review on the quantum walk algorithm. Concerning modeling quantum dynamics we can recall some remarkable applications in condensed matter. For instance, in topological insulators realization [25,24], for spintronic applications [23], and for quantum transport on 2-dimensional layer, as in graphene-like materials [20].

This thesis is specifically devoted to only one of the various DTQWs' applications: quantum simulation of quantum systems dynamics. In particular, let us remind the reader that, rather than a monograph, this manuscript appears as a collection of our main publications with an accurate introductory section to each of them. All our works deal with QWs propagating in one spatial dimension.

The thesis is organized in three main parts and four chapters. In the first part, the Chapter (1) is devoted to introduce DTQWs and report their main features, largely explored in the recent literature. In particular, we present the method to perform the continuous limit and its re-lated sufficient conditions. Moreover, we formally prove that not all families of homogeneous DTQWs (HDTQWs) admit such conditions. In Chapter (2), we extend the same analysis to the inhomogeneous DTQWs (IDTQWs) and we introduce the idea of synthetic gauge field. Above all, we prove that some special families of IDTQWs mimic the propagation in curved (1 + 1)spacetime of a Dirac fermion, coupled to a gauge electric potential. Let us recall that similar results have been recently obtained by our group in a (1+2) spacetime dimension (Arnault and Debbasch [5]). We conclude the first part with Chapter (3), introducing a randomized version of the IDTQW. We analytically provide detailed calculations demonstrating the loss of coherence and the emergent classical diffusive behavior of the walk.

In the second part and last chapter, Chapter (4), our analysis concerns the very new problem of equilibration in QWs. Generally, equilibration is the process through which a system, under a given dynamics, reaches a steady state. By thermalization we mean that this steady state can be interpreted as a thermal state at a given temperature. Thermalization is observed typically in nonlinear systems, but all QWs describe the motion of a single particle and are linear by definition. Moreover nonlinear QW-like implementations have been considered in different contexts. In the context of optics, the so-called Optical Galton Boards (OGB) [11] have been generalized to include nonlinear Kerr effects by Navarrete-Benlloch et al. [31] or in Shikano et al. [39], who proposed a DTQW with feed-forward quantum coin displaying high nonlinearities. In those models very complex and rich behaviors have been observed. The aim of this last chapter, and our main contribution, is to observe and characterize thermalization phenomena in a spatially-periodic version of the Nonlinear OGB.

INTRODUCTION

Bibliography INTRODUCTION [37] [15,26,46], genetic sequence location [45,28,33], optimal search strategies [29,4,36], diffusion and mobility in materials [25,13,44], exchange rate forecasts in economical sciences [30,6,24] and information spreading in complex networks [35,34,41]. Furthermore, they can successfully implement efficient algorithms, for example they can solve differential equations [3,20], optimization [4,5] and clustering problems [38,48]. Random walks spread to every domain of science for more than a century and are still an important source for researchers nowadays.

Even though in the second half of last century the interest in a quantum analogous to the classical stochastic process led to further investigation of quantum mechanics and quantum information, since the 1960s, numerous scientists [42,21,17,14,12,1,23,37] have extended the Brownian motion and stochastic calculus to particles which exhibited quantum effects. Schwinger was the first to demonstrate the importance of coherence effects in the evolution of a Brownian quantum particle. Just few years later, Fjeldso et al. [11] proposed the first discrete model with the intention of recovering a quantum version of the CRW. This first example represented a quantum planar rotor on a lattice which dynamics can be approximated by a random walk. Another important milestone was the work of Gudder [19] who studied systematically quantum Markov processes and established formally a connection with Feymann path integral formalism.

Interestingly some years later, Godoy and Fujita [16] and Grossing and Zeilinger [18] had independently the same intuition of developing a quantum analogous to the CRW in discrete time and discrete space. Godoy proposed a one-dimensional Markovian quantum walk displaying a diffusive behavior and Grossing and Zeilinger [18] formalized the first model of unitary Quantum Walk (QW), which dynamics was fully propagative. Let us remark that, although the latter was defined as a QCA this model is completely equivalent to the one proposed later by Aharonov et al. [2].

Thanks to its features, especially the unitarity, QWs were immediately considered a new and efficient tools for solving, in a wider range of applications, technical problems in a more convenient way than classical random walks.

From Classical to Quantum random walks

In explicating the main differences between DTQW and its classical counterpart, we briefly recall the definition of a discrete time random walk (DTRW) on an unrestricted line. Let us imagine a man moving along a line, taking, at random, steps to the left and to the right with equal probability. The step are of unit length so that his position can take on only the value n, where n 2 Z. We want to know with what probability he reaches, at a given point, a distance n from the origin after a given elapsed time. The traditional way to solve this problem is to allow the walker to take steps at time N (N 2 N) at which time he must jump either left or right, with equal probability. The DTRW can be written:

P (n, N + 1 | n 0 , N 0 ) = 1 2 [P (n + 1, N | n 0 , N 0 ) + P (n °1, N | n 0 , N 0 )] (1.1)
where P (n, N | n 0 , N 0 ) is the probability to find the particle at time N at position n, if at time N 0 it was at n 0 . Once solved this equation and fixing one initial condition, the probability distribution is known as the Bernoulli distribution. It gives the probability of a total of n heads in tossing an unbiased coin N times:

P (n, N | 0, 0) = µ 1 2 ∂ N N ! ∑µ N °n 2 ∂ ! µ N + n 2 ∂ ! ∏ °1 (1.2)
Therefore, a DTRW, on a line, is defined in terms of classical particle's probabilities, and it may show a very rich and complex dynamics (Revesz [39]). Notably the square root of variance ae CW is proportional to p n.

Now, if we imagine the particle to be quantum and completely isolated, we should conclude that there is at least one fundamental difference from the classical case: a quantum analogous of CRW is described by complex amplitudes and not by probabilities. Moreover, the quantum nature of the particle implicates the existence of entanglement states, i .e. non factorizable states, that have no analogous in classical physics. The presence of internal quantum states is another striking difference, and we will see, for instance, that it will recover a fundamental role in the definition of QW.

Quantum walks 11

Now, let us introduce more properly the theoretical background for the DTQW.

General Setup of a Discrete Time Quantum Walks

Let K be a finite dimensional Hilbert space and let H = K ≠ Z d . The Hilbert space H is that of a QW in Z d with k=dim(K) internal degrees of freedom. The space K is called the coin space and we denote the canonical basis of K by {|b ø i} ø2K ={1,...,k} . The projectors on the basis of {|b ø i} are noted by P ø = |b ø ihb ø |. We shall indicate the canonical basis of Z d by {|mi} m2Z d so that any vector ™ 2 H can be written as:

™ = X m2Z d ™ m |mi, (1.3)
where

™ m = X ø2K √ ø m |b ø i 2 K, (1.4)
and

||™ m || 2 < 1 where ||•|| 2 = ( P m2Z d |•| 2 )
1 2 . We specify, in the following, the elements of the corresponding orthonormal composite basis of H by |b ø , mi = |b ø i ≠ |mi 2 H. Let us now introduce the shift operator T defined by:

T = X m2Z d X ø2K P ø ≠ |m + S(ø)ihm|. ( 1.5)
The walker at the site m with internal degree of freedom b ø , represented by the vector |b ø , mi, is just sent by T to one of the neighboring sites depending on ø determined by the shift function S : K ! Z d . The action of T is given as follows:

T |b ø , mi = |b ø , m + S(ø)i.

(1.6)

Consider now B(K) the set of complex valued matrices acting on K and consider the map:

B(K) : √ ! B(K)[™].
(1.7)

The one step evolution of the quantum walk is given by the map U :

H ! H ™ = X m2Z d ™ m |mi ! U [™] = T B(K)[™]. (1.8)
In order to recover a one space dimensional QW with two internal state, we simply choose d =1 and K = C 2 , whose basis we signal by {b ø } ø2K and K = {+1, °1}. The shift function corresponds to:

S : {+1, °1} ! Z, (1.9)
where S(±1) = ±1.

We then choose the quantum coin B as an element of U (2):

B AE,µ,ª,≥ = e i AE √ e i ª cos µ e i ≥ sin µ °e°i≥ sin µ e °i ª cos µ ! (1.10) where the set of the four parameters (AE, µ, ª, ≥) lives in R 4 . Note that B AE,µ,ª,≥ 2 SU (2) if AE = kº, k 2 Z.

Symmetries and properties of the above quantum coin have been widely investigated by Chandrashekar et al. [9]. Let us observe for µ = º/4, AE = °ª = º/2 and ≥ = º/2, the above quantum coin recover the so called Hadamard coin, B H defined as below:

B H = 1 p 2 √ 1 1 1 °1! (1.11)
Also note that the global phase exp(i AE) does not play any role in homogeneous walks and notably in the time evolution of the probability density. However, it will play an important role in the next chapter when the quantum coin will depend on time and space points of the lattice. 1 , where √ 0,m is a complex regular function, at least twice differentiable, and so that ||√ 0,m || 2 < 1. The couple of constants (c L , c R ) 2 C 2 .

Let us introduce the notation

™ j = U • ... • U [™ 0 ], where ™ 0 = P m √ 0,m (c L |b L i + c R |b R i)
The finite difference equations after 1-fold iteration of the map U , for each amplitude √ L j ,m = hb L , m|™ j and √ R j ,m = hb R , m|™ j are given by: √ L j +1,m = e i AE (cos(µ)e i ª √ L j ,m°1 + sin(µ)e i ≥ √ R j ,m°1 ) (1.12)

√ R j +1,m = e i AE (°sin(µ)e °i ≥ √ L j ,m+1 + cos(µ)e °i ª √ R j ,m+1 ) (1.13)

The above equations usually define in the literature DTQW evolution. Nonetheless, there are other remarkable definitions of DTQW. For instance, let us recall the embedded QW Chandrashekar [8], a generic quantum walk model that uses a quantum coin embedded in the unitary shift operator T , where an additional degree of freedom is added. Let us state K = K 1 ≠K 2 where K 1 is the usual internal coin space and K 2 is an external coin space. The main role of the T is to move the particle in the superposition of position space at each time step, eliminating the need for a separate coin toss operation after every unitary step.

Another model definition is introduced in some of the publications included in this thesis, which simply consists in inverting the shift operator T and the action of the quantum coin B :

™ = X m2Z d ™ m |mi ! U [™] = B(K)T [™]
(1.14)

Then the amplitude, coding the probability of the particle to go towards the left √ L j ,m and towards the right √ R j ,m after one time step, is:

√ L j +1,m = e i AE (cos(µ)e i ª √ L j ,m°1 + sin(µ)e i ≥ √ R j ,m+1 ) (1.15)
√ R j +1,m = e i AE (°sin(µ)e °i ≥ √ L j ,m°1 + cos(µ)e °i ª √ R j ,m+1 ) (1.16) Observe that this definition of the walk is totally equivalent to the usual standard definition in case of homogeneous and static quantum coin and symmetric density time evolution. Let us conclude with a comment on the formalism which will be useful in the following chapters. Suppose we want to represent the entire history of a quantum walk observed through a stroboscope of period n. For all (n, j ) 2 N 2 , the collection W n j = (™ k,m ) k=n j,m2Z .

(1.17)

This collection represents the state of the quantum walk at time k = n j . For any given n, the 1 More generally ™ 0 = P m √ 0,m (c L |b L i+e i ¡ c R |b R i where e i ¡ is a complex phase which can be tuned to modify the symmetry of the initial condition.

collection

S n = (W n j ) j 2Z , (1.18) where the final state observed through a stroboscope of period n. The evolution equations for S n are those linking W n j +1 to W n j for all j . These can be deduced from the original evolution (1.12, 1.13) or (1.15, 1.16) from the walk, which also coincides with the evolution equations of S 1 .

Qualitative description

A general qualitative description of quantum walk dynamics has been largely reported by Venegas-Andraca [47]. We recall here the main features of a DTQW moving on unrestricted line. For simplicity let us consider here the simplest version of the DTQW, the Hadamard walk. The discrete time equations read:

√ L j +1,m = 1 p 2 (√ L j ,m°1 + √ R j ,m°1 ) (1.19) √ R j +1,m = 1 p 2 (√ L j ,m+1 °√R j ,m+1 ) (1.20)
where the quantum coin is the Hadamard Coin (or Hadamard gate), B H has been defined in (1.11).

The Hadamard walk has been extensively studied and it is known that the final distribution depends on the initial state of the particle. For instance, let us acknowledge an initial condition as ™ 0 = P m √ 0,m (c L |b L i + e i ¡ c R |b R i), where √ 0,m is normally distributed and centered about m = 0 and ¡ 2 R. From Fig. 1.1 we can notice that (i) the probability distribution is significantly different from the DTRW's one, which is Gaussian; (ii) the unitary operator treats differently the initial condition if (a) c L = 1, c R = 1, ¡ = 0 or (b) c L = 1, c R = °1, ¡ = 0. Between the two cases, there is a phase difference and the asymmetry arises from the constructive interference on one side, and from the destructive interference on the other side of the position space. In order to obtain a symmetric final state we can prepare a symmetric superposition state, as for example:

™ 0 = X m √ 0,m 2 (|b L i + i |b R i). (1.21)
where c L = 1, c R = 1, ¡ = º 2 . Let us remark that for a Hadamard walk on the line, it is proved by Venegas-Andraca [47] that after N steps, the probability distribution is spread on the interval [°N p ] and decreases quickly outside this region, as we can see in Fig. 1.1.

The approximation of the root of the variance, for a large number of steps, varies as following: .22) Let us now examine the U (2) quantum coin as in Eq. (1.10). It is possible to disregard the role of the particle's initial state, in order to obtain a symmetric amplitude distribution, just tuning the three parameters (µ, ª, ≥). In fact, let us consider P ± encoding the probability to find the particle in position m = 1 and in position m = °1. Starting from a symmetric initial state

ae j / s 1 °1 p 2 j. ( 1
= P m √ 0,m (c L |b L i + e i ¡ c R |b R i) T where (i) c L = 1, c R = 1, ¡ = 0 or (ii) c L = 1, c R = °1, ¡ = 0 and symmetric initial state c L = 1, c R = 1, ¡ = º/2.
The spatial component is given by √ 0 = P m ± 0,m ™ 0 , with c L = 1, c R = 1, ¡ = 0, obeying to the Eqs.(1.12, 1.13), the probabilities P ± read:

P QW ®
= (1 ± sin(2µ) cos(≥ °ª)) (1.23) This probability distribution will be equal and lead to a left-right symmetry in position if ≥ °ª = (2k + 1) º 2 , k 2 N2 . Otherwise the parameters ª and ≥ will introduce bias in the spatial probability distribution from the very first step. Concerning again the variance of the walk, the three Euler angles do not change qualitatively the dynamics in respect to the Hadamard case, which always remains ballistic. However, now the root of the variance ae is a function of all the three Euler angles: ae / K µ,ª,≥ j.

(1.24)

A simple case can be shown taking into account, for instance, all the angles equal to zero except µ and K µ / p 1 °sin µ [7]. Observe in Fig. 1.4 that the value of K µ decreases as µ increases. In particular in the limit of µ=º/2 the walker does not spread, but it is localized around the origin, therefore K µ = 0. Conversely, if µ= kº, k 2 Z, the particle spreads with K µ = 1. Intermediate values of the variance are dominated by the interference effects, and the probability distribution spreads, after N steps, over the interval (-N cos µ, N cos µ). As in the Hadamard walk, the probability distribution decreases rapidly beyond |N cos µ|. The analytical description of QWs dynamics has represented over the last few years an important field of investigation. Two methods have been extensively used: (i) the Fourier method and (ii) the combinatorial approach. In the former, the QW dynamics is studied in Fourier space in order to get a closed-form of the coin amplitude equations and, therefore, computing statistical properties of the walker. It was first introduced by Nayak and Vishwanath [32] and later by Kosik [27]. In the combinatorial approach, we shall compute the amplitude for a particular position m component by summing up the amplitudes of all the paths which begin in the given initial condition and up in the same position m. This approach is also called discrete path integral approach and was developed mainly by Machida and Konno [31].

θ = π/4 θ = π/3 θ = π/8 θ = π/16 |Ψ 100,m | 2 σ ∝ (1-sinθ)
Here, we briefly recall the first method for the family of QWs obeying to Eqs. (1.12) and (1.13). Note that this could easily be extended to Eqs. (1.15) and (1.16). The closed-form of the coin amplitude is straightforward if we consider Eqs. (1.12) and (1.13) in Fourier space. Then let us look at Tj,k = P m ™ j ,m e i km Fourier transformed state of ™ j ,m . The finite difference equations read:

√ √L j +1,k √R j +1,k ! = e i AE √ e i (ª°k) cos µ e i (≥°k) sin µ °e°i(≥°k) sin µ e °i (ª°k) cos µ !√ √L j ,k √R j ,k ! (1.25)
or more compactly Tj+1,k =U k Tj,k , where U k is the Fourier transform of the unitary operator driving the walk and it is local. In order to compute the solution after N steps, we need to compute (U k ) N and this is easily accomplished if we diagonalize the unitary matrix associated with the operator U k . A unitary matrix is always diagonalizable as a consequence of the spectral theorem. Therefore U k is diagonalizable:

V k = O °1U k O (1.26)
where O are unitary matrix and V k is diagonal. The eigenvalues

∏ 1 = (V k ) 11 and ∏ 2 = (V k ) 22 are

Connections between Quantum Walks and Relativistic Wave Equations

real. Moreover, the power U N k is simply given by U k = O(V k ) N O °1 or by:

U N k = O √ ∏ 1 N 0 0 ∏ 2 N ! O °1 (1.27)
To compute the solution, let us take into account as initial condition T0,k = P m ™ 0,m e i km , where ™ 0,m = ± m,0 |b L i. Computations of the exact expression for the Fourier-transformed amplitudes √L and √R are straightforwards. For instance, we report the solution for the amplitudes obeying to a Hadamard walk, that is for

B H = B °º/2,º/4,º/2,º/2 : √ L j ,m = √ 1 2º Z º °º °i e i k 2 p 1 + cos 2 k (e °i (! k j °km) )dk ! (1.28) √ R j ,m = µ 1 2º Z º °º 1 + cos(k) p 1 + cos 2 k (e °i (! k j °km) )dk ∂ (1.29)
where

! k =sin °1( sin(k) p 2 ) and ! k 2 [°º 2 , º 2 ]
. As we can see, the amplitudes for even m at odd j vanish and the same for the odd n at even j . Having an analytical expression of ™ j allows the study of the asymptotical behavior of the probability distribution |™ j ,m | 2 . Results of this study are extensively treated in [47], and confirm the previous qualitative description. We summarize here two important conclusions: (i) both amplitudes √ L and √ R are almost uniformly spread over the region [°j / p 2, j / p 2], in general over [°j K µ,ª,≥ , j K µ,ª,≥ ], and they decay rapidly to zero outside the region; (ii) the position probability distribution spreads as a function of j , in contrast to the standard deviation of an unrestricted CRW on a line, which is of order O( p j ).

Connections between Quantum Walks and Relativistic Wave Equations 1.2.1 Quantum walks and Feynman's Checkerboard

In Feynman's picture, a particle zig-zags at the speed-of-light across a spacetime lattice, flipping its chirality from left to right with an infinitesimal probability each time-step. The resulting dynamics, in the continuum limit, is the Dirac equation, with the flipping rate determined by the mass of the particle. (Strauch [43])

The Feynman's Checkerboard was first introduced in its original version by Feynman and Hibbs [10] as a simple path integral representation for the retarded Dirac propagator in 1+1 dimensions. In this model the particle's motion is restricted to be either forward (right) or backward (left) at the velocity of light. Feynman's path integral results from a particular finite differencing of the massive Dirac equation in (1+1):

( @ @ t ± @ @ x )√ L,R = i M√ R,L . (1.30)
where M is the real and constant mass of the fermion. On a square spacetime lattice (m, j ) 2 (Z£N) as in Fig. (1.5) with a mesh size ≤, the real variables (x, t ) are discretized so that t j = ≤ j and x m = ≤m. The finite difference equations, as prosed by

Chapter 1. HOMOGENEOUS QUANTUM WALKS

Feynman, read:

√ L j ,m = √ L j °1,m°1 + i M √ R j °1,m+1 (1.31) √ R j ,m = √ R j °1,m+1 + i M √ L j °1,m°1 (1.32)
where M = ≤M . These equations may be interpreted as following: the probability amplitude for a particle to be at ( j , m) moving toward the right, is equal to the amplitude encountered at ( j °1, m+1) moving toward the right plus i ≤M times the amplitude encountered at ( j °1, m°1) moving toward the left. Instead the probability amplitude for a particle to be at ( j , m) moving toward the left, is equal to the amplitude encountered at ( j °1, m °1) moving toward the left plus i ≤M times the amplitude encountered at ( j °1, m + 1) moving toward the right. Iterating the finite difference Eqs. (1.32), √ L,R may be expressed as a sum over paths leading to ( j , m).

The retarded propagator K j ,m is, therefore, the sum of number of paths, with N steps, leaving the point m 0 in direction ¬ 0 and end at point m" in direction ¬ and switching direction B times. Thus K j ,m can be expressed as:

K ¬¬ 0 = X B (© ¬¬ 0 ) j ,m;B (i ≤M ) B (1.33)
The convergence of Eq. (1.33) to the exact continuum propagator in the limit ≤ ! 0 is demonstrated by Jacobson and Schulman [22]. We can recognize that the probabilities (1.23) of a homogenous QW are related with what Feynman found for the Dirac particle of mass M . In fact, imposing µ= ≤M 2 , the probability distribution reads: .34) The relationship between the Feynman's Checkerboard and QWs is even more clear if we consider the discrete path integral approach introduced in modeling quantum walks, for instance, by Machida and Konno [31]. Thus, we dedicate the next section and the next chapter to understand the nature of this connection.

P QW ± º (1 ± ≤M ). ( 1

Homogeneous DTQWs and Weyl equation

Differently to what Feynman did, in this section we want to investigate the formal continuous limit of DTQWs with homogeneous and static quantum coin operator B AE,ª,≥,µ .

In order to investigate the continuous limit we first introduce a time step ¢t and a space step ¢x. We then introduce for any quantity a appearing in Eq. (1.15) and (1.16), a function ã defined on R + £ R such that the number a j ,m is the value taken by ã at the spacetime point (t j = j ¢t , x m = m¢x). We then suppose ™(t j , x m ) to be at least C 2 and that its characteristic length § to be much larger than the lattice parameter ¢x. Let us consider here the equations (1.15) and (1.16):

√ √ L (t j + ¢t , x m ) √ R (t j + ¢t , x m ) ! = B AE,ª,≥,µ √ √ L (t j , x m °¢x) √ R (t j , x m + ¢x) ! (1.35)
If the formal continuous limit exists, it will be obtained formally: (i) expanding each terms of the equations in ¢t and ¢x at fixed t j and x m and finally (ii) moving ¢t and ¢x to zero.

Let us now introduce a time-scale ø 2 R and a length-scale ∏ 2 R and an ≤ 2 R so that ø≤ ø 1 and ∏≤ ø 1. Define ¢t = ø≤ and ¢x = ∏≤ ± , where ± is a strictly positive real number, which traces the fact that ¢t and ¢x may tend to zero differently. The Taylor expansion of each spacetime dependent function in (1.35), up to order O(≤ 2 ), reads:

™ L,R (t j + ¢t j , x m ) = ™ L,R (t j , x m ) + ø≤@ t ™ L,R (t j , x m ) + O(≤ 2 ) (1.36) ™ L,R (t j , x m ± ¢x) = ™ L,R (t j , x m ) ± ∏≤ ± @ x ™ L,R (t j , x m ) + O(≤ 2± )
The Eq. (1.35) then leads us to:

√ √ L (t j , x m ) √ R (t j , x m ) ! + ø≤@ t √ √ L (t j , x m ) √ R (t j , x m ) ! + O(≤ 2 ) = B AE,ª,≥,µ (1 + ∏≤ ± @ x + O(≤ 2± )) √ √ L (t j , x m ) √ R (t j , x m ) ! (1.37)
Note that, in order to cancel the zeroth order contribution out, we have to impose that B AE,ª,≥,µ tends to the identity matrix I 2 of dimension 2 £ 2. Therefore the following equations must be satisfied:

e i (AE+ª) cos µ = 1 (1.38) e i (AE°ª) cos µ = 1 e i (≥+AE) sin µ = 0
The above relations imply µ = kº, AE = (k

+ k + + k °)º, ª =(k + °k°) º, (k, k °, k + ) 2 Z 3 .
The angle ≥ does not enter this constraint and is therefore an arbitrary real value. Finally, letting ≤ to zero, the Eqs. (1.37) admits formally two continuous limits: (i) one for ± > 1. In this case there is no propagation and the particle is localized at the initial condition; (ii) the most interesting case is for ± = 1 because all contributions are then of equal importance.

Let us now regard T = t /ø, X = x/∏, the equations of motion for the continuous limit of S 1 ≤ in this simple homogeneous case read:

( @ @ T °@ @ X )√ L = 0 (1.39) ( @ @ T + @ @ X )√ R = 0 (1.40)
Taken together, these two coupled first-order PDEs forms a system of relativistic wave equations describing massless spin-1/2 particles, otherwise known as Weyl equations. Let us remember that, in flat two-dimensional spacetime, the Clifford algebra can be represented by matrices acting on two-component spinors. This algebra admits two independent generators ∞ 0 and ∞ 1 , which can be represented by matrices obeying to the usual anti-commutation relation {∞ a , ∞ b } = 2¥ ab I, where ¥ is the Minkowski metric and I is the Identity matrix. Consider the (non unique) representation ∞ 0 = ae 1 and ∞ 1 = °ae1 ae 3 = i ae 2 , where ae 1 , ae 2 and ae 3 are the three Pauli matrices:

ae 1 = " 0 1 1 0 # , ae 2 = " 0 °i i 0 # , ae 3 = " 1 0 0 °1# . (1.41)
The Eqs. (1.39,1.40) can be recast in the following compact form:

i ∞ µ @ µ ™ = 0 (1.42)
where @ 0 = @ @ T , @ 1 = @ @ X .

Publication: "Massless Dirac Equation from Fibonacci Discrete-Time Quantum Walk"

This publication is a joint effort between the collaboration of the JSPS fellow Di Molfetta with the Shikano group, and the internship of two undergraduate Australian students, Lauchlan Honter and Ben Luo, from The University of Western Australia. We present two modified versions of the aperiodic quantum walk introduced by Ribeiro et al. [40], as a simple but non trivial case of the continuous limit with a homogeneous and periodic unitary step operator.

We introduce two models of DTQWs. The first, labeled FDTQW-I, is a generalized Hadamard coin B (µ j ) = B j following the infinite Fibonacci sequence, as below:

B j +1 = B j B j °1 (1.43)
where B 0 = B (AE) and B 1 = B (Ø)B (AE). We found that this sequence when applied to the coin operator is six-steps periodic. The formal continuous limit of the stroboscope of period six coincides with a ballistic equation that recover the Weyl equation in two spacetime dimension.

In the second case, FDTQW-II, the Fibonacci sequence is followed by the QW's step opera-

tor U = T B : U j +1 = U j U j °1 (1.44)
where U 0 = T B (AE) and U 1 = T B (Ø)T B (AE). In this case the sequence is no longer periodic but aperiodic as Ribeiro et al. [40] proved in his paper. In particular they showed that this model displays a sub ballistic behavior. What we prove is that if we truncate the Fibonacci sequence in FDTQW-II the behavior is no longer sub-balistic and display a standard propagative behavior.

There still remains a problem concerning the type of transition occurring when the periodicity FDTQW-II increases towards infinity, and the mechanism by which the dynamical behavior of the periodic extension changes into that of the infinite sequence [40].

Introduction

Discrete-time quantum walks (DTQWs) are defined as quantum-mechanical analogues of classical random walks. The concept of DTQWs was first considered by Feynman [1] and then introduced in greater generality by Refs. [2][3][4]. They have been realized experimentally in Refs. [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] and are important in many fields, from fundamental quantum physics [15,[21][22][23] to quantum algorithm [24][25][26] and condensed matter physics [27][28][29][30][31][32][33][34][35][36][37]. Previously, it has been shown that several DTQWs on a line admit a continuous limit identical to the propagation equations of a massive Dirac fermion [39][40][41][42][43] and those of massless Dirac fermion equations [21,39]. Furthermore, the relationship between DTQWs and artificial electric and gravitational fields has been shown [21,44]. Thus, DTQWs can be regarded as quantum dynamical simulators [39,40]. Additionally, it is well known that the classical random walk leads to a diffusive behavior characterized by the time evolution of the standard deviation, with σ (t) ∼ t 1/2 , while the standard DTQW leads to ballistic behavior, as σ (t) ∼ t. Further, the standard DTQW can be considerably enriched by generalizing the quantum coin operator and arranging it along different sequences. It has already been shown that quasi-periodic coin sequences induced by the Fibonacci sequence lead to sub-ballistic behavior, whereas random sequences lead to diffusive spreading [48][49][50]. Here, we consider two different Fibonacci DTQWs with periodic coin sequences. The first model (FDTQW-I) considers a time-dependent quantum coin following the Fibonacci sequence, while the second model (FDTQW-II) considers a modified version of the unitary operator first defined in Ref. [48], where the Fibonacci sequence is applied to the step operator. We show numerically and analytically that the continuous limit of these models reduces to a massless Dirac equation in (1 + 1) dimensions. by the chiral-dependent translation operator T , which is defined as:

✓ u m-1, j d m+1, j ◆ = T ✓ u m, j d m, j ◆ . (1)
Here, we introduce the simplest quantum coin, the generalized Hadamard coin, which is expressed as:

Ĉ(θ) = ✓ cos(θ) sin(θ) sin(θ) -cos(θ) ◆ , (2)
where θ ∈ [0, 2π ]. The one-step discrete time evolution is then given by:

✓ u m, j+1 d m, j+1 ◆ = T Ĉ(θ) ✓ u m, j d m, j ◆ . (3)
First, we consider FDTQW-I, which is the simplest case as only the QCO is defined as a Fibonacci series. Here,

Û j = T Ĉ j , j ∈ N 0 , Ĉ j+1 = Ĉ j Ĉ j-1 , (4)
with the initial conditions

Ĉ0 = Ĉ(α) , Ĉ1 = Ĉ(α) Ĉ(β), (5)
where j is the time step. On considering the DTQW acted upon by the Fibonacci coin series, we analytically find that the time evolution of the coin is cyclic with period 6. These coin operators then reduce to,

Ĉ0 = ✓ cos(α) sin(α) sin(α) -cos(α) ◆ , Ĉ1 = ✓ cos(α -β) -sin(α -β) sin(α -β) cos(α -β) ◆ , Ĉ2 = ✓ cos(2α -β) sin(2α -β) sin(2α -β) -cos(2α -β) ◆ , Ĉ3 = ✓ cos(α) sin(α) sin(α) -cos(α) ◆ , Ĉ4 = ✓ cos(α -β) sin(α -β) -sin(α -β) cos(α -β) ◆ , Ĉ5 = ✓ cos(β) sin(β) sin(β) -cos(β) ◆ . (6)
Here, the collection W n j = ( m,k ) m∈Z,k=nj is defined for j ∈ N and n ∈ {0, 1, 2, 3, 4, 5}, and represents the state of the walk at time k = nj. For any given n we define S n = (W n j ) j∈Z , where S n represents the entire history of the walk observed through a stroboscope of period n. Successive application of the 6 unitary operators to an initial state then gives the stroboscopic recursion equations for S 6 . The discrete-step equations for S 6 read

u m, j+6 = 3 X k=-3 A 2k (α, β)u m+2k, j + B 2k (α, β)d m+2k, j , (7)
d m, j+6 = 3 X k=-3 B -2k (-α, -β)u m+2k, j + A -2k (-α, -β)d m+2k, j . (8)
Here, the index k ∈ {-6, -4, -2, 0, 2, 4, 6} and the coefficients A k , B k ∈ R are explicitly given as

A -6 = c 2 α c β c 2 α-β c 2α-β (9)
A -4 = - 1 4 c α c 2 α-β (c α-2β + 3c 3α-2β -5c α + c 3α ) A -2 = 1 16 (-6c 2(α-β) + 4c 4(α-β) -c 2(α+β) -c 2(α-2β) + 2c 4α-2β -c 6α-2β + c 6α-4β -2c 4α -2c 2β + 6) A 0 = 1 4 c 2 α-β (-6c 2(α-β) + c 4α-2β -2c 2α + c 4α + c 2β + 5) A 2 = 1 8 c 2 α-β (6c 2(α-β) -3c 4α-2β -2c 2α + c 4α + c 2β -3) A 4 = 1 2 s 2α s β c 2 α-β c 2α-β A 6 = 0 B -6 = 1 2 s 2α c β c 2 α-β c 2α-β (10)
B -4 = 1 8 s 2α s β (s β -s 4α-3β ) + c β (3s 2α-β -s 4α-β + 3s 4α-3β -s 6α-3β ) B -2 = 1 8 ⇣ (c 2α -3)s 4α-4β -2s 4α c 2 α-β ⌘ 123 B 0 = 1 16 -s 2α-4β + 4s 4α-4β + s 6α-4β + 2c α (s α+2β - s α-2β -3s 3α-2β + s 5α-2β ) -2s 2α + 2s 4α B 2 = - 1 16 c α (s 3α-4β + 3s 5α-4β + 8s 3 α c 2α-2β + 4s α -2s 3α ) B 4 = -c 2 α s β c 2 α-β c 2α-β B 6 = 0
with c θ := cos(θ) and s θ := sin(θ).

Let us define the time and space variables, t j = j t and x m = m x, where t and x are the time and space steps, respectively. As t and x tend to zero, this allows us to take a Taylor expansion of the recursion relations for the DTQW and, hence, derive a pair of partial differential equations (PDEs). To take the continuous limit, we define

t = ϵ , x = ϵ γ , ( 11
)
where ϵ is an infinitesimal and γ > 0 is a scaling parameter. The difference between the two expressions is to account for the fact that t and x may tend to 0 differently. Then, taking a Taylor expansion about ϵ up to the leading orders of Eqs. ( 7) and ( 8), we obtain the following

u(x, t) + 6ϵ∂ t u(x, t) = u(x, t) + ϵ γ ( p 1 ∂ x u(x, t) + p 2 ∂ x d(x, t)) + O(ϵ 2 ), (12)
d(x, t) + 6ϵ∂ t d(x, t) = d(x, t) + ϵ γ ( p 2 ∂ x u(x, t) -p 1 ∂ x d(x, t)) + O(ϵ 2 ), (13)
where

p 1 = 3 X k=-3 2k 6 A 2k (α, β) = - 1 6 (c 4α-2β + 2c 2α + c 2β + 2), (14)
p 2 = 3 X k=-3 2k 6 B 2k (α, β) = - 2 3 (s 2α c 2 α-β ). (15)
Choosing scaling of γ = 1 and then taking the limit as ϵ → 0, we obtain the following pair of PDEs:

∂ t u(x, t) = p 1 ∂ x d(x, t) + p 2 ∂ x u(x, t), ∂ t d(x, t) = p 2 ∂ x d(x, t) -p 1 ∂ x u(x, t). ( 16
)
This set of equations can be then be recast, such that

I∂ t + P∂ x = 0 , P = ✓ p 1 p 2 p 2 -p 1 ◆ , ( 17
)
where

I
b u = 1 Z ✓ p 2 ω -p 1 b u + b d ◆ , ( 18
) b d = 1 Z ✓ - p 2 ω + p 1 b u + b d ◆ , (19)
with ω = q p 2 1 + p 2 2 and Z a normalized constant. Hence, Eq. ( 17) in the new basis reads

I∂ t + v(α, β)σ z ∂ x = 0, (20)
where

v 1 (α, β) = q 8c 2 α c 2α-2β + c 4α-4β + 4c 2α + 5 3 √ 2 , (21)
can be seen as the propagation velocity of the continuous limit distribution and σ z is the third Pauli matrix. We may also evaluate the standard deviation σ j of the probability distribution depicted in Fig. 1 as a function of time, by considering the exponent η(α, β) in σ j ∼ j η(α,β) . For the FDTQW-I case, we observe ballistic behavior for general (α, β), that is, η(α, β) = 1. Now consider Eq. ( 20) in the covariant form, which is expressed as:

i(γ 0 ∂ 0 + γ 1 ∂ 1 ) = 0, ( 22
)
where γ 0 = σ x and γ 1 = -iσ y are the usual gamma matrices, ∂ 0 = ∂ s , ∂ 1 = ∂ x , and the rescaled coordinate x = x/v 1 (α, β). This equation can now be interpreted as the massless Dirac equation in the 1+1 space-time dimension. In Fig. 2, we observe the density profile of the FDTQW-I case at a time step of j = 800, with a symmetric Gaussian initial condition that is sufficiently regular and large with respect to the lattice interval x. A truly ballistic propagation can be noted, as the continuous limit suggests. This result confirms that FDTQW-I can be used to simulate massless Dirac dynamics.

Discrete-time quantum walk with Fibonacci sequence quantum walk operator

The FDTQW-II previously defined in Refs. [48,51] can be expressed in the most general case as a unitary evolution, and the unitary evolution operator can then be defined using the following Fibonacci sequence with the initial conditions

Û j = Û j-1 Û j-2 , (23)
Û0 = T Ĉ(α), Û1 = T Ĉ(α) T Ĉ(β). (24)
Since each step operator contains increasing numbers of translation operators, the boundary size increases at an exponential rate. To account for this, we now parametrize time using the number of translation operators that have thus far been applied. It can then be noted that this DTQW has a quasi-periodic coin sequence. Thus, we define a new parameter for time r , such that s r = r t, where r = P j-1 n=0 F(n) and F(n) is the Fibonacci sequence with initial conditions, F(0) = 1 and F(1) = 2. This means how many translation operators a quantum walker is operated. This allows us to easily study the continuous limits of Eq. ( 3) for the operator expressed in Eq. (23). For the continuous limit to exist, we first require that (s r + τ t) → (s r ). It can be shown that this is true only when τ ∈ 6N 0 , that is, when τ is any positive integer multiple of 6, as the sequence of unitary operators then reduces to the identity operator. We assume the simplest periodic case, for which τ = 6, and the discrete-time equations read

u m,r +6 = 3 X k=-3 A 2k (α, β)u m+2k,r + B 2k (α, β)d m+2k,r , (25)
d m,r +6 = 3 X k=-3 B -2k (-α, -β)u m+2k,r + A -2k (-α, -β)d m+2k,r , (26)
where

A k ∈ R and B k ∈ R are A -6 = c 4 α c 2 β ( 27
)
A -4 = c 2 α s α (c 2 β s α + 2c α s 2β ) A -2 = - 1 8 s 2α -2s 2α + s 2(α-β) + 5s 2(α+β) A 0 = - 1 8 3 + c 4α -(1 + 3c 4α )c 2β -16c α s 3 α s 2β A 2 = c 2 β s 4 α -2c 3 α c β s α s β + c α s 3 α s 2β A 4 = c 2 α c 2 β s 2 α A 6 = 0 B -6 = c 3 α c 2 β s α (28) B -4 = c α c β ⇣ s 3 α c β + c α (1 -2c 2α )s β ⌘ B -2 = 1 8 s 4α (3c 2β -1) -4s 2 α (2c 2α + 1)s 2β B 0 = 1 8 4s 2 α (2c 2α + 1)s 2β + s 4α (1 -3c 2β ) B 2 = c α c β c α (2c 2α -1)s β -s 3 α c β B 4 = -s α c 3 α c 2 β B 6 = 0.
Taking the continuous limit of S 6 about ϵ (as in the previous section) and noting that the zeroth order terms cancel, we arrive at the first order term expressions for this system

∂ s u(x, s) = p ′ 1 ∂ x d(x, s) + p ′ 2 ∂ x u(x, s), ∂ s d(x, s) = p ′ 2 ∂ x d(x, s) -p ′ 1 ∂ x u(x, s), ( 29
)
and

p ′ 1 = 3 X k=-3 2k 6 A 2k (α, β) = - 1 3 (c 2 2α-β + 2c β c 2α-β ) (30) p ′ 2 = 3 X k=-3 2k 6 B 2k (α, β) = 1 3 s 2α-β (c 2α-β + 2c β ). (31)
As in the previous case, we can rewrite these equations in matrix form, such that,

∂ s + P ′ ∂ x = 0, P ′ = ✓ p ′ 1 p ′ 2 p ′ 2 -p ′ 1 ◆ , (32)
Again, we diagonalize the differential operator acting on , which spans itself on the new basis, (b

′ u , b ′ d ). The new basis components read b ′ u = 1 Z ′ ✓ p ′ 2 ω ′ -p ′ 1 b u + b d ◆ , (33) b
′ d = 1 Z ′ ✓ - p ′ 2 ω ′ + p ′ 1 b u + b d ◆ , (34)
with

ω ′ = q p ′ 1 2 + p ′ 2 2
and Z ′ a normalized constant. Eq. (32), expressed in terms of , then becomes

∂ s + v(α, β)σ z ∂ x = 0, v 2 (α, β) = 1 3 c 2α-β + 2c β . (35)
This PDE implies ballistic propagation of the quantum walker for fixed α and β. Similar results are obtained when other periodic sequences of finite length are considered, with each result yielding a similar differential equation, albeit with varying functions for v, as seen in Fig. 3. As in the previous section, the latter equation can be reformulated in the following covariant form: ,

i(γ 0 ∂ 0 + γ 1 ∂ 1 ) = 0, (36)
where γ 0 = σ x , γ 1 = -iσ y , ∂ 0 = ∂ s , ∂ 1 = ∂ x , while the rescaled coordinate is now x = x/v 2 (α, β).
In this rescaled time-space, the velocity is 1 and the FDTQW-II can again be interpreted as a massless Dirac equation in the (1 + 1) dimension.

Conclusion

By taking the continuous limit, we have characterized the propagation behavior of the DTQWs based on both the Fibonacci sequence of the coin operators and the periodic extension of these sequences. Our analysis shows that both of these quantum walks are ballistic and that the continuous limits can reduce to the (1 + 1)-dimensional massless Dirac equation. There still remain several unexplored topics, however. For example, we have not examined the continuous limit on the infinite sequence of step operators, as originally defined in Ref. [48]. The propagation properties of this model were numerically analyzed in [48,Fig.3]. According to Ref. [48], this walk is either ballistic or sub-ballistic depending on the initial values of α and β. Other unknown aspects include the type of transition that occurs when τ increases towards infinity, and the mechanism by which the behavior of the periodic extension changes into that of the infinite sequence.

Inhomogeneous Quantum Walks as synthetic gauge fields simulators 2.1.1 Quantum Simulation

Non-human animals can be reductively explained as automata. (Descartes, De homines, 1662) Let the computer itself be built of quantum mechanical elements which obey quantum mechanical laws" (Feynman, 1982) Since the first ideas on simulating quantum mechanics were formulated, it emerged that the simulation of the quantum system dynamics was a very challenging matter. This problem was unsolved until Feynman [20] [21] envisaged a solution: simulate quantum mechanics through an other quantum device. Unfortunately even though he foresaw this idea, he did not know how to realize it. Nowadays, quantum computers (i.e. an ensemble of well-defined qubits that can be initialized, measured and on which universal quantum gates can be performed Georgescu et al. [23]) are a reality and quantum computing and information are very active research fields. More important though is that we now acknowledge the possibility of implementing quantum simulation without an all-purpose, perfect quantum computer. Indeed, simpler quantum systems that are able to mimic the evolution of other quantum systems can be used for this task. In general, quantum simulations can be used to investigate physical problems that are either not solvable on classical computers or experimentally challenging and, in particular, they are useful to test new theoretical models. Applications can range from condensed matter physics [12,43,40,38] to high energy physics [24,48,11], from cosmology [39,26,2] to atom [25,54], to quantum chemistry [32,18] and, last but not least, to open quantum systems [34,53] and quantum chaos [START_REF] Weinstein | Experimental implementation of the quantum baker?s map[END_REF]27].

What is a synthetic gauge field?

Eichinvarianz (Hermann Weyl, 1929) To introduce the idea of gauge invariance and gauge field, let us introduce a closed system described by the action:

S = Z t 2 t 1 d t Z d 3 x L(√(x),@ µ √(x)), (2.1)
where L(√(x),@ µ √(x)) is the Lagrangian density and √(x) is a complex, differentiable field.

The equation of motion follows from Hamilton's principle of least action. In other words, the action should be stationary:

±S = ± Z t 2 t 1 d t Z d 3 xL(√(x), @ µ √(x)) = 0. (2.2)
As usual, the fields √ have to vanish at points t 1 and t 2 . As consequence of the least action principle, we derive the Euler-Lagrange equations, which allow us to write down the equations of motion for the fields √(x): @L @√(x) = @ µ @L @(@ µ √(x))

(2.3)

Now consider to multiply the √(x) by a global phase e i µ , where µ is a real constant. We can easily verify that such a transformation leave tshe Lagrangian density unaffected:

±L = 0. (2.4)
Evidently, if the Lagrangian density is unchanged under a global phase transformation, the equation of motions are invariant too. Let us now make a more systematic investigation of phase invariance. Consider an inhomogeneous phase e i µ(x) and the transformation:

√(x) ! √(x) 0 = e i µ(x) √(x). (2.5)
As consequence of this local phase rotations on the field, all terms depending on the derivative of the field transform as below:

@ µ √(x) 0 ! e i µ(x) (@ µ √(x) + i (@ µ µ(x))√(x)) (2.6)
It is clear that the Lagrangian density, and therefore the equations of motion, are no more invariant under local phase rotation. However, the invariance can be restored if we replace everywhere @ µ by

D µ = (@ µ + i A µ ), (2.7)
where A µ is introduced to take into account the non-homogeneity of the local phase rotation.

Let us call the new derivative (2.7) the gauge-covariant derivative and the field A µ , the electromagnetic gauge field. In fact, the Eq. (2.7) together with the Eq. (2.6) demands that A µ (x), in order to preserve global gauge invariance, transforms as:

A µ (x) ! A 0 µ (x) ¥ A µ (x) °@µ µ(x). (2.8) Therefore, D µ √(x) ! e i µ(x) D µ √(x).
(2.9)

For instance, if we explicitly introduce the Lagrangian density for a Dirac field:

L 0 = √(i ∞ µ @ µ °m)√, (2.10)
and we substitute @ µ by D µ , the gauge invariant Lagrangian density reads:

L = √(i ∞ µ D µ °m)√ = L 0 °Aµ √∞ µ √.
(

The last term is known as the electromagnetic current and come from the only requirement of the global gauge invariance. Historically, Weyl [START_REF] Weyl | Gravitation and electricity[END_REF] was the first to introduce the idea of gauge invariance as a geometrical principle to formalize the electromagnetism. Let us finally elucidate that a gauge field is usually called synthetic when it is defined as the result of the effective dynamics of the system. We will prove in the following sections that the effective dynamics of the IDTQWs coincides with the quantum propagation of a Dirac fermion in gauge fields, such as the electric or the gravitational field.

From Inhomogeneous Quantum Walks to synthetic gauge fields

The inhomogeneous QWs are usually characterized by a spatial dependence of the coin. Here we extend this definition to a dependence on each spacetime point of the lattice. Over the past decade many theoretical and experimental researches have been addressed to this particular family of walks.

Time dependent quantum coins have been introduced by Ribeiro et al. [46]. They defined a walk with several biased step operators applied aperiodically, and Banuls et al. [5] proposed a model with a time dependent coin as a control mechanism over a possible phase arising during the walk (as for example a consequence of an additional interaction). Albertini and D'Alessandro [1], [15] extended the analysis to a d-dimensional lattice and to the cycle. Recently Montero [41] unveiled some unexplored invariance in quantum walks with timedependent coin.

Space dependent quantum coins have been first introduced by Linden and Sharam [35] who investigated a spatial inhomogeneous quantum walk, where the coin operator depended periodically on the position; Shikano and Katsura [49] numerically and analytically showed the energy spectrum of such walks and Konno et al. [31], [30] have investigated the localization problems and the limit measures. Cedzich et al. [14] studied propagation and spectral properties of a QW in an electric field introduced as a space dependent phase in the quantum coin.

One of the main contribution of this thesis is to show that inhomogeneous QWs can be used to provide the gauge invariance with respect to certain gauge transformations, and therefore code for the corresponding gauge field, much in the same way as explained in the previous subsection for the continuum. In the continuum limit, the QW gauge field tend toward the Dirac equation gauge fields; as we will show for both the electric and the gravitational field in subsection 2.2.5 and as well as their attached papers.

Spacetime dependent quantum coin

The QWs defined here belong to a family larger than the ones introduced in Chapter (1), where the quantum coin is homogeneous in spacetime. In the inhomogeneous case, Eqs. (1.15) and (1.16) are transcribed in the following system of finite difference equations:

√ √ L j +1,m √ R j +1,m ! = B j ,m √ √ L j ,m°1 √ R j ,m+1 ! (2.12)
The quantum coin B j ,m belongs to U (2) and is given by: B (AE j ,m , µ j ,m , ª j ,m , ≥ j ,m ) = e i AE j ,m √ e i ª j ,m cos µ j ,m e i ≥ j ,m sin µ j ,m °e°i≥ j ,m sin µ j ,m e °i ª j ,m cos µ j ,m ! (2.13) As in equations (1.15) (1.16) and (1.12) (1.13) the probability ¶ j = P m |™ j ,m | 2 is finite and can be normalized to 1. It does not depend on j because of the unitarity of the dynamics. To investigate the continuous limit of a general stroboscopic equation S n , exactly as in sec.

(1.2.2), we introduce a time step ¢t and a space step ¢x, scaling with ≤:

¢x = ≤ ± ∏ (2.14) ¢t = ≤ø (2.15)
Now we should take into account the time and space dependency of the quantum coin and investigate formally the continuous limit of the walks, we should split each scalar field a in its homogeneous part a 0 and in its non-homogenous part a. In particular we consider that the spacetime dependent component of the field scales as O(≤) :

a = a 0 + ≤ ! a(t j , x m ) (2.16)
The exponent ! is a real and strictly positive number and is in general different from ±.

Note that in the Eq. (2.16), the spacetime dependency of the perturbation a(t j , x m ) permits to identify uniquely a family of QWs admitting a continuous limit, because the homogenous and the inhomogeneous part are in general always distinguishable. Furthermore, the inhomogeneities AE(t j , x m ), µ(t j , x m ), ª(t j , x m ) and ≥(t j , x m ), upcoming in subsection 2.2.5 and as well as their attached papers, are linked to the emergent electric and gravitational fields . The proof of the existence of these two emergent gauge fields will be the main goal of the following sections. Let us first and foremost introduce to the reader one of them: the synthetic gravitational field.

A synthetic gravitational gauge field 2.2.1 Simulating the effects of a gravitational gauge field

How to create synthetic gravity is a question often related to the reproduction of the gravitational force. This physical and engineering question can pertain both to zero-gravity environments and to the Earth. This task can be accomplished by several methods: for example by mechanical procedures [47], or by producing electric [28] and magnetic [29] forces. However, in this section, synthetic gravitational field has a quite different meaning and refers to mimicking the effect of gravity on the spacetime curvature of a specific particle motion. As we know, curved spacetime is not accessible in the labs, what we can look for is an analogue curved spacetime.

An interesting challenge is represented, for instance, by the simulation and study of QFTs in curved spacetime. Cold atoms have become a very useful tool for simulating gauge fields. Zohar et al. [START_REF] Zohar | Quantum simulations of gauge theories with ultracold atoms: Local gauge invariance from angular-momentum conservation[END_REF] and Boada et al. [7] were the first to prove that using ultra-cold atoms in optical lattices makes possible to mimic various types of statics and dynamical curved spacetime through simple manipulation of the optical setup.

Here we propose a new idea of simulators that reproduces the propagation of massless particle with two internal states on a curved (1 + 1) spacetime. In particular, the correction to the spinor propagation due to the gravitational deformation of the spacetime, is not real, but emerge as a geometric property of the QW. However, before we introduce the model, let us define several important elements of general relativity and curved field theory.

Curved spacetime and chiral field theory: an introduction

In the pre-relativistic context the framework of scientific researchers was the flat space of Euclidean geometry, grounded on the five Euclid's postulates. In 1827, in the Disquitiones generales circa superficies curvas, Gauss distinguished for the first time the inner or intrinsic properties of a surface from the outer or extrinsic ones. The former is measured by an observer living on the surface. The latter is derived from embedding the surface in a higher-dimensional space.

According to Gauss, the most fundamental inner property of a surface was the shortest path between two points. He realized, for instance, that a distance measured over a conic surface was the same as that one measured over a cylinder. These two surfaces, in fact, share the same intrinsic properties and in particular, their metric is flat. We cannot say the same for a sphere, as a sphere cannot be mapped onto a plane without distortions. However, since a cylinder or a cone are round in radial direction, we might be led to think they are curved surfaces. This is due to the fact that we consider them as 2-dimensional surfaces in a 3-dimensional space, and we intuitively compare the curvature of the lines, which are on the surfaces, with straight lines in the flat 3-dimensional space.

The interplay between the inner and the outer properties of surfaces is one of the key point of the Einstein's theory of general relativity. To be more clear we need first to introduce the formalism and the main assumptions of the theory. Then we will be prepared to define properly the Dirac equation in curved spacetime.

A formal general setup

"Riemannian geometry is characterized by the facts that the infinitesimal neighborhood of every point P has a Euclidian metric, and that the magnitudes of two line elements that belong to the infinitesimal neighborhoods of two finitely distant points P and Q are comparable. However, the concept of parallelism of these two line elements is missing; for finite regions the concept of direction does not exist. The theory put forward in the following is characterized by the introduction, in addition to the Riemann metric, of a direction, or of equality of direction, or of parallelism for finite distances. Correspondingly, new invariants and tensors will appear in addition to those of Riemann geometry. ( Einstein [19] ,1928) With no extra constructs, and by simply using the only intuitive idea of parallelism, Einstein introduced, as reported in the above cited manuscript, for the first time the concept of connection and curvature in field theory. This later entailed the formulation of a new theory called the vierbein field theory that has permitted to derive a gauge field theory approach to general relativity [START_REF] Yepez | Einstein's vierbein field theory of curved space[END_REF]. In particular, we will employ this approach to derive Dirac equation in curved spacetime.

To that end, let us first introduce the standard mathematical framework to treat the curved spacetime and the sufficient ingredients, with a particular attention to the notion of spin connection and Lorentz invariance. All we introduce in the following introductory section can be find in any standard textbook on curved QFT and we address the reader in particular for detailed computations to Boulanger et al. [8], Brill and Wheeler [9] or Birrell and Davies [6] for a general discussion on chiral matter in curved spacetime or for a synthetic but rigorous approach Yepez [START_REF] Yepez | Einstein's vierbein field theory of curved space[END_REF]. Moreover, en elegant differential approach is proposed by Lasenby et al. [33].

Mathematical Framework

Let us consider a Lorentzian manifold (M , g ), where M is a pseudo-Riemannian manifold and g is a Lorentzian metric, that is g relates at each point p 2 M a Lorentzian scalar product g p on the tangent space T p . The differential basis is given by the partial derivatives of the coordinates at p: ê(µ) = @ (µ) , (2.17) where the hat denote the unit vector, related to the basis vector ê. The subscript (µ) denotes a component of the basis vector, that serves as the coordinate system representation. All vectors A 2 T p can be spanned on this basis and read:

A = A µ ê(µ) . (2.18)
The cotangent basis is defined by the dual space of the T p , and the differential basis is:

ê(µ) = dx (µ) . (2.19)
Let us call the dual space T § p . All vectors in this space have components:

A = A µ ê(µ) = g µ∫ A ∫ ê(µ) . (2.20)
Now let us bring in a non-coordinate unit vector basis, an orthonormal basis that is independent from the coordinates. In a (1+3)-dimensional spacetime the basis is composed by four linear independent elements and let us call it tetrads (group of four from Greek). Alternatively in (1+1)-spacetime we call in the following this orthonormal basis dyads, and we denote each component by ê(a) . ê(a)

• ê(b) = ¥ ab , (2.21)
where ¥ ab is the Minkowski metric of flat spacetime. The non-coordinate unit vector bases spans all coordinate basis with a given local curvature (depending on x):

ê(µ) = e a µ (x)ê (a) . (2.22)
The matrix e a µ (x) is invertible and named vierbein field. These fields has to be at least C 1 and are defined such that:

g µ∫ (x)e µ a (x)e ∫ b (x) = ¥ ab . ( 2
(a) ! ê(a 0 ) = § a a 0 (x)ê (a) , (2.25)
where § a a 0 (x) is the position dependent transformation, so that:

§ a a 0 (x) § b b 0 (x)¥ ab = ¥ a 0 b 0 . (2.26)

Affine connection and Spin connection

To introduce the notion of connection, let us describe the Fig. 2.1. Consider two spacetime points x AE ¥ A and x AE + ±x AE ¥ B and denote by X AE (x) a vector at point A. If we perform a translation of ±X AE (x) in spacetime, the transformed vector, at point B , (X AE (x)) 0 , is given by X AE (x) + ±X AE (x). In general, the pseudo-Riemann manifold possesses a non-vanishing curvature and, thus, we should consider a distortion effect in respect to the vector X AE + ±X AE as in a flat spacetime. From the Fig. Now let us include the following finite difference:

r ∞ X AE (x) ¥ 1 ±x ∞ (X AE (x + ±x) °[X AE (x) + ±X AE (x)]) (2.27)
and comment that, if the curvature is zero, or ±X AE (x) = ±X AE (x), the above expression coincides, in the continuous limit, with the standard partial derivative at point x.

Otherwise, it is possible to demonstrate [START_REF] Yepez | Einstein's vierbein field theory of curved space[END_REF] that the above definition transcribe in the following covariant derivative:

r ∞ X AE = @ ∞ X AE + °AE Ø∞ X Ø (2.28)
In other terms, in coordinate-based differential geometry, the covariant derivative of a tensor is given by its partial derivative plus a correction terms, one for each index, involving an affine connection contracted with the tensor.

Exactly for the same principle, in the non-coordinate based differential geometry, the spin connection takes the place of the affine connection. For each index we get a correction factor, conventionally positive for upper index, negative for lower index:

r ∞ X a b = @ µ X a b + ! a µ c X c b °!c µ b X a c . (2.29)
Let us acknowledge that, using the vierbein matrices, we can write the affine connections in terms of spin connections and vice versa:

°∫ µ∏ = e ∫ a @ µ e a ∏ + e ∫ a e b ∏ ! a µ b ! a µ b = e a ∫ e ∏ b °∫ µ∏ °e∏ b @ µ e a ∏ .
(2.30)

Another essential ingredient is the tetrad postulate. It states that the covariant derivative of the vierbein field exists and vanishes, or in other words, the vierbein field is invariant under parallel transport:

r µ e a ∫ = 0. (2.31)
An important consequence of this postulate, called metric compatibility, is that the metric tensor g ∫∏ itself is also invariant under parallel transport:

r µ g ∫∏ = 0. (2.32)
Let us now introduce the last ingredient, that is the relation between the affine connection and the metric tensor. Permute the indices in the latter equation as follows:

r Ω g µ∫ = @ Ω g µ∫ °°∏ Ωµ g ∏∫ °°∏ Ω∫ g µ∏ = 0, (2.33)
r µ g µΩ = @ µ g ∫Ω °°∏ µ∫ g ∏Ω °°∏ µΩ g ∫∏ = 0, (2.34) r ∫ g Ωµ = @ ∫ g Ωµ °°∏ ∫Ω g ∏µ °°∏ ∫µ g Ω∏ = 0, (2.35)
Take (2.33)-(2.34)-(2.35), multiply by g aeΩ and contract the indices ae and µ, the affine connections read:

°µ µ∫ = 1 2 g µΩ @ ∫ g Ωµ . (2.36)
The latter equation can be transcribed in:

°µ µ∫ = 1 2 Tr (g ∏Ω @ ∫ g Ωµ ) = 1 p °g @ ∫ p °g (2.37)
where g ¥ °Det(g Ωµ ). Observe that from the above relation, given the Eqs. (2.30), we can derive simply the relation between the metric tensor and the spin connection.

Lorentz invariance in curved spacetime

To formalize a relativistic curved field theory demands the Dirac equation to be Lorentz invariant. In flat space, the Lorentz transformation acts in general on the physical space and on the spinor states. Notably the Dirac gamma matrices transform under an internal Lorentz unitary transformation U ( §) as follows:

U ( §) °1∞ µ U ( §) = § µ ae ∞ ae . (2.38)
and therefore the Dirac equation changes as follows:

[i ∞ µ @ µ °m]√(x) ! U ( §)[i ∞ µ @ µ °m]√( § °1x).
(2.39)

We redirect the reader to Peskin and Schroeder [45] for a detailed demonstration of the Lorentz invariance of Dirac equation in flat space. In curved spacetime we should take into account the correction due to the curvature of the space. Let us call it °µ. Following [START_REF] Yepez | Einstein's vierbein field theory of curved space[END_REF] Gamma matrices transform under Lorentz transformations as follows:

§ µ ae ∞ ae ( § °1) ∫ µ = ∞ ∫ (2.40)
and the Dirac equation in curved space:

∞ µ (@ µ + °µ)√(x) ! § µ ∏ e ∏ a ∞ a ( § °1) ∫ µ (@ µ + °0 µ )U ( §)√(U ( §) °1x) (2.41)
where ∞ µ = e µ a ∞ a . A tedious computation of the RHS of the above transformation, leads us to the condition, to be respected in order to preserve the Lorentz invariance:

°0 ∫ = U ( §)°∫U ( §) °1 °@∫ (U ( §))U ( §) °1
(2.42)

Chapter 2. INHOMOGENEOUS QUANTUM WALKS AND CONTINUOUS LIMITS

Finally the Dirac equation in curved spacetime reads:

i ∞ a e µ a (x)D µ √ °m√ = 0 (2.43)
where the covariant derivative

D µ = @ µ + °µ (2.44)
where °µ transforms according to Eq. (2.42). The last step is to determine explicitly the covariant derivative and notably the form of the correction °µ. To do this, we should demand the spinor field to be invariant under Lorentz transformations. The correction to the spinor field is found to be:

°AE = 1 2 e Ø k (@ AE e Øh )S hk . (2.45)
Yepez [START_REF] Yepez | Einstein's vierbein field theory of curved space[END_REF].

Let us remark that all we presented in this section is casted purely in terms of QW by Arrighi et al. [3].

Dirac equation in curved spacetime in (1+1) dimensions

We have now all the elements to derive the Dirac equation in a curved spacetime of arbitrary dimension. The (1 + 1) dimensional case is straightforward 1 . The problem in fact simplify greatly. The product ∞ a e µ a (x)D µ in Eq. (2.43) transcribes in:

1 2 ∞ a 1 p °g @ µ ( p °g e µ a ) (2.46)
thus the Dirac equation can be written as follows:

[i ∞ a e µ a @ µ + i 2 ∞ a 1 p °g @ µ ( p °g e µ a ) °m]√ = 0. (2.47)
Observe that we have derived the Dirac equation from considerations of invariance with respect to Lorentz transformations, in analogous way the electromagnetism is derived from a gauge invariance principle.

Publication: "Quantum walks as massless Dirac fermions in curved spacetime"

The intent of this publication is to show how a particular family of IDTQWs can mimic a massless fermion propagating in a curved metric. In order to show a specific case, a spherically symmetric solution to Einstein equation in vacuum such as the Schwarzschild metric, which represents a Schwarzschild black hole, has been chosen. This kind of black hole is no charged and its angular momentum associated to its mass is zero. Observe, finally that, even though we model the motion of a particle in a curved (1 + 1)°spacetime, we do not mimic two dimensional gravity. We notably study the case of a particle confined on the radial dimension of a (3 + 1)-spherical metric. Because the particle is massless, it follows the geodesics for all times and consequently it can be described by a (1 + 1)-curved Dirac equation. However, gravity in (1 + 1) has been investigated for instance by Brown et al. [10] and Morsink and Mann [42], [36]. Despite its simplicity, this theory also has a number of striking classical and semi-classical features, including a well-defined Newtonian limit ( [37]), FRW cosmologies, gravitational collapse, cosmological solutions and gravitational waves (Sikkema and Mann [50]).

Di Molfetta et al. [17], [16] investigated systematically the continuous limit of the QW's stroboscopic equation S 2 and they formally recovered the (1 + 1)-Dirac equation in curved spacetime. The curvature of spacetime in continuous limit is linked (through the metrics g ) to the quantum coin of the walk, B j m , in discrete spacetime.

Albeit recents, these results have immediately inspired others to investigate the matter further. For instance Arrighi et al. [4] recovered the massive Dirac equation in curved (1 + 1)spacetime, performing the continuous limit by a new method. They introduced in the standard definition of the QW time evolution operator, three extra unitary operations: encoding, grouping and decoding. These operations allowed the authors to build a map between usual QWs and usual QWs of increased dimension, called paired QWs. The continuous limit is consequently taken on the latter. Furthermore, the continuous limit of Paired QW also offers a new potential perspective: it can be seen as an emergent property of spacetime grouping. For instance, seeing QWs as a specific case of QCA , we can change the local interactions betweens the cells and the spatial distributions of the cells, i.e., the QCA interaction graph, to obtain different emergent dynamics in the grouping operations. This can lead us to connect the geometrical aspects of the interaction graph with the geometrical properties of gauge QFT, in a similar fashion to the gravitational field connected to the curvature of the spacetime. In this way, through the geometrical aspects of its interaction graph as such, gauge fields are naturally built into the very fabric of the associated QCA, and by changing the QCA structure we change the quantum field theory that has been described. This can a be a new potential possibility to explore in the future. Let finally remark that, in consequence to we said in this paragraph, the model proposed by Arrighi et al. can recover our model, by choosing a specific spacetime grouping or a QCA interaction graph.

Another approach was recently introduced by S.Succi et al. [52], who explored the connection between Dirac equation, QWs and Quantum lattice Boltzmann (QLB). However, their method is even more different from what has been described before. They showed that the discretization of the Dirac equation in flat and curved spacetime, through lattice Boltzmann numerical schemes, can be recast in the QW formalism.

The first quantum walk (QW) was built by Feynman [1] as a possible discretization of the standard, massive Dirac dynamics in flat spacetime. General discrete-time QWs have then been introduced in the physics literature by Refs. [2,3] and the continuous-time version first appeared in Ref. [4]. QWs are the simplest formal analogs of classical random walks and are important in many fields, ranging from fundamental quantum physics [5,6] to quantum algorithmics [7,8], solid-state physics [9][10][11][12], and biophysics [13,14].

QWs have been realized experimentally; for example, as transport of trapped ions [15,16], of photons in wave guide lattices [6] or optical networks [17] and of atoms in optical lattices [18]. QW experiments of two photons [19] have recently been performed, with the possibility of simulating Bose or Fermi statistics [20], and cavity QED QWs have also been proposed [21]. Following Feynman's idea, several authors have studied the continuous limit of general QWs. Most publications [22][23][24][25][26][27][28][29] only envisage QWs with constant coefficients. The continuous limit of QWs with time-and space-dependent coefficients has been considered only recently, in Refs. [30][31][32]. These references present several families of QWs, in both (1 + 1) or (1 + 2) spacetime dimensions, whose continuous limit is described by a flat spacetime Dirac equation with a generalized mass term and electromagnetic coupling. The electromagnetic field is generated by the spacetime dependence of the angles defining the walks and thus vanishes if these angles are constant.

The other gauge field which couples naturally to a Dirac spinor is evidently gravity. Yet, QWs whose continuous limit are described by Dirac equations in curved spacetime have remained elusive. This article considers a particular family of discrete-time QWs with nonconstant angles in (1 + 1) spacetime dimensions and associates to this family a continuous limit which is described by a massless Dirac equation in curved spacetime. This family was not investigated in Refs. [30][31][32] and the limit procedure used in these references cannot be applied to the family discussed here. The result presented in this article opens the way to possible laboratory experiments simulating coherent quantum propagation in relativistic gravitational fields. It also establishes a connection between general relativity and all the aforementioned fields where QWs are useful.

We consider QWs defined over discrete time and discrete one-dimensional (1D) space, driven by time-and spacedependent quantum coins acting on a two-dimensional (2D) Hilbert space H. The walks are defined by the following finite difference equations, valid for all (j,m) ∈ N × Z:

ψ L j +1,m ψ R j +1,m = B θ j,m ψ L j,m+1 ψ R j,m-1 , (1)
where

B(θ ) = -cos θ isin θ -i sin θ + cos θ . (2)
The index j labels instants and the index m labels spatial points. For each instant j and each spatial point m, the wave

function jm = ψ L jm b L + ψ R jm b R has two components ψ L jm
and ψ R jm on the spin basis (b L ,b R ) and these code for the probability amplitudes of the particle jumping towards the left or towards the right. Note that the spin basis is interpreted as being independent of j and m. The total probability

π j = m (|ψ L j,m | 2 + |ψ R j,m | 2
) is independent of j ; i.e., conserved by the walk. The set of angles {θ j,m ,(j,m) ∈ N × Z} defines the walks and is arbitrary.

Consider now, for all (n,j ) ∈ N 2 , the collection W n j = ( k,m ) k=nj,m∈Z . This collection represents the state of the QW at "time" k = nj . For any given n, the collection S n = (W n j ) j ∈N thus represents the entire history of the QW observed through a stroboscope of "period" n. The evolution equations for S n are those linking W n j +1 to W n j for all j . These can be deduced from the original evolution equations (1) of the walk, which also coincide with the evolution equations of S 1 . In particular, the evolution equations of S 2 read

ψ L j +2,m = c j +1,m c j,m+1 ψ L j,m+2 -is j,m+1 ψ R j,m + s j +1,m s j,m-1 ψ L j,m + ic j,m-1 ψ R j,m-2 , (3) ψ R j +2,m = s j +1,m ic j,m+1 ψ L j,m+2 + s j,m+1 ψ R j,m -c j +1,m is j,m-1 ψ L j,m -c j,m-1 ψ R j,m-2 , (4)
where c jm = cos(θ jm ) and s jm = sin(θ jm ).

To investigate the continuous limit of S n , we first introduce a time step t and a space step x. We then consider that jm and θ jm are the values taken by a two-component wave function and by a function θ at the spacetime point (t j = j t,x m = m x). We finally suppose that and θ are at least twice differentiable with respect to both space and time variables for all sufficiently small values of t and x. The formal continuous limit of S n is defined as the couple of differential equations obtained from the discretetime evolution equations defining S n by letting both t and

x tend to zero. Let us therefore introduce a time-scale T , a length-scale L, an infinitesimal ϵ and write t = ϵT and x = ϵL. The continuous limit of S n can then be investigated by Taylor expanding in powers of ϵ the discrete equations defining S n . For the limit to exist, all zeroth-order terms must identically cancel each other and the differential equation describing the limit is then obtained by equating to zero the nonidentically vanishing, lowest-order contribution. Consider first S 1 , which is identical to the original walk. It is rather obvious that zeroth-order terms cancel each other only if the operator B defining the walk tends to unity as ϵ tends to 0 (see Ref. [30] for a detailed discussion of this point). The operator B defining the family of walks considered in this article does not depend on ϵ and is different from unity for all values of θ. Thus, S 1 does not admit a continuous limit for the family of walks defined by Eq. ( 1). But S 2 , on the other hand, does. Indeed, a straightforward computation delivers the following equation obeyed by the wave function :

T + (cos θ)P X = Q , (5)
where the operators P and Q are represented, in the base (b L ,b R ), by the matrices

P = -cos θ isin θ -i sin θ cos θ (6)
and

Q = -θ X sin 2θ 2 i 2 [θ T -θ X (cos 2θ)] i 2 [θ T + θ X (cos 2θ)] θ X sin 2θ 2 . (7)
In Eqs. ( 5) and ( 7), the subscript T (X) indicates a derivative with respect to the dimensionless variable T = t/T (X = x/L).

The operator P is self-adjoint and its eigenvalues are -1 and +1. Two eigenvectors associated with these eigenvalues are

b -= i cos θ 2 b L -sin θ 2 b R , (8)
b + = i sin θ 2 b L + cos θ 2 b R . (9)
The family (b -,b + ) forms an orthonormal basis of the twodimensional-spin Hilbert space. Let us now rewrite Eq. ( 5) in components, but in this new orthonormal basis. A tedious but straightforward computation leads to

ψ - T -(cos θ)ψ - X + θ X 2 (sin θ )ψ -= 0, (10)
ψ + T + (cos θ)ψ + X - θ X 2 (sin θ )ψ + = 0,
where ψ -and ψ + are the components of in the new basis. This form of the equations makes it easy to check that the continuous dynamics conserves the total probability π

(T ) = dX| (T ,X)| 2 = dX[|ψ -(T ,X)| 2 + |ψ + (T ,X)| 2 ]
, as it should.

Suppose now, to make the discussion definite, that cos θ is strictly positive and introduce in spacetime {(T ,X)} the Lorentzian, possibly curved metric G defined by its covariant components (G µν ) = diag(1, -1/ cos 2 θ ), where (µ,ν) ∈ {T ,X} 2 . This metric defines the canonical, scalar "volume" element D G X = √ -GdX = dX/cos θ in physical 1D X space, where G is the determinant of the metric components G µν . Dirac spinors are normalized to unity with respect to D G X, whereas is normalized to unity with respect to dX. We thus introduce = √ cos θ and rewrite the equations of motion (10) in terms of . We obtain

γ a e µ a ∂ µ + 1 2 1 √ -G ∂ µ √ -Ge µ a = 0, (11)
where µ ∈ {T ,X}, a ∈ {0,1}. The usual 2D gamma matrices are

γ 0 = 0 1 1 0 , γ 1 = 0 1 -1 0 , (12)
and the e µ a are the components of the diad (orthonormal basis) e 0 = e T and e 1 = (cos θ )e X on the original coordinate basis (e T ,e X ). Equation ( 11) is the standard [33] equation of motion for a massless Dirac spinor in (1

+ 1)-dimensional spacetime with metric G. The spin basis is (b -,b + ).
This result shows that QWs can be used to model quantum transport in any 2D gravitational field. Indeed, any 2D Lorentzian metric can be put under the above diagonal form by a suitable choice of coordinates. The single angle θ(t,x) is thus enough to describe any 2D gravitational field. Let us stress, however, that gravity is very different in 2D and in 4D since, in particular, all 2D spacetimes are conformally flat [34]. But Eq. ( 11) can also be used to model quantum transport in higher-dimensional spacetimes by QWs on the line. As an example, we now construct a QW on the line which mimics the radial motion of a Dirac spinor in a spherically symmetric 4D black hole.

A Schwarzschild black hole is a spherically symmetric solution of Einstein equation in vacuo. The corresponding 4D metric reads, in dimensionless Lemaître coordinates (τ,ρ,θ,φ) [35],

ds 2 = dτ 2 - r g r dρ 2 -r 2 d , (13)
where r(τ,ρ) = r

1/3 g [ 3 2 (ρ -τ )] 2/3 , d = dθ 2 + (sin 2 θ)dφ 2 .
The event horizon is located at r = r g , i.e., ρ = τ + (2/3)r g , and the singularity is located at r = 0, i.e., ρ = τ . The exterior of the black hole is the domain r > r g .The range of variations for the Lemaître coordinates is τ > 0, ρ > τ [i.e., r(τ,ρ) > 0], 0 6 θ 6 π, 0 6 φ < 2π.

Because of the spherical symmetry, a point mass which starts its motion radially will go on moving radially. Radial motion can be studied by introducing the 2D metric g, also singular at r = 0, with covariant components g τ τ = 1, g ρρ = -r g /r, g τρ = g ρτ = 0. The null geodesics of g are defined by dτ = ±[r g /r(τ,ρ)] 1/2 dρ. Note that the 2D projection of the horizon on the (τ,ρ) plane coincides with a null geodesics of g. We now identify the dimensionless time T with the time coordinate τ and the dimensionless space variable X with λρ, where λ is an arbitrary strictly positive real number (see Fig. 1). The "radius" r can then be expressed as a function of T and X:

r(T ,X) = 3 2 X λ -T 2/3 r 1/3 g , (14)
and the components of g in the coordinate basis associated with T and X are g T T = 1, g XX = -r g /(λ 2 r), g T X = g XT = 0. Note that the condition ρ > τ transcribes into X > λT .

Let D be the domain where -g XX > 1. This domain is characterized, in (T ,X) coordinates, by the condition X 6 λT + 2 3λ 2 r g . In D the metric g can be identified with the metric G. This identification defines a angle θ which depends on T and X by

(cos θ)(T ,X) = λ r(T ,X) r g (15)
and, by extension, a QW in D.

The condition defining D can be rewritten as r 6 r g /λ 2 . The domain D thus includes, for all λ, the singularity located at r = 0. For λ > 1, r g /λ 2 < r g and D is then entirely located inside the horizon. For λ = 1, D coincides with the interior of the horizon, and D extends outside the horizon for λ < 1.

The spatial density | (T ,X)| 2 of the walk defined by Eqs. ( 1), (2), and ( 15) is plotted in Fig. 1. All graphs have been obtained with r g = 150 and ϵ = 0.5. The initial condition is (0,X) = √ n 0 (X)(b L + ib R ) with an initial Gaussian density n 0 of variance X 0 = 2.5 centered on X 0 = 50.5. In Figs. 1(a) and 1(b) (in which the right limit of D coincides with the horizon), both branches of the QW starts their evolution inside the horizon and end up on the singularity. In Fig. 1(c) the QW starts exactly on the horizon; the right branch follows it while the left branch ends up on the singularity. In Fig. 1(d) the QW starts outside the horizon; its right branch propagates away from the horizon while the left branch still ends up on the singularity.

Figure 1 clearly shows that the QW closely follows null geodesics of the metric, i.e., it behaves as a massless fermion in the gravitational field of the black hole. The agreement between the geodesics and the density profile of the walk is all the more remarkable given that all graphs correspond to ϵ = 0.5, which lies well outside of the continuous limit ϵ ≪ 1 envisaged above. Note, however, that the right branch of the QW lags slightly behind the null geodesic when approaching the r = 0 singularity [see top of Fig. 1(b)]. The numerical results thus suggest that the main result of this article can somehow be extended beyond the continuous limit. This interesting point will be investigated in future presentations.

In summary, we have considered a particular family of one-dimensional QWs which does not admit a continuous limit by the procedure used in Refs. [30][31][32]. The procedure used here keeps one time step out of two, and gives a continuous limit to walks from this family. We have computed this limit and proven that it coincides with the propagation of a massless Dirac fermion in an arbitrary gravitational field. Note, however, that Fig. 1, where all time steps are retained, shows that the density clearly follows, at all times, the null geodesic predicted by the continuous limit. We have also constructed explicitly a QW which mimics the propagation of a fermion outside and inside the event horizon of a Schwarzschild black hole and illustrated this construction by numerical simulations.

Let us now discuss the above results. Keeping only one time step of the QW out of two to build a continuous limit might appear unnatural and might even look like a purely mathematical trick. Let us explain now in detail why this is not so. Consider first, as an instructive example, the sequence of numbers u j defined by u 0 = 1 and

u j +1 = σ exp(iωT )u j , (16)
for j ∈ N; here, ω is an arbitrary real number, T is a time-scale and σ = ±1 and does not depend on j . A direct computation shows that u j = σ j exp(iωt j ) with t j = j T . Suppose σ = +1. The sequence u j is then a simple circular function of the time t j . On the contrary, if σ = -1, the sequence u j is then a circular function of the time t j combined with an extra phase shift of π at every time step. The sequence u j admits a continuous limit if σ = +1. But, it does not if σ = -1, because of this extra phase shift. In particular, if σ = -1, the oscillating behavior of the sequence u j cannot be recovered by simply taking the continuous limit of the evolution equation (16). The best way to recover this oscillating behavior is then to consider the new sequence v k built out of u j by keeping only one time step out of two. Indeed, this new sequence obeys the discrete evolution equation v k+1 = exp(2iωT )v k . This equation admits a continuous limit, described by the ordinary differential equation dv dt = 2iωv, which clearly reveals the oscillating behavior in v and, thus, in u. Naturally, all information on the extra phase shift of π is lost in this procedure. Note, however, that this extra phase shift does not influence the sequence of the squared moduli |u j | 2 .

Let us now compare the preceding example with the QWs examined in Refs. [30][31][32] and in this article. The spinor plays a role similar to u and the evolution equation of the walk [Eq. ( 1) for the QWs considered in this article] has the same status as Eq. ( 16). The QWs studied in Refs. [30][31][32] are equivalent to the sequence u j obtained with σ = +1. They thus admit a standard continuous limit, which is fully described in these earlier publications. This limit coincides with the propagation in flat spacetime of a Dirac fermion possibly coupled to an electric field. On the other hand, the QWs defined by Eqs. ( 1) and ( 2) correspond to the sequence u j with σ = -1 and they do not admit a continuous limit. But the derived walks built by keeping only one time step out of two of the original QWs do admit a continuous limit because the squared matrix B 2 = 1, just as (-1) 2 = 1. Keeping one time step out of two is thus not a contrived, unphysical procedure. On the contrary, it is dictated by the very definition of the QWs we study in this article and it is the only one which delivers a continuous limit for these walks. The obtained continuous limit coincides with the physically interesting situation of a massless Dirac fermion propagating in curved spacetime. Since the retained limit procedure is itself dictated by the QWs studied in this article, the geometry appearing in their continuous limit is an intrinsic property of the walks themselves.

The explicit construction of a QW mimicking propagation in and around a black hole shows that it is possible, at least in principle, to simulate by laboratory experiments the propagation of quantum systems in interesting relativistic gravitational fields.

The work presented in this article should naturally be extended in several directions. One should first investigate systematically all QWs on the line defined with a quantum coin acting on a 2D Hilbert space and, for each family of walks, try and determine if and how a continuous limit can be obtained. One should also extend the main result of this article to QWs defined on physical space of higher dimension and/or defined by quantum coins acting on a higher-dimensional Hilbert space. A particular goal would be to obtain walks whose limits are described by a Dirac fermion coupled to both a gravitational and an electromagnetic field.

As noted earlier, the geometry appearing in the continuous limit of the QWs studied in this article is an intrinsic property of these walks and is not to be confused with the geometry the "real", "physical" spacetime which might be used to realize the walk experimentally. Another extension of this work would therefore be to consider walks defined on curved physical spaces (graphs) and to investigate how the geometry of the underlying space (graph) couples to the intrinsic geometry of the walk. This is not a purely academic problem, since QWs can model photon transport in networks of algae, which may have a nontrivial geometry.

The main result of this article also suggests that concepts from general relativity and differential geometry may play a key role in understanding the behavior of QWs on graphs and their use in quantum algorithmics. This role should also be investigated thoroughly. The following publication represents an extension of the previous one. In the previous section we have seen that QWs can model fermion transport on a curved spacetime and that the gravitational field was characterized, in the continuous limit, by the metric G(t , x) = °1, °cos °2 µ(t , x) ¢

A synthetic electric gauge

where µ(t , x) is a spacetime dependent parameter of the quantum coin defined in (2.13) and the other parameters are constant. Furthermore, in the following publication, we explore all the parameters' space of the quantum coin. The walk is therefore characterized by three time-and space-dependent Euler angles (µ, ª, ≥) and by a global, also time-and space-dependent phase AE.

Exploring systematically the QW's stroboscopic equation S 1 and S 2 , we prove that this class of QWs exhibits an exact discrete U (1) gauge invariance. Let us consider the local phase transformation

™ j m ! ™ j m e i ¡ j m , (2.48)
where {¡ j m ,(j , m) 2 N £ Z}.

The finite difference equations of the QWs are left invariant if we replace the parameter of the coin {AE j ,m , ª j ,m , ≥ j ,m , µ j ,m } by:

AE j ,m ! AE j ,m + ae j ,m 2 ª j ,m ! ª j ,m + ± j ,m ≥ j ,m ! ≥ j ,m °±j,m
(2.49)

µ j ,m ! µ j ,m
where

ae j ,m = ¡ j ,m+1 + ¡ j ,m°1 °2¡ j +1,m (2.50)
± j ,m = ¡ j ,m+1 °¡j,m°1 2 . (2.51)
Observe that this discrete gauge invariance transcribes, in the continuous limit, into the standard continuous U (1) gauge invariance of electromagnetism presented in sec. (2.1.2). The electric field generated by the gauge potential in the continuous limit:

A(t , x) = (V (t , x), A(t , x)) = °Ǣ(t , x), °ā(t , x) ¢ , ( 2.52)
is, therefore, closely related to the U (1) discrete gauge invariance of the QW.

Introduction

QWs are simple formal analogues of classical random walks. They were first considered by Feynman [1] as possible discretizations of the free Dirac dynamics in flat space-time. They have been introduced in the physics literature by Refs. [2,3] and the continuous-time version first appeared in Ref. [4]. They have been realized experimentally in Refs. [5][6][7][8][9][10][11][12] and are important in many fields, ranging from fundamental quantum physics [12,13] to quantum algorithmics [14,15], solid state physics [16][17][18][19] and biophysics [20,21]. Following Feynman's idea, several authors have studied the continuous limit of various QWs. The first publications [22-24,1,25-28] only addressed QWs with constant coefficients and recent work has extended the discussion to QWs with time-and space-dependent coefficients [29][30][31][32], in both (1+1) and (1+2) space-time dimensions. In particular, a new method was developed in Refs. [29][30][31][32] to investigate the continuous limit of QWs with nonconstant coefficients. This method delivers interesting results, not only for standard QWs, but also for 'derived' QWs obtained from original QWs by keeping only one time-step out of two [32]. So far, this new method has only been applied to particular families of walks. This article presents the systematic application of this method to all QWs in (1 + 1) space-time dimensions. The main conclusions are: (i) all families of walks do not admit a continuous limit (ii) when the limit exists, it coincides, in all cases but one, with the dynamics of a Dirac fermion coupled to an artificial electric field and/or relativistic gravitational field. These theoretical conclusions are illustrated by numerical simulations. Connections with previous results as well as other topics like transport in graphene are discussed and possible experimental realizations are also mentioned.

⇤

Fundamentals

We consider quantum walks defined over discrete time and discrete one-dimensional space, driven by time-and spacedependent quantum coins acting on a two-dimensional Hilbert space H . The walks are defined by the following finite difference equations, valid for all (j, m) 2 N ⇥ Z:

" L j+1,m R j+1,m # = B ✓ j,m , ⇠ j,m , ⇣ j,m , ↵ j,m " L j,m+1 R j,m 1 # , (1)
where

B(✓ , ⇠ , ⇣ , ↵) = e i↵
 e i⇠ cos ✓ e i⇣ sin ✓ e i⇣ sin ✓ e i⇠ cos ✓ .

(

) 2
This operator is in U( 2), and is in SU( 2) only for ↵ = p⇡ , p 2 Z, and ✓, ⇠ and ⇣ are then called the three Euler angles of B. The index j labels instants and the index m labels spatial points. The two functions L and R can be interpreted as the components of a wave function on a certain orthonormal basis (b L , b R ) independent of j and m. These two components code for the probability amplitudes of the particle jumping towards the left or towards the right. The total

probability ⇡ j = P m | L j,m | 2 + | R j,m | 2 is independent of j,
that is, it is conserved by the walk. The set of angles ✓ j,m , ⇠ j,m , ⇣ j,m , ↵ j,m , (j, m) 2 N ⇥ Z defines the walks and is at this stage arbitrary.

Consider now, for all (n, j) 2 N 2 , the collection W n j = ( k,m ) k=nj, m2Z . This collection represents the state of the quantum walk at 'time' k = nj. For any given n, the collection S n = (W n j ) j2Z thus represents the entire history quantum walk observed through a stroboscope of 'period' n. The evolution equations for S n are those linking W n j+1 to W n j for all j. These can be deduced from the original evolution Eqs. ( 1) and (2) of the walk, which also coincide with the evolution equations of S 1 . For example, the evolution equations of S 2 read

L j+2,m = A L j,m L j,m+2 + B L j,m L j,m + C L j,m R j,m + D L j,m R j,m 2 R j+2,m = A R j,m L j,m+2 + B R j,m L j,m + C R j,m R j,m + D R j,m R j,m 2 , (3)
where

A L j,m = c j+1,m c j,m+1 e i(↵ j+1,m +⇠ j+1,m +↵ j,m+1 +⇠ j,m+1 ) B L j,m = s j+1,m s j,m 1 e i(↵ j+1,m +⇣ j+1,m +↵ j,m 1 ⇣ j,m 1 ) C L j,m = c j+1,m s j,m+1 e i(↵ j+1,m +⇠ j+1,m +↵ j,m+1 +⇣ j,m+1 ) D L j,m = s j+1,m c j,m 1 e i(↵ j+1,m +⇣ j+1,m +↵ j,m 1 ⇠ j,m 1 ) , (4)
A R j,m = s j+1,m c j,m+1 e i(↵ j+1,m ⇣ j+1,m +↵ j,m+1 +⇠ j,m+1 ) B R j,m = s j,m 1 c j+1,m e i(↵ j+1,m ⇠ j+1,m +↵ j,m 1 ⇣ j,m 1 ) C R j,m = s j+1,m s j,m+1 e i(↵ j+1,m ⇣ j+1,m +↵ j,m+1 +⇣ j,m+1 ) D R j,m = c j,m 1 c j+1,m e i(↵ j+1,m ⇠ j+1,m +↵ j,m 1 ⇠ j,m 1 ) , (5)
with c jm = cos(✓ j,m ) and s jm = sin(✓ j,m ).

The QWs defined by (1) admit a remarkable exact U(1) gauge invariance. Consider indeed an arbitrary set of numbers { jm , (j, m) 2 N ⇥ Z}, and write jm = 0 jm exp(i jm ). A straightforward computation shows that 0 obeys

" 0 L j+1,m 0 R j+1,m # = B ✓ 0 j,m , ⇠ 0 j,m , ⇣ 0 j,m , ↵ 0 j,m " 0 L j,m+1 0 R j,m 1 # , (6)
with

↵ 0 j,m = ↵ j,m + j,m 2 ⇠ 0 j,m = ⇠ j,m + j,m ⇣ 0 j,m = ⇣ j,m j,m (7)
✓ 0 j,m = ✓ j,m and j,m = j,m+1 + j,m 1 2 j+1,m (8) j,m = j,m+1 j,m 1 2 . ( 9
)
It can also be shown that the S 2 -type QWs admit the same discrete invariance. As detailed in Sections 4 and 5, this discrete gauge invariance transcribes in the continuous limit into the standard continuous U(1) gauge invariance of Maxwell electromagnetism.

To investigate the continuous limit of a collection S n , we first introduce a time step 1t and a space step 1x. We then consider that jm , ✓ jm , ⇠ jm , ⇣ jm and ↵ jm are the values taken by a two-component wave function and by four functions ✓ , ⇠ , ⇣ and ↵ at the space-time point (t j = j1t, x m = m1x). Thus, Eq. ( 1) transcribes into

 L (t j + 1t, x m ) R (t j + 1t, x m ) = B ✓(t j , x m ), ⇠ (t j , x m ), ⇣ (t j , x m ), ↵(t j , x m )  L (t j , x m + 1x) R (t j , x m 1x) . ( 10
)
We finally suppose that and ✓ are at least twice differentiable with respect to both space and time variables for all sufficiently small values of 1t and 1x. The formal continuous limit of S n is defined as the couple of partial differential equations (PDEs) obtained from the discrete-time evolution equations defining S n by letting both 1t and 1x tend to zero.

How to determine the continuous limit

Let us introduce a time scale T , a length-scale L, an infinitesimal ✏ and write

1t = ✏T 1x = ✏ L, (11)
where > 0 allows 1t and 1x to tend to zero differently. We also allow the angles defining the walk to depend on ✏ and characterize the ✏-dependence of these angles near ✏ = 0 by the following scaling laws:

✓ ✏ (t, x) = ✓ 0 (t, x) + ✓(t, x)✏ ! ⇠ ✏ (t, x) = ⇠ 0 (t, x) + ⇠(t, x)✏ ⇣ ✏ (t, x) = ⇣ 0 (t, x) + ⇣ (t, x)✏ (12) ↵ ✏ (t, x) = ↵ 0 (t, x) + ↵(t, x)✏ ⌘ ,
where the four exponents !, , and ⌘ are all positive. We also suppose that all functions are at least C 2 in t and x. The above relations define 1-jets of quantum walks. We finally denote by B ✏ (t, x) the matrix B(✓ ✏ (t, x), ⇠ ✏ (t, x), ⇣ ✏ (t, x), ↵ ✏ (t, x)).

Expand now the original discrete equations obeyed by a jet S n

✏ around ✏ = 0. A necessary and sufficient condition for the expansion to be self-consistent at order 0 in ✏ is that B n 0 (t, x) = 1 for all t and x (note from Eq. ( 1) that this condition is self-evident for n = 1). This transcribes into a constraint for the zeroth-order angles ✓ 0 , ⇠ 0 , ⇣ 0 , ↵ 0 . Suppose this constraint is satisfied. The differential equations defining the continuous limit are obtained from the expansion by stating that the next lowest order contribution in ✏ identically vanishes. If one excepts zeroth-order terms, the terms of lowest order in the expansion scale as ✏, ✏ , ✏ ! , ✏ , ✏ , ✏ ⌘ (see for example the similar expansions performed on particular simple quantum walks and presented in Refs. [30,31]). The richest and most interesting scaling is thus

= ! = = = ⌘ = 1,
because this makes all the above terms of the same order and, thus, delivers a differential equation with a maximum number of contributions. This scaling will be retained in the remainder of this article. Note that Eqs. ( 11) and ( 12) have actually very different meanings. Indeed, (11) states that the relative variations of between j + 1, m and j, m ± 1 should be small, while (12) states that the angles defining the walk do not deviate much from their zeroth-order values.

We will now present in detail the continuous limit of the jets S n ✏ for both n = 1 and n = 2.

Limit of S 1 ✏

Zeroth-order values of the angles

The constraint on the zeroth-order angles reads

sin ✓ 0 = 0 e i(↵ 0 +⇠ 0 ) cos ✓ 0 = 1 e i(↵ 0 ⇠ 0 ) cos ✓ 0 = 1. ( 13
)
The above relations imply ✓ 0 = k⇡ , ↵ 0 = (k

+ k + + k )⇡, ⇠ 0 = (k + k )⇡, (k, k + , k ) 2 Z 3 .
The angle ⇣ 0 does not enter this constraint and is therefore an arbitrary function of t and x. For a given value of ✏, there is thus no meaningful distinction between ⇣ 0 and ⇣ . We will therefore from here on denote ⇣ 0 by ⇣ in all equations, if only to simplify the notation.

Equations of motion

Let now T = t/T , X = x/L, x ± = (T ± X )/2 and @ ± = @ x ± = @ T ± @ X . The variables x ± are null coordinates in the flat 2D space-time. With these notations, the equations of motion for the continuous limit of S 1 ✏ read

@ L i( ↵ + ⇠ ) L = + ✓ e i(✓ 0 +↵ 0 +⇣ ) R @ + R i( ↵ ⇠ ) R = ✓ e i(✓ 0 +↵ 0 ⇣ ) L , (14)
where ✓ 0 and ↵ 0 are arbitrary multiples of ⇡ (see the constraint above) and ⇣ is an arbitrary real function of T and X . Taken together, these two coupled first-order PDEs form a Dirac equation in (1 + 1) dimensions. Let us recall that, in flat two-dimensional space-times, the Clifford algebra can be represented by 2 ⇥ 2 matrices acting on two-component spinors. This algebra admits two independent generators 0 and 1 , which can be represented by 2 ⇥ 2 matrices obeying the usual anti-commutation relation

{ a , b } = 2⌘ ab I, ( 15
)
where ⌘ is the Minkowski metric and I is the identity (unit) matrix. Consider the representation 0 = 1 and 1 = 1 3 = i 2 , where 1 , 2 and 3 are the three Pauli matrices:

1 =  0 1 1 0 , 2 =  0 i i 0 , 3 =  1 0 0 1 . ( 16
)
Eq. ( 14) can be recast in the following compact form:

(i 0 D 0 + i 1 D 1 M) = 0, (17)
where

D µ = @ µ iA µ , @ 0 = @ T , @ 1 = @ X , A 0 = ↵, A 1 = ⇠ , M = diag (m , m + ) and m ⌥ = ±i ✓ exp {i(✓ 0 + ↵ 0 ± ⇣ )}.
This equation describes the propagation in flat space-time of a Dirac spinor coupled to the Maxwell potential A (the corresponding electric field is E X = @ T ⇠ @ X ↵). The discrete gauge invariance presented in Section 2 degenerates accordingly into the standard local U(1) invariance associated with electromagnetism. Indeed, suppose that the numbers jm (see Section 2) are the values taken by a function at space-time points (t j = j1t, x m = m1x). Expanding Eqs. ( 7) and ( 9) at first order in ✏ delivers

↵ 0 = ↵ ✏ @ @T ⇠ 0 = ⇠ + ✏ @ @X ⇣ 0 = ⇣ ✏ @ @X (18)
✓ 0 = ✓.

The first two equations imply

A 0 0 = 1 ✏ ↵ 0 ↵ 0 = 1 ✏ ✓ ↵ ✏ @ @T ↵ 0 ◆ = A 0 @ @T A 0 1 = 1 ✏ ⇠ 0 ⇠ 0 = 1 ✏ ✓ ⇠ 0 ✏ @ @X ⇠ ◆ = A X @ @X , (19)
which are simply the standard gauge transformation for the potential A. The fourth relation implies that the mass tensor M is gauge invariant. Since the continuous limit equation of motion ( 14) depends only on ⇣ 0 (as opposed to ⇣ ), the third equation is not relevant to the continuous limit investigated in this section.

The angles ✓ 0 and ↵ 0 are both multiples of ⇡. Both masses are therefore complex conjugates of each other. They are real, and therefore identical, if ⇣ 0 is an uneven multiple of ⇡/2. They are both real and positive, equal to

| ✓ |, if ⇣ 0 = ✓ 0 + ↵ 0 + ¯ + (2p + 1)⇡ /2
, where exp(i ¯ ) = sgn ✓ is the sign of ✓ and p is an arbitrary integer. Note that, even in this case, the mass may depend on both T and X .

Limit of S 2 ✏

Zeroth-order values of the angles

The constraint on the zeroth-order angles now reads cos ⇠ 0 sin(2✓ 0 ) = 0 e 2i↵ 0 e 2i⇠ 0 cos 2 ✓ 0 sin 2 ✓ 0 = 1 e 2i↵ 0 e 2i⇠ 0 cos 2 ✓ 0 sin 2 ✓ 0 = 1. (20) As for n = 1, ⇣ 0 does enter this constraint; it is therefore an arbitrary function of t and x, which we denote simply by ⇣ (see the discussion at the end of Section 4.1).

The first relation implies that cos ⇠ 0 = 0 (case 1) or sin 2✓ 0 = 0 (case 2). The first case corresponds to ⇠ 0 = (2k + 1)⇡ /2, k 2 Z. The second and third relations then transcribe into the single constraint ↵ 0 = (2k 0 + 1)⇡ /2, with k 0 2 Z. Note that ✓ 0 can then be an arbitrary function of t and x, as ⇣ 0 . This function will be simply denoted by ✓, just as ⇣ denotes ⇣ 0 .

In contrast, the second case corresponds to ✓ 0 = k⇡ /2, k 2 Z. If k = 2p + 1, p 2 Z (case 2.1), the last two constraint relations deliver ↵ 0 = (2k 0 + 1)⇡ /2, with k 0 2 Z. The angle ⇠ 0 is then arbitrary and will be denoted simply by ⇠ . If k = 2p (case 2.2), the last two constraint relations deliver

↵ 0 = k 0 ⇡/2, ⇠ 0 = ↵ 0 + k 00 ⇡, (k 0 , k 00 ) 2 Z 2 .
Cases 1, 2.1 and 2.2 partly overlap. Indeed, jets obeying ✓ 0 = k⇡ , ⇠ 0 = (2k + 1)⇡ /2 and ↵ 0 = (2k 0 + 1)⇡ /2 can be field under both case 1 and case 2. These are the only jets which can be field under both cases.

Let us now give the equations of motion of the continuous limit in cases 1, 2.1 and 2.2.

Equations of motion

: case 1 2 @ T (cos 2 ✓)@ X L 2i ↵ + (cos 2 ✓) ⇠ L e +i(⇣ ⇠ 0 ) (sin 2✓ )@ X R = (sin 2✓ )(@ X ✓) + i(sin 2 ✓)(@ + ⇣ ) L + e +i(⇠ 0 +⇣ ) ✓ i(@ ⇣ ) sin 2✓ 2 i ⇠(sin 2✓ ) + (@ T ✓) (@ X ✓)(cos 2✓ ) ◆ R ( 21
)
and

2 @ T + (cos 2 ✓)@ X R 2i ↵ (cos 2 ✓) ⇠ R e i(⇣ ⇠ 0 ) (sin 2✓ )@ X L = +(sin 2✓ )(@ X ✓) i(sin 2 ✓)(@ ⇣ ) R + e i(⇠ 0 +⇣ ) ✓ i(@ + ⇣ ) sin 2✓ 2 i ⇠(sin 2✓ ) (@ T ✓) (@ X ✓)(cos 2✓ ) ◆ L . (22)
These equations can be put into the more compact form

@ T + (cos ✓)P@ X = Q , (23)
where

= L b L + R b R , P = ✓ cos ✓ e +i(⇣ ⇠ 0 ) sin ✓ e i(⇣ ⇠ 0 ) sin ✓ cos ✓ ◆ ( 24
)
and

Q = ✓ Q L L Q L R Q R L Q R R ◆ (25)
with

Q L L = i ↵ + (cos 2 ✓) ⇠ sin 2✓ 2 (@ X ✓) + i 2 (sin 2 ✓)(@ + ⇣ ) Q R R = i ↵ (cos 2 ✓) ⇠ + sin 2✓ 2 (@ X ✓) i 2 (sin 2 ✓)(@ ⇣ ) (26)
Q L R = e +i(⇠ 0 +⇣ ) 2 ✓ i(@ ⇣ ) sin 2✓ 2 i ⇠(sin 2✓ ) + (@ T ✓) (@ X ✓)(cos 2✓ ) ◆ Q R L = e i(⇠ 0 +⇣ ) 2 ✓ i(@ + ⇣ ) sin 2✓ 2 i ⇠(sin 2✓ ) (@ T ✓) (@ X ✓)(cos 2✓ ) ◆ .
The operator P is self-adjoint and its eigenvalues are 1 and +1. Two eigenvectors associated with these eigenvalues

are b = ✓ cos ✓ 2 ◆ b L + e i( ⇣ +⇠ 0 ) ✓ sin ✓ 2 ◆ b R , ( 27
) b + = ✓ sin ✓ 2 ◆ e i(⇣ +⇠ 0 ) b L + ✓ cos ✓ 2 ◆ b R . (28)
The family (b , b + ) forms an orthonormal basis of the two-dimensional spin Hilbert space, alternate to the original basis

(b L , b R ). Let = b + + b + . Eq. (23) transcribes into @ T (cos ✓)@ X i ↵ i cos ✓ ⇠ + i 2 ((cos ✓ 1)@ + ⇣ ) + @ X ✓ 2 (sin ✓) = 0 @ T + + (cos ✓)@ X + i ↵ + + i cos ✓ ⇠ + i 2 ((cos ✓ 1)@ ⇣ ) + @ X ✓ 2 (sin ✓) + = 0. ( 29
)
Suppose now, to make the discussion definite, that cos ✓ is strictly positive and introduce in space-time {(T , X )} the Lorentzian, possibly curved metric G defined by its covariant components

(G µ⌫ ) = 1 0 0 1 cos 2 ✓ ! , (30)
where (µ, ⌫) 2 {T , X } 2 . This metric defines the canonical, scalar 'volume' element

D G X = p G dX = dX / cos ✓ in physical 1D X -space,
where G is the determinant of the metric components G µ⌫ . Dirac spinors are normalized to unity with respect to D G X , whereas is normalized to unity with respect to dX . We thus introduce = p cos ✓ and rewrite the equations of motion ( 29) in terms of . We obtain

a  e µ a D µ + 1 2 1 p G @ µ ⇣ p Ge µ a ⌘ = 0, (31)
where µ 2 {T , X } , a 2 {0, 1} and D µ = @ µ iA µ with

A T = ↵ + 1 cos ✓ 2 @ X ⇣ (32)
and

A X = ⇠ 1 cos ✓ 2 cos ✓ @ T ⇣ . (33)
The usual 2D gamma matrices are

0 = ✓ 0 1 1 0 ◆ , 1 = ✓ 0 1 1 0 ◆ , (34)
and the e µ a are the components of the dyad (orthonormal basis) e 0 = e T and e 1 = cos ✓ e X on the original coordinate basis (e T , e X ). Eq. ( 31) is the standard [33] equation of motion for a massless Dirac spinor propagating in (1 + 1)-dimensional space-time under the combined influence of the gravitational field G and the electric field E deriving from A. Since the Dirac field is massless, its components are not coupled and evolve independently of each other. Each component follows a null geodesic of the gravitational field, and the electric field only modifies the energy along a given geodesic. Numerical simulations of a QW propagating radially in the gravitational field of a Schwarzschild black hole are presented in Section 6.3.

Let us conclude this section by commenting briefly on how the discrete gauge invariance presented in Section 2 transcribes in the present context. The continuous limit equation ( 18) is of course valid. Combining these with (32), (33) and keeping only the lowest order terms in ✏ leads to the standard gauge transformation A 0 0 = A 0 @ /@T and A 0 1 = A 1 @ /@X. Just as it was the case in Section 4, the transformation law for ⇣ does not contribute to the continuous gauge transformation, but it is not for the same reason. In Section 4, the potential A itself does not depend on ⇣ . Here, the potential A does depend on ⇣ , but the gauge transformation for ⇣ generates terms of order ✏ in the gauge transformation for A, and these terms vanish as ✏ tends to zero. In the present context, the final, fourth equation in (18) reflects the fact that the gravitational field does not depend on the choice of gauge for the phase of the spinor .

Equations of motion: case 2.1

The equations of motion of the continuous limit read

2@ T L 2i ↵ L = +i(@ + ⇣ ) L + 2 ✓e +i⇣ (cos ⇠) R 2@ T R 2i ↵ R = i(@ ⇣ ) R 2 ✓e i⇣ (cos ⇠) L , (35)
where ⇠ and ⇣ are arbitrary functions of T and X . These equations are not PDEs, but ordinary differential equations (ODEs) in ± . Thus, there is for example no propagation in this case. Technically, this comes from the fact that ✓ 0 is here constrained to be an uneven multiple of ⇡/2.

Equations of motion: case 2.2

The equations of motion read

2@ L 2i ↵ + ⇠ L = +2 ✓ e i(2↵ 0 +⇣ ) (cos ⇠ 0 ) R 2@ + R 2i ↵ ⇠ R = 2 ✓ e i(2↵ 0 ⇣ ) (cos ⇠ 0 ) L , (36)
where ↵ 0 is a multiple of ⇡/2, ⇠ 0 ↵ 0 is multiple of ⇡ and ⇣ is an arbitrary function of T and X . Eq. ( 36) can be recast in the following compact form:

(i 0 D 0 + i 1 D 1 M) = 0 (37)
where

D µ = @ µ iA µ , @ 0 = @ T , @ 1 = @ X , A 0 = ↵, A 1 = ⇠ , M = diag (m , m + ) and m ⌥ = ±i ✓ exp {i(2↵ 0 ± ⇣ )} (cos ⇠ 0 ).
This equation describes the propagation in flat space-time of a Dirac spinor coupled to the potential A and with a mass tensor M.

Numerical simulations

Basics

In order to ascertain the validity of the continuous limits that were derived above, we wish to compare numerical solutions of the QW defined by the finite difference Eqs. ( 1) and ( 2) with the corresponding Dirac-type PDEs defined by Eqs. ( 31) and (17).

While the numerical integration of the QW finite difference equations poses no particular problem, controlling the error on numerical solutions of PDEs is a more involved matter. This hurdle can be avoided in the special case where the mass terms cancels, because one can then compare the numerical solutions of the QW finite difference equations with the numerical solutions of the ODEs defining the characteristics of the massless Dirac PDE (see Section 6.3).

We have chosen to use Fourier pseudo-spectral methods [34], for their precision and ease of implementation. We therefore restrict ourselves to 2⇡ -periodic boundary conditions. A generic field (x) is thus evaluated on the n collocation points x j = 2⇡ j/n, with j = 0, n 1 as j = (x j ). The discrete Fourier transforms are standardly defined as

(x j ) = n/2 1 X k= n/2 exp(ikx j ) ˆ k ˆ k = 1 n n 1 X j=0 (x j ) exp( ikx j ). (38)
These sums can be evaluated in only n log(n) operations by using Fast Fourier Transforms (FFTs). Spatial derivatives of fields are evaluated in spectral space by multiplying by i k and products are evaluated in physical space.

The original QW Eqs. ( 1) and (2) can also be simply cast in this setting, as the translation operator j ! j±1 is represented in Fourier space by ˆ k ! ˆ k exp(±i k2⇡ /n). In this setting, the continuous limit is automatically taken when n is increased.

QWs in constant a uniform electric field

As explained in Section 2, the QWs and the Dirac equation exhibit a U(1) gauge invariance. All choices of gauge thus correspond to the same physics. Within a pseudo-spectral code, the right gauge to work with a 1D constant uniform electric field E is A 0 = 0, A 1 = E T ; in particular, the other 'natural choice' A 0 = E X, A 1 breaks the spatial periodicity condition.

In all QW simulations, the retained choice of gauge has been implemented by choosing the following numerical values:

↵(T , X ) = 0 ⇠(T , X ) = 1.1 T ⇣ (T , X ) = ⇡ 2 ✓(T , X ) = 0.24.
We used initial data consisting of Gaussian wave packets of positive energy solutions to the free Dirac equation. The Gaussian widths X are such that they are well resolved within the used resolutions so that spectral convergence is ensured.

As discussed above (see end of Section 6.1) the QW and its Dirac continuous limit can be jointly simulated within the same pseudo-spectral algorithm. This allows for a very simple, direct evaluation of the discrepancy between the QW and the corresponding solution of the Dirac equation. This discrepancy can be measured by the relative difference N rel between the density of the QW and the density of the solution of the Dirac equation. Fig. 1 shows that such a typical relative difference scales as ✏, as expected. Indeed, for a single time-step, the discrepancy is theoretically of order ✏ 2 . Thus, after a fixed time T = O(✏ 1 ), the discrepancy is of order ✏ 1 ✏ 2 = ✏.

This result confirms that QWs can be used to simulate Dirac dynamics in constant electric field, as was done for example in Refs. [35,36]. Both QW and Dirac dynamics are very rich, as exemplified by Fig. 2, which compare with Fig. 2 of Ref. [35] and Fig. 3 of Ref. [36].

Note that, as X increases, the spatial dispersion of the wave packet also increases which makes the time evolution of the density more complex. The solution which is initially a positive energy planar wave starts to oscillate between positive and negative modes under the action of the constant electric field displaying high-frequency Zitterbewegung in Fig. 2(c)-(d).

Offering new results of the Dirac dynamics in the presence of an electric field is not the purpose of this article. Let us conclude this section by offering instead a brief historical overview of the very large literature already existing on the topic.

In 1929, Klein studied a relativistic scalar particle moving in an external step function potential. He found a paradox that, in the case of a strong potential, the reflected flux is larger than the incident flux although the total flux is conserved [37]. Sauter studied this problem for a Dirac spin 1/2 particle by considering a potential corresponding to a electric field with constant value E 0 on a given interval. He found an expression for the transmission coefficient of the wave through the electric potential barrier from the negative energy state to positive energy states [38]. This remarkable phenomenon was related, in 1936, to positron-electron pair creations by Heisenberg and his student Hans Euler [39].

Of course, in order to deal with anti-particles a massive reinterpretation of the Dirac equation theory is necessary [40], leading to modern field theory and quantum electrodynamics. The modern formula for pair creation in a constant external electric field was delivered, in 1951, by Schwinger [41]. It involves the same dominant exponential term exp( ⇡ m 2 e c 3 heE ) that was derived, 20 years before, by Sauter. A detailed review of these historical developments is given in the first sections of Ref. [42].

QWs in Schwarzschild black hole

A Schwarzschild black hole is a spherically symmetric solution of the Einstein equation in vacuo. The corresponding 4D metric reads, in dimensionless Lemaître coordinates (⌧ , ⇢, ✓ , ) [43]:

ds 2 = d⌧ 2 r g r d⇢ 2 r 2 d⌦, (39)
where r(⌧ , ⇢) = r

1/3 g ⇥ 3 2 (⇢ ⌧ ) ⇤ 2/3 , d⌦ = d✓ 2 +(sin 2 ✓)d 2 .
The event horizon is located at r = r g , that is, ⇢ = ⌧ +(2/3)r g , and the singularity is located at r = 0 i.e. ⇢ = ⌧ . The exterior of the black hole is the domain r > r g . The range of variations for the Lemaître coordinates is ⌧

0, ⇢ ⌧ , that is, r(⌧ , ⇢) 0), 0  ✓  ⇡, 0  < 2⇡ .
Because of the spherical symmetry, a point mass which starts its motion radially will go on moving radially. Radial motion can be studied by introducing the 2D metric g, also singular at r = 0, with covariant components g ⌧ ⌧ = 1, g ⇢⇢ = r g /r, g ⌧ ⇢ = g ⇢⌧ = 0. The null geodesics of g are defined by d⌧ = ± r g /r(⌧ , ⇢)

1/2 d⇢. Note that the 2D projection of the horizon on the (⌧ , ⇢)-plane coincides with a null geodesics of g. We now identify the dimensionless time T with the time coordinate ⌧ and the dimensionless space variable X with ⇢,

where is an arbitrary strictly positive real number (see Fig. 3). The 'radius' r can then be expressed as a function of T and X :

r(T , X ) =  3 2 ✓ X T ◆ 2/3 r 1/3 g , (40)
and the components of g in the coordinate basis associated to T and X are g TT = 1, g XX = r g /( 2 r), g TX = g XT = 0. Note that the condition ⇢ ⌧ transcribes into X T .

Let D be the domain where g XX 1. This domain is characterized, in (T , X ) coordinates, by the condition

X  T + 2 3 2 r g . (41)
In D the metric g can be identified with the metric G (see Eq. ( 30)). This identification defines an angle ✓ which depends on T and X by

(cos ✓ )(T , X ) = s r(T , X ) r g . (42)
A QW in D can be defined by complementing this choice of theta by a choice of the other three angles. All simulations were done with

↵(T , X ) = 0 ⇠(T , X ) = ⇡ ⇣ (T , X ) = ⇡ 2 . ( 43
)
This QW has already been considered in Ref. [32].

The condition defining D can be rewritten as r  r g / 2 , The domain D thus includes, for all , the singularity located at r = 0. For > 1, r g / 2 < r g and D is then entirely located inside the horizon. For = 1, D coincides with the interior of the horizon, and D extends outside the horizon for < 1. Ref. [32] offers plots of the spatial density | (T , X )| 2 for several initial conditions. These plots confirm that the QW follows to a great accuracy the radial null geodesics of the Schwarzschild metric, except perhaps as the QW approaches the singularity. This phenomenon is explored in detail by Fig. 3. The plots reveal the existence of interesting 'interferences' near the singularity (see (a1) and (a2)), which seem to disappear as ✏ tends to zero.

Conclusion

We have revisited the continuous limit of discrete time QWs on the line, keeping every step or only one step out of two. We have identified all families of walks which admit a continuous limit and obtained the associated PDEs. In all cases but one, the PDE describes the propagation of a Dirac fermion coupled to an electric field and, possibly, to a general relativistic gravitational field. We have also illustrated these conclusions by new numerical simulations.

Let us now briefly discuss the above results.

As mentioned in the introduction, all above literal computations are based on a new method first introduced in Refs. [30][31][32]29]. New to this article is all the material presented in Section 5. The Dirac equation obtained in Section 4 has already been presented in Refs. [30][31][32]29], but without the important discussion of the U(1) gauge invariance. The discrete gauge invariance presented in Section 2 is also new. Let us mention in this context that QWs coupled to a uniform and constant electric field have also been considered theoretically and experimentally in Refs. [44,45]. These so-called 'electric walks' are particular cases of the walks considered in Refs. [30][31][32]29] and in Section 4 of the present article. In Ref. [44], the constant and uniform electric field is put by hand on the equations of motion of the walks. In contrast, the approach developed in the present article makes it clear that the electric field is simply a manifestation of the time-and space-dependence of the angles defining the walks. This approach also allows for a straightforward generalization to nonconstant and/or nonuniform electric fields (Sections 4 and 5), and to gravitational fields 5. The electric and gravitational fields coupled top the QWs thus clearly appear as synthetic gauge fields [46].

Experiments with time-and space-dependent coins are now possible thanks to at least two recently developed techniques [47,48]. For example, one can use integrated waveguide circuits [47] where the quantum coin and step operator at a given point and time corresponds to a precise beamsplitter that can be individually set. These techniques could be used to implement experiments where arbitrary electric and/or gravitational fields are simulated by properly choosing the timeand space-dependence of the angles ↵, ⇠ and ✓.

The work presented in this article should be extended in several directions. One should first determine how the new method works, and what it delivers, when one keeps only one step out of n for arbitrary n > 2. Extensions to higher dimensional space and/or to higher dimensional Hilbert space are also desirable. In particular, the fact that some QWs on the line can be interpreted as the propagation of charged massless Dirac fermions suggests that QWs could be useful in modeling charge transport in graphene [49,50], both natural and artificial [51,52]. Let us note that the inherent discreteness would give QWs a strong computational advantage over the more traditional models based on PDEs. Finally, determining systematically the continuous limit of nonlinear QWs [53] and of walks in random media [54] should also prove interesting.

Quantum Decoherence: An introduction

There is no room for a real collapse in the purely unitary models of measurements [START_REF] Neumann | Mathematische grundlagen der quantenmechanik[END_REF] The study of quantum decoherence, a milestone of quantum theory, is a relatively young subject. The Copenhagen interpretation has been predominant in the first half of 20th century based on a strict distinction between the classical macroscopic world and the microscopic quantum realm [30]. The border between those two worlds has been unclear for a long time in that period and the debate about the nature of transition from one to the other pertained more to the philosophical domain than to the physical or mathematical one. The main problem concerned the duality inherent in the quantum theory's axioms. According to them there are two ways for the state of a quantum system to change: unitary evolution, which is deterministic and reversible, and measurement, which is probabilisitic and not reversible. Each of them obeys to two different dynamical laws, and the discontinuity in-between is marked by a sudden collapse of the wave function on one of the eigenstates of the system. All this led people to suspect that quantum theory was not complete [7].

The interaction between a single quantum system and both the classical environment and the frontier quantum/classical physics is still an hard puzzling problem. In a seminal paper Schrödinger [35] illustrated this mechanism through the paradox of an imaginary cat which exists in a superposition of dead and alive states. He coined the term Verschränkung (the entanglement) as the result of the interaction between the macroscopic apparatus and the microscopic system. In contemporary physics, decoherence is still based on this latter idea.

The physical system S is entangled with the environment E that encode the informon about Chapter 3. QUANTUM WALKS, DECOHERENCE AND RANDOM SYNTHETIC GAUGE FIELD some "monitored" states of S. The cat states are defined as the superposition of these states. Schrödinger's cat states and quantum decoherence have been widely investigated in several theoretical and experimental seminal works [27,39,33,15,10].

The first mathematical treatment of quantum decoherence has been proposed by von Neumann [28], who postulated that the measurement bring about the quantum system into a statistical mixture in some "preferred" basis, corresponding to the eigenvalues of the measured observable. Once the system is in this "mixture", the information about the system can be described by classical probability. The foundations of the theoretical contemporary approach were laid by H.D. Zeh and W. H. Zurek. Between 1970s and 1980s they contributed to build a complete theory on the emergence of classicality in the framework of open quantum systems ( [38]).

Following [42], the idea that the openness of quantum systems might have anything to do with the transition from quantum to classical was ignored for a very long time. Even though Gottfried [17] anticipated some later developments and Zeh [40,41] ventured saying that it was not possible to isolate a macroscopic quantum system, the understanding of "how the environment distills the classical essence from quantum systems" [42] has been reached much more recently.

The decoherence mechanism has been widely investigated in recent years, and within the outline of physical systems it has been the focus of various experimental setups. However, far to be simply a fundamental problem, decoherence represents today the key obstacle to the realization of quantum computers [24]. In order to minimize the effects of decoherence on a qubit, we should minimize the interaction of these systems with the environment, while keeping it sufficiently open to permit quantum operations. In order to prolong the characteristic time of decoherence, quantum physics researchers have developed several techniques and schemes, such as the active quantum error correction [25,14] or the environment engineering [37,11]. Furthermore, new controllable quantum systems have been built to better understand decoherence and quantum entanglement. Among these new actors, we should cite QA and among them the QWs. Let us recall that the interest upon these new quantum systems is twofold: (i) they are versatile tools to reproduce a wide variety of quantum complex dynamics and (ii) they play today a fundamental role in the development of quantum computing.

Quantum Walks and decoherence

An overview

We have seen in the first two chapters that QWs belong to a very rich and complex class of coherent quantum dynamics with a wide variety of applications in quantum computing, quantum transports and quantum simulations. Now let us affirme that QWs dynamics are unitary and completely reversible, and that they show ballistic behaviors. However, the loss of quantum coherence due to a breakdown of unitarity, turns into diffusive behaviors or in some particular cases into anomalous diffusion [32]. Non unitarity may be caused by a byproduct of some other operations (such as a measurement/projection operation) or by the unwanted interaction with the environment. Aharonov et al. [1] were the first to prove that performing orthogonal projections on spin state causes a breakdown of unitarity in QWs dynamics, and to observe that the system behaved asymptotically as a classical system. In particular, he derived asymptotically the expression of the second moment computed on the spatial probability distribution and showed that in contrast to unitary quantum walks, quantum random walks exhibit asymptotically diffusive spreading. After him, a vast literature has numerically and theoretically shown that QWs, in case of loss of quantum coherence, behave asymptotically like classical random walks [21]. However, even when displaying the same spreading behavior asymptotically, the dynamic mechanisms underlying quantum decoherence in QWs may be various. For this reason the variance of the walk can be used as an immediate method to characterize dynamical properties of a decoherent QWs, but it is totally useless to discriminate the different mechanisms of decoherence. Nonetheless, a detailed study on the time evolution of the density profile of decoherent QWs has been developed by Dür et al. [16].

Over the last few years, some important physical implementations have led us to a deeper understanding of decoherent QWs. Karski et al. [20] realized a QW with a single laser-cooled cesium atoms trapped in the potential wells of a 1D optical lattice. Here the coherence of the system was gradually suppressed by coupling internal and external degrees of freedom by diabatic transport. Broome et al. [8] introduced in a single-photon QWs dynamics a controlled amount of decoherence. In the coherent case photons propagated through six calcite beamdisplacers perfectly aligned and the tunable decoherence was introduced by an intentional misalignment of QW steps, which introduced a dephasing in the evolution operator of the system. The suppression of interference induced the QWs to fully loose coherence and to behave classically. In another experience, Schreiber et al. [34] introduced decoherence by activating tunable fluctuations in the vertical and horizontal polarization of a coherent diode laser reproducing the QWs unitary dynamics. In such a way they randomize the coin state of the system causing transition to diffusion behavior. The earliest experimental analysis on the different decoherence mechanisms was proposed by Alberti et al. [4]. They proposed how Markovian and non Markovian decoherence act on coin and position state space, and came up with new methods to discriminate between spin and spatial decoherence models by implementing a Hadamard QW in neutral atoms experiment.

Let us observe that the reason to introduce intentionally decoherence in unitary QWs1 is twofold. On one hand, it allows us to test the robustness of the quantum system in interaction with the environment, to measure the characteristic time of decoherence and to understand the underlying model of decoherence. This can be very useful for experimentalists. On the other, we can qualitatively and quantitatively characterize the transient dynamics between classical and quantum realms by a tunable toy model which is easily reproducible in physical experiences.

Theoretically speaking, there are two ways to model decoherence in QWs dynamics. One is to replace the unitary standard evolution by a nonunitary one, introducing projection operator that are responsible for the damping of the off-diagonal terms of the density matrix [21,9,36,13]. The other way [2,3,19,29], which is developed in the present chapter, consists in randomizing the parameters of the walks in coin or position state space, and then in averaging them (we will present this model in all the details in sec. (3.3)). In many cases both ways appear to be equivalent and in those cases we should be able to move from a formalism to the other. Before presenting both the methods, let us present some useful tools to describe qualitatively the QWs dynamics in presence of decoherence.

A qualitative picture

In modeling quantum decoherence, density matrices play an important role and they formally provide a method to represent the quantum statistics of the system. Let us recall that the density matrix is the matrix representation in a particular basis of the density operator Ω, corresponding to the pure state |™i, defined as:

Ω = |™ih™|. (3.1)
Suppose that |™i = P i c i |√ i i, the density matrix written in basis {|√ i i} then transcribes in:

Ω = |™ih™| = X i j c i c § j |√ i ih√ j |. (3.2)
All the elements of the density matrices for which i 6 = j encode the coherence of the system between the different components |√ i i and are often called interference terms or off-diagonal terms. These elements tend to vanish when the system behaves classically, and can thus be considered as a measure of the quantumness of the system. In one-dimensional QWs, presented in this chapter, randomizing the coin state destroys the phase coherence -the interference terms -among the components of the state vector. However, pure states are not sufficient to describe a quantum system in general, in particular when it is not a closed system anymore. In the interaction with the environment, the system looses (or disperses in the environment) information available to the observer. In such a situation we can describe the state as a mixed state, i.e. the set of all possible pure states √ i with associated classical probabilities p i . This is therefore a classical ensemble. The mixed state density matrix encodes the ignorance about the state of the system and reads:

Ω = X i p i |√ i ih√ i |. (3.3)
The measure of purity of the system is important to characterize the transition from pure to mixed state and to denote the information loss. This measure of the system's "mixedness" can be computed as Tr ( Ω2 ) and is bounded from below by 1/d , where d is the dimensionality of the Hilbert space and from above by 1, in case of pure state.

Another formal way to quantify the loss of information encoded in the initial pure quantum state is the von Neumann entropy, introduced by the mathematician in 1927. This measure can be seen as the classical notion of entropy in statistical mechanics and reads:

S( Ω) = °Tr ( Ω log 2 Ω) = °X i ≤ i log 2 ≤ i , (3.4)
where ≤ i denotes the ith eigenvalue of the density matrix. As in classical mechanics, the loss of information is necessarily associated with an increase of entropy. Let us remark that in case of unitary QWs, we can measure the total entropy as well as the entanglement entropy between coin state space and position state space. Therefore the entanglement entropy associated to one subsystem is derived from the total entropy tracing over the degree of freedom of the other subsystem:

S r ( Ωr ) = °Tr ( Ωr log 2 Ωr ).

(3.5)

In particular, as QWs are bipartite system consisting of position-coin states, for any ™ of the system it is possible to find a new bi-orthonormal basis in which the reduced density operators Ωr are both diagonal and share the same eigenvalues. This is possible thanks to the Schmidt decomposition 2 and the basis is called the Schmidt basis. As a consequence of the Schmidt decomposition theorem, the entanglement entropy of the two subsystems are equal.

If we consider the entropy of one subsystem, we can use it as a measure of composite system's entanglement. Let us look at for example the trace on the position states:

S r ( Ωr C ) = °X AE ≤ r AE log 2 ≤ r AE , AE = {c, c 0 }, (3.7)
where the reduced density matrix Ω r C = P m=m 0 Ω c,c 0 m,m 0 . In an analogous way we can define Ω r P , tracing on the coin states. Let us remember that the entanglement entropy always reaches a limiting value depending on the initial state of the coin state. In Fig. 3.1 we compute the entanglement entropy of a coherent QW, with a quantum coin B µ , for different values of µ. We can observe that the entanglement entropy S r ( Ωr C ) (or equivalently S r ( Ωr P )) admits an absolute maximum for µ=º/4 when the largest part of the probability is uniformly spread. It admits a minimum for µ=º/2 and µ=0 when the particle propagates without interference.

If the system is subject to decoherence the entanglement entropies of the subsystems are in general no longer identical and the total entropy and the entropy of the subsystems increase.

An alternative way to investigate the breakdown of unitarity is represented by the Wigner function, introduced recently by Hinarejos et al. [18] and independently by Alberti et al. [4] for discrete time QWs. This formalism is advantageous in studying several quantum mechanical aspects of QWs and their dynamical behaviors in the quantum-to-classical transition. In contrast to the density matrix approach, this method permits to investigate various mechanisms underlying the loss of coherence, which are more understandable in phase space than in position space [26]. In order to define the Wigner function let us introduce a simple QW in a one dimensional lattice with two internal states represented by the density operator Ω. The Wigner function W c,c 0 m (k), as introduced in [4] reads :

W c,c 0 m (k) = 1 º Z º 2 °º 2 dk 0 e °2i k 0 x hk °k0 , c| Ω|k + k 0 , c 0 i, (3.8)
where the moments k 2 R. The standard probability in position and in momentum space is simply recovered by computing the marginals of the Wigner function:

|™ m | 2 = X c,c 0 Z 1 °1 dkW c,c 0 m (k) |™(k)| 2 = X c,c 0 X m W c,c 0 m (k).
(3.9)

2 Theorem: (Schmidt decomposition) Suppose |™i is a pure state a composite system AB. Then there exist orthonormal states |i A i for system A, and orthonormal states |i B i of system B such that:

|√i = X i ∏ i |i A i|i B i (3.6)
where the coefficients ∏ i are non-negative real numbers satisfying The main difference between the Wigner function and the standard probability distribution is that it can take negative values, and we can interpret this non-positivity as a measure of non-classicality of the system. For instance in quantum optics, this is often interpreted as a signature of non-classical states of light, caused by a quantum interference phenomenon.

P i ∏ 2 i =
More specifically in the study of QWs dynamics, the Wigner function, which permits to characterize position-momentum quantum correlation and quantum interference in phase space, leads to a clearer comprehension of the differences between spatial and spin loss of coherence and their short-time behaviors for the walk [4].

Projections cause spin and spatial decoherence

A recent approach to quantum decoherence, widely used in quantum information, has been introduced by Preskill [31] who generalized the von Neumann projective measurement. His formulation facilitated a general description of the evolution of the density matrix, including the evolution of pure states to mixed states, exactly as unitary transformations supported a description of coherent quantum evolution. This section will be devoted to introduce one of the main ways used to describe quantum decoherence in QWs.

Consider the initial state of the quantum walk represented by the density matrix Ω j ,m,m 0 and U the one-step unitary operator (1.8), composed by a coin operator C followed by a translation operator T . The unitary evolution of Ω j ,m,m 0 after M steps will be described by:

Ω j +M ,m,m 0 = U M Ω j ,m,m 0U † M .
(3.10)

The Eq. (3.10) is equivalent to the Eq. (1.8) in case of pure states and can more generally describe the time evolution of any quantum states, such as mixed states. If we now want to describe decoherence, we should include in the evolution a discrete operation which violates the unitarity of the walk. In general this can be accomplished applying a projector operator to the evolution (3.10):

Ω j +M ,m,m 0 = (PU ) M Ω j ,m,m 0 (U † P † ) M . (3.11)
Observe that here we can assume, for simplicity, that the operator Q = PU is invariant under spatial translations. We should relax this assumption if we want to consider, for instance, inhomogeneous quantum coin operators. In addition to this, the nonunitary operator Q acts at each time step independently and following Kendon [22], it can be unpacked in a deterministic and reversible part and in a Markovian part. The latter contributes to the damping of the interference terms at constant rate p. For one time step the Eq. (3.11) reads in general:

Ω j +1,m,m 0 = (1 °p)U Ω j ,m,m 0U † + p X j Q j Ω j ,m,m 0 Q † . (3.12)
The above dynamical equation describes the decoherent evolution of a system that is projected by Q j at each time step with probability p. Alberti et al. [4] consider more in detail different nonunitary operations separately acting on the spin space and position space with probability p c and p s respectively. Kendon and Tregenna [23] proved that the action of decoherence only in position space (p c =0) or only in coin space (p s = 0) is largely different if we consider the short time behavior of the evolution of spatial density probability. In particular they showed that if decoherence acts on coin state, it causes the formation of a cusp around m = 0 in the transition from quantum-to-classic behavior. In contrast, in the case of an initial loss of coherence in position space, the distribution takes the shape of a top-hat profile. Lopez and Paz [26] and Alberti et al. [4] proved respectively, on cycle and on one dimensional finite lattice, that these dynamics displayed qualitatively different behaviors by investigating systematically QW dynamics in Wigner representation. Furthermore, the latter group of researchers demonstrated that spin decoherence acts on spatial coherences, suppressing with an exponential law the spatial correlation of the walk. Nevertheless, in some special cases, such as when the initial condition is not localized in momentum space, a spatial coherence can persist in a fully spin decoherent system. On the other hand the spatial decoherence destroys directly the off-diagonal terms Ω m,m 0 6 =m with a rate 1 °ps per time step and the decoherence in position space implies decoherence in coin state because the position is conditioned by the state of the coin [23].

Let us note that the model of decoherence presented above can be analyzed in momentum space. Straightforward computations allow us to derive analytical properties of the walk, as in Brun et al. [9] who found the exact analytical expression of first and second order moments and computed the diffusion coefficient of the asymptotic transport.

Publication : "Discrete-time Quantum Walks in random artificial Gauge Fields."

As we have anticipated previously, another formal way to model decoherence in QWs can be randomizing the parameters of the step operator and then averaging the properties of the walk. Several authors have developed this technique in different specific frameworks, introducing random fluctuations in the coin operator [2,3,19,29] or in the properties of the lattice [12].

In the model that we propose in the following publication, the parameters of the coin state are random in time according to a uniform distribution. The asymptotic properties of the averaged transport are derived and we briefly discuss the formal continuous limit of the dynamical averaged evolution of the density matrix. However, even if similar results have been already found in the literature, here we introduce for the first time a formal connection with the fermion propagation in random synthetic gauge field and discuss the implications.

I. INTRODUCTION

Discrete time quantum walks (DTQWs) are simple formal analogues of classical random walks. They were first considered by Feynmann in [1], and then introduced in greater generality in [2] and [3]. They have been realized experimentally [4][5][6][7][8][9][10] and are important in many fields, ranging from fundamental quantum physics [10,11] to quantum algorithmics [12,13], solid state physics [14][15][16][17] and biophysics [18,19].

It has been shown [20][21][22] recently that several DTQWs on the line admit a continuous limit identical to the propagation of a Dirac fermion in artificial electric and gravitational fields. These DTQWs are thus simple discrete models of quantum propagation in artificial gauge fields. Here, we consider artificial gauge fields which depend randomly on time and investigate analytically and numerically how this randomness influences quantum propagation. The analysis presented in this article is based on a direct analytical computation of the exact evolution equation obeyed by the average density operator. This presents several advantages. First, the average dynamics is thus known exactly, without the noise inherent in any numerical evaluation of averages. Second, knowing the exact average equations of motion makes it possible to study the average dynamics analytically. Finally, simulating directly the exact analytical equations of the average dynamics o↵ers a significant gain in computation time over alternative methods where the average evolution is determined by simulating successively a large number of realizations of the random DTQWs.

Random DTQWs have already been studied by several authors (see for example [23][24][25][26][27][28]) , but the influence of random gauge fields has never been the object of specific analytical computations. In particular, exact expressions of the asymptotic density profiles as functions of the randomness caracteristics have never been computed. Our main results are (i) DTQWs interacting with artificial gauge fields which are random in time decohere and behave asymptotically like classical random walks (ii) the asymptotic density profiles of the DTQWs are Gaussian and we give exact analytical expressions of the asymptotic di↵usion coe cients as functions of the noise amplitude which generates the randomness. We also support all results by direct numerical simulations of the average dynamics and finally discuss the continuous limits of the DTQWs interacting with random artificial gauge fields.

II. A FAMILY OF DTQWS COUPLED TO ARTIFICIAL ELECTRIC AND GRAVITATIONAL FIELDS

A. Wave-function evolution

In physical space

We consider discrete time quantum walks in one space dimension driven by a time-dependent quantum coin acting on a two-dimensional Hilbert space H. The walks are defined by the following finite di↵erence equations, valid for all (j, m) 2 N ⇥ Z:

 L j+1,m R j+1,m = B (✓ j , ⇠ j )  L j,m+1 R j,m 1 , (1)
where

B(✓, ⇠) =  e i⇠ cos ✓ isin ✓ i sin ✓ e i⇠ cos ✓ . (2)
The operator represented by the matrix B is in SU (2) and ✓ and ⇠ are two of the three Euler angles. The index j labels instants and takes all positive integer values. The index m labels spatial points. We choose to work on the circle and impose periodic boundary conditions. We thus introduce a strictly positive integer M and restrict m to all integer values between M and +M i.e. m 2 Z M . Results pertaining to DTQWs on the infinite line can be recovered by letting M tend to infinity.

For each instant j and each spatial point m, the wave function jm = L jm b L + R jm b R = a jm b a , a 2 {L, R}, has two components L jm and R jm on the spin basis (b L , b R ) and these code for the probability amplitudes of the particle jumping towards the left or towards the right. Note that the spin basis is interpreted as being independent of j and m. For a given initial condition, the set of angles {✓ j , ⇠ j , j 2 N} completely defines the walks and is arbitrary.

It has been proven in [20][21][22] that walks from this family are models of Dirac fermions coupled to artificial electric and gravitational fields. Details can be found in these references and in the first appendix to the present article.

In Fourier space

A practical tool to study quantum walks on the discrete circle is the discrete Fourier transform (DFT). Let (A m ) m2ZM be an arbitray sequence of complex numbers defined on the discrete circle. The DFT of this sequence is the sequence ( Âkn ) n2ZM defined by

Âkn = +M X m= M A m exp (ik n m) (3)
with k n = 2n⇡/(2M + 1), n 2 Z M . The original sequence can be recovered from its DFT by the relation:

A m = 1 2M + 1 +M X n= M Âkn exp ( ik n m) . (4)
For infinite M i.e. DTQWs on the infinite line, the DFT of an infinite sequence (A m ) m2Z becomes a function

Â(k) = X m2Z A m exp (ikm) (5)
defined for k 2 ( ⇡, ⇡) and the inverse relation reads:

A m = 1 2⇡ Z ⇡ ⇡ Â(k) exp ( ikm) dk. (6)
In Fourier space on the infinite line, the evolution equation (1) transcribes into

" ˆ L j+1 (k) ˆ R j+1 (k) # = C(✓ j , ⇠ j , k) " ˆ L j (k) ˆ R j (k) # (7)
where

C(✓ j , ⇠ j , k) =  e i⇠ cos ✓e ik i sin ✓e +ik i sin ✓e ik e i⇠ cos ✓e +ik (8)
for all k 2 ( ⇡, ⇡).

B. Density operator evolution

In physical space

The walks can also be described using the density operator ⇢ = ⇤ ⌦ . We introduce the basis v

1 = b L ⌦ b L , v 2 = b L ⌦ b R , v 3 = b R ⌦ b L , v 4 = b R ⌦ b R
⇢ LL j+1,m,m 0 ⇢ LR j+1,m,m 0 ⇢ RL j+1,m,m 0 ⇢ RR j+1,m,m 0 3 7 7 5 = Q (✓ j , ⇠ j ) 2 6 6 4 ⇢ LL j,m+1,m 0 +1 ⇢ LR j,m+1,m 0 1 ⇢ RL j,m 1,m 0 +1 ⇢ RR j,m 1,m 0 1 3 7 7 5 (9)
where

Q(✓, ⇠) = 2
with c = cos ✓ and s = sin ✓. The probability to find the walk at time j at point m is N jm = ⇢ LL j,m,m + ⇢ RR j,m,m and the sum P m N jm is independent of j i.e. it is conserved by the walk. Contrary to equation (1), equation ( 9) can be used to describe walks with initial conditions which are not pure states. Equation ( 9) is thus more general than (1).

In Fourier space

Consider now, for any instant j, the double DFT of the density operator ⇢ j,m,m 0 , which we denote by ⇢j (k, k 0 ) or, alternately, ⇢j (K, p) where K = (k + k 0 )/2 is conjugate to m + m 0 and p = (k 0 k)/2 is conjugate to m 0 m. For DTQWs on the infinite line, the range of both K and p is ( ⇡, +⇡). The DFT of the density operator obeys ⇢j+1 (K, p) = R (✓ j , ⇠ j , K, p) ⇢j (K, p) with

R(✓, ⇠, K, p) = 2 6
Note that the operator R governing the evolution of ⇢ is unitary. This can be checked by a straightforward computation and it is a direct consequence of the unitarity of the operator B.

III. RANDOMIZING THE FIELDS AND AVERAGING THE DYNAMICS A. Randomizing the fields

The Hadamard walk corresponds to ⇠ = ⇠ H = ⇡/2 and ✓ = ✓ H = ⇡/4; since these angles are constant, the Hadamard walk describes propagation in the absence of electric and gravitational field [21,22]. We now consider situations where one of the angles ⇠ and ✓ does depend on time and fluctuates around its Hadamard value. More precisely, we consider two cases. Case 1 corresponds to ✓ = ✓ H = ⇡/4 and ⇠ chosen randomly at each time-step with uniform probability law in the interval (⇡/2 /2, ⇡/2 + /2), where 2 (0, 2⇡) is a fixed i.e. j-independent positive real number. As proven in [20][21][22] and detailed in the first appendix to the present article, a time-dependent ✓ is equivalent to a space-time metric whose purely spatial part depends on time, and such a metric represents a time-dependent relativistic gravitational field. Case 2 corresponds to ⇠ = ⇠ H = ⇡/2 and ✓ chosen randomly at each time-step with uniform probability law in the interval (⇡/4 /2, ⇡/4 + /2). As proven in [22], a time-dependent ⇠ is equivalent to a time-dependent 'vector' potential, which represents a time-dependent electric feld.

Thus, in each case, a realization of the random gauge field is determined by a sequence ! = (! 1 , ! 2 , ...) of independent random variables, where ! j represents the value of the random angle ✓ or ⇠ at time j. If one follows the walk till time N , the relevant random sequence is the N -uple

! N = (! 1 , ! 2 , ..., ! N ).
For each value of and each instant j, ! j is uniformly distributed in the interval I = (! H /2, ! H + /2) centered on the Hadamard value ! H . The probability density of ! j in this interval is thus simply p (! j ) = 1/ and is independent of both ! j and j. The probability density for

! N = (! 1 , ! 2 , ..., ! N ) in I N is therefore P (! N ) = ⇧ N j=1 p (! j ) = 1/ N and is independent of ! N .

B. Averaging the dynamics

At fixed initial condition ⇢ 0 and for each time N , the density operator ⇢ N at time N depends on the realization ! N of the random angle up to time N . At fixed initial condition, the easiest way to compute statistical averages over ! N is to first compute the statistical average ⇢N of the density operator over

! N : ⇢N = Z I N ⇢ N (! N )P (! N )d! N = Z I N ⇢ N (! 1 , ..., ! N )p (! 1 )...p (! N )d! 1 ...d! N = Z I N ⇢ N (! 1 , ..., ! N ) 1 N d! 1 ...d! N . (12)
Let us work in Fourier space. One can then write, for any realization ! N = (! 1 , ..., ! N ) of the random angle up to time N :

⇢N (! 1 , ..., ! N ) = R(! N )⇢ N 1 (! 1 , ..., ! N 1 ) = R(! N )...R(! 1 )⇢ 0 , (13)
where the variables K and p have been omitted for clarity reasons. Since the ! j 's are statistically independent of each other and are identically distributed, one obtains:

⇢N = RN ⇢0 ( 14
)
where R is the statistical average of the evolution operator R over the random angle ! = ✓ or ⇠ (the other angle being fixed to its Hadamard value):

R(K, p, ) = Z I R(!, K, p)p (!)d! R(K, p, ) = Z I R(!, K, p) 1 d!. ( 15
)
The average evolution R is thus a function of (K, p) and of the noise parameter and can be computed analytically from (11). It determines the evolution of the average density operator completely and, therefore, the average transport. Since everything that follows pertains only to the average transport, we simplify the notation by droping the bar on the letter ⇢ and the density operator of the averaged transport will now be designated simply by ⇢.

A direct computation from (11)

It proves convenient for all subsequent computations to change basis in ⇢ space and introduce the new vectors

u 1 = v 1 + v 4 , u 2 = v 1 v 4 , u 3 = v 2 + v 3 , u 4 = v 2 v 3 .
In this new basis, the components of Re (K, p, ) and Rg (K, p, ) read:

Re (K, p, ) = 2 6 6 4 cos(2K) i sin(2K) 0 0 0 0 sinc( /2) cos(2p) i sinc( /2) sin(2p) i sinc( /2) sin(2K) sinc( /2) cos(2K) (1 sinc( )) 2 cos(2p) i (sinc( ) 1) 2 sin(2p) 0 0 i (1+sinc( )) 2 sin(2p) (sinc( )+1) 2 cos(2p) 3 7 7 5 , (18)
Rg (K, p, ) = We choose as initial condition the pure state defined by j=0,m=0 = (b L + ib R )/ p 2 and j=0,m = 0 if m 6 = 0. This state corresponds to the density operator

⇢ j=0,m=0,m 0 =0 = (b L ⌦ b L + b R ⌦ b R + i(b R ⌦ b L b L ⌦ b R )) /2 = (u 1 iu
4 )/2 and ⇢ j=0,m,m 0 = 0 if m 6 = 0 or m 0 6 = 0. In Fourier space, ⇢j=0 (K, p) = (u 1 iu 4 )/2 for all K and p. For any realization of the noise i.e. for any given value of !, the initial pure state evolves by the DTQW into a pure state. But the average evolutions descibed by Re and Rg both transform the initial pure state into a superposition. However, the average transport is symmetrical around the origin, as is the classical Hadamard walk generated from the same initial condition.

Let us finally stress that, contrary to the operator R governing the unaveraged transport, the operators R e/g governing the averaged transport are not unitary. This loss of unitarity generates qualitative di↵erences between the unaveraged and the averaged transport. In particular, the averaged transport looses quantum coherence and is asymptotically di↵usive. These two important consequences of the averaging process are analyzed in the remaining sections of this article.

IV. QUALITATIVE DESCRIPTION OF AVERAGE TRANSPORT

Typical density profiles of the average transport are shown in Figure 1 for random gravitational and random electric fields. For small enough values of the noise parameter , the average transport behaves at short times like the Hadamard walk and is ballistic. Ballistic behavior then gradually disappears and is replaced by di↵usive behavior. For larger values of , ballistic behavior is replaced, even at short times, by di↵usive behavior. Note that the Gaussian-like form of the asympotic density profiles presents a central dip when the DTQWs interact with random gravitational fields, but presents a central cusp when the DTQWs interact with random electric fields. Asymptotically, DTQWs in electric and gravitational fields which are random in time thus behave like classical random walks. This means that the randomness in the fields prompts the DTQWs to loose coherence. This can be confirmed by considering the spin coherence defined by

! "! #! $! %! &! '! (! )! "! !$ "! !# "! !" "! ! * + + !,!-& !,!-) !,"-" !,
C j = max m,m 0 | ⇢ LR j,m,m 0 | . ( 20
)
Figure 2 displays the typical time-evolution of the spin coherence for various values of the noise parameter . These results confirm that the average transport loses spin coherence and that a higher value of leads to a quicker loss of spin coherence.

A brief comment on spatial coherence is in order. The retained initial condition vanishes everywhere except at m = m 0 = 0. If one prefers, the Fourier transform of the initial density operator is flat in both K and p space. There is thus initially no spatial coherence. As time increases, the Fourier transform ⇢(K, p) of the density operator ⇢ becomes non flat in both K and p (see for example the asymptotic form (21) of ⇢). In other words, each K-mode acquires spatial coherence. But ⇢(p) = P K ⇢(K, p) remains flat in p (data not shown) i.e. there is no total gain of spatial coherence.

The entanglement of the averaged dynamics can also be quantified by the Shannon entropy S r of the reduced density operator ⇢ r in spin space. To be precise [29][30][31], ⇢ r = P m ⇢ m,m 0 =m and the Shannon entropy S r = tr(⇢ r log(⇢ r )). The time-evolution of S r is prensented is Figure 3, together with the entanglement entropy of the pure Hadamard walk with the same initial condition, which admits 0.872 as asymptotic value [32]. The increase in S r signals the loss of coherence and the figure confirms that this decoherence by noise gets more e↵ective as increases.

The scaling of the decoherence time for small values of can be evaluated by the following reasonning. As previously explained, the operator Re/g completely controls the average dynamics. For = 0, there is no noise and the DTQW never decoheres i.e. the decoherence time is infinite. The first non-vanishing terms in the expansion of Re/g around = 0 are of second order in . Thus, per time step, the e↵ect of the noise on the DTQW is of order 2 for small enough values of . The typical decoherence time therefore scales as 2 for small values of . The next section, together with the appendices, provides an analytical investigation of how coherence is lost. In particular, the asymptotic form of the density operator is computed exactly. The corresponding density is Gaussian, which confirms that the DTQW behaves asymptotically like a classical random walk. Also, the asymptotic density operator is proportionnal to u

1 = v 1 + v 4 = b L ⌦ b L + b R ⌦ b R .
This proves that the spin coherence, which measures the amplitude of the b L ⌦ b R component, vanishes asymptotically, in accordance with Figure 2.

V. QUANTITATIVE DESCRIPTION OF THE ASYMPTOTIC REGIME A. Central limit theorem

The average dynamics is entirely determined by the eigenvalues e/g r and corresponding eigenvectors w e/g r , r = 1, 2, 3, 4, of the operators Re/g . As evident from Figure 1, the density profiles of the average transport become larger and smoother with time. This suggest that the asymptotic dynamics can be understood by computing the eigenvalues and eigenvectors only for values of K much smaller than unity. The detailled analysis, though very instructive, is too involved to merit inclusion in the main body of this article and it is therefore presented in the Appendix. The main conclusion can be stated as follows.

Theorem. Let K j = K ⇤ / p j where K ⇤ is an arbitrary but j-independent wave number. The average density operator in Fourier space admits as the time j tends to infinity the following approximate asymptotic expression:

⇢e/g j (K j , p) ⇠ 1 2 exp ⇣ ↵ e/g (p, )jK 2 j ⌘ u 1 ( 21
)
where

↵ e (p, ) = 2 3 + (sinc( )) 2 + 2 (sinc( /2)) 2 (1 + sinc( )) + 4 cos(2p) ⇣ sinc( ) + (sinc( /2)) 2 ⌘ 3 + (sinc( )) 2 2 (sinc( /2)) 2 (1 + sinc( )) + 4 cos(2p) ⇣ sinc( ) (sinc( /2)) 2 ⌘ ( 22
)
and

↵ g (p, ) = 2 1 + (sinc( )) 2 1 (sinc( )) 2 . ( 23
)
This result is a central limit theorem which proves that the asymptotic density operator is approximately Gaussian in K-space, with a typical width (in K-space) which decreases as j 1/2 , as in classical random walks and non quantum di↵usions. Note that ↵ g is actually independent of p.

One of the consequences of ( 21) is that the projection of ⇢j=J (K J , p) on the subspace spanned by (u 2 , u 3 , u 4 ) tends to zero. Remembering the expressions of the u i in termes of b L and b R , this means that ⇢LR , ⇢LR and ⇢LL ⇢RR all tend to zero as J tends to infinity. The component along u 1 coincides with ⇢LL + ⇢RR and determines the asymtotic density of the averaged walk after summation over p and Fourier transform over K.

B. Asymptotic mean-square displacement

Let us now explicitly compute the asymptotic expression of the mean-square displacement m 2 e/g in the special case of a random DTQW on the infinite line. Switching back the original spatial variables m and m 0 involves a double integration over K and p. The 2D measure to be used in this integration is dkdk 0 = 2dKdp. The density N e/g jm at time j and point m is the trace over m 0 = m of the component of the density operator along the basis vector u 1 . Expression (C10) for ⇢e/g is only valid for K ⌧ 1 (see the Appendix). But the functions ↵ e/g (p, ) are always non vanishing. The width K(j, p) of ⇢e/g j (K, p) in K thus scales as 1/ p j and tends to zero as j tends to infinity. Thus, for large enough j, the density and mean square displacement are given by:

N e/g jm = 1 4⇡ 2 Z ⇡ p= ⇡ dp Z ⇡ K= ⇡ dK exp ⇣ ↵ e/g (p, )jK 2 ⌘ exp ( iKm) ( 2 4 )
and

m 2 e/g (j, ) = 1 4⇡ 2 X m2Z m 2 Z ⇡ p= ⇡ dp Z ⇡ K= ⇡ dK exp ⇣ ↵ e/g (p, )jK 2 ⌘ exp ( iKm) . ( 25
)
Since the width K(j, p) of ⇢j,K,p scales as 1/ p j, one can also replace all discrete summations over m by integrals over the real line, because K(j, p) ⇥ x = 1/ p j ⇥ 1 ⌧ 1 for large enough j. Indeed, a simple computation confirms that the integrated density R R dmN e/g jm (with N e/g jm given by ( 24)) is equal to unity at all times j. Replacing in (25) the discrete summation over m by an integral delivers

m 2 e/g (j, ) = j ⇡ Z ⇡ ⇡ ↵ e/g (p, )dp. ( 26
)
The computation of m 2 g (j, ) is trivial because ↵ g (p, ) does not depend on p. One finds

m 2 g (j, ) = 2D g ( )j (27)
with

D g ( ) = 2 1 + (sinc( )) 2
1 (sinc( ))

2 .

(

) 28
The exact expression for m 2 e (j, ) is more involved. A direct computation leads to:

m 2 e (j, ) = 2D e ( )j (29)
with

D e ( ) = 2 sinc( ) (sinc( /2)) 2 0 @ ⇣ sinc( ) + (sinc( /2)) 2 ⌘ + 2 (sinc( /2)) 2 ⇣ (sinc( )) 2 + 2sinc( ) 3 ⌘ s( ) 1 A (30)
with

s( ) =  ⇣ 3 + (sinc( )) 2 2 (sinc( /2)) 2 (1 + sinc( )) ⌘ 2 16 ⇣ sinc( ) (sinc( /2)) 2 ⌘ 2 1/2 . ( 31
)
In both electric and gravitational case, the asymptotic mean square displacement in physical space grows linearly in time, as for classical random walks and non quantum di↵usions. The functions D e and D g are the asymptotic di↵usion coe cients of the average transport. Both functions are strictly decreasing on (0, 2⇡). Thus, decoherence occurs more rapidly as increases (see Section IV), but the asymptotic di↵usion coe cients decrease with . We also note that D g ( ) < D e ( ) for all 2 (0, 2⇡).

Figure 4 shows the time-evolution of the relative di↵erence between the di↵usion coe cients computed from ( 28), (28) and the mean square displacement computed from numerical simulations for various values of . This figures clearly supports the analytical computation presented in this Section. 4: Left (right) figure: time-evolution of the relative di↵erence between the gravitational (eletric) di↵usion coe cients computed from numerical simulations and the exact analytical expressions.

! "! #! $! %! &! !& &! ! ' ( ( !()(!*% !()(&*& !()(&*# j ! "! #! $! "! !$ "! !# "! !" "! ! "! " "! # % & & !&'&$ !&'&( !&'&) j FIG.

VI. CONTINUOUS LIMIT

The formal continuous limit of the original, unaveraged evolution equations ( 9) and ( 1) has already been considered in [21,22] and coincides with the Dirac equation obeyed by a fermion minimally coupled to an electric field and/or a relativistic gravitational field. Let us now determine the formal continuous limit of the averaged evolution equations specified by the operators Re and Rg .

As shown and discussed in [21,22] for the unaveraged evolution equations, the object which admits a continuous limit for ✓ = ✓ H or ⇠ = ⇠ H is not the original walk, but the walk derived from it by keeping only one time step out of two [? ]. We thus search for the continuous limit of the following discrete equations:

⇢j+2 (K, p) = ⇣ Re/g (K, p, ) ⌘ 2 ⇢j (K, p). ( 32
)
To be specific, we restrain j to uneven positive integer values and decide to work on the inifinite line, so that K and p take all values in ( ⇡, +⇡). We now supppose that, for all uneven j = 2r + 1, ⇢ j=2r+1,m,m 0 (resp. ⇢j=2r+1 (K, p)) is the value taken by a certain function ⇢ (resp. ⇢) at 'time' t r = r and positions x m = m and x m 0 = m 0 (resp. momenta K and p). Roughly speaking, the continuous limit refers to situations where the function ⇢ (resp. ⇢) varies only little during one time step t r+1 t r = 1. A necessary and su cient condition for this to be realized is that Re/g (K, p, ) 2 be close to unity. Direct inspection reveals that this transcribes into ⌧ 1, K ⌧ 1 and p ⌧ 1. The last two conditions mean that ⇢ has caracteristic spatial variation scales much larger than the distance m + 1 m = 1 between adjacent grid points and the first condition states that the noise amplitude is small. Note that K, p and are a priori independent inifnitesimal quantities. In particular, there is no reason why K and p should be of the same order of magnitude.

The formal continuous limit is then obtained by expanding Re/g (K, p, ) 2 around K = 0, p = 0, = 0 and by replacing ⇢j+2 ⇢j by @ t ⇢. One thus gets equations of the form :

@ t ⇢(t, K, p) = ⇣ S e/g (K, p, ) 1 ⌘ ⇢(t, K, p) ( 3 3 )
where, for example,

S g (K, p, ) = 2 6 6 6 4
1 4K 2 2iK 2iK(1 p 2 2 )(1 2 6 ) 2 Kp(1 2 6 ) 2iK(1 p 2 2 )(1 2 3 ) (1 4K 2 )(1 p 2 2 )(1 2
3 )

p 2 (1 2 6 ) 2 ip(1 2 6 ) 2iK(1 2 6 ) 4K 2 (1 2 6 ) ( 1 4K 2 )(1 p 2 2 )(1 2 3 ) ip(1 4K 2 )(1 2 3 ) 2Kp(1 2 6 ) ip(1 4K 2 )(1 2 6 ) ip 1 4p 2 3 7 7 7 5 (34)
at second order in all three independent infinitesimals K, p and . These equations can be translated into physical space by remembering that iK and ip are the Fourier representations of @ X and @ y where X = (x + x 0 )/2 and y = x 0 x. The analysis presented in Sections III and IV above has been carried out with an initial condition which spreads over the whole K-and p-ranges. The resulting density operator does localize in time around K = 0, but it never localizes around p = 0 and remains spread in p-space. The continuous limit thus cannot be used to recover the results of Section IV. As can be checked directly from (34), the continuous limit equations nevertheless predict di↵usive behavior if K is much lower than both p and . A systematic study of the continuous limit dynamics for various scaling laws obeyed by K, p and falls outside the scope of this article and will be presented elsewhere.

VII. CONCLUSION

We have studied two families of DTQWs which can be considered as simple models of quantum transport of a Dirac fermion in random electric or gravitational fields. We have proven analytically and confirmed numerically that randomness of the fields in time leads on average to decoherence of the walks. The asymptotic average transport is thus di↵usive and we have computed exactly the di↵usion coe cients. We have also obtained and discussed the continuous limit of the model.

A few words about the loss of coherence in DTQWs may prove useful at this point. Pure, deterministic DTQWs are standard quantum systems in the sense that their time-evolution is unitary. They thus never loose coherence nor do they exhibit di↵usive behavior. As with any quantum system, the loss of coherence in DTQWs is induced by the so-called interaction with an environment. There are essentially two ways to model this interaction. The first one is to start from the unitary evolution of the density operator and to modify this unitary evolution into a non-unitary one by introducing so-called projector or measurement operators [23][24][25]33]. The second way of introducing decoherence is the one followed in this article. It consists in introducing some randomness in the parameters of the DTQW and in averaging over this randomness [26-28, 34, 35]. Contrary to the unaveraged density operator, the averaged density operator then follows a non-unitary evolution and this breakdown of unitarity induces the loss of coherence and the asymptotic di↵usive behavior displayed by the averaged transport.

The results of this article constitute/are an addition to the already extensive literature dealing with the asymptotic behavior of DTQWs and CTQWs. Standard deterministic QWs are famous for typically exhibiting asymptotic ballistic behavior. But di↵usive and anomalous di↵usive asymptotic behavior have also been observed [36,37], as well as localization [35,38,39] and soliton-like structures [34].

Let us conclude by listing a few natural extensions of this work. The random artificial gauge fields considered in this article have two main characteristics: they depend only on time and the associated mean fields vanish [40]. One should therefore extend the analysis presented above to situations where the mean fields do not vanish and where the artificial gauge fields depend not only on time, but also on position. In particular, the continuous limit equation derived in Section VI is markedly di↵erent from both the Caldeira-Leggett [41,42] and the relativistic Kolmogorov equation describing relativistic stochastic processes [43][44][45]. Indeed, because the random fields depend only on time, the dynamics considered in this article does not couple di↵erent (K, p)-modes, but these are coupled in both the Caldeira-Leggett and the relativistic Kolmogorov equation. Considering DTQWs coupled to artificial gauge fields which also depend randomly on position should therefore lead to master equations closer to the the Caldeira-Legget and the Kolmogorov models. Moreover, cases where both electric and gravitational fields vary randomly are certainly worth investigating.

Finally, at least some DTQWs in two spatial dimensions can be considered as models of quantum transport in electromagnetic fields [46]. The analysis presented in this article should therefore be repeated in higher dimensions to include random magnetic fields [47] and evaluate their e↵ects on spintronics.

p-space where the condition | K |⌧| p | does not apply actually shrinks to zero with time. For dynamics taking place on a finite circle (finite value of M ), computations are a little more involved but can nevertheless be carried out. We feel a detailled analysis of the problem for finite values of M does not bring any valuable insight on interesting physics or mathematics, and we thus restrict the analytical discussion of the asymptotic dunamics to DTQWs on the infinite line, where expressions (B4) and (B5) su ce.

A direct computation shows that the corrections to the eigenvectors are first order in K. By convention, we fix to unity the value of the first component of w e/g 1 (K, p, ) in the basis (u 1 , u 2 , u 3 , u 4 ). One thus gets for example

w g 1 (K, p, ) = {1, 2iKsinc( ) 2 1 sinc( ) 2 , 2iKsinc( ) 1 sinc( ) 2 , 2Ksinc( ) tan(p) 1 sinc( ) 2 }. (B6)
The expression of w e 1 is substantially more complicated and need not be reproduced here.

as is the initial condition {q 1 , ..., q N , p 1 , ..., p N } = {q, p} is the vector of the conjugate real variable q and p. The system is then described by the canonical equations:

@ @ t q = @H @p @ @ t p = °@H @q , (4.3)
where H is the Hamiltonian of the system. For a such system the Eq. ( 4.2) can be simplified in the following form:

@ @ t P (X, t ) + V(X) • rP (X, t ) = 0, (4.4)
because the conjugate variables in X satisfy the divergence relation r • V(X) = 0. Now let us suppose that Eq. (4.1) admits a set of integrals K i . We can write for each integral:

dK i d t = @K i @X @X @t = rK • V(X) = 0 (4.5)
The Eq. (4.4) with the Eq. (4.5) allows us to compute, if it exists, the stationary probability:

P st (X) = 1 Z e °Ø P i K i (4.6)
where Z is a normalization factor. The solution (4.6) is known as absolute equilibrium. In thermodynamics it represents a thermodynamic equilibrium, called Gibbs distribution, and Ø corresponds to the inverse of the thermodynamic temperature of the system. The evolution from the initial condition to the absolute equilibrium is called equilibration, or thermalization if it does not depend on the initial conditions.

Thermalization and absolute equilibria in Galerkin truncated PDEs

Absolute equilibria are typically observed, for instance, in real fluid at very high Reynolds numbers, where the energy transfers from macroscopic scales to molecular thermal energy [9]. In non-dissipative flow, artificial microscopic systems can cause an effective dissipation, leading to thermalization. This effective dissipative effect is a consequence of a Galerkin truncation that consists in taking the Fourier transform of the PDEs and retaining only a finite number of the Fourier coefficients as dynamical variables. This operation can be represented by a Galerkin projector which suppresses all the modes k greater than a given k max , i .e. the truncation wavenumber. The truncated systems recover the original system only if the convergence of the numerical scheme is ensured. In other words the spectrum has to fall down rapidly at a wave number much smaller than the truncation wavenumber.

The thermalization phenomenon of the truncated system due to the effective viscosity was discovered by Lee [17] and studied in detail by Kraichnan [12] in 3D. Later Cichowlas et al. [7] analyzed an ideal non-dissipative flow obeying to the truncated Euler equations, grasping the underlying mechanism of thermalization in this system. Recently Krstulovic and Brachet [13] [14] extended this analysis to the magneto-hydrodynamics (MHD) and the Gross-Pitaevskii (GP) equations. These equations together with the truncated compressible and incompressible Euler equations are all conservative truncated systems. The property that this 4.1. Absolute equilibrium in conservative systems 101 class of systems relaxes, with a very rich behavior, toward the statistical equilibrium is a very general fact.

From microcanonical to grand canonical ensemble

Let us examine the Eqs. (4.3), and suppose that the system admits the conserved quantities K i and the equilibrium point (4.6). Following Landau and Lifshitz [16], we should be able to define the microcanonical ensemble by the probability d w of finding the system in states with K i = K 0 i (the overscript denotes the initial time t = 0):

d w / e S Y i ±(K i °K 0 i ), (4.7)
where S = l og ≠ is the entropy with the number of accessible microstates. To obtain the microcanonical statistical states we integrate Eqs. (4.3) until the system reaches the thermodynamical equilibrium. However, the associated temperature Ø °1 is not always accessible. This is because, we cannot always obtain an explicit form for d w and S. Following Krstulovic et al. [15], to control the temperature we should use the grand-canonical distribution probability given by the Boltzmann weight:

P st = 1 Z e °ØF , (4.8)
where F = H °Pi µ i K i and Z is a normalization function, Ø °1 is the temperature and µ i are the Lagrange multiplier associated to each conserved quantities K i . Let us mention that, as the microcanonical ensemble is controlled by the variables K 0 i , the grand canonical ensemble is characterized by the conjugate µ i .

The key point that allows us to use the grand canonical ensemble, is a well known equivalence under very general circumstances between the microcanonical and grand canonical ensemble. In fact, in the thermodynamic limit, that is if the number of Fourier modes N ! 1, both statistical ensembles are expected to be equivalent [16].

Let us now observe that computing the stationary distribution (4.8) is not a trivial problem, especially when the Hamiltonian in Eqs. (4.3) is not quadratic and therefore the statistical distribution is not Gaussian. Nonetheless Krstulovic et al. [15] have introduced an algorithm that generates absolute equilibrium for all Hamiltonian dynamics of the form (4.3) admitting a stable equilibrium point X § . The idea of this algorithm is to construct a stochastic process that converges to a stationary solution with the equilibrium distribution given by (4.8).

The stochastic process can be defined by the following Langevin equations:

@ @ t q = °∫ @H @q + p 2¥∫ª 1 (t) @ @ t p = °∫ @H @p + p 2¥∫ª 2 (t), (4.9)
where hª i (t ), ª j (t 0 )i = ± i j ±(t °t 0 ) define a white Gaussian forcing term, with hª i (t )i = 0. The coefficients (¥, ∫) 2 R + £ R + . The probability distribution P of the stochastic process, defined by Eqs. (4.9), can be shown to obey the following Fokker-Planck equation: @ @ t P = ∫ ∑ @ @q ( @ @q HP) + @ @p ( @ @p HP) ∏ + ¥∫r 2 P, (

where r 2 = @ @p 2 + @ @q 2 . The latter equation is the Fokker-Planck equation of the transition probability P (p, q), a special case of the differential Chapman-Kolmogorov equation. It describes mathematically a diffusion process. The terms in square brackets are known as the drift terms and the Laplacian corresponds to the diffusion part. Let us comment that when the diffusion term is zero, the noise in (4.9) is zero. In this case the (4.10) reduces simply to the Liouville's equation, which describes a fully deterministic motion.

This Fokker-Planck equations shares the same stationary probability (4.8) 1 and we should now directly control the µ i to fix the mean values hK i i to K 0 i and to reproduce the same absolute equilibrium reached by (4.3).

In conclusion we should briefly discuss the case of complex variables. Let us introduce w = p + i q and w § = p °i q, and the Hamiltonian is the same in Eqs. (4.3). The canonical equation for the complex variables w and w § read:

@ @ t w = °i @H @w § @ @ t w § = i @H @w . (4.11)
The Langevin process (4.9) transcribe in:

@ @ t w = °∫ @H @w § + p 2¥∫•(t ) @ @ t w § = °∫ @H @w + p 2¥∫• § (t ) (4.12)
where •(t ) = ª 1 (t ) + i ª 2 (t ) is now a complex Gaussian white noise defined by h• i (t )• j (t 0 )i = ± i j ±(t °t 0 ) and h• i (t )i = 0 The Fokker-Planck associated to latter stochastic process reads: @ @ t P = °@ @w ( @ @w § P ) °@ @w § ( @ @w P ) + ¥∫ @ 2 @w@w § P (

Let us just finally observe that in case the system is described by fields we have to replace derivative by functional derivative and delta Kronecker by a Dirac delta. Nevertheless, in this case we have to fix an ultraviolet cut-off, for instance Galerkin truncate the Fourier modes of the field, to prevent the infinite energy solutions.

Nonlinear QW-like models and thermalization

Thermalization in closed quantum systems

Thermalization in closed quantum systems is a very recent research topic and concerns mainly the understanding of complex quantum many-body systems (QMBs). When, in 1902, the sem- 1 Let us notice that the Fokker-Planck presented here is different from the one Krstulovic et al. [15] found because the latter was associated with a different stochastic process. However, both shares the same stationary probability.

Nonlinear QW-like models and thermalization 103

inal work "Elementary principles in statical mechanics" by Gibbs was published, quantum mechanics was not yet formalized. In contrast to classical statistical mechanics, it is not clear yet, in quantum theory, in what way interacting QMBs can equilibrate and thermalize.

Lately some results have been achieved on how thermal states of local systems emerge from complex quantum dynamics [20,21,22,8]. Nonetheless, many questions are still open: how typical values emerge for macroscopic local observables with suppressed quantum fluctuation [21,8]? How does the system behave close to phase transitions when the correlation length is comparable to the size of the system itself and thermalization is not expected? How does the macroscopic picture of thermodynamics and its second law arise from the microscopic description?

Discrete QCA are naturally good candidates to investigate thermalization because of the discrete structure of the matter. Meyer [18], for instance, was the first to introduce QCA to describe quantum many-body interactions, in particular quantum gas. Differently from QCA , DTQWs are one-particle presenting no interaction with other systems, and obey to reversible and linear dynamics.

Nonetheless, different Nonlinear QW-like (NQWs) models on the one-dimensional unrestricted line have been considered recently by Shikano et al. [25] and Navarrete-Benlloch et al. [19]. These models display complex behaviors in physical space and chaotic dynamics phases. The main aim of this chapter is to prove that at least one of them thermalizes. In particular preliminary works on several of these QW-like (but non linear) automaton on a cycle, showed that they mix and rich a limiting stationary distribution. Otherwise, usual QWs do not mix [2].

Let me mention an interesting line of research by Romanelli [23] who fist studied thermal equilibria in QW from a radically different point of view. In fact the QWs they study is perfectly linear. Yet, what they look at is coin degree of freedom versus position degree of freedom, and find that the position acts like a thermal bath on the coin. They have close forms Liouville equations for the coin states and in the thermodynamical limit they define a temperature to control the relaxation process. Their work started with the Hadamard Walk [23], but was extended to a generalized QW Romanelli and Segundo [24]. Let us notice that they compute the temperature as function of the initial condition. Usually in the literature this process is called equilibration and not thermalization, that in general does not depend on the initial condition, like in the model we introduce in sec. 4.3.

The section will be organized as follows: in sub-section (4.2.2) we define a QW on a Ncycle, recalling the concept of mixing time; in sub-section (4.2.3) we introduce two different examples of NQW-like model. Lastly, we present, in section (4.3), the publication where we proved that a Nonlinear Optical Galton Board on a cycle, a NQW-like model, converges in longterm to an absolute equilibrium (interpreted by a classical thermalization mechanism) and displays a final steady probability distribution independent of initial conditions.

QWs on N-cycle and limiting distribution

QWs on general graphs have been extensively investigated [26] and represent a cornerstone in contemporary quantum computing. They were first introduced by D. Aharonov et al. [1] and their limiting distribution have been then fully examined by Chisaki et al. [6] and Konno [11], respectively for the discrete time and continuous time quantum walks. Graphs can be very different, with their intrinsic topological properties, and the QWs behavior can depend strongly on these properties. Here we specialize on a 2-regular graph.

Let us recall that a general connected graph G is defined by N connected vertices V and M edges E . The graph is undirected if the edges are not oriented. Classical random walks on graphs have been largely explored and their comprehension appears fundamental in the understanding of Markov chains in computational models. As in ergodic Markov chain, CRWs on graph tend to a limiting distribution and this distribution is unique. The number of steps, that the system performs before its distribution becomes infinitesimally close to its limiting distribution, is called mixing time. More formally, if P (t ) is the probability distribution of a CRW on a connected undirected graph, the mixing time ø ≤ is defined as:

ø ≤ = min t {t |t ∏ T ) ||P u (t ) °º|| < ≤} (4.14)
In other words ø ≤ represents the first time when the distance between the stationary distribution º and P (t ) is less than ≤. Not surprisingly, quantum walks do not reach any stationary distribution because the unitarity prevents the walk to reach a steady state.

However, it is proved that the temporal average of the probability distribution does converge to a limiting distribution º § under special assumptions2 . The mixing time is therefore defined as follows:

M ≤ = min{T |8t ∏ T ) ||P QW (t ) °º § || < ≤} (4.15)
In the particular case of N -cycle, with N odd, the characteristic mixing time of an homogeneous QW grows linearly, N log N , in contrast to O(N 2 ) in CRW.

A Nonlinear Quantum Walk-like model on N-cycle

In order to introduce a NQW-like model on N-cycle, let us maintain that the formal definition of QW, introduced in the first chapter, extends simply on cycle considering that now the canonical basis of position |mi 2 Z N .

Let us assume that the unitary step operator U depends, in some way to be defined, on the wave vector itself:

U [™] = T B[™] (4.16)
where T is the coin state dependent translation operator and B is the coin operator. The dependence of the coin operator on the state vector can take several forms. For instance in Shikano et al. [25] the quantum coin depends at each space time point of the lattice on a sitedependent rate function which includes the nearest-neighbor interactions:

g j ,m = |√ L j ,m°1 | + i |√ R j ,m+1 |, (4.17)
and the quantum coin operator transcribes in:

B[g j ,m ] = 0 @ g j ,m °q1 °|g j ,m | 2 q 1 °|g j ,m | 2 g § j ,m 1
A . (4.18)

Publication: "Nonlinear Optical Galton Board: thermalization and continuous limit" 105

This feed-forward DTQW leads, on an unrestricted line, to an anomalous diffusion behavior and a q-Gaussian limiting distribution. Analytical computations on the associated Markov model show a good agreement with numerical data. This model has been then generalized on the N-cycle by (Di Molfetta et al., 2014) 3 in the following way:

B[g j ,m ] = 0 @ g AE,h j ,m °q1 °|g AE,h j ,m | 2 q 1 °|g AE,h j ,m | 2 g AE,h § j ,m 1 A (4.19)
where the new rate function reads

g AE,h j ,m = |√ L j ,m°h | AE + i |√ R j ,m+h | AE .
The resulting behavior appears very complex and depends strongly on the parameter (AE, h) (data not shown). Let us remark that the parameter h amplifies the nonlocal strength of the interaction. We have shown numerically that the behavior of these new walks, even in the case presented by Shikano et al. [25], but on the cycle, displays a chaotic dynamic and tends in long time to a uniform probability distribution (data not shown). Such a behavior is not typical for this model, but seems to share the main features with another NQW-like model introduced by Navarrete-Benlloch et al. [19]. The proposed scheme describes a classical light radiation that acquires intensity-dependent phase while traveling in a nonlinear medium as, for instance, in Kerr medium (like an optical fiber). In QWs formalism this can be represented by the following quantum coin:

B[√ L,R j ,m ] = 1 p 2 √ exp (i g |√ L | 2 ) exp (i g |√ R | 2 ) exp (i g |√ L | 2 ) °exp (i g |√ R | 2 ) ! (4.20)
The authors in Navarrete-Benlloch et al. [19] have found that the dynamics on a one-dimensional finite lattice appears very complex in the presence of formations of soliton-like structures, inelastic shocks and chaotic behaviors. In the next section we have reproduced exactly the same model but on N-cycle. We have shown numerically and analytically, at least in the continuous limit, that the walk tends to a limiting uniform distribution and thermalizes reaching the absolute equilibrium of the system.

Publication: "Nonlinear Optical Galton Board: thermalization and continuous limit"

One of the main ways to reproduce the QW dynamics is suggested by the seminal work of Sir Francis Galton in 1877 on random walks. In order to study the classical random walk he proposed a simple device stylized in Fig. 4.1, that is now called Galton Board or Galton Quincunx.

As we can see in Fig. 4.1 balls roll down on a vertical board and are scattered in a regular grid by pins, reproducing a random walk and leading to a binomial distribution. In spite of its simplicity this model is still an inspiring source for researchers [5].

Now, let us imagine that waves replace balls: the main striking difference is that waves, in contrast to balls, do not run on a single path but traverse all possible paths simultaneously. Bouwmeester et al. [4] implement an optical version in a linear resonator, in which the frequency levels mimic the rolling balls, and the birefringent crystals play the same role of the scattering pins in the classical version of the quincunx. The main result of this experience, then realized by Knight et al. [10] in an optical bidirectional ring, was that the final probability distribution appeared very different from the Gaussian shape of the classical Galton board and presented the same features of standard QWs distribution. Indeed, QWs can be classically reproduced and optically implemented; in particular the Optical Galton Board (OGB) represents a natural way in which, the quantum probability amplitude for each component of the spinor, is replaced with the spectral intensities of the radiation. passes through a beam splitter (BS) and it splits in two different paths, upper and lower. In the lower path the frequency of the light is increased of ¢! by a tuned electro-optical modulator. Then a system of mirrors in the cavity reflects back to the starting point both the components of radiation, and the OGB scheme is implemented again. Note that the BS here plays the role of the quantum coin B and the tuned electro-optical modulator shifts the frequency similarly to the shift operator and T translates the QWs on the physical lattice.

Nonlinear Optical Galton Board (NLOGB) was instead introduced by Navarrete-Benlloch et al. [19] and analyzed in QWs formalism. Differently from the linear case, in NLOGB the radiation acquires phase because of the nonlinearity of the medium. The scheme presented above remains exactly the same but it is possible to add optical fibers, for instance, before the modulators.

We present, in the next section, the NLOGB in long-term regime and we perform analytical computations to derive the formal continuous limit. We then provide numerical calculations to prove that the system thermalizes and demonstrate that the solution converges to a Gibbs distribution.

Nonlinear Optical Galton Board: thermalization and continuous limit

Giuseppe Di Molfetta and Fabrice Debbasch The nonlinear optical Galton board (NLOGB), a quantum walk like (but nonlinear) discrete time quantum automaton, is shown to admit a complex evolution leading to long time thermalized states. The continuous limit of the Galton Board is derived and shown to be a nonlinear Dirac equation (NLDE). The (Galerkin truncated) NLDE evolution is shown to thermalize toward states qualitatively similar to those of the NLOGB. The NLDE conserved quantities are derived and used to construct a stochastic di erential equation converging to grand canonical distributions that are shown to reproduce the (micro canonical) NLDE thermalized statistics. Both the NLOGB and the Galerkin-truncated NLDE are thus demonstrated to exhibit spontaneous thermalization.

I. INTRODUCTION

At the fundamental level, quantum theory is linear. Yet, non linear models are often useful to take into account interaction in an e ective manner. Two examples are the so-called nonlinear optical Galton board (NLOGB) [1] and the non-linear wave equations describing the dynamics of Bose-Einstein condensates (BEC). Though the NLOGB is discrete and wave equations are by definition continuous, these models have much in common. Indeed, the NLOGB is essentially a non-linear quantum walk (QW), and the formal continuous limits of linear QWs are wave equations [2][3][4][5]. Typical such wave equations are the Dirac or the Schrödinger equation which non-linear version, called the Gross-Pitaevskii equation (GPE), is used to model BEC [6]. Also, QW descriptions of BEC have been proposed in [7,8].

Finally, numerical solutions of continuous wave equations are actually solutions of discrete systems approximating the continuous equations.

The NLOGB can be seen as a discrete model of nonlinear waves similar to those which propagate in BEC. One can therefore expect the NLOGB to display properties similar to those of the standard nonlinear model of BEC: the GPE. One such property which has until now never been explored on the NLOGB nor, more generally, in the context of QWs and quantum automata, is the so-called spontaneous thermalization.

In the context of (nonlinear) BEC, microcanonical equilibrium states are well-known to result from long-time integration of the so-called truncated (or Galerkin-projected) Gross-Pitaevskii equation (GPE) and involve a condensation mechanism [6,[9][10][11]. Furthermore, such thermalization is also known to happen in discretized (rather than spectrally-truncated) GPE [12]. Classical Galerkin-truncated systems have been studied since the early 50's in fluid mechanics. In this context, the (time reversible) Euler equation describing spatially-periodic classical ideal fluids is known to admit, when spectrally truncated at wavenumber k max , absolute equilibrium solutions with Gaussian statistics and equipartition of kinetic energy among all Fourier modes [13][14][15][16]. Furthermore, the dynamics of convergence toward equilibrium involves a direct energy cascade toward small-scales [17,18].

The aim of the present work is to study thermalization phenomena in a spatially-periodic version of the NLOGB and relate it to the thermalization of its (Galerkin-truncated) continuous limit.

The paper is organized as follows. Section II is devoted to the definition of the NLOGB model and its numerical solution. The main result of this section is to display and characterize the complex behavior of the log-time regime. Section III is devoted to the behavior of the continuous limit (sect. III A), its conserved quantities (sect.III B) and the long-time behavior and thermalization of its Galerkin-truncated version (sect. III C). Finally section IV is our conclusion. Technical details are given in appendices.

II. NONLINEAR DISCRETE TIME QUANTUM WALK

A. Fundamentals

We consider discrete time spatially periodic quantum walks defined by the following equations: The index m = 0, ..., N ≠ 1 labels points on the circle and the index j oe N labels instants. The complex wave function or spinor =  ≠ b ≠ +  + b + is defined by its two components  ± on a certain timeand space-independent basis (b ≠ , b + ). The parameter g fixes the importance of the non-linearity. For g = 0, equations (1) coincide with the evolution equations of the standard Hadamard walk. The probability

 ≠ j+1,m = 1 Ô 2 [e ig| ≠ j,m+1 | 2  ≠ j,m+1 +e ig| + j,m+1 | 2  + j,m+1 ] (1)  + j+1,m = 1 Ô 2 [e ig| ≠ j,m≠1 | 2  ≠ j,m≠1 ≠ e ig| + j,m≠1 | 2  + j,m≠1 ]
j = q m 1 |Â ≠ j,m | 2 + |Â + j,m | 2 2 = q
m j,m is independent of j i.e. it is conserved by the walk. We will henceforth denote it by .

B. Asymptotic behavior of the DTQWs

As displayed in Fig. 1a, the family of DTQWs defined by equations (1) exhibits a very complex dynamics, much richer than the dynamics of the Hadamard walk shown for comparison in Fig. 1b. Of particular interest is the asymptotic behavior of the family. Let j,m be the probability of finding the quantum walk at time j at point m and let H j (p) be the distribution of the the probability at time j i.e. H j (p)dp is the number of values of the space-coordinate m for which the probability j,m falls between p and p + dp. Direct numerical simulation (DNS) shows that H j (p) tends towards a stationary distribution H ú (p) which depends only on the parameter g and not on the initial condition. Figure The existence of H ú (p) is typical of non-linear chaotic systems. These systems also exhibit a great sensitivity towards initial conditions, and this sensitivity is confirmed by DNS of the NLOGB. Indeed, starting a DNS of the Hadamard walk with a symmetric initial condition delivers a numerical solution which is symmetrical at all times, whereas using the same initial condition in a DNS of the NLOGB delivers a numerical solution which is not symmetric (see Fig. Note that this sensibility on initial conditions does not impact the determination of the stationary asymptotic distribution H ú (p), since this distribution is the same for all initial conditions.

III. NONLINEAR DIRAC EQUATION

A. A non linear Dirac equation as continuous limit of the DTQWs

The asymptotic aspects of the NLOGB dynamics can be understood by investigating the continuous limit of these walks. The method employed is the same as in [2,19,20] and detailed computations are given in Appendix A. The formal continuous limit of the NLQWs read:

3 IˆT ≠ ‡ 3 ˆX ≠ 3ig 4 M( , † ) 4 = 0 (2) with M( , † ) = † M , (3)
M = I + ‡ 2 3 ( 4
)
where

I is the identity, ‡ 2 = A 0 ≠i i 0 B , ‡ 3 = A 1 0 0 ≠1 B ( 5
)
are the second and third Pauli matrices. The continuous limit of the NLOGB is thus described by a non-linear Dirac equation (NLDE). The non-linearity is confined to the mass term, which depends quadratically on the spinor . Note that (spatially Twodimensional) NLDE have also been used to describe experimental BEC on 2D hexagonal lattice [21][22][23].

The NLDE ( 2) is formally equivalent to Nambu-Jona-Lasinio-like equations (NJLE) (Nambu and Jona-Lasinio, [24]) in 1+1 dimension, which describe a non linear interaction between fermions with chiral symmetry. The constant g corresponds to a non linear coupling constant and if g = 0, (2) degenerates into the Weyl equation.

As detailed in Appendix B, the validity of the continuous limit is best confirmed by using Fourier pseudo-spectral methods [25], which are precise and rather easy to implement. In particular, Fig. (5) displays for di erent values of g the relative di erence between the solution of equations ( 1) and ( 2) as a function of the ' parameter which controls the continuous limit.

Fig. (3) shows the typical profile of asymptotic probability density (T, X) and the stationary distribution H(p) of this density, as obtained form a Galerkin-truncated simulation of the NLDE, dealiased in a way that ensures conservation laws in the truncated system, see Appendix B. Both plots are strikingly similar to the corresponding plots presented in Fig. 2.a and 2.b obtained by numerically integrating the NLOGB. In other words, the NLOGB and the Galerkin-truncated NLDE seem to have very similar asymptotic behavior. We will now analyze in detail the asymptotic behavior of the Galerkin-truncated NLDE. We will first identify the conserved currents for the NLDE (Section III B) and then show that the asymptotic statistics Galerkin-truncated NLDE is identical to the so-called grand canonical statistics (III C).

B. Lagrangian formulation and conserved quantities

The NLDE derives from the following Lagrangian density:

L( , † ) = i 2 # " µ (ˆµ ) ≠ (ˆµ )" µ $ ≠ g 2 ! N " 2 (6)
with

N = " 0 + 1 Ô 3 " 5 , ( 7
)
" 0 = ‡ 1 = A 0 1 1 0 B , " 1 = i ‡ 2 , " 5 = i " 0 " 1 , = † " 0 and ˆ0 = ˆT , ˆ1 = ˆX .
There are two conserved currents and these generate three integrals of motion (conserved quantities). The first current is simply the 2-current J µ = " µ associated to the U (1) invariance of the NLDE. The corresponding integral of motion is the total probability = s † dX of finding the fermion somewhere in space.

The other current is associated to the space-time translation invariance of the NLDE and is the stressenergy tensor

T µ‹ ( , † ) = i 2 [ " µ (ˆ‹ ) ≠ (ˆ‹ )" µ ] ≠ ÷ µ‹ L (8)
where ÷ µ‹ = diag(1, ≠1). The associated conserved quantities are the energy E and the momentum P , which are defined by

E[ , † ] = ⁄ T 00 ! (X), † (X) " dX (9)
and

P [ , † ] = ⁄ T 01 ! (X), † (X) " dX (10)
with

T 00 ( , † ) = ≠ i 2 # " 1 (ˆX ) ≠ (ˆX )" 1 $ ≠ g 2 ! N " 2 (11)
and

T 01 ( , † ) = ≠ i 2 # " 0 (ˆX ) ≠ (ˆX )" 0 $ . ( 12
)

C. Thermalization in the Galerkin-truncated NLDE

If one studies the NLDE on the circle, it is natural to write at all times the spinor (T, X) as a spatial Fourier series and to replace the NLDE by an evolution equation obeyed by the time-dependent Fourier coefficients ˆ (T, k). In performing a Galerkin truncation [26], one retains only a finite number of these coefficients as dynamical variables, say ˆ (T, k) with k = ≠ N 2 , ..., N 2 ≠1, and replaces the exact NLDE dynamics by a new dynamics which, at small k, approximates at least formally the original NLDE dynamics. By Fourier transforming the ˆ (T, k), k = ≠ N 2 , ..., N 2 ≠ 1, back to original physical space (i.e. the circle), one obtains a set of N spinors m (T ), m = 0..., N ≠ 1, which are to be interpreted as the values (T, X m ) taken by the spinor field (T, X) at point X m = 2fim N (see Appendix B). The spinors (T, X m ) are on the same footing as the ˆ (T, k), k = ≠ N 2 , ..., N 2 ≠ 1, and can be viewed as the dynamical variables of the Galerkintruncated NLDE. We now denote by ˜ (T ) the collection { m (T ) = (T, X m ), m = 0, ...N ≠ 1}. All integrals over space of quantities involving the Dirac field can be replaced by Riemann sums. Thus, the total probability, the energy and the momentum can now be viewed as functions of the collection ! ˜ (T ), ˜ ú (T ) " . These functions will still be denoted by , P and E and are conserves by the Galerkin truncated dynamics, see Appendix B.

We now introduce two Lagrange multipliers µ and ◊, define H ◊µ = E ≠ ◊P ≠ µ and consider the following stochastic di erential equations:

d m dT = ≠ ˆH◊µ ˆ ú m + ÷ › m (T ) d ú m dT = ≠ ˆH◊µ ˆ m + ÷ › ú m (T ) ( 13
)
where ÷ is a real coe cient and the › m (T )'s are complex independent Gaussian white noises [27] with correlation functions The density f of this stochastic process obeys the exact Fokker-Planck equation [10,28,29]

< › m (T )› ú n (T Õ ) >= " mn "(T ≠ T Õ ). (14)
ˆtf = ÿ m ; ≠ ˆ m ( ˆH ˆ ú m f ) ≠ ˆ ú m ( ˆH ˆ m f )+ ÷ 2 ˆ2f ˆ m ú m < . ( 15
)
The stationary solution f ÷◊µ of this equation is of the form

f ÷◊µ = 1 Z ÷◊µ exp 3 ≠ 2 ÷ 2 H ◊µ 4 (16)
which is the so-called grand canonical distribution with inverse temperature 2/÷ 2 . We have simulated the stochastic equations ( 13) to obtain numerically the stationary distribution f ÷◊µ and compare it with the asymptotic long-time distribution of the Galerkin-truncated NLDE. Typical results are displayed in Fig. 4 and confirm that the Galerkin-truncated NLDE and the stochastic equations (13 are described by very similar distributions (see Appendix B).

IV. DISCUSSION

We have considered the NLOGB confined to the circle and we have shown that the continuous limit of this NLOGB is a NLDE identical to the NJLE-model. Pseudo-spectral numerical simulations reveal that the asymptotic behavior of the NLOGB is similar to the asymptotic behavior of the Galerkin-truncated NLDE and we have shown that the associated asymptotic statistics is identical to the grand-canonical statistics. Thus, both the NLOGB and the Galerkin-truncated NLDE exhibit spontaneous thermalization.

Previous work on other non-linear quantum walk [30] suggests that this observed spontaneous asymptotic thermalization is not a particular feature of the systems studied in this article, but will also be encountered in other non linear quantum walks, whatever the dimensions of the underlying physical space or of the coin space may be. It is obvious that quantum walks which thermalize will explore space in a very di erent manner from walks which do not thermalize, and their importance for quantum computing should certainly be explored in depth. In a di erent direction, it would be interesting to exhibit and analyze spontaneous thermalization in QWs couple to synthetic gauge fields [2,20,31]. , for di ents '=2fi/n, (from right to left) n=2 6 ,2 7 ,2 8 ,2 9 . be evaluated in only N log(N ) operations by using Fast Fourier Transforms (FFTs). Spatial derivatives of fields are evaluated in spectral space by multiplying by ik and products are evaluated in physical space. The original QW equations can also be simply cast in this setting, as the translation operator m ae m±1 is represented in Fourier space by ˆ k ae ˆ k exp (±ik2fi/N ). In this setting, the continuous limit is automatically taken when N is increased. As we can observe in Fig. (5) the relative di erence scales as expected as ' for di erent values of Ê. However the pseudo-spectral code solving the NLPDEs generates a problem called aliasing [32]. In general the fields needs to be de-aliased by proper spectral truncation. Here, we used the so-called 2/3-rule in all our numerical schemes in the same way as done in reference [10]. De-aliasing is fundamentally important to preserve the conservation of the Galerkin truncated non linear dynamics as we can observe in Fig. (6). Indeed, although it is straightforward to show that Eq.( 2) can be written

ˆT m = ≠i ˆE ˆ ú m (B1) ˆT ú m = i ˆE ˆ m (B2)
and thus formally conserves the energy, it can be shown that exact conservation requires proper dealiasing (see appendix of ref. [10]).

As displayed in Fig. 4, the statistical distributions generated by the NLDE dynamics Eq.( 2) and by the stochastic equations ( 13) are really close and this can be justified on very general grounds.

First, by construction, the stochastic equations ( 13) generate the grand canonical distribution ( 16) that is controlled by the inverse temperature 2/÷ 2 and the Lagrange multipliers µ and ◊. On the other hand, as the spectrally-truncated dynamics (2) conserves , P and E, its long time behavior should be described by the so-called micro canonical distribution

f ≥ "(E ≠ E in )"( ≠ in )"(P ≠ P in ). ( B3
)
that is determined by the values (E in , in ,P in ) of the conserved quantities given by the initial condition in . As is well-known [33], under very general circumstances both grand canonical and micro canonical distribution yield similar statistical results (provided that the 2/÷ 2 and the Lagrange multipliers µ and ◊ have values that correspond to E in , in ,P in ). Fig. 4 indicates that, in this case, both distributions yield identical results for density fluctuations.

Part III

Conclusions and Perspectives

CHAPTER 5

CONCLUSIONS AND PERSPECTIVES

Conclusions

In this thesis we have studied several QWs in discrete space and discrete time with homogeneous and inhomogeneous quantum coin. In the first part, we proved that, under some conditions and employing the method presented for the first time in [5] and [6], these walks coincide with the propagation of a Dirac particle. In Chapter (1) we introduced the simplest case of the homogenous QW and formally derived its continuous limit, which recover the massless Dirac equation in flat (1+1)-spacetime.

In Chapter (2) we extended the previous analysis to the inhomogeneous QWs. In particular in [7] we proposed an explicit construction of DTQWs mimicking the propagation of a fermion in and around a (1+1)-dimensional black hole. This case has been validated by analytical calculation for the specific case of a 4D Schwarzschild non-charged black hole, in Lemaître coordinates, for a particle that starts its motion radially. By Direct Numerical Simulations (DNS) we proved that the agreement is surprisingly valid still when the initial condition is no longer smooth in respect to the lattice parameter. We have also proved (data not shown in this thesis) that the spacetime curvature, generated by the inhomogeneities of the QW, is still visible within a very small amount of time steps (for instance in 8 and 16 time steps).

The main goal obtained in [8] is that a large class of IDTQWs admit remarkably an exact discrete gauge invariance, which at the continuous limit coincides with the standard U (1) Maxwell gauge invariance. In continuous limit these walks recover the Dirac equation of a fermion propagating in a gauge electric field. We validate this result with accurate pseudospectral numerical simulations in the specific case of a constant and static electric field.

In Chapter (3) we proposed two families of DTQWs as models mimicking quantum transport in time-dependent random electric and gravitational field. These systems, in long time behavior or large decoherence rate, loose spin and spatial coherence and the probability density tends to a limiting Gaussian distribution. However, in the transient regime, the QWs present two characteristic density profiles with a cusp in the case of the random electric field and a dip in the case of gravitational field. We should remark also that the loss of coherence is largely slower in the case of a random electric field (data not shown in this thesis). For both kinds of system we explicitly computed the diffusion coefficient and the analytical expression of the second moment in physical space. Let us finally regard that because of the spatial translation invariance of the random fields, the averaged transport operator does not couple the different Fourier modes.

We were the first to prove to the best of our knowledge (and we thus reported in Chapter (4)) that a nonlinear QW-like model thermalizes and reaches an absolute equilibrium in long time behavior. To prove this, we adopted the NLGOB introduced by Navarrete-Benlloch et al. [11]. Several new results have been presented at this stage: (i) the continuous limit of this discrete model corresponds to the NJL equations introduced as dynamical model of ele-mentary particles based on analogy with superconductivity; (ii) this system relaxes toward the equilibrium with a rich and interesting transient. As in the Galerkin-truncated Euler-Voigt-AE equation, (Di Molfetta et al. [9]), in the Galerkin truncated NJL, the thermalized small scales act as thermostats generating a pseudo-dissipation at large scales. (iii)The same asymptotic behavior is observed in the NLGOB. Let us acknowledge that the algorithm introduced by Krstulovic et al. [10] to generate absolute equilibrium of spectrally truncated compressible flows in the hydrodynamics framework, has been adapted to the Galerkin truncated NJL. All these results are validated by precise numerical pseudo-spectral simulations for the spatial discretization and by a 4th order Runge-Kutta method for the time stepping.

Perspectives

All the works presented in this thesis should certainly be continued in several directions. One should first extend the main results of the first part of this thesis to QWs defined in physical space of higher dimension and/or defined in a higher-dimensional Hilbert space. Our group has started to work on the former direction, showing that certain QWs on 2D square lattices describe the transport of two-component spin 1/2 fermions coupled to arbitrary electromagnetic fields. In particular, Landau levels can be defined for these QWs (Arnault and Debbasch [1]). The next goal will be to exhibit families of 2D QWs on square lattices, that can be interpreted as the transport of two-component spin 1/2 fermions in arbitrary gravitational fields.

QWs defined in a higher-dimensional Hilbert space have been instead developed in a preliminary work carried out by Di Molfetta and Debbasch. We showed that several 1D QWs with more than two components (but still defined on square lattices) can be interpreted as quantum transport in non-abelian synthetic gauge fields. In particular we have proved 1 that QWs living in SU (2) ≠ SU (n) can simulate SU (n) Yang-Mills theory in flat and curved spacetime.

Another interesting perspective should be the extension of all the results above for QWs propagating on graph. In fact, on one hand graphs possess intrinsic geometrical properties, on the other QWs can be interpreted as fermions interacting with gauge fields, including gravitational fields. Now, gauge theories are based on geometry; for example, gravitational fields arise from the curvature of space-time. It is logical at this point to wonder what effects the interplay between the intrinsic geometry of a graph and the geometry of gauge fields has on coherent quantum transport.

Let us notice that a recent method of spacetime grouping introduced by Arrighi et al. [2] in the QWs framework may allow us to relate QWs to the interaction graph of a specific QCA . By changing the local interactions between the cells and the spatial distributions of the cells, we should obtain different emergent dynamics. Therefore, there could be two interesting extensions related to our works: (i) the former is to generally explore what possible QFTs can be described by different grouping operations; (ii) the latter, investigating the classical limit of these configurations by, for instance, a coarse graining procedure.

Perspectives 123

Chapter (3) should be instead expanded considering other types of noises, such as positionand time-dependent noise, and investigating if they could permit a closer connection with some specific class of relativistic stochastic processes (Debbasch et al. [3]). Moreover, previous analysis carried out by G. Di Molfetta and F. Debbasch have related decoherent QWs dynamics to classical Bounded Velocity Stochastic Processes [4] and this subject clearly demands further analysis.

In Chapter (4) we have just presented some preliminary new results on a specific nonlinear QW-like model and its continuous limit. We should first extend the same analysis to the model presented by Shikano et al. [12], notably in the generalized form proposed by Di Molfetta. Then, we have to move our analysis further by (i) proving that these quantum walks do thermalize and finally (ii) analyzing the thermalization in terms of gauge fields.

Let us recall, in the end, that we need to investigate all these QWs using some experimental setups in order to definitely validate our theoretical results. Note that eq. B.4 is a solution of the associated Liouville equation of the system B.2, therefore when is given by p.d.f. B.4, the velocity field remains in a statistically stationary state. More in detail, Cichowlas et al. [2] found that in long time behavior a clear spontaneous scale separation appears in energy spectrum. Large and small scales display a different behavior: at small scales the thermalized modes of energy spectrum rise as k 2 and act by a pseudodissipative effect (as an effective thermal bath) on large scale. The progressive thermalization extends to larges scales till a fully repartition of the energy. Note that the large scales behavior, before fully equilibration, is compatible with a K41 scaling, that is the typical Kolmogorov spectrum in case of viscous turbulent hydrodynamics ( [4,3]). In this temporal evolution of spectrum the wave number k th , defined as the mode triggering the process of thermalization plays an important role and quantifies the degree of thermalization of the whole system.

B.2 Eddy-damped quasi-normal Markovian theory (EDQNM)

The eddy-damped quasi-normal Markovian theory was introduced by Orszag [13] in order to treat analytically turbulence and in particular it has been used to study thermalization from Kolmogorov k °5/3 scaling to absolute equilibria. It was also largely used in the past in numerical simulations as it allows a large scale separation with a relatively low computational cost.

This theory takes into account some important assumptions to determine the evolution of the energy spectrum: (i) the isotropy of the flow, (ii) the Gaussian statistical distribution of the Fourier modes in the statistics of the flow. Note that the definition of quasi-normality derived from the fact that this assumption is relaxed on the third moment that does not vanish (as they would be for zero mean Gaussian variable). The last assumption is (iii) Markovianity, that is reached dropping out the dependence of the past introducing an eddy-damped term ¥ k which takes into consideration the nonlinear interaction. The resulting equation is:

£ @ t + 2∫k 2 § E (k, t ) = ZZ ¢ £ kpq (x y + z 3 )[k 2 pE(p, t )E (q, t )p 3 E (q, t )E (k, t )] d pd q pq (B.5)
where the characteristic time £ kpq is defined by:

£ kpq = 1 °exp(°(¥ k + ¥ p + ¥ q )t ) ¥ k + ¥ p + ¥ q (B.6)
In the EDQNM theory the eddy damped ¥ k is defined as:

¥ k = ∫k 2 + ∏ s Z k 0 s 2 E (s, t )d s (B.7)
This theory introduced by Orszag [13] has been then used by Bos and Bertoglio [1] in the inviscid case (∫ = 0) and introducing an effectivity viscosity, the eddy viscosity, used to explain the thermalization behavior and the dependence of E on the Fourier modes.

B.3 Self-truncation

As we have seen, the effect of Galerkin truncated Euler equation is responsible for the energy cascade toward small-scales and contains transient that mimics the irreversible viscous effects produced by the microword of small scale partially thermalized Fourier modes. Recently Krstulovic and Brachet [8] [7] have introduced in the case of superfluids the Galerkin-B.4. Publication A1: "Self-truncation and scaling in Euler-Voigt-AE and related fluid models" 137 truncated Gross Pitaevskii equation. In that case absolute equilibria exist and energy cascade is accompanied by an energy transfer to large scale. Furthemore it has been observed that partial thermalization is independent of the truncation wave number and that can be observed defining a self-truncation wave number k st . This can be explained by the fact that an increase of dispersion (that is controlled by a physical parameter called "healing length") produces a slowdown of the energy transfer at small scales inducing a bottleneck and the emergence of a self-truncation.

B.4 Publication A1: "Self-truncation and scaling in Euler-Voigt-AE and related fluid models"

In the following publication we proposed a generalization of the Euler-Voigt-AE introduced by Larios et al. [11] in a different context. This model was characterized by a penalization of the small scales and recovered the standard Euler equation in the limit AE = 0. The main aim of the article was to show self-truncation in inviscid generalized Euler-Voigt-AE model and investigated all the temporal regimes. We introduced Leith and EDQNM model to describe analytically the self-truncation and we showed that in long time behavior the spectrum is self-similar and k s t / t ¥ . The exponent ¥ was determined analytically.

Self-truncation and scaling in Euler-Voigt-↵ and related fluid models The values of the self-similar exponents of these models are found analytically.

I. INTRODUCTION

Classical Galerkin-truncated systems have been studied since the early 50's in fluid mechanics. In this context, the (time reversible) Euler equation describing spatiallyperiodic classical ideal fluids is known to admit, when spectrally truncated at wavenumber k max , absolute equilibrium solutions with Gaussian statistics and equipartition of kinetic energy among all Fourier modes [1][2][3][4]. Furthermore, the dynamics of convergence toward equilibrium involves a direct energy cascade toward smallscales and contains (long-lasting) transient that mimic (irreversible) viscous e↵ects that are produced by the "gas" of high-wavenumber partially-thermalized Fourier modes generating (pseudo) dissipative e↵ects [5][6][7].

In the the case of superfluids, the relevent equation is the so-called truncated (or Galerkin-projected) Gross-Pitaevskii equation (TGPE). In the TGPE case, absolute equilibrium can also be obtained by a direct energy cascade, in a way similar to that of the truncated Euler case, with final thermalization accompanied by vortex annihilation. Furthermore, increasing the amount of dispersion produces a slowdown of the energy transfer at small-scales inducing a bottleneck and a partial thermalization that is independent of the truncation wavenumber and takes place below a 'self-truncation' wavenumber k st (t) that is observed to slowly increase with time [8][9][10].

The purpose of the present paper is is to find and study such self-truncation phenomena in the simpler context of classical hydrodynamics of an ideal fluid. This is obtained by using equation of motion of the Euler type. To wit, we study here a simple generalization of the standard 3D Euler-Voigt-↵ model [11,12]. Compared to the Euler equations, this conservative model penalizes the formation of small scales. We show that this penalization is enough to produce a self truncation regime. Our main findings are that the self-truncation regime of this generalized model reproduces the behavior of the truncated Euler equation [5]. The long-time behavior of the energy spectrum appears to be self-similar.

To understand this self-similarity we further introduce two di↵erent models that are ↵-Voigt versions of the Eddy-Damped Quasi-Normal Markovian (EDQNM) model and the Leith model, respectively. Both models are shown to present behaviors that are similar to that of the Euler-Voigt-↵ model. The relative simplicity of these models allows us to determine the analytical values of the self-similar exponents.

The paper is organized as follows. Section II is devoted to our generalized model. Basic definitions of the Euler-Voigt-↵ model are given in section II A. Numerical methods and performed computations are detailed in section II B. Our results on the self truncation regime are described in section II C. The long-time behavior is studied in section II D. Related theoretical models are presented in section III. Section III A is devoted to the ↵V-EDQNM model and section III B to the ↵V-Leith model. Finally, our main results are summarized in sec-tion IV where we give our conclusions.

II. EULER-VOIGT-↵ MODEL

A. Definition of the model

The standard 3D Euler-Voigt-↵ model [11,12] is a partial di↵erential equation for the 3D velocity field u(x, y, z, t) that explicitly reads

(1 ↵ 2 r 2 ) @u @t = (u • r)u rp r • u = 0. (1)
The operator ↵ 2 r 2 @ @t (as we will see later) suppresses the formation of scales smaller than ↵. The associated wave number to this scales is denoted k ↵ = ↵ 1 . We refer to the operator in Eq.1 as the ↵-term.

Let us define the generalized 3D Euler-Voigt-↵ model:

(1 + ( ↵ 2 r 2 ) 2 ) @u @t = (u • r)u rp r • u = 0, (2)
where the power is an even integer. We refer to as the penalization exponent as its increase enhances the suppression of small scale generation. When ↵ = 0, the generalised model Eq. ( 2) reduces to the standard 3D incompressible Euler equations

@u @t + u • ru = rp, r • u = 0. (3)
We consider here spatially-periodic solutions defined in the domain ⌦ = [0, 2⇡] 3 . The kinetic energy spectrum E(k, t) associated to (3) is defined as the sum over spherical shells

E 0 (k, t) = 1 2 X k2Z 3 k 1/2<|k|<k+1/2 |b u(k, t)| 2 , (4)
and the energy

E 0 = 1 2(2⇡) 3 Z ⌦ |u(x, t)| 2 d 3 x = 1 2 X k2Z 3 |b u(k, t)| 2 ,
is independent of time when u satisfies the 3D Euler equations (3). The conserved energy associated to the generalized Euler-Voigt-↵ model ( 2) is straightforward to obtain and reads in physical space

E ↵ = 1 2(2⇡) 3 Z ⌦ u • [1 + ( ↵ 2 r 2 ) 2 ]u d 3 x,
and in spectral space

E ↵ = 1 2 X k2Z 3 [1 + (↵k) ]|b u(k, t)| 2 . (5)
Consequently, the generalised energy spectrum is defined as

E ↵ (k, t) = 1 2 X k2Z 3 k 1/2<|k|<k+1/2 [1 + (↵k) ]|b u(k, t)| 2 . (6)
In the following, we refer to E ↵ (k, t) as the energy spectrum and E(k, t) as the kinetic energy spectrum. Equations (2) also conserve the generalised helicity

H ↵ = 1 2 X [1 + (↵k) ] b u(k, t) • b !( k, t). (7)
In this work we only consider flows with H ↵ = 0. Let us remark that the di↵erential operator multiplying the r.h.s. of our generalized 3D Euler-Voigt-↵ model ( 2) can be written in Fourier space as 1 + (↵k) = 1 + (k/k ↵ ) . The formal limit ! 1 of Eq. ( 2) thus corresponds to a standard spherical Galerkin truncation (û(k) = 0 for |k| > k max ) of the Euler equation ( 3) at

k max = k ↵ .
It is well known that the truncated Euler equation admits statistically stationary solutions given by the microcanonical distribution determined by the invariants [1,3]. These solutions are the so called absolute equilibrium and lead to equipartition of energy among Fourier modes. The Euler-Voigt-↵ model considered as a truncated system, also admits absolute equilibrium solutions. When fully thermalized, Fourier modes can be described by the canonical Gibbs distribution û(k) ⇠ Z 1 exp { E ↵ [û]}, where Z is the partition function [33]. As the invariant energy is quadratic, the absolute equilibrium is Gaussian and the Fourier modes are independent. This leads to the spectra

E ↵ (k) = 3E ↵ k 2 k 3 max , or E(k) = 3E ↵ k 2 k 3 max (1 + ↵ k ) . (8)
The large scale behavior of the kinetic energy spectrum thus depends on the value of the penalization exponent as E(k) ⇠ k 2 . Therefore, in thermal equilibrium the small scales of u (i.e. k k ↵ ) are penalized when > 0. Note that, in Fourier space, the di↵erential operator in (2) can be defined for real values of 0. Choosing = 11/3 yields an absolute equilibrium E(k) ⇠ k 5/3 and thus a fully thermalized field following Kolmogorov scaling. In the same vein, choosing in two-dimensions = 2/3 also yields Kolmogorov scaling, this time for the equipartition of enstrophy. This can represent an interesting alternative to the fractal decimation method that was used in reference [14].

B. Numerical method

The generalized 3D Euler-Voigt-↵ equations (2) are solved numerically using standard [15] pseudo-spectral methods with resolution N . Time marching is performed using a second-order Runge-Kutta scheme and the solutions are spherically dealiased by suppressing, at each time step, the modes for which the wave-vector exceeds two-thirds of the maximum wave-number N/2 (thus a 2048 3 run is truncated at |k| > k max = 682). We consider here solutions of Eq. ( 2) that correspond to the so-called Taylor-Green (TG) [16] (2⇡-periodic) initial data u(x, y, z, 0) = u TG (x, y, z), with u TG = (sin(x) cos(y) cos(z), cos(x) sin(y) cos(z), 0). ( 9) The simulations reported in this paper were performed using a special purpose symmetric parallel code developed from that described in [17][18][19][20]. The code uses the symmetries of the Taylor-Green initial data to speed-up computations and optimize memory usage. The workload for a timestep is (roughly) twice that of a general periodic code running at a quarter of the resolution. Specifically, at a given computational cost, the ratio of the largest to the smallest scale available to a computation with enforced Taylor-Green symmetries is enhanced by a factor of 4 in linear resolution. This leads to a factor of 32 savings in total computational time and memory usage. The code is based on FFTW and a hybrid MPI-OpenMP scheme derived from that described in [21]. At resolution 2048 3 we used 512 MPI processes, each process spawning 8 OpenMP threads.

When compared with standard Euler equation (3) runs that were performed in reference [19], the only computational advantage of the the generalized 3D Euler-Voigt-↵ model (2) stems from the much weaker Courant-Friedrichs-Lewy (CFL) condition on the time step t which is conditioned by k ↵ = ↵ 1 rather than k max . We have performed a number of high-resolution runs that are summarised in Table I was enough to insure both stability and energy conservation up to 0.3% in the worst case.

C. Self truncation

A first indication on the dynamics of the generalized Euler-Voigt-↵ model (2) with TaylorGreen initial data ( 9) is given by the behavior of the Energy spectra for a run at resolution 1024 3 with = 4 and k ↵ = 80 (run 6 Table I) that is displayed in Fig. 1. Di↵erent regimes are clearly observed. First the energy is transferred toward small scales as in the standard Euler equation evolution (t  5). Then the energy reaches the wavenumber k ↵ = 80 and the ↵-term in (2) starts to suppress the energy transfer for k > k ↵ . As a consequence, the energy piles up around this wavenumber, similarly to the truncated Euler case with a cut-o↵ ⇠ k ↵ . Note that a compatible Kolmogorov k 5/3 scaling is observed at large scales (t = 15), followed by a partially thermalized zone in k 2 extending up to k st . The wavenumber k st (t) then slowly growth from its initial value k ↵ until it eventually reaches the simulation cut-o↵ k max . For t ! 1 the system fully thermalises independently of the parameters (data not shown) and the spectrum is then described by the absolute equilibrium (8). Note that for k > k st the energy quickly decays and the partial thermalization regime is thus independent of the simulation cut-o↵ k max . We thus refer to k st as the self-truncation wavenumber. A further discussion and justification for its name will be given later (see below, paragraph following Eq.10).

Figure 2 displays the spectra of simulations at resolution 2048 3 for di↵erent values of and k ↵ taken at time t = 10.2. The self-truncation is apparent for the smaller k ↵ at all values of and, for larger 4, at all values of k ↵ . The three ranges delimited by the thermalization wavenumber k th and the self-truncation wavenumber k st are clearly visible (see Fig. high resolution. Formally, these zones correspond to the Kolomogorov regime (k ⌧ k th ), the thermalization range ( k th ⌧ k ⌧ k st ) and the exponential energy decay (k k st ). Note that, at this intermediate time (t = 10.2), the = 6 case appears to be already behaving somewhat like what is expected of the = 1 limit. Indeed, Fig. 2(c) is reminiscent of the previously studied truncated-Euler case (see Fig. 1 of reference [5]). Of course, in the truncated Euler case no third decreasing zone exists as, when = 1, E ↵ (k) = 0 for k > k ↵ . In order to make the comparison with the = 1 limit more quantitative, following [5], we have computed the thermalized energy E th (t) = P k th (t)<k E(k, t) and the e↵ective dissipation "(t) = dE th (t) dt . The time evolutions of k th , E th and " are presented on Fig. 3. The agreement with the truncated Euler data appears to be good (compare with Figs 2 and 4 in reference [5]). Note that these quantities do not appreciably depend on the value of .

In references [5,22], an e↵ective generalized Navier-Stokes model for the dissipative dynamics of modes k close to k th (t) was suggested for the original Euler case with fixed truncation at k = k max . In this case, the e↵ective viscosity of the model was given by ⌫ e↵ = p E th /k max . If we assume that the generalized largecase behaves similarly to the Euler case truncated at k st , we find that the dissipative wavenumber k d should be given by

k d ⇠ " 1/4 ( p E th /k st ) 3/4 (10)
The consistency of this estimation of the e↵ective dissipation with the results displayed in Fig. 3a-c 2) is seen to emulate the dynamics of the Euler equation, spectrally-truncated at k max = k st . For this reason, we call the regime where this behavior takes place the self-truncation regime.

In the self-truncation regime the energy spectrum E ↵ (k) sharply decreases for for k > k st (see Fig. 2). In order to quantify this behavior we have used the so-called analyticity strip (AS) method [23]. Let us recall that the idea is to monitor the 'width of the analyticity strip' ( 0) as a function of time, e↵ectively measuring a 'distance to the singularity' [24]. Using spectral methods [15], (t) is obtained directly from the high-wavenumber exponential fall o↵ of the spatial Fourier transform of the solution [25]. The AS method has been used to numerically study putative Euler singularities (see e.g. reference [19] for implementation details). The common procedure 11): Horizontal lines correspond to kmax = 2 (dashed black line) and k↵ = 2 (dot-dashed black line). Exponential law (dot black line) = 2.7 exp( t/0.56) from reference [19].

is to perform a least-square fit at each time t on the logarithm of the energy spectrum E ↵ (k, t), using the functional form ln E ↵ (k, t) = ln C(t) n(t) ln k 2k (t). (11) Energy spectra are fitted on the intervals 2 < k < k ⇤ for t < 2 and on the interval 300 < k < min(k ⇤ , k max ) for t > 2, with k ⇤ = inf {k | E(k) < 10 32 } denotes the beginning of round o↵ noise. The fit presented in Fig. 4.a in good agreement with the data. The time evolution of the fit parameter is displayed in Fig. 4.b, where it is compared with the exponential in time law (t) = 2.7e t/0.56 followed by the Euler equation (see reference [19]). The horizontal lines show the value of the length 2/k max (dashed) and 2/k ↵ (dotted-dashed). It is apparent that the model follows the Euler dynamics as long as k ↵ 1. The measure of the fit parameters is reliable as long as (t) remains larger than a few mesh sizes ( k max 1), a condition required for the smallest scales to be accurately resolved and spectral convergence ensured. Thus the dimensionless quantity k max is a measure of spectral convergence. Therefore, the self-truncation solutions are solution of the full partial di↵erential equation (2) and not of the spectrally truncated system.

From the mathematical point of view, Eq.( 2) can be considered as ordinary di↵erential equation as the right hand side can be shown to be a bounded continuos Lipschitz map between Banach spaces. In this case, existence and regularity can be proved using the same techniques than in (finite dimensional) ordinary di↵erential equations [26].

D. Long-time behavior of kst

The self-truncation observed in Fig. 2 is accompanied by a very slow growth of k st until it reaches the simulation cut-o↵ k max . Such a behavior was also observed in the dispersive self-truncation of the truncated Gross-Pitaevskii equation in two [10] and three [8,9] dimensions. In the Gross-Pitaevski case, the self-truncation wavenumber k st was shown to obey a power-law scaling

k st (t) ⇠ t ⌘ , t 1. ( 12
)
Fig. 5 displays the temporal evolution of k st for di↵erent values of and k ↵ . The wavenumber k st is determined by the weighted average k st = p (5/3) P k E ↵ (k, t)k 2 /E ↵ . With this definition, in the case of the absolute equilibrium (8) we obtain k st = k max . Note that k st is defined for all times, even before the self-truncation starts to take place. It is in fact proportional to the inverse of the Taylor microscale [25].

Figure 5.a compares k st (t) for di↵erent values of the penalization exponent with k ↵ = 20 (lower three curves) and k ↵ = 80 (upper three curves). Note that for the short times when k st < k ↵ the temporal evolution of k st is, as expected, independent of . Furthermore, for the very short times, the penalization of small scales introduced by the ↵-term in the Euler-Voigt-↵ model is negligible and thus k st is also independent of k ↵ as the dynamics is given by the (standard) Euler equation (black curve with crosses).

A behavior compatible with a power law is observed at long times. Note that the exponent seems to depend on the value of thought the data do not allow a clear determination of the exponent ⌘. The power-law behavior is contaminated by the initial dynamic as it is actually expected to have the form k st (t) ⇠ (t t 0 ) ⌘ , where t 0 is the time when self-truncation starts. Simulations where t t 0 and k st ⌧ k max are di cult to reach with the present choice of parameters and resolution.

In order to explore such power-law behavior a series of runs (19 21, see Table I) have been performed in resolution 512 3 and at small value of k ↵ = 4. With this choice of parameter, self-truncation starts at low wavenumbers (⇠ k ↵ ) and thus no Kolmogorov scaling can be observed. However such a choice allows very long temporal integrations (up to t = 2300) and to clearly observe the powerlaw behavior of k st , as apparent in the inset of Fig. 5. The value of the exponent clearly depends on the penalization exponent . The values measured for ⌘ are displayed in Table II. The power-law observed for k st strongly sug- gests to look for self-similar behavior of the energy spectrum E ↵ (k, t), where the only temporal dependence of the spectrum is given by k st (t). A self-similar form of the energy spectrum compatible with the conservation law ( 5) is given by

E ↵ (k, t) = E 0 k st (t) ✓ k k st (t) ◆ (13)
Where E 0 is a constant with dimension of energy. The function (z) is expected to behave as (z) ⇠ z 2 for z ⌧ 1 and exponentially decay for z 1. Figure 6 displays k st (t)E ↵ (k, t) as a function of k/k st (t) where the tendency to converge towards a self-similar distribution is confirmed for long times. The assumption of self-similarity implicitly supposes that the scale ↵ is very small. Such hypothesis is valid only for long times such that k ↵ ⌧ k st (t) ⌧ k max . Discrepancies with the full self-similar form are certainly due to the finite values of the infra-red and ultra-violet cut-o↵ of the simulation. In order to develop further this idea, in the next sections we introduce two theoretical models that allow to obtain both a clear numerical support of self-similarity and an analytic expression for the self-truncation exponent ⌘.

III. THEORETICAL MODELS A. Eddy-damped Quasi-Markovian Euler-Voigt-↵ model

A popular model of turbulence is the so called Eddy-Damped Quasi-Markovian (EDQNM) closure [4]. It was derived in the 60-70's and is based on statistical closure of the velocity correlations in Fourier space plus some adhoc modeling of the dissipative time scales. It furnishes an integro-di↵erential equation for the spectrum E(k, t). EDQNM has been proved to be a powerful theoretical and numerical tool in the las 30 years as it allows to achieve a very large scale separation. It was also shown by Bos and Bertoglio [22] that EDQNM reproduces well the dynamics of the truncated Euler equation, including the k 5/3 and k 2 scalings together with the relaxation to equilibrium. It was also used in [27] to give an analytic prediction of the e↵ective viscosity acting on the large scales of truncated Euler flows.

The extension of EDQNM model to the Euler-Voigt-↵ case is straightforward. First, following Orszag derivation [4] but using equations (2), we directly find @E(k, t)

@t = 1 1 + ↵ k T NL (k, t) (14)
where the nonlinear transfer T NL is modeled as

T NL (k, t) = Z Z 4 ⇥ kpq (xy + z 3 )
 k 2 pE(p, t)E(q, t) 1 + ↵ k p 3 E(q, t)E(k, t) 1 + ↵ p dp dq pq .

In ( 15), 4 represents a strip in (p, q) space such that the three wavevectors k, p, q form a triangle. x, y, z, are the cosine of the angles opposite to k, p, q. ⇥ kpq is a characteristic time defined as

⇥ kpq = 1 exp ( (⌘ k + ⌘ p + ⌘ q )t) ⌘ k + ⌘ p + ⌘ q . ( 16
)
The standard EDQNM equations are recovered by setting ↵ = 0. Note that absolute equilibrium ( 8) is a stationary solution of ( 14) that satisfies detailed balance (no flux solution). The eddy damped inverse time ⌘ k is defined as

⌘ k = 0 v u u u t k Z 0 s 2 E(s, t) 1 + ↵ s ds. (17)
In the limit ↵ ! 0 the standard eddy damped inverse time is recovered [28]. Note that for ! 1, the standard eddy time is also recovered for k < k ↵ = 1/↵ whereas ⌘ k = 0 for k > k ↵ , consistently with the dynamics of the truncated Euler equation, as for k > k max the dynamics is frozen. The constant defines a time scale and we use the standard value = 0.36. The truncation is imposed by omitting all interactions involving waves numbers larger than k max in (15). We refer to (14)(15)(16)(17) as ↵V-EDQNM model.

A number of simulations of the ↵-EDQNM equation has been performed and give a behavior that is similar to that of the DNS of the full Euler-Voigt-↵ model, including comparable time scales. Figure 7.a displays the energy spectrum for = 2, 4, 6 and k ↵ = 10 4 for a simulation with k max = 43918 (and 14.2 points per octave). Fig. 7.a shows that the three zones observed in the Euler-Voigt-↵ model (see Fig. 2), are also apparent in ↵V-EDQNM model spectra, with scaling laws extending for more than two decades. In the same spirit than in the Euler-Voigt-↵ DNS runs with a smaller value of k ↵ = 5 have been performed, allowing a clear determination of the self truncation exponent ⌘. The values are presented on Table II. Analogously to the Euler-Voigt-↵ DNS, we look for a self-similar behavior of the energy spectrum. The collapse is manifest in Fig. 7.b where the self-similar form (13) is displayed for = 4. The ↵V-EDQNM model allows for directly looking for such self similar behavior and find the exponent ⌘ by counting powers. Indeed, introducing the self-similar form and the variable z = k/k st (t) in ( 13) in the eddy damped inverse time (17) we find ⌘ z (t) = E the left and right hand sides of ( 14) we obtain k st k 2 st ⇠ k 3 st . Finally, the scaling k st (t) ⇠ t ⌘ leads to

⌘ = 1 3 1 + (19)
The exponent can be computed using the ansatz (z) = z Let us now introduce another model that shares the same dynamics properties than Euler-Voight-↵. It is a spectral di↵usion model that generalizes the so-called Leith model [29]. The original Leith model it is a phenomenological non-linear (local) spectral di↵usion equation for the energy spectrum that admits the stationary solutions corresponding to an absolute equilibrium and Kolmogorov scaling. When forcing and dissipation is added to the model, a steady state containing mixture of constant flux and thermal equilibrium was observed in [30]. It also known to posses self-similar solutions [31].

The simplest generalization of the Leith model to take into account the ↵-term that conserves the total energy E ↵ (k) is given by

@E @t = 1 (1 + ↵ k ) @F @k , (20)
where F (k) is a spectral (non-linear) flux. Following Leith's original derivation, we assume that the spectral flux is defined in terms of a di↵usion coe cient D(k) and a potential Q(k) such that

F (k) = k 2 D @Q @k (21)
Assuming locality in the flux, by dimensional analysis we obtain

D = k 9/2 m (E/k 2 ) n (1 + ↵ k ) p (22)
Q = k m (E/k 2 ) 3/2 n (1 + ↵ k ) q (23)
The dimensionless coe cient set the global time-scale and n, m, q, p are free parameters to be determined. The first contraint is given by the no-flux solution or the absolute equilibrium (8). Imposing that F (k) = 0 for the absolute equilibrium E ↵ (k) ⇠ k 2 leads to m = 0 and n = 3/2 q. The second contraint is for the Kolmogorovlike solution E ↵ ⇠ k 5/3 . Imposing that for such solutions the flux is given by F = ✏/(1 + ↵ k ) r , we obtain p = 3/2 q r, with r a free parameter. Such a flux can be easily interpreted in the limit of large . Indeed, for k ⌧ k ↵ the ↵ term is negligible, and the flux becomes constant. The Kolmogorov phenomenology is thus recovered. On the other hand for k k ↵ the flux vanishes as expected. Taking into account the previous constraints, we obtain a family of a of di↵usive ↵V-Leith models indexed by the parameter r:

@E ↵ @t = 2q 3 @ @k  k 13/2 (1 + ↵ k ) r @ @k h E 3/2 ↵ k 3 i . (24)
The standard Leith model is recovered by setting = 0 and rescaling the time.

The behavior of the Euler-Voight-↵ is indeed reproduced by the ↵V-Leith model (24). Figure 8.a displays the energy spectrum for di↵erent values of and r = 2. The Kolmogorov k 5/3 , the equilibrium k 2 regimes and the fast decay for large k is manifest. The absence of a dissipative zone is apparent in Fig. 8.a, when compared both with the ↵V-EDQNM case (see Fig. 7.a) and the Euler Voigt-↵ model (see Fig. 2). This is certainly due to the locality in Fourier space of the ↵V-Leith model.

As in the previous models, we look for self-similar solution of (24). Introducing (13) in to (24) we obtain in

This prediction coincides with the one of EDQNM for r = 2 + / (see Eq.( 19)). The self-similarity behavior of E ↵ (k, t) is apparent in Fig. 8.b, where the self-similar for is displayed for = 4 and r = 2. The inset shows the temporal evolution of k st (t) for di↵erent values of . A power-law growth is manifest. The measured values of the exponent ⌘ are presented inTable II in good agreement with the prediction (25).

The self-similar analysis leads to a non-linear secondorder ordinary di↵erential equation for (z) (see (13)). This equation cannot be solved analytically but an asymptotic analysis predicts (z) ⇠ z 2 for z ⌧ 1 and (z) ⇠ (cte z) 3/2 for z & 1. Data is compatible with this result (not shown). Note that, unlike the Euler-Voigt-↵ and ↵V-EDQNM models, the Leith model presents a sharp cut-o↵ instead of an exponential decay.

IV. CONCLUSION

In summary the Euler-Voigt-↵ model allowed us to show that its self-truncation regime reproduces the behavior of the truncated Euler equation [5]. We also found evidence for self-similarity in the long-time behavior of the energy spectrum. Introducing two di↵erent simplified models, the ↵V-EDQNM model and the ↵V-Leith model, we were able to show that they present behaviors similar to that of the Euler-Voigt-↵ model. We were able to determine the analytical values of the self-similar exponents of the simplified models.

In the present work we have used only integer values for . As was noted in Sect.II A, choosing = 11/3 yields an absolute equilibrium E(k) ⇠ k 5/3 and, in two-dimensions, = 2/3 also yields Kolmogorov scaling. This can represent an interesting alternative to the fractal decimation method that was used in reference [14].

The present work can be naturally extended to the 2D and 3D Ideal MHD equations. In this context, it was recently shown that dynamo action can be triggered by turbulence in absolute equilibrium [32]. Thermalization with = 11/3 would allow for a more realistic velocity spectrum, mimicking the infinite Reynolds limit.

Figure 1 . 1 :

 11 Figure 1.1: Spread of probability density of a Hadamard walk on n grid points with asymmetric initial state ™ 0= P m √ 0,m (c L |b L i + e i ¡ c R |b R i) T where (i) c L = 1, c R = 1, ¡ = 0 or (ii) c L = 1, c R = °1, ¡ = 0 and symmetric initial state c L = 1, c R = 1, ¡ = º/2. The spatial component is given by √ 0 = P m ± 0,m

Figure 1 . 2 :

 12 Figure 1.2: Square root of variance of a Hadamard Walk with symmetric initial state and localized position in the origin. The slope is (1 °1 p 2 ) 1/2 .

Figure 1 . 4 :

 14 Figure 1.4: Value of K µ for a given µ.

Figure 1 . 5 :

 15 Figure 1.5: Feynman's Checkerboard representation in the discrete space time lattice (x, t ).

2

  Discrete-time quantum walk with Fibonacci sequence coin Let us consider the two-dimensional spin state m, j ∈ C 2 , spanned by the orthonormal basis (b u , b d ), and defined by its discrete one-dimensional position m ∈ Z and discrete time j ∈ N 0 . The standard DTQW's time evolution is given by the application of the quantum coin operator (QCO) Ĉ on m, j = u m, j b u + d m, j b d = ✓ u m, j d m, j ◆ , followed

is the 2 ×

 2 2 identity matrix. To diagonalize the operator acting on , we perform a change of basis from (b u ,b d ) to the new basis, (b u , b d ), with = ub u + db d . The new basis components are

Fig. 1 , 123 Fig. 2

 11232 Fig. 1 Velocity contour plot for the Fibonacci sequence on the QCO, with α and β ∈ [0, π 2 ]

Fig. 3

 3 Fig. 3 Contour plot for velocity in the case in which Fibonacci sequence is on the step operator for the three periodic Fibonacci sequence

Figure 2 . 1 :

 21 Figure 2.1: Distortion in a vector field in the transport from a point i to a point j of the curved spacetime.

5 FIG. 1

 51 FIG. 1. (Color online) Density of the QW vs null geodesics (solid curves) of the 2D Schwarzschild metric for various values of λ [see text and Eq. (14)]. The singularity is represented by the dotted and dashed line on the left and the horizon (which is a null geodesic) is represented by the dashed line. (a), (b) The two branches of the QW which starts inside the horizon end up on the singularity. The (red) solid line represents the limit of the definition domain D of the QW. (c) One branch of a QW which starts on the horizon stays on the horizon while the other branch ends up on the singularity. (d) One branch of the QW which starts outside the horizon propagates away from the black hole, and the other branch ends up on the singularity.

field 51 2. 3 A

 513 synthetic electric gauge field 2.3.1 Publication: "Quantum Walks in artificial electric and gravitational gauge fields."

Fig. 1 .

 1 Fig. 1. Relative difference between the density of the QW and the density of the solution of the Dirac equation N rel = p < (N QW N D ) 2 >/ < (N D ) > plotted at T = 100 versus ✏ = 1X in Log-Log representation. We have plotted the relative difference in the same conditions as in Fig. 2(a) but for five different resolutions (i.e. value of ✏): from the left n = 2 8 , 2 9 , 2 10 , 2 11 , 2 12 , 2 13 . The solid black line represents the expected ✏ 1 scaling.

Fig. 2 .

 2 Fig. 2. Quantum simulation of Eqs. (1) and (2) representing a Dirac particle in a constant and static electric background (see Section 4.2 and Eq. (39)). The initial condition is a Gaussian wave packet of positive energy solutions with width X = 0.005 in (a), 0.01 in (b), 0.03 in (c), 0.08 in (d) and resolution n = 2 9 .

Fig. 3 .

 3 Fig. 3. (Color online) Time evolution density of the QW vs. null geodesics (solid curves) of the 2D Schwarzschild metric with = 1. The initial condition is (0, X ) = p N 0 (X)(b L + Ib R ) with an initial Gaussian density N 0 of X = 0.5 centered on X 0 = 50.5. The singularity is represented by the dotted and dashed line on the left and the horizon (which is a null geodesic) is represented by the dashed line. The two branches of the QW which starts inside the horizon end up on the singularity. The (red) solid line represents the limit of the definition domain D of the QW. In (a1) note that the right branch of the QW lags slightly behind the null geodesic when approaching the r = 0 singularity. The agreement between the geodesics and the density profile of the walk gets better as we increase the resolution of the simulation: 200 gridpoints in (a1), 800 in (a2) and 1600 in (a3).

1

 1 Quantum Decoherence: An introduction . . . . . . . . . . . . . . . . . . . . . . . . 69 3.2 Quantum Walks and decoherence . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 3.2.1 An overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 3.2.2 A qualitative picture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 3.2.3 Projections cause spin and spatial decoherence . . . . . . . . . . . . . . . . . 74 3.3 Publication : "Discrete-time Quantum Walks in random artificial Gauge Fields." 75

FIG. 1 :

 1 FIG. 1: (Color online) (left) Probability profile of the average transport in a random ⇠-field (a) and ✓-field (c) vs grid point m at j = 160 (a) and j = 70 (c) for di↵erent values of the noise parameter . (right) Probability profile of the average transport in a random ⇠-field (b) and ✓-field (d) vs grid point m for = 0.5 (d) and = 0.8 (d) at di↵erent time steps. Square marker represents fully decoherent regime, solid line the intermediate regime and dashed line indicates a fully coherent state.

1 FIG. 3 :

 13 FIG.2: Log-lin plot of time evolution of the spin coherence Cj in a random ✓-field for various values of the noise parameter .

Figure 4 . 1 :

 41 Figure 4.1: Classical Galton Board.

Figure 4 . 2 :

 42 Figure 4.2: (Left) The optical ring cavity for the implementation of the Optical Galton Board. The EOMs are the electro-optical modulators and the BS is the beam splitter. The solid gray rectangle is the partially reflecting mirror serving as input and output port and the black ones are the fully reflecting mirror of the cavity.(Adapted from Knight et al. [10]) (Right) On the top the probability distribution of a classical Galton board, and on the bottom the density profile of frequencies in the wave-mechanical case within the Landau-Zener crossings. © 1 is a control parameter proportional to time. (Adapted from Bouwmeester et al. [4])

Figure 1 .

 1 Figure 1. (Color online) Time evolution of the probability density j,m for (a) NLOGB (g=10 fi) and (b) an Hadamard DTQW (g=0) with a symmetric initial condition 0,m = " 0,m 2 (b≠ + ib+). Resolution N = 256.

  Fig. (2.a) displays how the probability j,m typically depends on m at fixed large values of j and Fig. (2.b) displays H ú (p).

  1.a). This symmetry breaking becomes greater with the time j (see Fig.(2.d) and depends on the resolution of the DNS and the strength of non linearities. In particular Fig.2.d show that the symmetry breaking starts from the round-o noise and we have verified that adding a non-symmetric noise to the initial condition produces a translation of the start point (data not shown) confirming that the symmetry breaking is due entirely to the round-o noise.

Figure 2 .

 2 Figure 2. (a) Long time behavior. Log-Lin plot of the density j,m at time j = 5 ◊ 10 5 (yellow cercles) for the NLOGB with a symmetric initial condition 0,m = " 0,m 2 (b≠ + ib+) (black point) for g=10 fi. Resolution N = 64. (b) Time evolution of the probability density j,m in same conditions that Fig. 1.b for short time behavior. (c) Log-Log histogram H(p) of the density j,m at time j = 5 ◊ 10 5 for the NLOGB with a symmetric initial condition 0,m = " 0,m 2 (b≠ + ib+) and g=10 fi. Resolution N = 64. (d) Asymmetry measure j versus time for di erent values of the resolution. j = q N/2≠1 m=0 j,m ≠ q N ≠1 m=N/2 j,m.

Figure 3 . 2 "

 32 Figure 3. (Color online) (a) Log-Lin plot of the density (T, X) at time T = 5 ◊ 10 5 (red square) obeying the Dirac equation with a symmetric gaussian initial condition (0, X) = f (X) 2 (b≠ + ib+) (black point) for g =10fi. Resolution N = 64. The gaussian shape f (X) = 1 2fi ‡ exp ! ≠X 2 / Ô 2 ‡ 2 " where ‡ =10 x (b) Histogram H(p) of density .

Figure 4 .

 4 Figure 4. Histograms H(p) of the thermalized state density (X) for the NLDE and for the stochastic equations(13) for g = 10fi. The conserved quantities and the noise coe cient are: top E= -19.8, P = 0, = 1 and ÷ =1; bottom E= -17.12, P = 0, = 1, T0= 250, T1=312 and ÷ =1.6.

Figure 5 .

 5 Figure 5. Log-Log plot of the relative di erence "j at time j = 400, defined as Ô <(fl QW ≠fl D ) 2 > <fl D >

Figure 6 .

 6 Figure 6. (Color On line) Time evolution of energy (top) E(T ), momentum (center) P (T ) and particle number (T ) simulated by a pseudo-spectral code for the spatial part and a 4th-order Runge Kutta for the time step. The resolution N = 128. The blue courbes represents the non dealiased code.

FIG. 1 .

 1 FIG. 1. (Color online) Temporal evolution (indicated by arrows) of the energy spectrum E↵(k) for = 4 and k↵ = 80. Resolution 1024 3 (kmax = 342). The dashed lines respectively display the Kolmogorov k 5/3 and the equipartition k 2 scaling. The self-truncation wavenumber is indicated by the small vertical arrow.

FIG. 2 .

 2 FIG. 2. (Color online) Energy spectrum E↵(k) versus k for = 2 and k↵ = 50, 100 and 200 at t = 10.2. (b) Same conditions than in (a) but = 4. c) Same conditions than in (a) but = 6. The dashed-black lines corresponds to Kolmogorov scaling E↵(k) ⇠ k 5/3 and energy equipartition E↵(k) ⇠ k 2 . The thermalization wave number k th and selftruncation kst are indicated by arrows in (c).

FIG. 3 .

 3 FIG. 3. (Color online) Temporal evolution of (a) k th , (b) E th , (d) ✏ = @E th @t , (e) k th /k d with k d estimated based on the self-truncation wavenumber kst, see Eq.(10).

  requires that k d ⇠ k th . The ratio k th /k d is displayed on figure Fig.3.d. It indeed seen to be of order unity and reasonably constant in time. Thus the large-dynamics of Eq. (

  FIG. 4. (Color online) (a) Exponential decay of E↵(k) at t = 7.8. b) Time Evolution of energy spectrum fit parameter of Eq.(11): Horizontal lines correspond to kmax = 2 (dashed black line) and k↵ = 2 (dot-dashed black line). Exponential law (dot black line) = 2.7 exp( t/0.56) from reference[19].

FIG. 5 .

 5 FIG. 5. (Color online) Time evolution of kst(t) at resolution 1024 3 (see Table I) . a) Di↵erent values of for k↵ = 20. The crosses show the Euler evolution of kst (↵ = 0). b) = 4 and di↵erent values of k↵. The horizontal dashed lines represent the values of k↵/kmax.

FIG. 6 .

 6 FIG. 6. (Color online) a) Temporal evolution of the selfsimilar function (k/kst(t)) = E↵(k, t)kst(t) (see Eq.13) for = 4 and k↵ = 4. Data from direct numerical resolution of (2) at resolution 512 3 . The inset shows the temporal evolution of kst(t)/kmax for di↵erent values of .

5 ± 6 ⇥

 56 10 3 0.25 ± 3 ⇥ 10 3 0.07 ± 5 ⇥ 10 3 ↵V-EDQNM 0.33 ± 5 ⇥ 10 5 0.11 ± 9 ⇥ 10 5 0.085 ± 1 ⇥ 10 4 ↵V-Leith 0.33 ± 10 6 0.15 ± 2 ⇥ 10 4 0.09 ± 9 ⇥ 10 6

FIG. 7 .

 7 FIG. 7. (Color online) a) Energy spectra E↵(k) at t = 3.5 of the ↵V-EDQNM model, for di↵erent values of , obtained with k↵ = 10 4 and kmax = 43918 (corresponding to resolution 131072, with 14.2 points per octave). b) Temporal evolution of the self-similar function (k/kst(t)) = E↵(k, t)kst(t) for = 4 and k↵ = 5 and kmax = 692 (28.4 points par octave). The inset b.1 shows the temporal evolution of kst(t) for different values of . The inset b.2 shows the ↵V-EDQNM (red dashed line) and ↵V-Leith (solid green lines) theoretical predictions (19) and (25) for the self-truncation exponent ⌘ and numerical data from table II.

  2 exp [ z] and it reads = min [ , 5/2]. The theoretical prediction confronted with the data are presented in good agreement in the inset b.2 of Fig.7. B. ↵-Voight Leiht model.

FIG. 8 .

 8 FIG. 8. (Color online) a) Energy spectra E↵(k) of the ↵V-Leith model (24) for di↵erent values of obtained with k↵ = 400 and kmax = 4000. b) Temporal evolution of the self similar function (k/kst(t)) = E↵(k, t)kst(t) for = 4 and k↵ = 2. The inset shows the temporal evolution of kst(t) for di↵erent values of .

1 Quantum walks 1.1.1 Introduction Classical

  B. Schumacher and R. F. Werner. Reversible quantum cellular automata. arXiv preprint quantph/0405174, 2004. 1.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9 1.1.2 From Classical to Quantum random walks . . . . . . . . . . . . . . . . . . . . 10 1.1.3 General Setup of a Discrete Time Quantum Walks . . . . . . . . . . . . . . . . 11 1.1.4 Qualitative description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.1.5 Quantitative description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.1.5.1 Fourier methods for the Hadamard walk . . . . . . . . . . . . . . . . 16 1.2 Connections between Quantum Walks and Relativistic Wave Equations . . . . . . 17

1.2.1 Quantum walks and Feynman's Checkerboard . . . . . . . . . . . . . . . . . . 17

1.2.2 Homogeneous DTQWs and Weyl equation . . . . . . . . . . . . . . . . . . . . 19 1.2.3 Publication: "Massless Dirac Equation from Fibonacci Discrete-Time Quantum Walk" . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 1.random walks (CRWs) are employed to model phenomena as chemical reactions

  1/2 .

	0.18
	0.16
	0.14
	0.12
	0.1
	0.08
	0.06
	0.04
	0.02
	0	-1	-0.8	-0.6	-0.4	-0.2	0	0.2	0.4	0.6	0.8	1
							m/j

  2.1 we can geometrically observe that the curvature of the manifold is measured by the quantity [X AE +±X AE ]°[X AE +±X AE ] = ±X AE °±X AE . The affine connection represents the discrepancy between the variation ±X AE and the vector ±X AE . It represents a multiplicative factor between X Ø itself and the finite displacement ±x ∞ . Let us denote this factor by °AE Ø∞ X Ø .

  Corresponding author. Tel.: +33 144277286; fax: +33 144277287.

E-mail address: fabrice.debbasch@gmail.com (F. Debbasch). 0378-4371/$ -see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.physa.2013.11.036

  1 known as Schmidt coefficients.

		Chapter 3. QUANTUM WALKS, DECOHERENCE AND RANDOM SYNTHETIC GAUGE FIELD
	6
		θ=0, π/2
		θ=π/36
	5	θ=π/18
		θ=π/9
		θ=π/4
	4	θ=π/3 θ=2π/5
		θ=2π/11
	3						6
							5
	2					100	4	DTQW Hadamard B -π/2,π/4,π/2,π/2
						H	3
	1						2	B -π/2,0,π/2,π/2	B -π/2,π/2,π/2,π/2
	0 0	10	20	30	40	50	60 0 1	70 0.5	80 θ	1	90	1.5	100
						j

s Figure 3.1: Variation of entropy of measurement (or entanglement) (S r ) j VS time for different values of µ. Inside the figure: the variation of entropy VS µ for a given number of steps.

  6 4 c 2 e 2iK ics e +i⇠ e 2ip +ics e i⇠ e +2ip s 2 e 2iK ics e +i⇠ e 2iK c 2 e +2i⇠ e 2ip s 2 e +2ip +ics e +i⇠ e 2iK +ics e i⇠ e 2iK s 2 e 2ip c 2 e 2i⇠ e +2ip ics e i⇠ e 2iK s 2 e 2iK +ics e +i⇠ e 2ip ics e i⇠ e +2ip c 2 e 2iK

	3
	7 7 5 .

  LERMA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ. Paris 6, UMR 8112, F-75014, Paris France

	Marc Brachet Laboratoire de Physique Statistique de l'Ecole Normale Supérieure / PSL Research University, associé au CNRS et aux Universités Pierre-et-Marie-Curie Paris 06 et Paris Diderot, 24 Rue Lhomond, 75231 Paris, France
	(Dated: June 13, 2015)

  Giuseppe Di MolfettaLERMA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universits, UPMC Univ. Paris 6, UMR 8112, F-75014, Paris France limit is shown to correspond to Galerkin truncation of the Euler equation. Direct numerical simulations (DNS) of the model are performed with resolutions up to 2048 3 and Taylor-Green initial data. DNS performed at large demonstrate that this simple classical hydrodynamical model presents a self-truncation behavior, similar to that previously observed for the Gross-Pitaeveskii equation in Krstulovic and Brachet [Phys. Rev. Lett. 106, 115303 (2011)]. The self-truncation regime of the generalized model is shown to reproduce the behavior of the truncated Euler equation demonstrated in Cichowlas et al. [Phys. Rev. Lett. 95, 264502 (2005)]. The long-time growth of the self-truncation wavenumber kst appears to be self-similar. Two related ↵-Voigt versions of the EDQNM model and the Leith model are introduced. These simplified theoretical models are shown to reasonably reproduce intermediate time DNS results.

		Giorgio Krstlulovic
		Laboratoire Lagrange, UMR7293, Université de Nice Sophia-Antipolis,
		CNRS, Observatoire de la Côte dAzur, BP 4229, 06304 Nice Cedex 4, France
		Marc Brachet
	Laboratoire de Physique Statistique de l'Ecole Normale Supérieure / PSL Research University,
		associé au CNRS et aux Universités Pierre-et-Marie-Curie Paris 06 et Paris Diderot,
		24 Rue Lhomond, 75231 Paris, France
		(Dated: April 13, 2015)
	A generalization of the 3D Euler-Voigt-↵ model is obtained by introducing derivatives of arbitrary
	order	(instead of 2) in the Helmholtz operator. The	! 1

  . The CFL condition t k ↵ ⇠ 0.1

	Run Res.	k↵ tmax Run Res.	k↵ tmax
	Euler 1024	30	11 2048 2 100 15
	1	1024 2 20 50	12 2048 2 200 15
	2	1024 2 40 50	13 2048 4 50 15
	3	1024 2 80 50	14 2048 4 100 15
	4	1024 4 20 50	15 2048 4 200 15
	5	1024 4 40 50	16 2048 6 50 15
	6	1024 4 80 50	17 2048 6 100 15
	7	1024 6 20 50	18 2048 6 200 15
	8	1024 6 40 50	19 512 2 4 2300
	9	1024 6 80 50	20 512 4 4 2300
	10 2048 2 50 15	21 512 6 4 2300
	TABLE I. List of runs of the generalized 3D Euler-Voigt-↵
	model (2) with TaylorGreen initial data (9) and maximum
	integration time tmax.

TABLE II .

 II Values of the exponent ⌘ of the self-truncation wavenumber kst(t) ⇠ t ⌘ (see (12)) obtained from direct numerical simulation of the Euler-Voigt-↵ model (2), ↵V-EDQNM(14)(15)(16)(17) and ↵V-Leith model (r = 2).(24).

Note that in case of Eqs.(1.15, 1.16), the probability to find the particle in position m = ®1 starting from a symmetric initial state, reads P ® = (1 ® sin(2µ) cos(≥ + ª)). This probability distribution will be equal and lead to a left-right symmetry in position if ≥ + ª = (2k + 1) º 2 , k 2 N.

Let us notice that Dirac equation in (1+1) curved spacetime has been derived by Sinha and Roychoudhury[51],Morsink and Mann[42], Carter and McLenaghan[13], Birrell and Davies[6], Fulling[22], Parker[44]

Note that a non-unitary version of QWs was introduced by Attal et al.[5][6]. These QWs, called Open Quantum Walks (OQWs), are formulated directly as quantum Markov chains. OQWs are very useful to model dissipative quantum systems dynamics.

The unitary operator of the QW on the Cayley graph of an Abelian group must have all distinct eigenvalues. See Theorem

3.5 in D.[1]. In particular in the non-degenerate case, verified only in the case of odd N , the distribution is uniform. For even N, the limiting distribution displays more interesting stationary profile distributions then in odd case and they have been investigated by Bednarska et al.[3].

Di Molfetta Giuseppe, Yutaka Shikano, Fabrice Debbasch, Marc Brachet, Work in progress, project with the financial support of JSPS summer program (SP14203).

G. Di Molfetta, F. Debbasch Non-Abelian Quantum Walks, preprint,

Remerciements

Acknowledgments The authors acknowledge useful discussion with Fabrice Debbasch and Jingbo Wang. G.D.M. expresses thanks to the JSPS summer program (SP14203) for financial support, while L.H. and B.L. thank the Department of Education of the Australian Government for financial support through the AsiaBound Grants Program. T.W. thanks the JSPS KAKENHI 25400188 and IMS Joint Study Program for financial support, while Y.S. thanks the DAIKO Foundation, also for financial support.

Chapter 1. HOMOGENEOUS QUANTUM WALKS

Bibliography

Acknowledgements: We acknowledge useful scientific discussions with Claude Bardos, Uriel Frisch, Annick Pouquet, Samriddhi S. Ray and Edriss Titi. The computations were carried out at IDRIS and the Mésocentre SIGAMM hosted at the Observatoire de la Côote d'Azur.

Appendix A: Interpretation in terms of artificial gauge fields It has been proven in [20][21][22] that quantum walks in (1 + 1) dimensional space-times can be viewed as modeling the transport of a Dirac fermion in artificial electric and gravitational fields generated by the time-dependance of the angles ✓ and ⇠. We recall here some basic conclusions obtained in [20][21][22] and also o↵er new developments useful in interpreting the results of the present article.

The DTQWs defined by (1) are part of a larger family whose dynamics reads:

where B(✓, ⇠, ⇣, ↵) = e i↵  e i⇠ cos ✓ e i⇣ sin ✓ e i⇣ sin ✓ e i⇠ cos ✓ .

(A2)

The walks in this larger family are characterized by three time-and space-dependent Euler angles (✓, ⇠, ⇣) and by a global, also time-and space-dependent phase ↵. They have been shown to model the transport of Dirac fermions in artificial electric and relativistic gravitational fields generated by the time-dependence of the three Euler angles and of the global phase. In a (1 + 1) dimensional space-time, an electric field derives from a 2-potential A j,m = (V j,m , A ,mj ) and a relativistic gravitational field is represented by 2D metrics G j,m . The walks considered in this article correspond to

where ⇠j and ✓j are random variables which depend on the time j and ↵ = 3⇡/2. According to [22], these walks model the transport of a Dirac fermion in an electric field generated by the 2-potential

and in a gravitational field caracterized by the metric

Since relativistic gravitational fields are represented by space-time metrics [48], making the angle ✓ a time-dependent random variable is equivalent to imposing a time-dependent random gravitational field. To better understand the electric aspects of the problem, let us recall that the DTQWs defined by (A1) exhibit the following exact discrete gauge invariance [22]:

where j,m = j,m+1 + j,m 1 2 j+1,m j,m = j,m+1 j,m 1

(A7)

and is an arbitrary time-and space-dependent phase shift. Let us now define a new quantity E j,m by

where the actions of the operators D s and D s on an arbitrary time-and space-dependent quantity u j,m are (D s u) j,m = u j,m+1 u j,m 1 2 (A9) and (D t u) j,m = 2u j+1,m u j,m+1 u j,m 1 2 .

(A10)

The operators D s and D t are discrete counterparts of space-and time-derivatives. It is straightforward to check that the quantity E j,m is gauge invariant and coincides, in the continuous limit, with the standard electric field E(t, x), defined by E(t, x) = @ x V + @ t A. The quantity E j.m is thus a bona fide electric field in discrete space-time. For the DTQWs considered in this article, this electric field depends only on the time j and is related to the angle ⇠ by E j = ⇠j+1 ⇠j . Making this angle a time-dependent random variable is thus equivalent to imposing a random electric field.

Appendix B: Aymptotic computation of the eigenvalues and eigenvectors of the averaged transport operators Let us here compute the eigenvalues e/g r and eigenvectors w e/g r , r = 1, 2, 3, 4 only for values of K much smaller than unity. We do not perform an expansion in p because the initial condition is uniform in p and the average evolution does not localize the density operator around p = 0. Indeed, the initial condition is localized at x 0 = x i.e. does not exhibit any spatial correlation and the dynamics does not create spatial correlations.

The second order expansions of the operators Re and Rg in K read:

and

For K = 0, these two matrices are both block diagonal and we write Re/g 2 (K = 0, p, ) = diag(1, M e/g (p, )), where M e/g (p, ) are 3 ⇥ 3 matrices acting in the space spanned by (u 2 , u 3 , u 4 ). The matrices Re/g 2 (K = 0, p, ) share u 1 as common eigenvector, which we identify as w e/g 1 (K = 0, p, ); the associated eigenvalue is e/g 1 (K = 0, p, ) = 1. The other eigenvectors and eigenvalues, at zeroth order in K, are those of M e/g (p, ). These eigenvalues can be computed analytically by solving the third-order characteristic polynomials associated to these matrices. The explicit expressions of these eignevalues are quite involved and need not be replicated here. What is important is how the moduli of these eigenvalues compares to unity. Direct inspection reveals that the moduli of all three e r (0, p, ), r = 2, 3, 4 are strictly inferior to unity if is not vanishing. The same goes for all three eigenvalues in the gravitational case, except for one of them which reaches ±1 independantly of for p = ±⇡ and is also equal to +1 for p = 0; the eigenspaces corresponding to g 4 (±⇡, ) and g 4 (0, ) are identical and generated by u 4 , which we choose as w g 4 (p = ±⇡, ) = w g 4 (0, ). For other values of p, the eigenvalue g 4 (p, ) and the eigenvector w g 4 (p, ) are defined by continuity. All other eigenvectors need not be specified for what follows.

Let us now turn to non vanishing values of K. The characteristic polynomials of Re/g 2 (K, p, ) contain terms of order 2 and 4 in K; at lowest order in K, the corrections to the eigenvalues e/g j (K = 0, p, ) thus scale generically as K 2 . Let be the variable of the characteristic polynomials. At second order in K, the K-dependent correction to each of the zeroth order eigenvalues e/g r (K = 0, p, ) can be found by expanding the characteristic polynomial of Re/g 2 (K, p, ) at first order in ⇣ e/g r (K = 0, p, ) ⌘ and by keeping only the terms scaling as K 2 . This gives rational expressions for the corrections to the eigenvalues; these rational expressions can be further simpified by a final expansion around K = 0 if p is treated as a finite, non infinitesimal quantity i.e. | K |⌧| p |. One then finds:

and

(B5)

Note that ↵ g is actually independent of p. Note also that the condition | K |⌧| p | does not hinder asymptotic computations, at least on the infinite line. Indeed, as time increases, the density operator becomes more and more localized around K = 0, but it does not localize in p-space [? ]. If one works on the infinite line, both K and p are continuous variables and the localization of the density operator around K = 0 implies that the size of the region in Let us now use the above results to determine the time evolution of the average density operator in both cases under consideration. The first step is to express the initial condition, ⇢j=0 (K, p) = (u 1 u 4 )/2 for all (K, p), as a linear combination of the eigenvectors w e/g r (K, p, ). We thus write, for a = 1, 2, 3, 4

Now, for all r,

since e/g r (K = 0, p, ) = 1. It follows that, for small enough K, the contributions to (C5) proportional to ⇣ e/g r (K, p, ) ⌘ J are much smaller than the contribution proportionnal to ( e/g 1 (K, p, )) J for all values of p and such that | e/g r (K = 0, p, ) |< 1. According to the above discussion, this is realized for all r 6 = 1 and for all values of p and , except in case 2 (random gravitational field) for r = 4, p = ±⇡ or p = 0 and all values of . What happens at p = ±⇡ has no incidence on the computation of the density operator in physical space. Indeed, for finite values of M , the maximum value p max of | p | is p max = (2M/(2M + 1))⇡ < ⇡. Thus ±⇡ is only reached in the limiting case of infinite M i.e. for quantum walks in the infinite line. However, ±p max = ±⇡ then only appear as upper and lower bounds for integrals over p, and the values taken by ⇢(J, K, p) at points ±⇡ does not modify the values of the integrals. Moreover, all current computations are only valid for | p |⌧| K | and are thus sl a priori invalid for p = 0. What happens around p = 0 has however no relevance to asymptotic computations on the infinite line because, as time increases, the density operator becomes more and more localized around K = 0 (see discussion below (B5)).

For large enough J and small enough K, the double sum in (C5) thus simplifies into:

At lowest order in K, 1 ↵ e/g (K, p, )K 2 J = 1 ↵ e/g (K, p, )JK 2 . We will now restrict the discussion to scales K and times J obeying JK 2 K i.e. JK 1. Note that the maximum spatial spread of ⇢ at time J is L max (J) = 2J, so that the minimum value of K for which ⇢ takes non negligible values at time J is of order K min (N ) = 1/J. The condition JK 1 thus restricts the discussion to length scales much smaller than L max (N ). In particular, consider the time-dependent scale K J = K ⇤ / p J, where K ⇤ is an arbitrary time-independent wave-vector. The wave-vector

K J for su ciently large J. Thus, the possible di↵usive behavior of the averaged transport is encompassed by the present discussion.

With the above assumption, equation (C7) implies the following approximate but very simple expression for the long time (large J) density operator in Fourier space:

In particular, for K J = K ⇤ / p J (where K ⇤ is an arbitrary but J-independent wave number) and large enough J,

Absolute equilibrium in conservative systems

A general introduction

In classical statistical mechanics all dynamical systems in equilibrium or out of equilibrium can be described by a Liouville equation, a milestone in many-body theory. In the following we try to introduce the usual notation and definition that the reader can find in a standard textbook of statistical mechanics (e.g. landau2013course).

The Liouville equation describes the time evolution of the probability of finding the system in a given region of the phase space and expresses the conservation of the probability current in this space. Let us introduce a general dynamical system defined by:

The 2N-dimensional vector X, which represents the coordinates of the system, evolves in the 2N-dimensional phase space with a probability distribution P (X, t ). The probability P (X, t )d 2N X obeys to the following PDE:

which is called the Liouville equation, derived for the first time by the french mathematician Joseph Liouville in 1838. Now let us consider that the system is Hamiltonian and that X = Consider for all (n, j) oe N 2 , the collection W n j = ( k,m ) k=nj,moeZ . This collection represents the state of the NLOGB at 'time' k = nj. For any given n, the collection S n = (W n j ) joeN thus represents the entire history of the NLOGB observed through a stroboscope of 'period' n. The evolution equations for S n are those linking W n j+1 to W n j for all j. The method employed here to obtain the continuous limit of a generic S n was introduced in [2,20]. One first introduces a time-scale • , a length-scale ⁄, an infinitesimal ' and interpret the space-index m as referring to position x m = m'⁄ = m x and the time index j as referring to the instant t j = j'• = j t. The formal continuous limit is obtained expanding the equations defining S n in Taylor series around ' = 0 and by letting ' tend to zero. For the limit to exist, all zeroth order terms of the Taylor expansion must identically cancel each other and the di erential equation describing the limit is then obtained by equating to zero the non identically vanishing, lowest order contribution.

The original NLOGB S 1 does not admit a continuous limit because the zeroth order terms do not cancel each other identically. The equations defining S 2 read:

where

and

These equations admit a formal continuous limit, which reads:

where

and T = t/• and X = x/⁄. The operator P is self-adjoint and its eigenvalues are ≠1 and +1. Two eigenvectors associated to these eigenvalues are

and

The family (B ≠ , B + ) forms an orthonormal basis of the two dimensional spin Hilbert space. In this new basis, equation (A3) reads:

where

Appendix B: Numerical Methods

We restrict ourself to 2fi-periodic boundary conditions. A generic field (X) is thus evaluated on the N collocation points X m = 2fim/N , with m = 0, N ≠ 1 as m = (X m ). The discrete Fourier transforms are standardly defined as

Awards and Fellowship

Selected by CNRS for the JSPS Summer Program 2014, Fellowship "JSPS", CNRS -ENS Paris -Institute for Molecular Science Okazaki, Supervisor: Yutaka Shikano. The equations simulated by our numerical code are the finite difference equations of the QWs, the linear and nonlinearDirac equations and the stochastic hyperbolic equations. In all these cases, except when the parameters of the DTQW are non-homogeneous in space, the discrete dynamical systems associated to those equations, have been simulated by the pseudo-spectral methods. In all other cases we have employed the Direct Numerical Simulations (DNS).

Workshop

A.1 Spectral Methods

A.1.1 Fundamentals

Let us introduce the periodic fields f that verify: f (x + L) = f (x) where L is the periodicity of the lattice. In order to simplify, let us choose L = 2º. Then the continuous function f (x) can be expressed on the Fourier basis:

We define the scalar product as:

where the Fourier coefficients read:

Generally, we can access to the knowledge of f through its value at the discrete finite set of points:

with x j = j ¢ x and j = 0, 1, 2, ..., N °1, so that the discrete space interval ¢x = 2º N . This set is called "physical space". The Fourier coefficients are therefore approximated by:

with x n = 2ºn N , n = 0, ..., N °1. The latter equation define what we call the Discrete Fourier Transform (DFT). We then project f over a basis of N functions (e.g . the trigonometric function), and finally define f N , the approximated function of f , just inverting the last DFT:

The points k = °N 2 , ..., N 2 °1 represent the discrete "spectral space". Now let us observe that computing the DFT and its inverse requires O(N 2 ) operations. This has represented a real problem, until Cooley and Tukey [1], introduced in 1965 the Fast Fourier Transform, FFT. This algorithm permits to reduce the number of operations to O(N log 2 N ). Another advantage of the FFT, and in general of the Fourier transforms, concerns the factorization of the convolutions in 2 and 3 dimensions, but it will be not useful in this thesis where all physical systems are described in (1 + 1) dimensions.

In conclusion, let us keep in mind that in this thesis, all the field vector ™(x) and its Fourier transformed are approximated in 1-dimensional space by:

(

A.2 Convergence in spectral methods

In order to approximate a solution f of a PDE, we can employ finite difference methods of order p. In these methods the precision of the approximation depends on the order p of the Taylor development, and it is in general of O(¢x p ). In spectral methods, instead, the coefficients fN are computed over the N points fixed by our code resolution. In other words, we can increase the precision increasing the numerical resolution and notably it is not fixed a priori, such as in the finite difference methods. For example, for an interval ¢x / O(1/N ), the error in the pseudo-spectral code is of order O[(1/N ) N ]. The error decreases faster than any polynome of order N . This is what we call an exponential convergence.

A.3 Approximate a PDE by spectral method

Now, let us suppose that we need to solve a PDE in 1-dimensional space and that it takes the following form:

where D is a linear differential operator, T is a positive number and ™ is defined in [0, 2º]£]0, T ]. The approximation of ™ is defined as the solution of the approximated PDE. We should find a function ™ N (x, t ) in [0, 2º]£]0, T ], so that:

A.3. Approximate a PDE by spectral method 131

where ™ N is defined in [0, 2º]£]0, T ].

In order to solve this problem we should compute the Fourier coefficient ™ N so that D(™ N ) vanishes. In this thesis we have employed two important spectral methods: the Galerkin method and the pseudo-spectral method.

A.3.1 Galerkin method

In order to illustrate the idea of this method, let us introduce a simple example from fluidodynamics and present the inviscid Burgers equation:

with the following periodic conditions:

The DFT of ™, therefore, obeys to the following equation:

If the linear terms contain only wave number |k| ∑ N /2, the nonlinear term of order 2 contains wave number |k| ∑ N . The Galerkin method consists in retaining only the terms with wave number |k| ∑ N /2. The above equation, after simplifications, reads:

The sufficient number of operations to evaluate the convolution is of order O(N 2 ), that represents a high computational cost. For this reason we compute the nonlinear term in physical space (by just N operation) and we translate the solution from the spectral to physical space and vice versa by FFT (O(N log N )). We will see, in the next section, that the pseudo-spectral method is based on this idea.

A.3.2 Pseudo-spectral method

The pseudo-spectral methods consist in computing the approximation of a function over a determined basis, interpolating over a finite set of grid points (or collocation points). Therefore, the PDE will be exactly solved at the collocation points. We have employed this idea in computing the convolution in the physical space: we define as basis the DFT of the trigonometric functions (Dirac ±), which corresponds exactly to our collocation points. For the nonlinear term in A.12, we use two inverse FFT in order to compute ™ and @ @ x ™ in the physical space. Then, we calculate the product and we use again a FFT in order to come back to the spectral space:

X q=0 e °i 2º l (n°p°q) N Tp @ ™ @ x q (A. 16)

In other words, the approximated Fourier coefficients in the above equation are composed by all the exact coefficients along with other terms for which the function e i nx is indistinguishable from the e i nx [N ] (modulo N). This phenomenon is know as aliasing. Note, in conclusion, that extending these results to a generic nonlinear Dirac equation in (1 + 1) is straightforward.

A.3.3 De-aliasing

Let us study a lattice of N grid points where the wave number k varies in the interval [°N /2, N /2]. Therefore the values k 1 + k 2 vary between °N and N and, because of the aliasing problem, are replicated in the same interval. In order to solve this problem, we can truncate the spectrum by a cut-off |k max |. For a nonlinearity of order 2, it is possible to prove that 2k max °N < °kmax or equivalently k max < N 3 . In other words, we should eliminate all the values of k greater than N /3 and smaller than °N /3. In this way, all replicated values vanish at each time step. Let us remark that solving the aliasing problem is necessary to preserve correctly the integration by parts. In particular, the violation of this property leads to the non-conservation of the integrals, when they exist, in the considered EDP. More generally for a nonlinearity of order d we can de-aliase by the following rule:

A.3.4 Time-stepping

In the previous section we have introduced the pseudo-spectral methods to discretize the spatial part of the initial PDE and obtain an ODE system. Typically to discretize the PDE in an ODE, we should employ a method which is explicit for the nonlinear term and an implicit one for the linear part. In particular in this thesis we will use an Euler scheme (implicit or explicit) and a Runge-Kutta method of 4th order. Let us recall these three methods separately. Consider the partial differential equations of the form:

where L is a linear operator and G(™) is a nonlinear operator. Suppose that the initial condition is known ™(0, x) = ™ 0 (x). Discretization in time is given by all the t j 2 R + = j d t where j 2 N.

• Implicit Euler (order O(d t))

• Runge-Kutta (order O(d t 4 ))

The RK-schemes is the more accurate and we will use it when we need to verify conserved quantities. Otherwise we use the following combined scheme:

.

A.4 Discrete Fourier Transform

Let us now briefly describe the Fourier representation of the DTQWs. We know that QWs can be numerically implemented in two equivalent ways: in direct physical and in Fourier space. Generally the latter is preferable because the transformed unitary step operator becomes local and the translation operator is diagonal in Fourier basis. We therefore systematically prefer numerical simulation in Fourier space, when it is possible.

Consider the QWs, with the space homogeneous quantum coin, in physical space:

Let us restrict the domain of ™(x m , t j ) to a 2º-periodic boundary conditions and for any integer N >0 consider the set of points:

where m = 0,...,N °1 referred to as nodes or grid points; now define the lattice parameter ¢x, equal to 2º N . The discrete Fourier coefficients of the complex-valued function ™(x m , t j ) with respects to the grids points are:

where °N 2 ∑ k ∑ N 2 °1. Because of the orthogonality, over the interval [0, 2º], of the set of functions ¡ k (x m ) = e °i kx m , we can write the following inversion formula:

Tk (t j )e i kx m .

(A.34)

The above polynomial represents the Discrete Fourier transform (DFT) of ™(t j , x m ).

Let us mention that this DFT can be accomplished by the Fast Fourier Transform algorithm introduced by Cooley and Tukey [1]. In general FFT causes N to be a power of 2 and, if the data are fully complex, it requires O(N log 2 N ) instead of O(N 2 ) real operations required by straightforward sum. Moreover FFT incurs less errors, due to round-off, than the direct summation method Cooley et al. [2]. The way by which translation is accomplished in spectral space is given by the following operation:

e ±i k¢x √k (t j )e i kx m (A. 35) and therefore, the equations A.31 read:

™(x m , t j + ¢t ) = C (AE(t j ), µ(t j ), ª(t j ), ≥(t j ))

T k Tk (t j )e i kx m , (A. 36) where the trasformed translation operator is diagonal and reads as: Let us remark in conclusion that in this setting the continuous limit is automatically taken when N is increased. This allowed us to evaluate with high precision the discrepancy between the QW and the corresponding solution of the Dirac equation. where the projector P AEØ∞ = k Ø P AE∞ + k ∞ P AEØ with P AEØ = ± AEØ °kAE k Ø /k 2 . P AEØ is the projector into divergence-less function ( û(k)•k=0) and P G is the Galerkin projector that ensures the û(k) ¥ 0, 8|k| > k max . Conservation of energy E and Helicity H was already shown by Lee [12] and Kraichnan [5]. Thermalization was proved by Orszag [13] and Cichowlas et al. [2] showed that B.2 displays a Kolmogorov spectrum. Both conserved quantities, the energy and the helicity, help us to characterize the transition of the system to the fully relaxation. For instance, let us keep in mind the case of vanishing helicity, the only conserved quantity is the Energy and it reads:

TRUNCATED EULER-VOIGT-AE

and the absolute equilibrium is simply described by the Boltzmann weight: