Analog quantum simulation of the spinor-four Dirac equation with an artificial gauge field

TL;DR

This paper demonstrates how ultracold atoms in optical lattices can be engineered to simulate the full spinor-four Dirac equation, including effects of an artificial magnetic gauge field, advancing quantum simulation capabilities.

Contribution

It introduces a method to emulate the complete Dirac equation with artificial gauge fields using modulated ultracold atom systems in optical lattices.

Findings

Successfully simulates the spinor-four Dirac equation

Creates an artificial magnetic gauge field

Replicates relativistic quantum particle behavior

Abstract

A two-dimensional spatially and temporally modulated Wannier-Stark system of ultracold atoms in optical lattices is shown to mimic the behavior of a Dirac particle. Suitable additional modulations generate an artificial gauge field which simulates a magnetic field and imposes the use of the full spinor-four Dirac equation.

Equations

$$H_D=c\mathbf{\alpha }\cdot \left(\mathbf{p}-q\mathbf{A}\right)+\beta m_D{c}^2$$

$$\alpha_{i}=\left(\begin{array}[]{cc}0&\sigma_{i}\\ \sigma_{i}&0\end{array}\right)\qquad\beta=\left(\begin{array}[]{cc}\boldsymbol{1}&0\\ 0&-\boldsymbol{1}\end{array}\right).$$

$$H_0=\frac{p_x^2+p_z^2}{2m^{\ast }}+V_L(x,z)+F_{x}x+F_{z}z$$

$$V_L(x,z)=-V_1\left[\cos (2\pi x)+\cos (2\pi z)\right]$$

$${\Phi }_{n,m}(x,z)={\phi }_n^{(x)}(x){\phi }_m^{(z)}(z)$$

$$E_{n,m}=E_0+n{\omega }_B^{(x)}+m{\omega }_B^{(z)}$$

$${\phi }_n^{(x)}(x)={\phi }_0^{(x)}(x-n),$$

$$i{\partial }_t\left[\Psi \right]=\left\{c{\alpha }_x\left(p_x-A_x\right)+c{\alpha }_z\left(p_z-A_z\right)+\beta m{c}^2\right\}\left[\Psi \right]$$

$$\bar{V}(x,z,t)=$$

$$V_z\cos \left(\pi x\right)\cos \left(2\pi z\right)\sin ({\omega }_B^{(z)}t).$$

$$\psi (x,z,t)$$

$$\times {\phi }_n^{(x)}(x){\phi }_m^{(z)}(z)$$

$$i{d}_t{c}_{n,m}=\sum _{j,k\in \mathbb{Z}}$$

$$\times \exp \left[-i\left(j{\omega }_B^{(x)}+k{\omega }_B^{(z)}\right)t\right]c_{n+j,m+k}\}$$

$$⟨{\phi }_n^{(x)}|\cos \left(2\pi x\right)|{\phi }_{n\pm 1}^{(x)}⟩=⟨{\phi }_0^{(x)}|\cos \left(2\pi x\right)|{\phi }_1^{(x)}⟩\equiv {\Omega }_x$$

$$⟨{\phi }_n^{(x)}|\cos \left(\pi x\right)|{\phi }_n^{(x)}⟩⟨{\phi }_m^{(z)}|\cos \left(2\pi z\right)|{\phi }_{m\pm 1}^{(z)}⟩$$

$$={\left(-1\right)}^n⟨{\phi }_0^{(x)}|\cos \left(\pi x\right)|{\phi }_0^{(x)}⟩⟨{\phi }_0^{(z)}|\cos \left(2\pi z\right)|{\phi }_1^{(z)}⟩$$

$$={\left(-1\right)}^n{\Omega }_z$$

$$i{d}_t{c}_{n,m}=$$

$$+(-1{)}^n{T}_z\left(c_{n,m+1}-c_{n,m-1}\right)$$

$$i{d}_t\left[\Psi (x,z,t)\right]=H_S\left[\Psi (x,z,t)\right]$$

$$i{d}_t{c}_{n,m}=$$

$$+T_x\left(c_{n+1,m}-c_{n-1,m}\right)$$

$${\bar{V}}_A(x,z,t)=$$

$$V_z^{A}x\cos \left(\pi x\right)\cos \left(2\pi z\right)\cos ({\omega }_B^{(z)}t).$$

$$i{d}_t{c}_{n,m}$$

$$+T_x\left(c_{n+1,m}-c_{n-1,m}\right)$$

$$+T_x^{A}m\left(c_{n+1,m}+c_{n-1,m}\right)$$

$$+(-1{)}^{n}n{T}_z^A\left(c_{n,m+1}+c_{n,m-1}\right)$$

$$(E-m{c}^2)\varphi$$

$$(E+m{c}^2)\chi$$

$$(E^2-m^2{c}^4)\varphi =$$

$${\varphi }_{\pm }=(1,\pm i{)}^{⊺}{\psi }_{\pm }(x,z)$$

$$\left[\frac{p^2}{2\mu }+\frac{\mu }{2}(x^2+z^2)+L_y\right]{\psi }_{\pm }=(\tilde{E}\mp 1){\psi }_{\pm }$$

$$\tilde{E}=(E^2-m^2{c}^4)/(B{c}^2)$$

$$E=\sqrt{m^2{c}^4+2B{c}^{2}n}$$

Full text

I Introduction

The Dirac equation, unifying quantum mechanics and special relativity, is a major achievement in physics. It has vast implications in several fields, e.g. particle physics where it describes spin-1/2 leptons and their interactions with photons, condensed matter where it is used as a model for several types of quasiparticles, and also in the fast progressing field of topological insulators.

Although the Dirac equation represents a paradigm for modern field theory, this physics remained relatively elusive, restricted to high-energy situations or to exotic materials. Recent developments both in condensed matter and ultracold-atom systems have generated a burst of interest, in particular, the concept of “quantum simulator” opened a new window in the study of Dirac physics (Park et al., 2019; Garreau and Zehnlé, 2017; Lopez-Gonzalez et al., 2014; Qu et al., 2013; Tarruell et al., 2012; Mazza et al., 2012; Zhang et al., 2012; Gerritsma et al., 2011; Lamata et al., 2011; Gerritsma et al., 2010; Longhi, 2010; Dreisow et al., 2010; Witthaut et al., 2011). Inspired by an original insight of Feynman (Feynman, 1982), an analog quantum simulator (as opposed to quantum-computer simulations) is a “simple” and controllable system that can mimic (some aspects of) the behavior of a more complex or less accessible ones. The flexibility of ultracold-atom systems also prompted for innovative ideas like the generation of the so-called “artificial” gauge fields acting on (neutral) atoms that mimic the effects of an electromagnetic field (Dalibard et al., 2011), allowing for the study of quantum magnetism (Aidelsburger et al., 2013), engineered spin-orbit coupling, and topological systems (Galitski and Spielman, 2013). The mixing of these ideas with the physics of ultracold atoms in optical lattices has proven extremely fruitful (Bloch et al., 2008; Dalibard et al., 2011; Bloch et al., 2014; Georgescu et al., 2014; Goldman et al., 2016; Gross and Bloch, 2017; Garreau, 2017).

The Dirac equation is the most complete formulation describing a relativistic charged fermion of spin 1/2. It leads to a spinor-4 description which automatically includes the spin and the electron antiparticle, the positron. In the rest frame of the free particle, which always exists for a massive particle, the spinor-4 components can be separated in two spinor-2 obeying equivalent equations. In the presence of electromagnetic fields, however, the two spinor-2 for a massive particle generally obey different equations, and a full spinor-4 description is necessary. The aim of the present work is the elaboration of a minimal model analog to a 2D Dirac equation in the presence of an artificial gauge field (related to a vector potential $\boldsymbol{A}$ with non-zero rotational). In this case the spinor-4 description is mandatory, and we show that the main characteristics of a Dirac particle are obtained.

In a previous work (Garreau and Zehnlé, 2017) we introduced a general model allowing to mimicking a spinor-4 Dirac Hamiltonian in 1D without magnetic field. In the present work we combine time- and space-tailored optical potentials acting on independent atoms in order to obtain a 2D Dirac Hamiltonian of the form

[TABLE]

where $c$ is the velocity of light, $m_{D}$ and $q$ are the Dirac particle’s mass and charge, $\boldsymbol{p}$ is the momentum and $\boldsymbol{A}$ the vector potential. The $4\times 4$ Dirac matrices constructed from Pauli matrices $\sigma_{i}$ ($i=x,y,z$) are:

[TABLE]

and $\boldsymbol{1}$ is the $2\times 2$ unit matrix. With this we i) generalize this approach to the 2D case, ii) introduce an artificial gauge field, and iii) demonstrate its ability to simulate known properties of the Dirac equation. To the best of our knowledge, such a complete Dirac simulator has not been proposed in the literature, and opens new ways for a deeper exploration of the Dirac physics.

II Spinor-4 Dirac quantum simulator in 2D

We first construct a modulated 2D tilted optical lattice model that mimics a 2D spinor-4 Dirac equation with no field. The model builds on the general approach introduced for the 1D case in Ref. (Garreau and Zehnlé, 2017). A 2D “tilted” (or Wannier-Stark) lattice in the $x,z$ plane 111As it will be seen in what follows, this choice makes the writing in terms of the conventional Pauli matrices possible. We will later introduce an artificial magnetic field in the $y$ direction. is described by the 2D (dimensionless) Hamiltonian

[TABLE]

with $\boldsymbol{p}=p_{x}\mathbf{x}+p_{z}\mathbf{z}$ the momentum in 2D ($\mathbf{x},\mathbf{z}$ are unit vectors in the directions $x,z$), $\boldsymbol{F}=F_{x}\mathbf{x}+F_{z}\mathbf{z}$ a constant force and $V_{L}$ a square lattice formed by orthogonal standing waves

[TABLE]

where space is measured in units of the step $\mathsf{a}$ 222We use sans serif characters to represent dimensionfull quantities, except when no ambiguity is possible, e.g. $\hbar$. of the square lattice $V_{L}$, $\mathsf{a}=\lambda_{L}/2$ if the lattice is formed by counter-propagating beams of wavelength $\lambda_{L}=2\pi/\mathrm{k}_{L}$. Time is measured in units of $\hbar/\mathsf{E}_{R}$ where $\mathsf{E}_{R}=\hbar^{2}\mathrm{k}_{L}^{2}/2\mathrm{M}$ is the atom’s recoil energy ($\mathrm{M}$ is the atom’s mass). With these choices, $m^{*}=\pi^{2}/2$ , $\hbar=1$ and $p_{j}=-i\partial_{j}$ ($j=x,z)$. Since the Hamiltonian Eq. (2) is separable, its eigenstates can be factorized in terms of localized Wannier-Stark (WS) functions $\varphi_{n}^{(x,y)}$ (Wannier, 1960; Zener, 1934) which are exact eigenstates of a tilted one-dimensional lattice (Garreau and Zehnlé, 2017):

[TABLE]

where the integer index $n$ indicates the lattice site along the $x$ direction where the eigenfunction $\varphi_{n}^{(x)}(x)$ is centered (resp. $m$ in the $z$ direction). The energies of the system are then

[TABLE]

where $E_{0}$ is an energy offset with respect to the bottom of the central well $n=m=0$. The energy spacing in directions $x$ and $z$ define the so-called Bloch frequencies, $\omega_{B}^{(x)}=F_{x}$ ($=\mathsf{F}_{x}\mathsf{a}/\hbar$ in dimensionfull units) and $\omega_{B}^{(z)}=F_{z}$, where we intentionally choose $F_{z}\neq F_{x}$. This defines the lowest “Wannier-Stark ladder” of energy levels separated by integer multiples of $\omega_{B}^{(x,z)}$, Eq. (3). According to the potential parameters, there can be “excited” ladders, also obeying Eq. (3), but with higher values of $E_{0}$. As in Ref. (Garreau and Zehnlé, 2017), we assume here that the excited ladders are never populated and that the dynamics occurs only among the lowest-ladder WS states of each site.

The Hamiltonian Eq. (2) is invariant under a translation by an integer multiple $n$ of the lattice step $a=1$ (in dimensionless units) in the $x$ direction provided that the energy is also shifted by $nF_{x}$ (resp. $mF_{z}$ in the $z$ direction), hence the eigenstates are such that

[TABLE]

(resp. $\varphi_{m}^{(z)}(z)=\varphi_{0}^{(z)}(z-m)$ in the $z$ direction).

Our aim is to obtain an effective evolution equation for the system equivalent to a 2D Dirac equation Eq. (1) for a particle of mass $m$

[TABLE]

where $\partial_{t}\equiv\partial/\partial t$, $\left[\Psi\right]$ denotes a spinor-4, and $\boldsymbol{A}=A_{x}\mathbf{x}+A_{z}\mathbf{z}$. The effective speed of light is $c=\hbar\mathrm{c}/\mathsf{a}\mathsf{E}_{R}$, the gauge field $\boldsymbol{A}=q\mathsf{a}\mathbf{A}/\hbar$, the mass $m=\mathsf{E}_{R}\mathsf{a}^{2}m_{D}/\hbar^{2}$, and effective charge is set to unity.

Consider first the case without gauge field ($\boldsymbol{A}=0$). We can induce controlled dynamics (Garreau and Zehnlé, 2017) in the system by adding to $H_{0}$, Eq. (2), a resonant time-dependent perturbation of the form

[TABLE]

The utility of the term proportional to $\cos(\pi x)$ with spatial periodicity $2$ will appear below.

The general solution of the corresponding Schrödinger equation can be decomposed on the eigenbasis $\Phi_{n,m}(x,z)$ :

[TABLE]

where $c_{n,m}$ are complex amplitudes. Taking into account the Hamiltonian $H_{0}+\bar{V}$, one obtains the following set of coupled differential equations for the amplitudes

[TABLE]

where $d_{t}\equiv d/dt$ and we simplified the notation $\left|\varphi_{n}^{(x)}\right\rangle\otimes\left|\varphi_{m}^{(z)}\right\rangle\equiv\left|n,m\right\rangle$. To the dominant order, one can neglect non-resonant terms and keep nearest neighbor couplings only, i.e take only $j=\pm 1$ and $k=\pm 1$ in Eq. (7) (see Appendix A).  The time modulation at the Bloch frequencies $\omega_{B}^{(x)}$ (resp. $\omega_{B}^{(z)}$) resonantly couples first-neighbor WS states in the $x$ direction (resp. $z$ direction). Note that as we assume $\omega_{B}^{(x)}\neq\omega_{B}^{(z)}$ resonant couplings in both direction $x$ and $z$ can be controlled independently.

Taking into account the choice of perturbation $\bar{V}$ of Eq. (6), the couplings in Eq. (7) are proportional to the spatial overlap amplitudes

[TABLE]

[see. Eq. (4)] in the $x$ direction, and

[TABLE]

in the $z$ direction. The $(-1)^{n}$ contribution comes from the spatial modulation in $\cos(\pi x)$ in Eq. (6) which has a period of two lattice steps [see Eq. (23) in Appendix A] and introduces a distinction between even and odd sites in the $x$ direction.

From Eq. (7), one obtains a set of differential equations with nearest-neighbor couplings.

[TABLE]

where $T_{x}=-iV_{x}\Omega_{x}/2$ and $T_{z}=-iV_{z}\Omega_{z}/2$.

Assuming that the wave packet is large and smooth at the scale of the lattice step one can take the continuous limit of Eq. (8), and associate the discrete amplitudes $c_{n,m}(t)$ to two continuous functions. Since the site parity dependence in Eq. (8) must be taken into account, we introduce the functions $s_{ee}(x,z,t)$ which is the envelope of $c_{n,m}$ for $n,m$ even, $s_{eo}(x,z,t)$ which is the envelope of $c_{n,m}$ for $n$ even and $m$ odd, and analogously for $s_{oe}$ and $s_{oo}$. These functions can be arranged as components of a spinor-4 $[\Psi]=\left(s_{ee},s_{oo},s_{eo},s_{oe}\right){}^{\intercal}$ describing 4 coupled sub-lattices, which, from Eq. (8), obey a Hamiltonian equation

[TABLE]

where $H_{S}$ can be easily shown to be analogous to a 2D Dirac Hamiltonian of the form $H_{S}=c_{S}\mathbf{\boldsymbol{\alpha}}\cdot\boldsymbol{p}$ with an effective velocity of light $c_{S}=V_{x}\Omega_{x}=V_{z}\Omega_{z}$ (these terms can be made equal by adequately tuning the modulation amplitudes $V_{x,z}$, but one can also create an anisotropic model with $c_{S_{x}}\neq c_{S_{z}}$, with the possibility of an effective violation of Lorentz invariance (Park et al., 2019)). The above system can be broken into two equivalent set of equations (Garreau and Zehnlé, 2017) which correspond to the well-known massless Weyl spinor-2 fermion. The twofold-degenerated dispersion relation deduced from Eq. (9) is $\omega(k)=\pm\left(\Omega_{x}^{2}k_{x}^{2}+\Omega_{z}^{2}k_{z}^{2}\right)^{1/2}$ and corresponds, as it could be expected, to a Dirac cone.

The Hamiltonian of a massive particle can be generated by adding a static perturbation $V_{0}\cos(\pi x)\cos(\pi z)$ to Eq. (6). Then, neglecting derivatives $\partial_{x,z}^{2}c_{m,n}$ or higher, Eq. (8) takes the form

[TABLE]

where $T_{0}=V_{0}\left\langle\varphi_{0}^{(z)}\right|\cos(\pi z)\left|\varphi_{0}^{(z)}\right\rangle\left\langle\varphi_{0}^{(x)}\right|\cos(\pi x)\left|\varphi_{0}^{(x)}\right\rangle$. In the continuous limit one obtains a 2D Dirac Hamiltonian for a finite-mass particle $H_{S}=c_{S}\mathbf{\boldsymbol{\alpha}\cdot p}+\beta mc^{2}$ with an effective rest energy $mc^{2}=T_{0}$. This term opens a gap of width $2mc^{2}$ in the the dispersion relation, separating “particle” and “antiparticle” states, exactly as in the Dirac equation.

The approximations used in construction of the above model are well understood and controlled. As the system is to a very good approximation a closed one (even experimentally), there is only a small broadening of the resonances, which make these resonances highly selective. Next-to-neighbor couplings are out of resonance and are thus negligible, and the same is true for intrawell couplings between different ladders. For particular values of the potential parameters $(V_{x,z},F_{x,z})$ “accidental” resonances can occur, but these exotic situations are not considered here. The accuracy of our approach can also be checked by comparing the results of the model to a simulation of the full Schrödinger equation, which takes into account all existing couplings, as we have done with the simpler, but conceptually equivalent, 1D model in Ref. (Garreau and Zehnlé, 2017). For the very same reason (the sharpness of the resonances) the system is not expected to be particularly sensitive to heating caused by experimental noise.

III Dirac equation with an artificial gauge field

In this section, we show how the analogous of a “non-trivial” (i.e. of non-zero rotational) vector potential $\boldsymbol{A}=A_{x}\mathbf{x}+A_{z}\mathbf{z}$ with $A_{x}\propto z$ and $A_{z}\propto x$, can introduced in the system. The Dirac Hamiltonian of Eq. (1) can be realized by adding to Eq. (6) an additional perturbation

[TABLE]

This term can be treated in the same way as in Sec. II and leads to

[TABLE]

where $T_{z}^{A}=\frac{1}{2}\Omega_{z}V_{z}^{A}$, $T_{x}^{A}=\frac{1}{2}\Omega_{x}V_{x}^{A}$. The additional terms on the last two lines in Eq. (12) are due to the potential Eq. (11) which generates the slopes proportional to $m$ (resp. $n$) in direction $z$ (resp. $x$) 333Experimentally, such linear terms can be obtained by superposing two laser fields of frequencies separated by $\delta\omega$, which generate a beat note varying in space as $\sin\left(\delta\omega x/c\right)\sim\delta\omega x/c$; if $c/\delta\omega$ is large enough compared to the lattice wavelength the linear approximation can be valid over a large number of neighbor sites where the wavepacket is concentrated..

In the continuous limit, and neglecting second order and higher derivatives of the spinor components, we obtain the Hamiltonian Eq. (1), with an artificial gauge field $\boldsymbol{A}=-2(zT_{x}^{A}\mathbf{x}+xT_{z}^{A}\mathbf{z})/c_{S}$ (cf. App. A). The symmetric gauge considered in the next section can be realized by tuning the modulation amplitudes in $\bar{V}_{A}$ so that $T_{x}^{A}=-T_{z}^{A}$, so that $\boldsymbol{A}=2T_{z}^{A}(z\mathbf{x}-x\mathbf{z})/c_{S}$, corresponding to a uniform magnetic field in the $y$ direction $\boldsymbol{B}=4T_{z}^{A}\mathbf{y}/c_{S}$ 444The Landau gauge can be obtained by setting $T_{z}^{A}=0$ (or $T_{x}^{A}=0$)..

IV Quantum simulation of Dirac physics

We briefly recall well-known theoretical results of Eq. (5) for a spinor-4 Dirac particle in the presence of a symmetric vector potential $\boldsymbol{A}=2T_{z}^{A}(z\mathbf{x}-x\mathbf{z})/c_{S}$. We will compare them with the numerical results of the discrete model Eq. (12) and show that they are in very good agreement, proving that the discrete model reproduces to a good level of accuracy the Dirac physics.

In the following, we write the spinor-4 as \left[\Psi\right]=\left(\begin{array}[]{c}\phi\\ \chi\end{array}\right) where $\phi,\chi$ are spinor-2s. With this convention, the Dirac equation can be decomposed in 2 coupled equations

[TABLE]

These equations are symmetric under the transposition E,$$\phi, $\chi$ $\longleftrightarrow$ -E,$$-\chi, $\phi$ so that the negative energy states are easily deduced from their positive energy counterparts.

Eliminating $\chi$ from Eqs. (14) gives $(E^{2}-m^{2}c^{4})\phi=c\left[\boldsymbol{\sigma}\cdot\left(\boldsymbol{p}-\boldsymbol{A}\right)\right]^{2}\phi$, and, after some straightforward algebra, one finds

[TABLE]

where $L_{y}=(zp_{x}-xp_{z})$ is the angular momentum component in the $y$ direction, with a “diamagnetic” term proportional to $x^{2}+z^{2}$ and “paramagnetic” terms of the type $\boldsymbol{L}\cdot\boldsymbol{B}$ and $\boldsymbol{\sigma}\cdot\boldsymbol{B}$. This equation has “spin up” $\phi_{+}$ and “spin down” $\phi_{-}$ spinor-2 555We label these solutions “spin up” and “spin down” (with quotation marks) because they are proportional to the eigenstates $(1,\pm i)^{\intercal}/\sqrt{2}$ of $\sigma_{y}$ (remember that $\boldsymbol{B}\parallel\mathbf{y}$). This is not to be confused with spinor-4-particle’s spin components. solutions

[TABLE]

where the functions $\psi_{\pm}(x,z)$ obey the following equation

[TABLE]

with

[TABLE]

and $\mu=B/2$, which strongly evokes a 2D harmonic oscillator of mass $\mu$ and natural frequency $\omega=1$ in a magnetic field.

Solutions of Eq. (17) are the well-known “Landau levels”, with spectrum

[TABLE]

where $n$ (not to be confused with the site index $n$, which does not appear in the present section) is a strictly positive natural number for “spin up” solution $\phi_{+}$ and a natural number for “spin down” $\phi_{-}$. For each $n$ the corresponding energy is infinitely degenerated.

The solutions of Eq. (17) are also well known, and together with Eq. (16) and Eq. (14) lead to solutions for the full eigenspinors $\left[\Psi\right]$. In polar coordinates $x=r\cos\varphi$, $z=r\sin\varphi$ and $\rho=r/\ell$, where $\ell=\sqrt{2/B}$ is an effective dimensionless “magnetic length”, an example of solution for the fundamental level $E=mc^{2}$ is

[TABLE]

where $\mathscr{N}$ is a normalization factor. There is a corresponding negative energy $E=-mc^{2}$ counterpart with $\phi\rightarrow\phi$ and $\chi\rightarrow-\chi$. Examples of eigenstates belonging to the $n=1$ excited state family $E=\left(m^{2}c^{4}+2Bc^{2}\right)^{1/2}$

[TABLE]

or

[TABLE]

where $\beta=-2c/\left[\ell(mc^{2}+E)\right]$ 666Note that these states are also eigenstates of the angular momentum operator $L_{y}+\varSigma_{y}$ , where \varSigma_{y}=\frac{1}{2}\left(\begin{array}[]{cc}\sigma_{y}&0\\ 0&\sigma_{y}\end{array}\right) is the spin operator. .

Numerical simulations of the discrete model reproduce the characteristics of the above solutions of the Dirac equation. According to Eq. (19), plateaus are obtained at integer values $n=\left(E^{2}-m^{2}c^{4}\right)/2Bc^{2}$. As shown in Fig. 1, the numerical solution of Eq. (12) displays such a behavior, with plateaus (red disks) appearing as expected at integer values. The simulation also shows eigenenergies which fall at non-integer values (blue circles); we checked that these are edge states due to finite size effects, as the numerical simulation is performed in a square box containing finite number of lattice sites (80$\times$80) 777We numerically asserted that a $150\times 150$ lattice gives essentially the same results..

Figure 2 shows a numerical example of the the four components of the eigenspinor $\left[\Psi\right]=\left[s_{ee},s_{oo},s_{eo},s_{oe}\right]^{\intercal}$ for the bulk state, corresponding to the analytical solution of the Dirac excited mode of Eq. (20); injecting the magnetic length $\ell=10$ in Eq. (20), we get a quantitative agreement with the simulation result. This confirms the ability of our model to mimic the behavior of a spinor-4 Dirac particle. Figure 3 is the false-colors representation of the border state shown in Fig. 1, clearly displaying its chiral edge-state nature which evokes Quantum Hall Effect states. Recently, the existence of chiral Quantum-Hall-Effect-like states in acoustic systems received great attention (Peri et al., 2019; Wen et al., 2019), opening prospects for the simulation of Dirac physics in such systems. However, in contrast with our spinor-4 simulator, these acoustic simulators are described by a Weyl equation for a spinor-2. Finally, Fig. 4 display the chiral vector fields of the true Dirac spinor-4 current of probability for the bulk and for the edge states (see App. B), further evidencing the ability of our system to reproduce spinor-4 Dirac physics.

V Conclusion

The present work demonstrates the quantum simulation a Dirac particle in the presence of a magnetic field, putting together quantum simulations of “exotic” dynamics (from the point of view of low-energy physics) and artificial gauge fields. This is intrinsically expected to lead to topological systems, as the Dirac cone is one of the main (and simplest) examples of topology in physics, and opens new ways to explore exciting new possibilities in this very active field. Other interesting prospects would be the inclusion of particle-particle interactions and/or disorder in this analog quantum simulations. These are intriguing and highly challenging tasks, beyond the scope of the present work.

Acknowledgements.

This work was supported by Agence Nationale de la Recherche through Research Grants K-BEC No. ANR-13-BS04-0001-01 and MANYLOK No. ANR-18-CE30-0017, the Labex CEMPI (Grant No. ANR-11-LABX-0007-01), the Ministry of Higher Education and Research, Hauts-de-France Council and European Regional Development Fund (ERDF) through the Contrat de Projets Etat-Region (CPER Photonics for Society, P4S).

Appendix A Derivation of the Dirac equation

This appendix discusses in more detail the developments leading to Eqs. (8) and (10) for the site amplitudes $c_{n,m}$ and explains how the equations for their continuous envelopes can be put in an form analogous to the Dirac equation. The derivation of the Dirac equation in presence of a gauge field of Sec. III follows the same steps and is also sketched below.

Injecting the general wave function $\psi(x,z,t)=\sum_{n,m}c_{n,m}(t)\exp\left(-i(n\omega_{B}^{(x)}+m\omega_{B}^{(z)})t\right)\varphi_{n}^{(x)}\varphi_{m}^{(z)}$ in the Schrödinger equation for the Hamiltonian $H=H_{0}+\widetilde{V}$ Eq. (6), one obtains the following set of coupled differential equations

[TABLE]

To the dominant order, neglecting non-resonant terms, and keeping nearest neighbor couplings 888This is justified by the rapid vanishing of overlapping integrals between WS states localized in different sites. only $j,k=\pm 1$ in Eq. (21) one gets

[TABLE]

The overlap integrals can be simplified using the Wannier-Stark state translation properties as in Sec. II. For instance:

[TABLE]

defining the coupling parameter $\Omega_{z}$ [the $(-1)^{n}$ factor come from the perturbation in $\cos(\pi x)$]. The same calculation can be done for overlaps $\left\langle n,m\right|\cos(2\pi x)\left|n\pm 1,m\right\rangle$ and gives

[TABLE]

These results lead to Eq. (8).

When a static perturbation $V_{0}\cos(\pi x)\cos(\pi z)$ is added to Eq. (6), we have, in the resonant approximation, an additional term corresponding to a “self-coupling” contribution

[TABLE]

which is to be added to the R.H.S. of Eq. (8). The overlap integral is then

[TABLE]

which in turn gives Eq. (10).

In taking the continuous limit, we have to account for the factors $(-1)^{n}$ and $(-1)^{n+m}$ in Eq. (8) or Eq. (10) which lead to parity-dependent smooth envelopes. Taking as an example $n,m$ even, Eq. (22) shows that $c_{n,m}$ is coupled to the next-neighbors amplitudes with $n$ and $m$ of opposite parities, i.e. $(n\pm 1,m)$ or $(n,m\pm 1)$

[TABLE]

with the coupling constants $T_{j}=-iV_{j}\Omega_{j}/2$ ($j=x,z$) defined as in Sec. II. Taking the continuous limit $c_{n+1,m}-c_{n-1,m}\simeq 2\partial_{x}c_{oe}$ and $c_{n,m+1}-c_{n,m-1}\simeq 2\partial_{z}c_{eo}$, one has

[TABLE]

where $p_{x}=-i\partial_{x}$, $p_{z}=-i\partial_{z}$. The calculation is analogous for the remaining smooth components $s_{oo},$ $s_{oe}$ and $s_{eo}$.

The introduction of the artificial gauge field in Sec. III follows the same lines as above. The additional modulation terms [see Eq. (11)] generate extra terms in the evolution equation for the amplitudes $c_{n,m}.$ For instance, the term proportional to $V_{x}^{A}z\cos\left(2\pi x\right)\cos(\omega_{B}^{(x)}t)$ gives in the R.H.S. of the evolution equation for amplitudes $c_{n,m}$ the following contribution:

[TABLE]

with overlap integrals that can then be written as

[TABLE]

where the constant $\bar{z}$ is a small offset which can be canceled by a $z$ translation. Equation (12) is then easily obtained, as well as its continuous limit for the envelopes of $c_{ab}$ ($a,b=e,o$), assuming they are smooth enough to allow neglecting second order derivatives (for instance $c_{n+1,m}+c_{n-1,m}\simeq 2s_{ee}(x,z,t)$ for $n$ odd and $m$ even, and so on) .

Appendix B Dirac current

The Dirac equation leads to a standard probability current continuity equation $\partial_{t}\left(\Psi^{\dagger}\Psi\right)+\mathbf{\boldsymbol{\nabla}}\cdot\boldsymbol{j}=0$ with the following current components

[TABLE]

where $\left[\Psi\right]=\left[\phi;\chi\right]$ as in Sec. IV. For instance, taking the of Eq. (20) the current is

[TABLE]

which is in excellent agreement with Fig. 4(a).