Unruh effect for interacting particles with ultracold atoms

The Unruh effect is a quantum relativistic effect where the accelerated observer perceives the vacuum as a thermal state. Here we propose the experimental realization of the Unruh effect for interacting ultracold fermions in optical lattices by a sudden quench resulting in vacuum acceleration with varying interactions strengths in the real temperature background. We observe the inversion of statistics for the low lying excitations in the Wightman function as a result of competition between the spacetime and BCS Bogoliubov transformations. This paper opens up new perspectives for simulators of quantum gravity.


Introduction
Since the seminal idea of Feynman [1] the capabilities and regimes of applications of quantum simulators [2,3] have rapidly grown in the last decades. In particular, systems naturally defined on a lattice like trapped ions [4,5] and neutral groundstate and Rydberg atoms in optical lattices and tweezers [6][7][8], are excellent platforms for simulating coursed-grained models. They are obviously well suited for simulating condensed matter phenomena where the lattice is intrinsically there. For instance, by engineering synthetic gauge fields by laser means [9] in real or synthetic lattices [10,11] one can experimentally realize Hofstadter model [12,13] and visualize the edge states of quantum Hall effect [14][15][16], as well as realize other topological models like the Haldane model [17]. Similarly by lattice shaking [18], one can achieve non-Abelian gauge fields [19,20] and observe classical magnetism [21,22] (for comprehensive reviews see e.g. [23]). Such single-particle simulators are a first step towards experiments that could help solving open questions in many-body phenomena through quantum simulation, as for instance the relation between the Hubbard model and high-T c superconductivity [24][25][26]. The combination of synthetic gauge fields with interactions could lead for instance to the realization of (quasi-)fractional quantum Hall states in ladders [27,28] or quantum spin liquid [29] states in triangular lattices with ultracold atoms.
Lattice based atomic simulators are also well suited for simulating high-energy physics where the lattice is adopted as regularization tool for studying, e.g., strongly coupled gauge theories [30]. Lattice gauge theories can be interpreted as the natural next step after synthetic gauge fields, where the phases representing the classical gauge fields are promoted to operators acting on the gauge degrees of freedom living on the links of the lattice. By adopting convenient gauge invariant truncations of such degrees of freedom, one can design simulators of both Abelian [31][32][33][34][35][36][37] and non-Abelian gauge theories [38][39][40] capable to probe confinement, string breaking, and to study the dynamics of charges as in the first proof-of-principle experimental realization of the Schwinger model with four ions [41,42]. Together with the new classical simulation approach based on tensor network, e.g. [43][44][45][46][47][48][49][50], such simulators offer a new way to study challenging open questions in high-energy physics like the dynamics and the phase diagram of strong interactions at finite temperature and densities, for review see [51][52][53]).
Recently, lattice based quantum simulators have been pushed forward also for the study of gravitational phenomena [54]. There is a natural and unresolved tension between quantum physics and general relativity that produces spectacular and counterintuitive effects. Indeed, the quantization requires the choice of a preferred time direction, which fixes the notions of vacuum state and particles, choice that it is at odd with general covariance and makes such notions depending on the observer [55,56]. Such tension becomes manifest in presence of an event horizon and leads to the appearance of phenomena like the Hawking radiation of a black hole [57] and the Unruh effect [58] for an accelerated observer. The latter sees the vacuum as a thermal state with a temperature proportional to acceleration essentially because in his/her reference frame (Rindler metric) there is an event horizon and he/she can access only part of the vacuum [59][60][61]. Simulating the Unruh effect and the Hawking radiation is interesting on one hand because these striking effects are hard to observe in nature, on the other hand because they are at the moment the only experimentally accessible phenomena that can bring us closer to quantum gravity.
Alternatively, one can simulate the motion of artificial Dirac fermions in artificial gravity background by tailoring Dirac Hamiltonian on the lattice properly. This strategy allows for a systematic quantum simulation approach. As shown in [54] (see also [81][82][83]) one can engineer the motion of Dirac fermions in optical metrics [84] by making the Fermi velocity position dependent, that is by engineering position-dependent tunneling rates on the lattice, which can be equivalently translated in random walk processes [85][86][87][88]. One can apply such procedure to any lattice formulation of Dirac Hamiltonian [89], as for instance artificial graphene obtained from ultracold fermions in a brick-wall lattice [90] (see [91,92] and [93] for realizations of hexagonal optical lattices). Note that one can obtain the motion of Dirac fermions in certain curved space-times by physically bending a graphene sheet [94,95]. In principle, one can observe signatures of the Unruh effect with some subtleties in such a system [96][97][98] as well as in more complicate scenarios [99].
Engineering artificial gravity through position-dependent tunneling rates for ultracold atoms in optical lattices offers several advantages, as we have shown recently in Ref. [ 3. It allows to go beyond single-particle phenomena and explore many-body physics in curved space-times.
Here we take a first step in the this largely unexplored world and consider the manifestation of the Unruh effect for interacting fermions in two spatial dimensions. The model we consider is essentially equivalent to the 2D version of Thirring model [101] (the thermal nature of Hawking radiation was verified for 1D version in [102]). We propose an experimental protocol for studying the interplay between relativistic invariant quartic interactions and the acceleration with ultracold atoms. In the experiment, the acceleration of the (interacting) vacuum can be obtained by quenching the Dirac Hamiltonian as proposed first in [100], while the interactions can be controlled for instance by Feshbach resonance [103]. Alternatively, interaction can be mediated by another bosonic atom, as proposed for the Thirring model in flat space in [104]. We show that when the system can be described thought BCS theory the Unruh effect occurs for the Cooper pairs: due to the "inversion of statics" [60] we observe a crossover between bosonic and fermionic thermal response in two spatial dimensions.
The paper is organized as follows. In Sec. 2 we review the derivation of the singleparticle and interacting naive Dirac Hamiltonian in the Rindler metric on a square lattice. The expert reader can directly to Sec. 3 that is the hearth of the paper where we discuss the theoretical and experimental aspects of the Unruh effect for interacting Dirac fermions. We consider the case of attractive interactions and derive the Wightman response functions in the mean field limit. We obtain it by determining the Bogoliubov transformation that relates the quasi-particle states in Rindler and Minkowski spacetime, once the normalization of the Rindler Hamiltonian relative to the Minkowski one is chosen such that the former equals the latter in the far horizon limit, i.e. on the boundaries of the cylindrical lattice. We find the expected thermal response function for the quasi-particles, in particular at low energy the effective description in terms of Cooper pairs translates in a Fermi-Dirac-like Planckian response characteristic of bosons. In Sec. 4 we take a step back and analyze in details the BCS theory for the naive Dirac Hamiltonian in Minkowski space. In particular, we discuss and explain its subtle dependence on the boundary conditions on finite-size systems (further explanation can be found in the Appedindix that contains also derivations relevant for Sec. 3). Finally, we draw our conclusions in Sec. 5, we resume our results and comment about the new perspectives and developments opened up by our work.

Derivation of discrete Dirac Hamiltonian
In this section, we present the derivation of the discrete Hamiltonian of interacting massless Dirac fermions in the Rindler Universe in (2+1) dimensions. For the sake of clarity and to be self contained, we first review the derivation of the non-interacting Hamiltonian [54,105], by starting from the general Lagrangian density for Dirac particles in (d+1) curved spacetimes. We use this section also to fix the notion used in the following sections.

Lagrangian density
The Lagrangian density of a massless non-interacting Dirac spinor ψ in a (d + 1) dimensional curved spacetime reads (see, for example [106][107][108]) Here D µ ψ is the covariant derivative of a spinor ψ, introduced to guarantee the invariance under the local Poincaré transformations, w ab µ (x) is the spin-connection which gauges the Lorentz group, and γ's are the gamma matrices in a curved spacetime where g µν is the metric tensor ds 2 = g µν dx µ dx ν . Note that we are using the mostly-positive metric signature, where the norm of timelike vectors is negative. Choosing the opposite signature would require to multiply all gamma matrices by the imaginary unit and, as a consequence, changing sign of the Lagrangian density. In (2+1) dimensions γ's are 2×2 matrices, and can be expressed in terms of Pauli matrices. A possible choice for the "flat" gamma matrices in Minkowski space is Irrespectively of the choice of the gamma matrices, in (2 + 1) dimensions the product of all γ's is proportional to the identity γ 0 γ 1 γ 2 = −1.

Spin-connection and vielbein
The relation between gamma matrices in curved spacetimes γ µ ≡ γ µ (x) and the flat matrices γ a is given be a vielbein field e a µ , i.e. γ µ (x) = e a µ (x)γ a . The vielbein field introduces locally flat Cartesian frame of reference. It is defined by where e µ a and e a µ are inverse matrices, and e µ a e b µ = δ b a . Note that the vielbein e a µ and the spinconnection w ab µ (x) fields cannot be independent, as Lorentz translations must relate orthogonal frames in different points. The requirement of a covariantly constant vielbein [109] imposes where the indices µ and ν are antisymmetrized. The spin-connection can be written explicitly using the Christoffel symbols Γ µ

Lagrangian density in the Rindler Universe
We consider (2 + 1) dimensional static metric of the Rindler Universe [60,61,109] which describes a flat Minkowski spacetime from the point of view of an observed moving with the constant acceleration a = 1/|x|. The Rindler universe is static. Indeed, with the above coordinate choice the Rindler metric (10) is manifestly time-independent, with the time-like Killing vector B = ∂ t . The manifest time-translational invariance allows to construct a conserved Hamiltonian function. We can obtain it in a simple manner from the Dirac Lagrangian in Rindler space we derive below.
From (10), we choose the nonzero components of the vielbein as which produce a particularly transparent form of curved gamma matrices Accordingly with the choice of a vielbien, we can construct the spin connection, which the only nonzero components are where the last equality comes from the antisymmetricity of w ab µ in its internal (flat) indices. The Dirac Lagrangian density L 0 written explicitly in the (2+1) Rindler Universe becomes

Dirac Hamiltonian
The Hamiltonian density is obtained by the Legendre transformation For the Dirac Lagrangian density L R 0 in the Rindler Universe we have Finally, we obtain the Rindler Hamiltonian by integrating the Hamiltonian density on a spacelike hypersurface H R 0 = dΣ H R 0 , where the volume element dΣ = √ −g dxdy includes the metric determinant. It is very convenient to write the Hamiltonian in the fully symmetric form, i.e., resulting in where we express the γ's in terms of the Pauli matrices (4).
Once written in the form of (18), we can properly discretize the Rindler Hamiltonian simply by replacing integration dxdy with summation d 2 m,n , where d is the lattice spacing, and by replacing the derivatives with finite differences which fulfill the Leibnitz rule for differentiation. Finally, the Rindler Hamiltonian on the lattice reads which has the form of a Hubbard Hamiltonian with non-Abelian non-uniform tunneling terms Note that the Rindler Hamiltonian is scale invariant, as it does not depend on the lattice spacing d (d 2 from the volume element cancels due to the replacement ψ m,n → d −1 ψ m,n , done to guarantee the normalization condition dxdy ψ † ψ = 1).
A very similar derivation can be performed to obtain the Dirac Hamiltonian in a flat Minkowski space. In fact, since the Rindler Hamiltonian in the symmetric form (18) does not contain spin connection terms, it suffices to replace the metric determinant with a constant, i.e. √ −g = c, which consequently leads to a discrete Hubbard Hamiltonian with a constant tunneling t (hereafter we choose t as an energy scale, i.e. we put t = 1).
The Minkowski Hamiltonian (23) should be recovered from the Rindler Hamiltonian (21) in the asymptotic limit of small acceleration a = lim |m|→∞ 1/|m| = 0. Therefore, we should choose a common energy scale at the boundary of the system. This is done by rescaling the tunneling values (22) where m = −M, . . . , M .

Lagrangian density for interacting particles
In order to describe interacting Dirac particles we need to include a nonlinear term in the Lagrangian density (14) The simplest choices for L int , which would guarantee Lorentz invariance and U(1)-current conservation are • density-density interaction [110] • current-current interaction [101,111] Note that in (3+1) spacetime dimensions, other possible choices would be L int = λ 2 ψ γ 5 ψ 2 or L int = λ 2 ψ γ 5 γ µ ψ ψ γ 5 γ µ ψ [112], but in (2+1) dimensions a γ 5 matrix is trivial -the product of all gamma matrices is proportional to the identity (5), reflecting the well-known fact that chirality is absent in even spatial dimensions.

Interacting Dirac Hamiltonian
The Legendre transformation of the full Lagrangian density (25) leads to the Hamiltonian density Similarly to the non-interacting case, one gets the full Hamiltonian by integration over a spacelike hypersurface The discretized version of (29) is straightforward, as neither choice of the Lagrangian interaction density (26) - (27) includes derivatives. For the Rindler Universe and densitydensity interactions (27), we obtain where we decompose the lattice spinor operator ψ m,n = [c ↑,m,n c ↓,m,n ], introduce the number operators n σ,m,n = c † σ,m,n c σ,m,n , n m,n = σ n σ,m,n , and define U m = |m|λ/d ≡ |m|U . Also, we can easily check that the other choice of Lorentz-invariant interactions (27) leads to an equivalent term as apart from a different coefficient multiplying n ↑,m,n n ↓,m,n , both interaction terms in (2+1) dimensions have the same form. Finally, let us compare (30) with the corresponding interaction Hamiltonian in the Minkowski space In order to make the two terms to coincide in the limit of zero acceleration we have to renormalize U m accordingly to such to make the value of the couplings to coincide at the lattice boundaries U M = U as previously done in the non-interacting case for the tunneling, t x,y M ∼ t = 1.

Rindler Universe Dirac Fermions with attractive interactions
In this section we consider a model of Dirac fermions with attractive interactions in the meanfield regime. The primary focus of this section is to analyze the power spectrum of the Rindler noise with increasing interaction strength, and show the difference between interacting thermal particles in a flat space and interacting accelerating particles at non-zero Unruh temperature T U . In order to compute it we first find quasi-particle basis in the Minkowski and in the Rindler space for appropriately normalized Hamiltonians defined in the previous section. Indeed, Rindler Hamiltonian has to coincide with the Minkowski one in the limit of zero acceleration.
On the lattice, this requirement translates in matching the coupling on the boundaries in x (we exploit the translational invariance along y by taking it periodic such that our lattice is a cylinder). Then, we find the unitary transformation between the two quasi-particle basis and determine the Wightman function on the Minkowski ground state in the Rindler universe, i.e., as measured after quenching the Hamiltonian from Minkowski to Rindler spacetime. Finally, we study the case in which the Minkowski background in which the observer accelerate is the thermal state of the interacting Dirac Hamiltonian. We conclude the section by discussing the experimental requirements for observing the Unruh effect in presence of interactions with ultracold fermions in optical lattices.

Self-consistent Minkowski Hamiltonian
Let us start with writing effective mean-field Hamiltonian in the Minkowski spacetime at T = 0 and express eigensolutions in terms of quasiparticle modes. The full Hamiltonian of interacting Dirac particles in the Minkowski spacetime is a sum of noninteracting (23) and interacting (32) terms or explicitly in (2+1) dimensions where U is an attractive interaction strength U = −|U |, and ψ m,n = [c ↑,m,n c ↓,m,n ], n σ,m,n = c † σ,m,n c σ,m,n , n m,n = σ n σ,m,n . We tackle the problem by approximating the interaction term with a mean field averages [113,114] where we consider all possible terms, i.e.
Notice that (39) does not conserve spin (or species) of particles and therefore, we expect it to be identically zero Λ ≡ 0. Consequently, the mean field Minkowski Hamiltonian for interacting particles is where H M 0 is a free particle term given by (23), H M n is an on-site potential energy term (where we include the chemical potential µ) and H M ∆ is a pairing term Before we proceed to the eigensolutions of (40), let us stress that ∆ and W are given in (37) - (38), where the averages are calculated in the ground state |0 M of H M mf . Consequently, ∆ and W are not external parameters and the eigenvalue problem should be solved self-consistently. Furthermore, since both Minkowski and Rindler metrics are translation invariant in y direction, we consider the eigenvalue problem on the cylinder (see Appendix A for the analytical solution on a torus in Minkowski space). After performing the Fourier transformation ψ m,n = 1/ N y ky e ikyn ψ m,ky , we write the Hamiltonian in a compact matrix form where a spinor Ψ ky is defined as and where P x is a discrete derivate in x direction P x ψ m,ky = −i t ψ m+1,ky − ψ m−1,ky and ky = 2 t sin k y . In the Hamiltonian matrix (45) we drop the on-site potential term (41), as we choose the half-filling condition which yields W = U/2 at µ = 0. Eventually, we can express the spinor Ψ ky in terms of the quasiparticle eigenmodes are column eigenvectors to E M ky,p and −E M ky,p , respectively (see Appendix A). The expression (46) is a canonical Bogoliubov transformation between interacting particles and quasiparticles which diagonalize the Hamiltonian (45).

The mean field Rindler Hamiltonian
In the Sec. 2 we derive a kinetic (21) and interacting (29) terms of the Dirac Hamiltonian in the Rindler Universe where primes denote the rescaled tunnelings (24) and interactions (33). Let us stress again that an accelerated observes moves in the physical vacuum |Ω , which is a ground state of the Minkowski Hamiltonian |Ω = |0 M . Although the Rindler Hamiltonian governs the dynamics of an accelerated observer, its ground state |0 R does not obviously need to coincide with the physical vacuum |0 R = |Ω . Therefore, from the point of an accelerated observed, its background is an excited state.
For that reason, the mean field averages for the interaction terms should now be calculated not in the ground state |0 R , but in the physical vacuum |Ω = |0 M . For example, the Rindler pairing function ∆ R reads where ξ m = |m|/M is a rescaled distance from the horizon. Apart from the position-dependent elements, the Rindler Hamiltonian (48) has the same form as the Minkowski Hamiltonian (35), and therefore we obtain its mean field counterpart practically automatically with where ξ is a rescaled position operator ξψ m,ky = ξ m ψ m,ky , ∆(ξ) = ξ∆ is a position-dependent pairing function, ky (ξ) = ξ ky a position-dependent tunneling energy (transverse to the horizon), and where we remind the reader that the prime in summation indicates that p in p runs over the positive eigenvalue, E R ky,p > 0. Now, let us focus on how the Unruh temperature T U influences the Rindler pairing ∆ R (in Sec. 4 we discuss how the pairing ∆ in the Minkowski space changes with the physical temperature T ). Since ∆ R (m) = ξ m ∆ and ξ m is proportional to the inverse of acceleration 1/a = |m|, we could be tempt to write that the ∆ R is proportional to the inverse of Unruh temperature T U = a/(2π). Nevertheless such conclusion would not be correct. The Unruh temperature is defined locally on x = constant-hypersurfaces of the Rindler Universe and therefore is different for inequivalent observers. The expression (49) tells us how the pairing function changes with |m| from the point of view of the Minkowski observer. Since the proper time slows down the closer we are to the horizon, then the interactions seem to be weaker. Therefore ∆ R (m) for the observer at |m| should be rescaled with the proper time ∆ R (m)/ξ m = ∆. Consequently, the pairing strength seen by an accelerated observer does not depend on the Unruh temperature. At the same time, an observer at |m| sees the pairing to be weaker ∆ R (m)/ξ m < ∆ when she/he looks towards the horizon and stronger ∆ R (m)/ξ m > ∆ when she/he looks opposite to the horizon.

Wightman function for an accelerated observer
The expressions (46) and (52) are canonical Bogoliubov transformations [113,114] between interacting fermionic particles and noninteracting quasiparticles in the Minkowski and Rindler spacetimes, respectively. As quasiparticles modes (46) -(52) diagonalize the mean-field Hamiltonian (43) -(50), we can write where p in the sum runs over positive eigenvalues and A refers to either Minkowski (M ) or Rindler (R). Note that in the last equality of (53) we use the fact that H A ±ky have the same spectra (see the discussion of the symmetries of H A ky in Appendix C). From (53) we see that the groundstate of H A mf is deprived of quasiparticle excitations. Therefore, it is anihilated by all quasiparticle operators A state that fulfills (54) is of a form where |0 where is a particle vacuum. Since quasiparticle operators mix particle creation and annihilation processes, we directly see that the ground state |0 A must contain particles, and that the groundstates of Minkowski and Rindler Hamiltonians are different as they have different quasiparticle excitations. Also, because quasiparticle modes (46) and (52) are different for the two (stationary and accelerated) observers, the act of creation of a particle is seen differently. Combining (46) and (52) together and using R † ky ,p R k y ,p = δ ky,k y δ p,p we obtain which is the spacetime Bogoliubov transformation [60,61] between two noninteracting quasiparticles from the point of view of different observers. It is straightforward to find out that in general β R ky ,p |0 M = 0, and so indeed |0 M = |0 R . In particular, we expect a different response from a particle detector for different observers. Let us write the Wightman function for an accelerated observer G m,n (t) = 0 M |c † σ,m,n (t)c σ,m,n |0 M , and its Fourier time transform which is the power spectrum of the Rindler noise.
Because the time evolution c † σ,m,n (t) is different for the stationary and accelerated observers, and since |0 M = |0 M , it is intuitive that a response function G m,n (ω) should also be different for the two reference frames. However, the most interesting part is far less intuitive: (i) the response function G m,n (ω) for an accelerated observed exhibits thermal behavior, (ii) in even spacetime dimensions thermal distribution of fermions (bosons) is Fermi-Dirac (Bose-Einstein), but in odd spacetime dimensions the statistic interchange.
In particular, in the continuous limit the power spectrum of a thermal noninteracting gas in (2+1) Minkowski spacetime is [60] G T M (ω) = whereas an accelerated observer sees where the Unruh temperature is T U = a/(2π).
One might wonder how to relate the fermionic and bosonic power spectra of a cold thermal gas (T → 0) in a flat space, to the ones seen by an accelerated observer in the limit of zero acceleration (T U → 0). In fact, it is easy to find out that in the zero temperature limit, the modulus of the Bose distribution is Fermi-Dirac therefore, except for a singular point at ω = 0, the power spectra match exactly in the zero temperature limit. Applying the spacetime Bogoliubov transformation (56) we get explicitly the Wightman function Since we expand the Dirac fields in the eigenmodes of the Rindler Hamiltonian, we expect that the Fourier transform of (62) should express the power spectrum as seen by the accelerated observer a = 1/|m|. However, this statement is true only if t is the proper time of the observer. In the Rindler Universe the proper time of an observer is ξ m t, and consequently, to compare the response of different observers, the power spectrum of noninteracting particles needs to be rescaled as ξ m G m (ω/ξ m ) [100]. Furthermore, in the interacting problem we have a band gap separating valance and conductance bands. The band gap 2∆ R (m) seen by a Minkowski observer is |m| dependent. To account for that, we need to compare shifted power spectra, i.e. ξ m G m (ω − ∆ R (m))/ξ m . Note that in the continuous limit, the Unruh temperature of noninteracting fermions is T U = 1/(2π|m|). The numerical results for the finite lattice system show that indeed T U decreases with increasing |m|, but the functional dependence T U (m) might be in principle different. We investigate the relation T U (m) by fitting the continuum limit power spectrum (59) to the numerical results. Eventually, we find that T U (m) = A/|m| + B with A ≈ 0.2 and B ≈ 0.01 reproduces quite well the data for 5 < |m| < 45, see Fig. 1. Thus, we find good agreement with the results of Bisognano-Wichmann theorem [115,116] (see also [117]) that holds in the infinite lattice limit. Indeed, the theorem, which follows from axiomatic field theory, implies that all Lorentz-invariant local field theories display the Unruh effect with T U = 1/(2π|m|). In the following we discuss in details the behavior of the Rindler noise in the different regimes of interest, extending the results obtained by us in [100] in the non-interacting Dirac fermions. As in [100] the numerical results are in the agreement with the continuum limit in the linear dispersion regime, i.e. ω/ξ < 1.

Power spectrum at T U ≈ 0 and T = 0
We start by discussing the limit of zero acceleration, T U → 0, in which we expect the Wightman function in the Rindler space to coincide with the Wightman function in the Minkowski space at T = 0. By using the analytical solution for the Minkowski system on a torus (see Appendix A), we find for the latter and consequently where is a positive quasiparticle eigenenergy of a system and ρ(ω) is the negative energy density of states.
The numerical density of states ρ(ω) for different values of ∆ is plotted on Fig. 2 (left  panel). For the noninteracting system the density of states can be well approximated with the continuous limit density, i.e. ρ(ω) ∝ |ω| for |ω| 1. For the interacting system we can approximate for |ω| 1, i.e. a nonzero value of the pairing function ∆ introduces a 2∆ band gap and a ∆ step jump in the density of states. As ∆ is much smaller than the half band width w = 2 √ 2 ≈ 2.83 of the noninteracting system, for |ω| ∆ the density of state is basically free fermionic, and only at |ω| ≈ ∆ encounters the step function jump. As we expect that G T =0 M (ω) replicates the power spectrum of the Rindler noise for T U = 0, we can interpret this Fermi-Dirac like behavior as a response of composite bosons (i.e. Cooper pairs), which should be fermionic (statistic inversion in even spatial dimensions). In other words the pairing influences power spectrum only near |ω| ≈ ∆, which is expected as Cooper pairs tend to form near the Dirac cones (corresponding to E k ≈ ∆) and their number increases with ∆ (see Appendix B).
Since we expect that (68) estimates the power spectrum of the Rindler noise for T U = 0 can infer that for T U 0 the qualitative behavior of the power spectrum should be for |ω| 1. It is trivial to realize that (69) in the limit T U → 0 recovers (68). Also, for ∆ ≈ 0 and ω/T U 1 we can drop the plus one in denominator (69) and therefore recover the power spectrum of noninteracting fermions (60).

Power spectrum at T U > 0 and T = 0
In this section we calculate the power spectra of the Rindler noise using the explicit formula In order to both minimize lattice artifacts and to mimic realistic measurement schemes of finite space resolution, as in [100] we present theG(ω) convoluted with a Gaussian filter of 2-site width. For |ω| 0 the fermionic power spectrum is recovered. Note that the curves were rescaled to account for different proper times of different observers and shifted by ∆ (due to the spectral gap).
The numerical results reproduce quite well the expected thermal behavior (69) with the Unruh temperature inversely proportional to the distance to the horizon T U ∝ 1/|m|. For |ω| 0 the power spectra are approximately linear (just like in a free fermionic case), while near |ω| ≈ 0 we observe a clear Fermi-Dirac profile (i.e. a bosonic response), which is more evident for ∆ ≈ 0.5, since the number of Cooper pairs is more significant. Note that, as expected, the Fermi Dirac plateau is nearly twice as higher in for ∆ ≈ 0.5 then in ∆ ≈ 0.25.
After rescaling (71) we can directly compare different power-spectrum curves. The results are plotted on Fig. 4. It turns out, that the postulated scaling is indeed valid near |ω| ≈ 0, see the discussion in the previous section.
Note that the chosen values of the pairing function ∆ on Fig. 3 are one order of magnitude greater than the order of the Unruh temperature. For ∆ ∼ T U we find that the power spectra are changed only marginally in comparison to the noninteracting ones.

Power spectrum at T U > 0 and T > 0
In this section we again consider the power spectrum of Dirac fermions in the Rindler Universe, although, now we consider as background the thermal state at T > 0 of the interacting Dirac Hamiltonian in Minkowski space. Such situation is relevant also from the experimental point of view as the ultracold fermions in optical lattices have typically a non-negligible temperature compared to the band width, thus order one in our units (t = 1). Thus, it is crucial to establish the interplay between the Unruh temperature T U and the physical temperature T such to determine the visibility of the Unruh effect for interacting Dirac fermions, as done in [100] for the noninteracting case. The Wightman function can be written as where and the power spectrum is In order to compare with the results of the previous section, we choose the interaction U = 3.55 and manipulate the physical temperature T in such a way to obtain ∆(T ) ≈ 0.5, 0.25 (as in Fig. 3).
The numerical results are presented on Fig. 5. As expected, we find that the thermal background does not affect the Unruh effect until ∆(T )/T ∼ 1. In particular, for T = 0.315 and ∆(T = 0.315) ≈ 0.25 (left panel) we observe a thermal positive frequency contribution from the quasiparticles above the Fermi sea. For T = 0.1 and ∆(T = 0.1) ≈ 0.5 (right panel) the power spectrum is very close to the T = 0 case.

Experimental realization and detection of Unruh effect
In [100] we presented a detailed proposal how to simulate the noninteracting Dirac Hamiltonian in Minkowski (35) and Rindler (48) spacetimes in the optical lattice setup. Since the tunneling matrices of a two component spinor ψ m,n is purely off-diagonal, then flipping its components at every second site ψ m,n → σ x ψ m,n diagonalize the tunneling and the two components do not mix. Therefore, the noninteracting naive Dirac Hamiltonian on a square lattice is equivalent to two independent copies of a π-flux Hamiltonian, and the optical lattice simulation can be done with one component only. In other words, one can exploit that the lattice is bipartite and that the unit cell has dimension 2 due to the π flux to encode the two components of Dirac spinor in two different sublattices (and reduce the doubling of Dirac points to 2). The key feature of our experimental proposal in [100] is that the tunneling is assisted in both x and y directions by Raman lasers that induce a synthetic magnetic field of flux π in the symmetric gauge. The intensity of the tunneling is controlled by the intensity of (one of) Raman lasers. Such a scheme allows both for the preparation of the Minkowski vacuum as groundstate of the corresponding non-interacting Dirac Hamiltonian with uniform tunneling rates, and for its "acceleration" via a sudden quench of the tunneling rates to a V -shape profile of the Dirac Hamiltonian in Rindler spacetime. The generalization of the scheme above to the case of interacting Dirac fermions requires the introduction of spatially tunable interactions. One possibility would be to realize the scheme in [100] with dipolar fermionic gases like erbium, where it has been very recently experimentally demonstrated that dipolar interactions are stable and tunable by Feshbach resonances [118]. In such a scheme, each of the two components of the spinor c ↑ , c ↓ are identified with the two spin species of erbium that occupy two different sublattices. The interaction term n ↑ n ↓ , thus becomes a nearest-neighbor density-density interaction provided by the magnetic dipole moment of the atoms. The magnetic field (or in alternative the light shift) inducing the Feshbach resonance has to be then tuned spatially such to provide the desired V -shape interaction profile. In alternative, at the price of doubling the number of Dirac points, we can consider two spin states of fermionic atoms with Feshbach resonance like potassium [119] loaded in a spin-independent π-flux square lattice and perceiving the same laser-induced tunneling term. The interactions are now on site. A third possibility would be to incorporate in the set-up in [100] an additional lattice that hosts bosonic atoms that mediate the interaction between fermions in the same spin state as in [104]. Again the interaction between the bosonic and fermionic species needs to possess a Feshbach resonance that allows to tune the scattering length appropriately.
In [100] we proposed an experimental scheme to measure G(ω) by using the one-particle excitation spectroscopy [120], which corresponds to a frequency-resolved transfer of the atoms to initially unoccupied auxiliary band of negligible width. In the weak-coupling limit the number of atoms transferred to the auxiliary band as a function of frequency detuning ω is proportional to G(ω). A similar scheme can be adopted here if the interactions in the auxiliary excited state are negligible.

Pairing function ∆(U, T )
In Sec. 3.2 we argue that the pairing strength∆ R = ξ −1 m ∆ R seen by an accelerated observer does not depend on the Unruh temperature, and is the same for all inequivalent observers in the Rindler Universe∆ R = ∆. In this section we shall consider how the Minkowski pairing function ∆ depends on the physical temperature T and the interaction strength |U |.
At the half filling for the noninteracting system U = 0, we have ∆ = 0 and the ground state of the Minkowski Dirac Hamiltonian (23) is a Dirac semimetal with two valence and conduction bands touching at Dirac points [121][122][123][124][125]. Similarly to a standard BCS theory [113,114], in a half-filled attractive Fermi Hubbard model on a square lattice even arbitrarily small attractive interactions U < 0 give rise to a nonzero ∆ pairing [126][127][128]. On the contrary, it is known that below a critical interaction strength |U c | = 0, Dirac fermions on a honeycomb lattice do not form Cooper pairs, since the density of states vanishes linearly at Dirac points [129][130][131][132][133]. At T = 0 and at the critical interaction |U c | the system undergoes the quantum phase transition between semimetal and a paired superconductor [134], although different methods, i.e. mean-field, variational and Monte Carlo give different estimates of the critical interaction |U c | ∼ 2 − 5. Note that honeycomb and square lattices are bipartite and therefore the particle-hole transformation allows us to relate attractive and repulsive Fermi Hubbard models at the half-filling [135].
We analyze ∆(U, T ) dependence for two types of boundary conditions: (i) open in x and periodic in y (cylinder), (ii) periodic in both x and y (torus). We find that at T = 0 the boundary conditions have strong affect the pairing properties of finite systems. In particular, we show that arbitrarily small U < 0 gives rise to pairing in a small finite system on a torus.

Solution on a cylinder
Here we study in details the behavior of BCS pairing for the Hamiltonian (35) on a cylindrical square lattice. For concreteness, we present and discuss the numerical results for a cylinder of size 50 × 50 lattice system (open in x direction).  In particular, we find that: (i) U c increases with the temperature T (left), (ii) the temperature phase transition is more evident when the interaction strength |U | increases (right). Near T = T c for U U c qualitatively recover the critial behavior known from the standard BCS theory, see the main text.
In Fig. 6 (left panel) we plot the pairing gap as a function of the interaction strength for several temperatures T . Our numerical results show that the pairing properties of the system described by the Hamiltonian (35) are similar to the Hubbard model on the honeycomb lattice, as one could expect. The quantitative differences might be due to the different dispersion relations away from the Dirac points. At T = 0 we recover the critical value |U c | ≈ 3 as obtained for an attractive π-flux model [124]. For T > 0, the critical interaction increases.
We plot the temperature dependence of the pairing function in Fig. 6 (right panel). As expected, we find that the pairing gap ∆(T ) diminishes with the increasing temperature, and at some point we reach a normal unpaired state. In finite systems, it is always the crossover. Nevertheless, we observe that with increasing U the behavior of ∆(T ) starts to resemble a phase transition known from the standard BCS theory. It is known that in the standard BCS theory, the near critical behavior is universal [136] and where A ≈ 3.07 T c , B ≈ 1.764 T c , A/B ≈ 1.74 are independent of material. We find that the numerical results on Fig. 6  One may question the validity of the meanfield approach for values of the interactions of the order of the band width. While we can expect (small) quantitative deviations with respect to more precise approaches like diagrammatic quantum MonteCarlo, the qualitative behavior is known to be well captured by the mean field approach.

Solution on a torus
Here we analyze a system described by the Hamiltonian (35) in the meanfield approximation on a torus. Interestingly, in this case we find that for finite systems, the pairing does happen for arbitrarily small |U | at T = 0. We comment more on this point at the end of this section. In this case the expression for ∆ can be computed analytically (see the Appendix A for details) and reads where is a positive quasiparticle eigenenergy of a system. From (77) we get where is finite in the thermodynamic limit N x , N y → ∞. The summation in (80) goes over all k in the Brillouin zone except for k ∈ {(0 0), (0 π), (π 0), (π π)}, since for ∆ = 0 the dispersion relation has four Dirac points which require a separate treatment as E k = 0.
First of all, let us analyze the equation (79) in thermodynamic limit.
which is finite for all T < ∞ and reproduces the numerical results for a system on a cylinder. Indeed, for T ≈ 0 we obtain Instead, for finite systems we have where now the second term cannot be dropped.
In particular, |U c (T )| ∝ T for small T and goes to zero in the limit T → 0+. Note that if one puts T = 0 in (79) before taking the limit ∆ → 0+ then the finite size result would be different.
Because lim T →0+ |U c (T )| = 0, the pairing of Dirac fermions on torus can take place for an arbitrarily small interaction strength |U |. This apparent discrepancy between the two finite-size solutions of the Dirac Hamiltonian (23) with different boundary conditions can be explained as following. From (77) at T = 0 the highest contribution to ∆ comes from the smallest positive eigenvalues E k . In particular, for small ∆ ≈ 0 the highest contribution comes from the Dirac cones K ∈ {(0 0), (0 π), (π 0), (π π)} (the density distribution of Cooper pairs is discussed in details in the Appendix B) which greatly affects small systems but is irrelevant in the thermodynamic limit, see Fig. 7. At the same time, this argument does not apply to the solution on the cylinder, as for even number of points N the two bands touch only in a thermodynamic limit, cf. Fig. 6.
Different than for the naive Dirac Hamiltonian on a square lattice (23), the boundary conditions do not play a significant role for the honeycomb lattice, where the critical interactions is nonzero |U c | > 0 also on a torus. Again, such behavior admits a simple explanation in terms of the band structure of the model. The graphene dispersion relation is minimalized only in the thermodynamic limit as the Dirac cones' coordinates K = ±(0, 4π/(3 √ 3)) are not integer multiple of π. Thus, as it happens for finite cylindrical lattices the Dirac points do not contribute to the computation of ∆ also in finite toroidal honeycomb lattices.

Conclusions and Outlook
In this paper we have investigated and proposed the experimental observation of the Unruh effect for interacting particles with ultracold fermions in optical lattices by a quantum quench. We have shown that achieving tunable Lorentz-preserving interactions with fermionic atoms is possible. Thus, it is possible to simulate an accelerated observer in an interacting background by simulating the corresponding Hamiltonian in Rindler space. Observing the Unruh effect reduces then to measuring the Wightman response function in Rindler spacetime for the interacting background, the ground state of the interacting Dirac Hamiltonian in Minkowski space.
We have studied the Wightman response function detected by an accelerated observer for attractive relativistic interactions in the meanfield approximation, for varying interactions and real temperature T of the background. In this approximation, the Unruh effect results from the interplay between the two different Bogoliubov transformations that relate the notion of particles for inertial and accelerated observers, and of particles and BCS quasi-particles, respectively. In the low-energy limit, in which the lattice system is with good approximation relativistic invariant, we have found that the Wightman function (precisely its power spectrum) displays the Planckian spectrum characteristic of the Unruh effect with a peculiarity. When the interactions grow, up to dominate over the Unruh temperature, there is a crossover between normal and "double" inversion of statistics determined by the (bosonic) Cooper pairs. Remarkably, in the low-energy limit our meanfield lattice calculations for interacting fermions not only give that the response is thermal, but also that the functional relation between the Unruh temperature T U and the proper acceleration a is the same as for a free theory, Such finding is in agreement with the predictions of the Bisognano-Wichmann theorem, which are valid under very general assumptions for any relativistic quantum field theory. On one hand, such fast convergence of lattice calculations to the correct relativistic behavior indicates that our lattice approach to the quantum simulation of quantum field theories in curved spacetime is promising. On the other hand, it can be read as a further evidence of the robustness of the Bisognano-Wichmann predictions recently observed in [137], where the equivalence between entanglement Hamiltonian and the Hamiltonian perceived by accelerated observers is used to access the entanglement spectrum of lattice models.
The present paper opens up interesting perspectives. It offers for instance a natural setup for testing quantum thermometry [138] and quantum thermodynamics [139] in curved spacetimes in presence of interactions (for a review on recent trends and developments in quantum thermodynamics see e.g. [140]).
It offers also a experimental playground for toy models of Lorentz-violating and trans-Planckian physics [141]. Indeed, as we argue above there is a tight relation between Unruh effect and basic principles of relativistic quantum field theory as Lorentz invariance and locality. This tight relation explains the universality of the thermal behavior of the Wightman function [142,143]. Conversely, deviations from such behavior can signal e.g. the breaking of Lorentz invariance [144], the deformation of the uncertainty principle [145], or even open a window on quantum gravity [146]. Ultracold atom simulators of quantum interacting matter in artificial curved spacetime may serve for analyzing/testing such scenarios (cf. with Lorentz violations in neutrino physics [147,148]).
Beyond the Unruh effect, another interesting direction for the quantum simulator we propose here is the study of dynamical chiral symmetry breaking in curved spacetime [149]. Indeed, by considering further atomic species we can in principle include more than one flavour and engineer the Gross-Neveu model [150] in an arbitrary optical (or more complicated) metric.
Last but not least, our study can be seen a further important step in the long journey to the simulation of self-gravitating quantum manybody physics.

A Analytic equation for ∆ on torus
In this section we seek for a fully periodic solution of the mean-field Minkowski Hamiltonian at the half filling. After expressing the field operators in the momentum space, we can write where ψ † k = c † ↑, k c † ↓, k , and where ∆ = |U | c ↑,m,n c ↓,m,n , g k = 2 t (sin(k x ) sin(k y )), and σ = (σ x σ y ), such that g k • σ = 2 t (sin(k x )σ x + sin(k y )σ y ). The eigenvalues of (86) read λ k = ± ∆ 2 + 4 sin 2 (k x ) + sin 2 (k y ) ≡ ±E k .
The two eigenvectors associated to the positive eigenvalue are with G k = 2t E k (sin(k x ) + i sin(k y )). Let us the write the field operator explicitly in terms of the eigenmodes where X
This simple observation tells us that the pairing is dominated by Dirac-like excitations near the Dirac cones. This is the reason why the bosonic response (the Fermi Dirac plateau) only appears in the low frequency regime (Fig. 3).
We plot ρ(k x = 0, k y ) in Fig. 8. The figure shows that up to ∆ ≈ 1, with good approximation the composite bosons (Cooper pairs) are formed around Dirac cones only. First of all, it is easy to find out that a unitary operator where 1 X is the N x × N x -identity matrix over the lattice coordinate m along x, transforms the Hamiltonian as U † 1 H A ky U 1 = −H A * -ky .
Thanks to this symmetry, we obtain in (46) and (52)  . Additionally, after making a gauge choice such that ∆ ∈ R, then we see that there is another symmetry of the Hamiltonian where At the half filling, we also have where All four possible combinations of U 1 , U 2 , U 3 , i.e. U 1 U 2 , U 1 U 3 , U 2 U 3 and U 1 U 2 U 3 are necessarily symmetries of H A ky . In particular, defining U 4 = U 1 U 3 we see therefore, we conclude H A ky and H A -ky have the same spectrum.