Dirac Equation For Cold Atoms In Artificial Curved Spacetimes

We argue that the Fermi-Hubbard Hamiltonian describing the physics of ultracold atoms on optical lattices in the presence of artificial non-Abelian gauge fields, is exactly equivalent to the gauge theory Hamiltonian describing Dirac fermions in the lattice. We show that it is possible to couple the Dirac fermions to an"artificial"gravitational field, i.e. to consider the Dirac physics in a curved spacetime. We identify the special class of spacetime metrics that admit a simple realization in terms of a Fermi-Hubbard model subjected to an artificial SU(2) field, corresponding to position dependent hopping matrices. As an example, we discuss in more detail the physics of the 2+1D Rindler metric, its possible experimental realization and detection.


Introduction
The study of ultracold quantum matter in artificially designed external gauge fields is one of the most rapidly developing areas of the physics of ultracold atoms [1]- [5]. Originally, these studies arose from investigations related to the response of superfluids, such as Bose-Einstein condensates (BEC), to rotation (cf [6,7]). On the one hand, rotation induces quantized vortices and/or vortex Abrikosov lattices [8]. On the other, its effects are equivalent to those of an artificial constant (Abelian) magnetic field. The latter analogy has led to the idea of realizing strongly correlated quantum liquids, such as the celebrated Laughlin liquid of the fractional quantum Hall effect (FQHE) (cf [9]) by means of rapid rotation [10,11]. Unfortunately, reaching the FQHE regime with rotation is experimentally very difficult; it has been achieved recently, but in a system of only 1 < N 10 atoms in rotating microtraps (lattice site potential wells) [12,13]. Several researchers have proposed, thus, alternative approaches involving, for instance, laserinduced gauge fields that employ dark states in three-level systems [14], or laser-induced gauge fields acting on atoms confined in an optical lattice (OL) [15]- [17]. Other approaches concerned rotating OLs [18,19], or interactions of lattice atoms embedded in a rotating BEC [20]. Interestingly, the proposals involving lasers can be relatively straightforwardly generalized to particles possessing internal 'color' states subjected to artificial non-Abelian gauge fields [21,22].
Over the last few years, there has been a large number of works reexamining these ideas and proposing experimental realization within the reach of the present state of the art. The NIST group employed an approach similar to [14] and used Raman (Bragg) transitions in sodium to realize experimentally the first nonzero constant vector potential (corresponding to zero 'artificial' magnetic field) [23]. Note that when the gauge symmetry is provided by Nature, 4 artificial non-Abelian fields in lattices is exactly the HLGT version of the Dirac Hamiltonian. That is to say, in the scheme we consider, all excitations of the fermion field in the lattice are Dirac-like, not only the low-lying ones. The advantage of this point of view is that it allows one to couple the simulated Dirac fields to external fields, or quantum fields in a straightforward way, which is not so clear in artificial graphene-like systems.
In particular, we will show that it is possible to couple Dirac fermions to an artificial gravitational field, that is, to consider the Dirac equation in curved space. We will identify and focus on a special class of spacetime metrics that admit a simple formulation of the Dirac lattice Hamiltonian in terms of a Fermi-Hubbard model subjected to an artificial SU(2) field [61], corresponding to tunneling matrices with position-dependent overall hopping rate.
We will not consider here the fermion doubling problem as it is inessential for the discussion of gravitational effects 6 . Feasible solutions of such a problem in OL simulations of (3 + 1) fermions have been very recently proposed in [66] and [31].
The motivation for simulating the Dirac equation in curved spacetime is at least twofold. The propagation of fermions in curved space is not experimentally accessible in highenergy physics/cosmology as the (piece of the) Universe that we are nowadays able to probe is practically flat [78]. The lattice realization proposed in the paper is an appropriate description of the continuum physics for processes occurring over scales of several lattice spaces. The realization with cold atoms in OL provides us with 'laptop experiments', paving the way for observing exotic effects such as the thermalization theorem, also known as the Fulling-Davies-Unruh effect 7 [79]- [81] (for a review on the subject, see [82]). Roughly speaking, it states that an accelerated observer perceives the Minkowski vacuum as a thermal bath. The origin of the above phenomenon, as well as of Hawking radiation [83], is the existence of a nontrivial Bogoliubov transformation between the Minkowski vacuum and a spacetime with a horizon. The effect of the latter, in quantum information language, is to 'trace-out' part of the system, giving rise to a thermal reduced density matrix. Another exotic gravitational effect to observe might be the curved spacetime version of Zitterbewegung [84], to compare with its recent observation achieved in flat space [54].
On other hand, simulating the Dirac equation in curved spacetime gives rise to the possibility of modeling the analogues of graphene ripples on a flat, square lattice. Such a possibility could help in disentangling the action of ripples (which admit natural 'geometric' interpretations [85] or gauge field interpretations, e.g. [86]) on the carrier density from other contributions.
Different approaches to quantum simulation include the exact simulation of the dynamics of strongly correlated systems by unitary gates [87], or the the creation of interesting strongly correlated states that are ground states of interesting Hamiltonians [88].
Let us stress that the models we propose in this paper are experimentally challenging, but feasible. Indeed, non-Abelian OLs have already been realized experimentally and reported in [33] where the spin-orbit coupling is simulated.
The paper is organized as follows. In section 2, we demonstrate how the lattice Dirac Hamiltonian in flat spacetime [89] can be obtained as the discretization of the continuum 6 The mixture of the two Dirac points due to the gravitational background potentially induces a gauge field coupling to the composite fermions (for a graphene-like lattice, see for example [42], [75]- [77]). However, the field contribution is relevant in the presence of conical singularities, disclination or dislocation in the graphene language, while it is negligible when the metric is smooth. 7 Also, the analogue of the Unruh effect in superfluid helium has been discussed in [49, section 6.15].
Dirac Hamiltonian H . This Hamiltonian is identified then as the Hamiltonian of the SU(2) Fermi-Hubbard model considered for instance in [61], where the hopping terms are given by the Pauli matrices σ i , i = x, y, times a constant hopping rate J . Although the discretization can be done at this stage in many different ways, we obtain it by coarse graining the 'symmetric' formulation of H in the fermion field ψ and its conjugate ψ † . This strategy turns out to be very convenient in order to find the lattice version of the Dirac Hamiltonian in curved spacetime, discussed in section 3. In order to favor the accessibility to a broad audience, we provide a physically based introduction to the Dirac equation in curved spacetime in section 3.1, minimizing as much as possible the technical details and fixing the notation. As the first exercise, after reviewing the form of the continuum Hamiltonian on a generic manifold (admitting a 'time' isometry), we successfully apply the above strategy to the interesting case of the (2 + 1) Rindler Universe [90] in section 3.2. The resulting lattice Hamiltonian differs from the Hubbard Hamiltonian of the flat case in the following: the hopping rates exhibit a linear dependence with position. Such a simple form is due to the cancellation of the spin-connection, once ψ and ψ † are treated in the same manner: as a consequence, the hopping matrices need not be locally rotated. We find that this simplification happens in fact for any static spacetime. For details of the derivation, see appendix B. We discuss under which conditions and how the fermion propagation in such spacetimes can be engineered and detected on OLs in sections 4 and 5 that are the heart of this work. There, several different experimental ways of implementing the Hubbard model of interest are presented. We propose the density of states as the relevant observable to capture the Dirac physics. It is a measurable quantity both in graphene [91]- [94] and in OL [95]- [102] experiments. We compute the theoretical value of the density of states analytically in the continuum limit using perturbation theory.
We present our concluding remarks and discuss further developments in section 6.

Dirac Hamiltonian on a lattice
In this section we show, by discretization of the spatial coordinates, that Dirac's Hamiltonian in (2 + 1) dimensions is a hopping Hamiltonian with non-Abelian tunneling matrices. In order to make the subject accessible to a nonspecialized reader, we first derive the Dirac Hamiltonian in a pedagogical way starting with the Dirac equation. For more details, we refer the reader to his favorite quantum field theory textbook (see for instance [103].) As originally shown by Dirac, in any spacetime dimensions the appropriate Lorentz covariant and local equation of motion for a spin-1/2 field is linear in the spacetime derivative. For a massless fermion in d spatial dimensions (in this paper we discuss the simulation of massless fermions only), it takes the simple form where a = 0, . . . , d, a = 0 indicates the time direction, ψ is a two-component spinor and the γ a are the Dirac gamma matrices. From the mathematical point of view, the gamma matrices close a Clifford algebra and their commutators γ ab ≡ 1 2 [γ a , γ b ] form an irreducible representation of the (d + 1) Lorentz group SO (1, d). Such properties are encoded in the key relation for the anticommutator of the γ a , where η ab is the (d + 1) Minkowski metric, η ab = Diag(−1, 1, . . . , 1). In particular, the above relation implies that γ a ≡ η ab γ b = (γ a ) −1 and that the indices of the gammas behave as 6 spacetime indices, raised and lowered with the metric η ab . Using this property it is easy to verify that ψ satisfies the Klein-Gordon equation for a massless particle, implying the relativistic dispersion relation E = | p|. In what follows, we adopt the mostly positive convention, = 1. Thus, γ 0 is anti-Hermitian and squares to minus the identity, while the spatial gammas square to one and are Hermitian. In (2 + 1) dimensions, the case we are interested in, the γ a can be expressed in terms of the Pauli matrices. A possible choice is γ 0 = iσ z , γ 1 = σ y and γ 2 = −σ x , which is consistent with (2) and with the condition γ 0 γ 1 γ 2 = −1 that is specific to 2 + 1.
As the Dirac equation is linear, we can exploit the analogy to Schrödinger's equation to obtain the Hamiltonian density H. The time evolution of the field ψ is from which we can easily read off the single-particle Hamiltonian. As was first pointed out in [89], in (2 + 1) dimensions and on a lattice of spacing the discretized version of equation (3) is with ψ(x) = ψ(m , n ) = ψ m,n . Integrating to obtain the Hamiltonian operator H = dx dy ψ † Hψ, we obtain on the lattice It is easy to see that (5) is nothing more than a Fermi-Hubbard Hamiltonian with non-Abelian hopping matrices where the interactions and the effect of the trap have been neglected. In the notation of [61], we have U x = iγ 0 γ 1 = iσ 1 and U y = iγ 0 γ 2 = iσ 2 (which implies that γ 0 = iσ 3 ) and the Abelian flux = 0.
For further reference, we note that the lattice Hamiltonian of (5) can be alternatively obtained by computing with the substitution of spatial derivatives with finite differences over one lattice spacing , ∂ x ψ → ψ m+1,n −ψ m,n and ∂ y ψ → ψ m,n+1 −ψ m,n .

Dirac Hamiltonian in curved spacetime
We will apply the same strategy as in the previous section to derive the lattice formulation of the Dirac Hamiltonian in curved spacetime. In order to streamline the presentation, we first review the formulation of the Dirac equation on a generic spacetime and single out the form of the Hamiltonian density. As a second step, we show how such density can be integrated to give a consistent and conserved Hamiltonian charge in the case of a Rindler Universe, for both the continuum and lattice versions. The formulation of the continuum Hamiltonian and its lattice companion in the generic case is analyzed in the appendix.

Dirac equation in curved spacetime
In this section, we review the formulation of the Dirac equation on a curved spacetime.
Although some technical complications are unavoidable, we will make use of the analogy between coupling the fermions to gravity and to a gauge field. Again, the guiding principle that originally led Fock to the formulation of the equation is covariance. Indeed, in the flat case the Dirac equation is appropriate because γ a ∂ a ψ transforms as a spinor under global Poincaré transformations. However, coupling the field to gravity is equivalent to asking for a modified equation transforming as a spinor under local Poincaré transformations. Roughly speaking, it is like considering the gauge theory of the Poincaré group: two gauge fields enter the game. One is the spin-connection w ab (x) = w ab µ (x) dx µ that gauges the Lorentz group, i.e. rotations and boosts. The index µ (Greek letters) runs over the (d + 1) spacetime coordinates {t, x 1 , . . . , x d } and for reasons that will be clear soon, is called the curved index, whereas the indices a, b (Latin letters) should be now regarded as internal indices and are usually named as flat. From the practical point of view, the introduction of w ab µ amounts to replacement of the ordinary derivative with the covariant derivative ∂ µ ψ → D µ ψ ≡ ∂ µ + 1 4 w ab µ γ ab ψ, and w ab µ has to transform as an ordinary (but noncompact) non-Abelian gauge group to compensate for the effect of a local infinitesimal Lorentz transformation (x) ab γ ab on the spinor.
The other is the vielbein e a (x) = e a µ (x) dx µ that defines the local translations e a ≡ e µ a ∂ µ (where e µ a = (e a µ ) −1 is a (d + 1) square matrix, according to (7)). The vielbein can be regarded as the analogue of the Cartesian reference frame, and tells us how the orthonormal basis is oriented in each point of the spacetime. Mathematically, this means that e a µ (x)η ab e b ν (x) = g µν (x), or equivalently where g µν (x) is the spacetime metric at the point x, which is used to raise and lower the curved indices as η ab does on the flat ones (as η ab = η ab = (η ab ) −1 , it follows that g µν = (g µν ) −1 ). Roughly speaking, it can be said that vielbein translates curved indices in flat indices and vice versa. Hence, to assign the vielbein is equivalent to fixing the notion of distance and, upon reparameterization, also the contrary is true. Actually, by assigning the orthonormal frame in each point, the choice of e a µ determines the local Lorentz transformation relating the frames in two different points. This principle is known as parallel transport and implies that the spinconnection is not an independent field but exists as a unique choice compatible with the metric g µν . Mathematically, it translates into the requirement of a covariantly constant e a µ (see for example [104]) which allows us to compute w ab µ explicitly. The square brackets indicate that the indices µ and ν are antisymmetrized.
In view of the above considerations, it should be natural to conclude that the correct covariantized Dirac equation is of the form where γ µ are the curved spacetime gamma matrices, γ µ (x) = e a µ (x)γ a .
Let us now focus on (2 + 1) dimensions. By separating the time component from the spatial as γ µ = (γ µ ) −1 . Hence, we have identified the Hamiltonian density for a Dirac fermion on a curved manifold. In the next section, we will show how, after a proper integration over the space, the corresponding Hamiltonian operator can be discretized on the lattice for the paradigmatic example of the Rindler Universe, in full analogy to the flat case. The implementation of the resulting model on an OL will be discussed in section 4.

An example: the Dirac Hamiltonian in a Rindler Universe
As an example, we will derive the Dirac Hamiltonian in a particularly simple spacetime. Let us start by recalling that the Hamiltonian operator in curved spacetime is not necessarily a conserved and well-defined quantity, as the notion corresponding to time translational invariance in flat space is not generic. When it is the case, i.e. technically when a time-like isometry (Killing vector) exists, H is the integral of the Hamiltonian density on a time-like hypersurface. In (2 + 1) dimensions we have The differential is the volume element on a time-like slice and includes the determinant of the metric, d 2 = √ −g dx dy. In a more down to earth fashion, it can be proven that in the presence of generalized time translation there exists (in any spacetime dimension) a choice of coordinates such that the metric does not depend on time, ∂ t g µν = 0 and, as a consequence, the measure and the Hamiltonian density are separately independent of time. In these coordinates the symmetry generated by H is proportional to ∂ t , with the 'clock pace' depending on the position.
Maybe the simplest although nontrivial spacetime enjoying such a property is given by the Rindler Universe. In Rindler space, the metric takes the form ds 2 = −(ax) 2 dt 2 + dx 2 + dy 2 .
We are interested in this metric for two reasons. Firstly, it is the metric seen by an observer in constant acceleration in flat spacetime. For instance, it could have implications for earth-dwelling detectors observing cosmic background radiation [105]. Secondly, it is the nearhorizon metric of a Schwarzschild black hole. The Rindler metric suggests the choice of driebein e 0 = |ax| dt, e 1 = dx, e 2 = dy. Using (9), we may compute the spin connection, whose only nonvanishing component is w 01 t = a x |x| . The Dirac equation (11) greatly simplifies in this spacetime, where we have used that γ 0 γ 1 γ 2 = −1, which holds only in (2 + 1) dimensions. In what follows, we adopt the gamma matrices representation choice In order to carry out the discretization analogous to how it is done for the case of a gauge theory, we note that the Hamiltonian can be written in terms of the Hamiltonian density (14) as In this symmetric form, the spin-connection term disappears and it turns up again on integrating by parts. The discretized version of the Hamiltonian (15) is simply obtained by the substitution ∂ x ψ → ψ m+1,n −ψ m,n and ∂ y ψ → ψ m,n+1 −ψ m,n , with x = m and y = n. One readily obtains Therefore, a lattice with hopping matrices, J i U i = imσ i , i = x, y, growing linearly in the x-direction gives an appropriate description of a free massless fermion in Rindler spacetime. Such a Hamiltonian can be, in principle, implemented in an OL. This problem will be tackled in the next section.

Dirac Hamiltonian in curved spaces and optical lattices
In view of the explicit realization in an OL, we focus on the classes of spacetime where the massless fermion propagation can be described by the Hubbard model of the form As discussed in full detail in appendix B, we find that the lattice Hamiltonian reduces to this simple form only if the metric is static. In (2 + 1) dimensions it is equivalent to saying that it can be written (in a certain coordinate system) in a diagonal form as With this choice of coordinates the hopping rate is simply given by J mn = e (x m ,y n ) . This is not the only requirement to be satisfied in order to reproduce such propagation on an OL, however. Indeed, it is important to note that, even if the Dirac Hamiltonian written in the symmetric form (17) is the same for any choice of the function f in (18), the corresponding Hamiltonian system is distinct for each metric as the canonical momentum is and depends explicitly on f . Now, in OL experiments, the canonical momentum is fixed by the anticommutation relation to be simply iψ † that implies that f = e . Thus the cold fermions in the OL simulate the propagation of massless fermions in a metric of the form 8 ds 2 = −e 2 dt 2 + dx 2 + dy 2 .
The Rindler metric (13) corresponds to = log(ax). 8 It is worth noting that for a generic function the metric is curved and not Weyl-invariant, cf [106].
Before moving to the explicit implementation of the Hamiltonian (17) in OLs, let us briefly discuss to what extent it is a good description of the dynamics of massless fermions in a spacetime given by (19). There are two kinds of limitations. On the one hand, due to finite size of the OL we are able to cover only a finite portion of spacetime. On the other hand, the lattice approximation is valid when the metric, or the function , is sufficiently smooth and slowly varying over one lattice space . It is worth noting that the second problem can be circumvented by using techniques from lattice gauge theory (the continuum limit is obtained by extrapolation).
These two limitations should not obscure the physical content, however, as OLs with up to 300 × 300 sites can be achieved. This implies that the overall variation of the metric over the lattice can be of the order of one. Further considerations in the case of Rindler space are given in section 6.

Experimental realization
In this section, we discuss briefly how one can achieve the appropriate site-dependent hopping rate J mn in an OL. We propose two different techniques.
Recently it was proposed in [31,107] that the intensity of the hopping can be tailored almost at will by considering bichromatic spin-independent superlattices that trap the hyperfine states of alkali bosonic or fermionic atoms. The split of the hyperfine levels is controlled by a magnetic field. The hopping between neighboring Zeeman sublevels of the F (lower) hyperfine manifold, i.e. our 'electrons', is induced via adiabatic elimination of an intermediate F = F ± 1 (upper) hyperfine manifold coupled to F via an off-resonant Raman transition. For details of the scheme for fermionic 40 K atoms, see the original proposal [31]. 40 K have F = 9/2 and allow us, in principle, to simulate any pseudo-spin F 9/2, employing splitting of the Zeeman sublevels in a magnetic field and optical pumping to the relevant sublevels. For the purposes of the present paper it is sufficient to have F = 1/2 or alternatively to use from the very beginning atoms with F = 1/2 in the ground state manifold, such as 6 Li.
To be more concrete, consider a 6 Li Fermi gas loaded on a 2D square OL of size L × L, where the relevant information about our quantum simulator is encoded in the Zeeman sublevels of the hyperfine manifold F = 1/2. Laser-assisted tunneling methods allow us to design arbitrary operators U rν dressing the hopping between lattice sites r → r + ν, where r = mx + nŷ, m, n, ∈ {1 . . . L}, and ν ∈ {x,ŷ}. Usually, such schemes rely on Raman couplings to auxiliary states trapped in the links of the original lattice and belonging to a different hyperfine manifold. Here, following [31,107], we use bichromatic spin-independent superlattices, and use the secondary minima of F = 3/2 as bus states to mediate the hopping. Note that individual addressing of each hopping rate is granted by the Zeeman splitting within the hyperfine manifolds, and the different detunings of the Raman lasers. These detunings can be quite large, so that the lifetimes of atoms on the lattice (limited by photon absorption and spontaneous emission) can be quite large, of the order of τ l ∼ 1 s. By making the Raman laser intensity/detunings and/or Zeeman splitting spatially dependent, one obtains the desired spatially dependent hopping rates, which is necessary for the realization of curved spacetimes, equation (17). On top of that, it is possible to use Feshbach resonances to turn off the atom-atom scattering and make the system essentially noninteracting.
An alternative, and perhaps even a simpler, method to realize the Fermi-Hubbard model of the form of (17) can be achieved by taking into account the finite waist of the lasers used for the generation of the hopping terms. In general, this is an undesirable feature, and it can usually be neglected. Indeed, typically Gaussian laser beams of waist w are used, characterized by a Rayleigh length -the distance along the direction of propagation from the waist to the place where the area of the cross section is doubled-z R = πw 2 /λ with λ denoting the wavelength. That is to say, within a volume of w 2 × z r , the ideal planar wave is a good approximation, at least around the center of the beam. For a lattice with L = 30, its linear size is 30 × λ/2, so that focusing a laser on the whole lattice w 15λ leads to z R 700λ, so the plane wave description can even be used for an array of a few hundreds of 2D 30 × 30 lattices.
There are, however, no technical obstacles to focusing the lasers on much smaller spots, smaller than Lλ/2, keeping z r still quite large. Actually, such a spatial modulation of the intensity due to the waist of 'real' lasers can be used to induce hopping terms that depend nontrivially on the position [29]. In general, the hopping rate, i.e. the modulus of the hopping term, is proportional to the intensity of the laser producing it. For instance, taking the paradigmatic example of the hopping induced by the Raman transition in Jaksch and Zoller's setup [15], an OL implementing the Hamiltonian of (17) can be achieved for Raman lasers propagating all in the same direction, once we consider only the radial waist and neglect the waist along the beam 9 . In this case, the shape of the laser intensity will correspond to the e factor of the metric.
Comparing the two methods explained above, the former has the advantage that, in principle, any shape of J mn , i.e. any metric of the form (19), can be engineered, but at the price of dealing with quite an involved experimental apparatus, while the latter, although it allows for a restricted choice of J mn (for instance a Gaussian shape), is almost free. The desired hopping rate profile is obtained by reverse engineering of the laser waist.

Density of states at the Fermi level as a simple observable
Let us turn to a discussion of the possible detection schemes of Dirac physics in curved spacetime. A simple observable that characterizes such physics and that contains information about the effects of the beam waist is the density of states [85,108,109]. This quantity is routinely measured in graphene by using scanning tunneling spectroscopy [93,94] and electron transmission spectroscopy [91,92]. For the ideal graphene at zero temperature, the density of states is zero at the Fermi level and is proportional to |E| (once we have fixed E F = 0) as the charge carriers are described by massless fermions propagating in flat space. In the presence of deformations of the graphene sheet, its deviation from free behavior can be analytically computed at first order from the propagator of the Dirac equation. Following [85], by modeling such deformations as perturbations of the Minkowski metric, it is possible to compute the Green function treating the correction to the free equation as an interacting term V . In fact, in this section we will first reproduce the computation of [85] for a metric of the form (19), instead of the spatial deformation (as in (A.1)) considered there.
Before entering into the details of the calculation, a comment on which kinds of corrections to the density of states are expected is in order. In the case of a gauge field deformation, for example by applying a magnetic field to the Dirac system, the main correction to the density is due to the shift of the Dirac points. This implies a nonzero correction to the density of states at the Fermi level. Such a correction is given by the occupation of the corresponding Landau level, as discussed for instance in [49, section 6.8]. For the gravitational deformation we are interested in, the situation is different. Indeed, for the class of metric (19), there is no shift in the position of the Dirac points, as is made evident in the symmetric form of the Hamiltonian (17) where the spin-connection cancels out, but a position-dependent effective Fermi velocity. Hence, at lowest order a correction linear in the energy is expected, with the coefficient depending on the specific metric.
Our final goal is to determine the local density of states, defined by in terms of the Feynman propagator (see the appendix) using the relation where r = (x, y).
In order to findŜ F (w, r, r) we start by using the defining equation for the fermion propagator where the x indicates a point of the spacetime, to a metric of the form (19). By retaining terms linear in , the above equation can be written as where is the effective 'external potential'. As we are interested inŜ F (w, r, r ) and due to the time translation invariance it is convenient to perform the Fourier transformation in time of (23), The above equation can be solved consistently within the first order approximation bŷ whereŜ 0 F (w, r, r ) = 13 is the free fermion propagator and space translation invariance holds. By using the Fourier transformation of (r ) equation (26) can be explicitly computed by performing the integration in r . The relevant contribution to the trace turns out to be linear in w. Explicitly with The above integral is logarithmically divergent, but its imaginary part is not. It is easy to show that this is the only part contributing to the density. Indeed, as (p) is the Fourier transformation of a real function, (p) * = (−p), and (w, p) is even in p, (w, p) = (w, −p), one finds that which immediately implies that δρ(w) = sign(w) 1 π d 2 r d 2 p (2π) 2 e ip · r (p) Im (w, p).
The explicit expression for Im (w, p) is where and p ≡ |p|. For more details see appendix D.
Hence, the density of states always receives a correction proportional to (r) itself:

Example: density of states for a Gaussian beam
We are now able to compute the correction to the density of state when the hopping J has a Gaussian shape due to the finite laser beam waist. Under the assumption that the Raman lasers propagate along the y direction, this implies that (x, y) = −( x a ) 2 , where a is of the order of 10 2 lattice spacings [2]. As (p) = + (2π ) 2 a 2 δ ( p x )δ( p y ), the second term of (34) turns out to be zero and the correction to the density of state simply reduces to 14 providing a clear experimental signature. The local density of states gets a quadratic correction in x. Incidentally, the same cancellation happens in the case of exponential behavior of J , i.e. for a linear (r). Consequently, the simple relation (35) applies also to this case.

Experimental detection
A recent review of the detection methods that can be applied for investigating Dirac physics with ultracold fermions in non-Abelian gauge fields is contained in an article authored by one of us [95]. Here we just summarize this discussion with a particular focus on density of states. Let us start by observing that for a noninteracting Fermi gas at T = 0, the total number of fermions where µ is the chemical potential, is equal at T = 0 to the Fermi energy E F . We see that ρ(E F ) = dN F /dE F so that measuring the variance of N F with E F allows us to determine ρ(E F ). If the system is confined additionally in a slowly varying harmonic potential V (r), a local chemical potential can be introduced: µ(r) = E F − V (r), and the corresponding local density of states, related to the local density by where (.) is the Heaviside (step) function. In this case, we obtain dn(r)/dµ(r) = ρ(µ(r)), i.e. a similar formula as the Streda formula used in [95] for the detection of Hall conductivity. The determination of density of states can be achieved by • Measurements of the total number of fermions as a function of the chemical potential. Here, the best currently available methods are: direct in situ individual atom detection [96]- [98] or quantum spin polarization spectroscopy [99,100].
• Measurements of the (coarse-grained) local density of fermions as a function of the local chemical potential. Again, the best currently available methods are direct in situ individual atom detection [96]- [98] or quantum spin polarization spectroscopy with spatial resolution [101].
• Measurements of the frequency-momentum resolved single-particle excitation spectrum, such as those being carried out in Bragg (Raman) scattering spectroscopy (for a state-ofthe-art report, see [102]). The spectrum in such processes is proportional to the density of initial states of the scattering process.
Of course, many other methods, such as atomic ARPES, noise interferometry or even absorbtion and/or phase contrast imaging, can give at least indirect information about ρ(E). All of these methods are well developed in experiment with ultracold atoms (see [95] and references therein).

Conclusions and outlook
In this paper, we discussed the simulation of the Dirac equation in artificial curved spacetime with cold atoms. We showed that, using state-of-the-art techniques, it is possible to simulate relativistic fermion dynamics in curved spacetimes with a flat 2D square lattice for an interesting class of (2 + 1) metrics. Moreover, we pointed out the relation between a certain class of Hubbard models and the Dirac Hamiltonian in curved backgrounds, which can be employed to make analytic computations in the continuum limit of the former. We proposed to characterize the nature of Dirac fermions on the lattice by measuring the density of states at the Fermi level. This observable, on the one hand, can be analytically computed in perturbation theory in terms of the Dirac propagator and, on the other hand, is accessible to measurements. The present study opens the way to the direct observation of elusive effects such as Rindler noise. Because we deal with the odd dimensional (2 + 1) Rindler system, the Dirac thermal noise, measured by an ideal point-like De-Witt detector as a consequence of the local acceleration, is expected to be 'anomalous' (see chapter 8 of [82]), i.e. it should follow the Bose-Einstein distribution. This issue is currently under investigation.

Acknowledgments
AC thanks J M Pons and J Garriga for useful discussions and J G Russo for providing insights into the thermalization theorem and for pointing out [82]. Now we consider a different situation where the (2 + 1) metric is spatially nontrivial. Such a case is relevant in describing the properties of a graphene sheet with ripples. The most generic spatial deformation (at least in some patch) can always be written as ds 2 = −dt 2 + e 2 (x,y) (dx 2 + dy 2 ).
The driebein are e 0 = dt, e 1 = e (x,y) dx, e 2 = e (x,y) dy, i = x, y, and the spin-connection can be chosen to be nontrivial in the spatial part only: It follows that the curvature is For instance, the slices of the metric (A.1) at constant time will be spheres or hyperboloids for a positive or a negative quadratic form of x and y, respectively. Applying (11) to this case, we find We are tempted to interpret the above Hamiltonian as that of fermions coupled to a 'geometric' non-Abelian vector potential A ≡ (∂ y , −∂ x )σ z . Adopting the the gamma matrices' representation of the previous section, we can write it as where the presence of A indicates that the rotation in the x y-plane, i.e. the SO(2) subgroup of the Lorentz group, is promoted to a local symmetry in the background described by the metric (A.1). Such identification is related to the treatment of the conical defects in the graphene sheets (dislocations and disclinations) as sources of magnetic fluxes. However, this interpretation here is misleading as the gauge group does not commute with the spacetime symmetry, as σ z anticommutes with σ x and σ y . Taking in account this fact, the Hamiltonian density can be more appropriately written as where the symmetric role of x and y is evident. Now we are ready to compute the total Hamiltonian. By rewritting (12) we obtain Using the the same manipulations as in the previous section, we can recast it into a form where the spin-connection is not present, It follows that the discretized version of H , as in Rindler spacetime, takes the form of a SU(2) Fermi-Hubbard model with the modulus of hopping depending on the position At first sight it seems very surprising that the discretized Hamiltonian of massless fermions in Rindler geometry coincides with the one of fermions propagating in metric of the form of (A.1). Indeed, by taking = ln(ax) the expression (A.9) reduces to (16). Such an apparent contradiction disappears upon closer inspection. At the end, from the point of view of the Dirac Hamiltonian for both metrics, what has changed with respect to the flat case is the effective speed of light, or equivalently the hopping rate, which becomes position (and direction) dependent. Although the origin of the position-dependent hopping rate is different in the two cases (it comes from the Hamiltonian density in the Rindler case whereas it is due to the invariant measure in the other), the effect is the same. Roughly speaking, there are two possible ways of modifying the effective speed of light in one spacetime direction, let us say x: one way is to change g tt , while the other way is to change g x x by an inverse factor.
Nevertheless, the eigenfunctions and the spectra of the two Schrödinger problems remain different as the Hamiltonian densities in the two cases are.

Appendix B. The generic case
In order to treat the generic case let us retrace a few steps and analyze the formal expression for the Hamiltonian (12). First of all, we show that it corresponds to the Legendre transformation of the relativistic Lagrangian: once we have chosen the coordinates to have that the time-like Killing vector is K = ∂ t , which implies that ∂ t g µν = 0. Indeed, by defining H = δL δ∂ t ψ ∂ t ψ − L, the expression (12) is recovered: Now, it is instructive to check explicitly that the Hamiltonian above is a Hermitian operator due to the existence of a time-like isometry. By using that (γ µ ) † = γ 0 γ µ γ 0 and noting that ψ γ µ w ab µ γ ab ψ † = −ψ w ab µ γ ab γ µ ψ, we obtain that In order to compare the above expression with H it is convenient to integrate by parts and rewrite it as The cancellation follows from the identity as ∂ µ ln √ −g = ν µν and ∂ t γ t = γ t ∂ t ln √ −g = 0 if and only if ∂ t g µν = 0. At this point, we can use expression (B.3) for H † to obtain a Hamiltonian symmetrical in ψ and ψ † . By writing H ≡ 1 2 (H + H † ), we find Let us characterize the term ≡ − 1 4ψ w ab µ {γ ab , γ µ }, which can be regarded as the obstruction to write the lattice Hamiltonian simply as as found for the special cases discussed in the previous sections. Using the relation between the spin-connection and the driebein and the anticommutator {γ ab , γ µ } = 2e µc γ abc in (2 + 1) dimensions it reduces to {γ ab , γ µ } = − 1 3 e µc abc and it follows that We will express the number of states as the number of poles in the propagator when averaged over the set of eigenstates of H , {|n }, which we demand be complete and normalizable (as is the case for any sound Hamiltonian operator). Furthermore, we assume that the position operator eigenstates |r can be completed to give an orthonormal basis |r |i where i encodes all internal degrees of freedom, such as the spin. Using the theorem of residues and integrating on a rectangular contour around the real axis of width 2 > 0, we can write where d 2 w p = d 2 p (2π) 2 2|p| is the Lorentz invariant measure and the normalization of the momentum states |p, ± is fixed accordingly. By taking S = d 2 r as the volume of the system and by going to polar coordinates, one finds that ρ(E) = S π Im +∞ 0 d p 2π 2 pE p 2 − E 2 − i sign(E) . (C.5) The last integral can be solved in many ways, for example by changing the variable to z = p 2 and regularizing the integral with cutoff 2 0 dz 2π As the imaginary part is independent of the cutoff, the final result, as expected, is which is in agreement with the result quoted in [110]. It is worthwhile to note that ρ(E) can be written in terms of the Feynman propagator. This fact can be derived in a more general setting. By definition, 1 E−H +i is the Fourier transformation in time of the retarded propagator defined by the equation with boundary conditions G + (t − t ) = 0 for t < t . On the right-hand side of (C.8), the identity in spin space has been written explicitly to remind the reader that G + (t − t ), in general, acts as a matrix on the internal degrees of freedom. After multiplying on the left by γ t and on the right by γ t and taking the expectation value with eigenstates of the position operator, the above equation gives −iγ µ ∂ µ G + (t − t , r, r )γ t = δ(t − t )δ d (r − r ) ⊗ I spin , (C.9) where we use γ t γ t = 1. Hence, we conclude that G + (t − t , r, r )γ t is related to the Feynman propagator as it solves the same equation. The precise relation can be derived by taking the boundary conditions into account. This can be explicitly checked by Fourier-transforming to momentum space. Indeed, . (C.10) The above relation between the retarded propagator and the Feynman propagator can be extended using perturbation theory to the interactive second quantized formalism. To conclude, let us remark that it is better to use the local definition of density (21) because it is easy to make it generally covariant in order to apply it in a curved gravitational background. In this way, the generalized notion of inner product is properly taken into account due to (23).