Simulating quantum field theory in curved spacetime with quantum many-body systems

Run-Qiu Yang,1 Hui Liu,2 Shining Zhu,2 Le Luo,3 and Rong-Gen Cai4, ∗ Center for Joint Quantum Studies and Department of Physics, School of Science, Tianjin University, Yaguan Road 135, Jinnan District, 300350 Tianjin, P. R. China National Laboratory of Solid State Microstructures and School of Physics, Collaborative Innovation Center of Advanced Microstructures, Nanjing University, Nanjing, Jiangsu 210093, China School of Physics and Astronomy, Sun Yat-sen University, 519082, Zhuhai, China CAS Key Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China


I. INTRODUCTION
Quantum field theory in curved spacetime is a semiclassical approximation of quantum gravity theory, where the curved background spacetime is treated classically, while the matter fields in the curved spacetime are quantized. Although a fully successful quantum gravity theory is still not yet available, such a semiclassical approximation framework has offered us a large amount of interesting new phenomena, such as the Hawking radiation of black hole, particle production in an expanding universe etc., see Refs [1,2] for some review articles. Since in general these phenomena are extremal weak, they are extremely difficult to be observed in the real gravity situations. Analogues of black-hole or other phenomena in curved spacetimes in laboratory offer us new perspectives on quantum effects in curved spacetimes, which might help us deeply understand the nature of gravity. Following original work of Unruh [3,4], which studies Hawking radiation in sonic analog of a black hole, large amount of systems have been proposed and explored, such as surface wave in water flows [5], Bose-Einstein condensate (BEC) [6][7][8], optic systems [9][10][11][12][13], ultracold atoms in optical lattices [14] and so on. See Refs. [15][16][17] for reviews and references therein.
In spite of impressive progresses have been made in theoretically and experimentally on various analogue gravity systems, it is still interesting to seek some analogue models which are more "pure" in theory and more easily controlled in experiment. In condensed matter physics, there exist three basic quantum many-body models, the hopping model, Hubbard model and isotropic XY model [18], * cairg@itp.ac.cn which are of wide applications in many fields. In this paper we find that these models are also of interesting applications in quantum field theory in curved spacetime. The hopping rate in natural materials is constant, there are no enough motivations for physicists in condensed matter physics and materials to study the "site-dependent" hopping cases. Here we do see that the "site-dependent" hopping cases occur in the simulation of quantum field theories in curved spacetime by quantum many-body systems.
So far most of analogue gravity models have focused on the simulation of Hawking radiation or other types of spontaneous particle creation, such as Unruh effect, particle creation in the universe, dynamical Casimir effect and so on (see, e.g., Refs. [19][20][21]). Let us notice that over the past decades, a remarkable progress in gravity and relevant fields is the proposal of AdS/CFT correspondence [22][23][24], which says that a quantum gravity in the anti-de Sitter (AdS) spacetime is dual to a conformal field theory (CFT) living in the AdS boundary. The connection between geometry in the bulk and entanglement entropy in the boundary is also suggested in [25,26]. Recently based on the AdS/CFT correspondence, quantum scrambling has been suggested as a powerful tool for characterizing chaos in black holes [27,28], and Refs. [28][29][30] conjectured that black hole has the fastest scrambling and is a most quantum chaotic system in nature.
One of remarkable features of the AdS/CFT correspondence is the strong/weak duality: a weak gravity theory in the AdS bulk is equivalent to a strong coupled CFT in the AdS boundary. Although there exist many pieces of evidence to show the correspondence is true, it is extremely difficult, if it is not impossible, to prove the AdS/CFT correspondence. The analogue gravity models provide the possibility to test experimentally the AdS/CFT correspondence.
In this paper we will show that there exists a one-toone correspondence between quantum field theories in an arbitrary two dimensional spacetime and site-dependent bosonic hopping model, free Hubbard model or isotropic XY model in quantum many-body systems. As some applications of our analogue gravity model, we will study Hawking radiation of black hole and its entanglement, and show that black holes are most chaotic systems and the fastest scramblers in nature, predictions of the AdS/CFT correspondence. We also will use a concrete example to show how to use picture of back hole physics to learn something about quantum many-body systems.

II. QUANTUM FIELDS IN CURVED SPACETIME
We consider a 2-dimensional background spacetime with signature (+, −). In the static case, the metric can always be given in the Schwarzschild coordinates {t, x} as, In the most cases, we are interested in the static black hole spacetime with a single non-degenerated horizon, i.e., f (x) > 0 for x > x h and there is only a point at g h is the surface gravity of the horizon, which gives the Hawking temperature T H = g h /(2π) of the black hole. The metric (1) in the coordinates {t, x} is singular at the horizon. To overcome this shortage, we can define an infalling Eddington-Finkelstein coordinate by the coordinates transformation, The metric (1) in the infalling Eddington-Finkelstein coordinates {v, x} becomes In this case the metric has no longer the coordinate singularity at the horizon. Let us first consider a scalar field in the 2-dimensional curved spacetime. The Klein-Gorden equation of a complex scalar field in the metric (3) reads By introducing the variable ϕ Eq. (4) can be rewritten into two coupled 1st order equations Now make variable transformation φ = w √ f and we can rewrite the above equations into the following forms [31] ∂ In the massless limit m → 0, the above two equations decouple and there is only one independent evolutional equation, A similar result can also be obtained for Dirac field. The Dirac equation with the general vielbein e µ a and metric g µν can be written as [14,32] iγ a e µ a ∂ µ ψ + Here g is the determinate of metric g µν . The γ-matrices in the two-dimensional case are chosen such that γ a = (σ z , iσ y ). Choose the vielbein to be and take the decomposition into account, we find that there are two independent equations In the massless limit m → 0, there is only one evolutional equation remained, which is the same as Eq. (8).

A. Theory model
Now let us discretize the system. The spatial position is discretized as x = x n = nd with n ∈ N and d λ 0 , where λ 0 is the effective average wavelength in the system. The functions in the fixed spacetime are then transformed into discrete forms as follows, The spatial derivatives in Eq. (8) are approximated by central differences. Upon a variable transformation w n = (−i) n e −iµvw n , Eq. (8) can be rewritten into the following form with Here µ is an arbitrary constant. We will see later that it can be interpreted as the chemical potential in quantum many-body systems. Due to the discretization, discrete form is a well approximation for continuous fields if fields are slowly varying, i.e., Eq. (11) is valid in the low energy limit. Now let us quantize these fields themselves. This can be done by promoting fieldw n into operator. For bosonic field, we use the replacementw n →â n / √ d and introduce bosonic commutators such that The evolutional equation for the field operator then reads, Considering the evolutional equation in Heisenberg picture i∂ vân = [â n , H], Eq. (13) implies a following Hamiltonian This Hamiltonian describes a bosonic hopping model and can be treated as a limit case of a certain of different wellstudied quantum systems. For example, in condensed matter systems, it is the Bose-Hubbard model [33][34][35][36] with site-dependent hopping amplitude and zero on-site self-interaction. For the Dirac field, we can do the similar thing. Take the replacementw n →ĉ n / √ d and introduce anticommutators such that we can obtain a following Hamiltonian form This is just the free Hubbard model with site-dependent hopping. This model has been widely studied and can be realized in various different platforms, see Refs [37][38][39] , for instance. The Hamiltonian (15) can also be rewritten into other well-studied model in condensed matter physics: the isotropic XY model [40,41]. To do that, let us introduce the following operators according to Jordan-Wigner transformation and σ z n = 1 − 2ĉ † nĉn with the periodic/anti-periodic boundary condition. Upon a constant, Hamiltonian (15) can be rewritten as Now introduce the Pauli matrices then the above Hamiltonian reads, This is nothing, but the isotropic XY model with sitedependent hopping.

B. Experimental simulation
The Bose-Hubbard model in Eq. (14) can be realized in laboratory with various systems for implementing quantum simulation, such as optical lattices, superconducting qubits, and trapped ions etc.. Here we just concentrate to a simple case, which consists of a linear chain of ions in a linear Paul trap. In a linear trap, ions are arranged in a Coulomb chain. Assuming x as one of the transverse directions and z the trap axis, the Hamiltonian of the chain with N ions is where ω α , α = x, y, z are the trapping frequencies in each direction, V C is the Coulomb energy, while V L is the coupling between different axial modes, and t i,j are the hopping energies that are induced by a pair of Raman laser. For a linear trap ω x,y ω z , the ions form a chain along the z axis and occupy equilibrium positions. Phonons in the z direction can be described approximately by [42] Note that t i,j can be precisely adjusted to site(mode)dependent by varying the phase and the detuning of the Raman beams. By this scheme, we use the phonon modes of trapped ions to realize the Bose-Hubbard model with zero on-site energy. To simulate n-site Hubbard modes, we need to trap N ions in a linear trap and use N-1 pairs of laser to drive photon transitions between the N axial modes.

A. Hawking radiation and its entanglements
In this section, we will use the above quantum manybody model to study quantum aspects of black holes in gravity. Let us consider the bosonic hopping model as an example. For convenience in numerical computations, let us specify the function f (x) = α tanh x and d = 0.1. In this case, there is a horizon at x = x h = 0 with the Hawking temperature T H = α/(4π). It is worth noting that κ n = 0 at the horizon though f (x) is zero at the horizon. Without lose of generality, we set µ = 0 as the total particle number is conserved.
To study the black hole evaporation, we set a particle in the inner region of the black hole by initial state |Ψ(0) = |e n0 and choose n 0 = −2/d as an example. It describes an initial particle which is localized at the n 0th site. Based on the picture of "pair creation" in Hawking radiation, "particle-antiparticle pairs" can be created around the horizon. The antiparticle (negative energy) falls into the black hole and annihilates with this particle inside the black hole, the particle outside the horizon is materialized and escapes into infinity. Note that the pair creation/annihilation is a virtual processe, and the really materialized result is that the original particle inside the black hole disappears but an identical particle appears outside the horizon. This leads to an equivalent picture to understand Hawking radiation via quantum tunneling: the particle inside the horizon escapes to outside by quantum tunneling. According to Refs. [43][44][45], neglecting the back reaction of the radiation, the probability of finding this particle outside the horizon and its energy should obey the following blackbody spectrum, In Fig. 1 we show the numerical results about the probability of finding a particle of energy E n in the outer region. Here E n is the positive eigenvalue of Hamiltonian for the outer subregion. The evolutional time is chosen so that the radiation does not touch the cut-off boundary. For the details of numerical calculation, one may refer to the supplemental materials. We see that numerical results show that P (E) satisfies the blackbody spectrum Eq. (23) approximately with the temperature T = T H . Note that the numerical results for smaller energy deviate from the blackbody spectrum (23), because our finite size cut-off cannot cover the low energy region with E O(2πα/(Ld)) and so leads to the deviation. In addition, we also compute the entanglement entropy between the inner region and outer region, which is given by and the reduced density matrix for outer region is given by It shows that entanglement between the inner and outer regions increases during the Hawking radiation. Because there is only one particle in the black hole, the evaporation will stop in a short time and so the entanglement entropy saturates.

B. Quantum chaos and fastest scrambling
In this subsection, let us exhibit how to use our analogue model to study some new features of quantum field theory in curved spacetime: quantum chaos and fastest scrambling of black holes, appearing from the AdS/CFT correspondence. To supply an asymptotic AdS 2 black hole background, we consider f (x) = x 2 (1 − x h /x) as an example.
To describe the quantum chaos, it was proposed recently that the "out-time order correlation" (OTOC) may serve as a useful characteristic of quantum-chaotic behavior. For two local operatorsŴ (t) andV (t) in the Heisenberg picture, their OTOC is typically defined as HereŴ (0) andV (0) can be same or different, · stands for average in an initial state. Ref. [30] shows that, with a few general assumptions on the underlying field model and in thermal equilibrium state, the growth of a general OTOC C(t) satisfies where λ L is the Lyapunov exponent and satisfies following "chaos bound", Here T is the temperature of the system. The exponential growth (25) will be broken after the "scrambling time" Here N f stands for the degrees of freedom of the system. It is conjectured in [28][29][30] that black hole is a most chaotic system and has the fastest scrambling, i.e., it saturates the bounds (26) and (27). Now let us employ our model to check if it can exhibit the exponential growth of OTOC and gives us a positive Lyapunov exponent. As an example, we numerically study the following OTOC, HereN n0 is a local operator associated to particle number operator at the n 0 -th site, Here l 0 has the length scale and stands for the width of distribution ofN n0 . The time evolutional op-eratorN n is given by Heisenberg pictureN n (v) = exp(−iHv)N n exp(iHv) . The reason we use the Eq. (29) to define the local operatorN n0 rather thanN n0 = a † n0ân0 is that Eq. (29) is a well-defined smooth local operator in the continuous limit d → 0. Instead, N n0 =â † n0ân0 will become a δ-like function in continuous, which is singular. The initial state is a thermal state with the temperature same as the temperature of black hole Here Z is the normalized factor which insures Tr(ρ) = 1 and the summation contains all the positive energy modes of outside Hamiltonian H out (as the negative modes are assumed to fall into black hole). H out is obtained by only extracting the sites outside the horizon in Eqs. (15), (16) and (18). The time-evolution of C(v) is obtained numerically. The results are shown in Fig. 2. For convenience, we defineC(v) = C(v)/C(0.02), which does not change the slop of ln C(v). For the numerical details, one can refer to the appendix. We can observe that C(v) exponentially growes approximately in the early time. The slop of the fitting curve is found to be 2πT H approximately. The chaos bound (26) is saturated approximately.
In the pure gravity theory, the effective degree of freedom will be proportional to G −2 [30], where G is the Newton's constant. Here we neglect the backreaction of matters on geometry, which means the limit of G → 0. Thus, in principle, the OTOC will increase forever, i.e., t * → ∞. However, as we here use the lattice model, the operators and their commutators are bounded and so exponential growth will stop at a finite time. We study how the C(v) depends on the discrete distance d. The results show that the time scale of exponential growth will increase if we decrease d but fix the horizon radius x h and distribution width l 0 ofN n0 . This suggests that the time scale of exponential growth will become infinity in the continuous limit d → 0, as expected. Strictly speaking, to claim a system to be chaotic, either classically or quantum mechanically, the positive Lyapunov exponent is necessary but not sufficient. The positive Lyapunov exponent only indicates the sensitivity to the initial perturbations, which is the necessary condition of chaos. For example, in the classical case, we also require that the trajectory is dense in a neighborhood of phase space (i.e., ergodic). However, the linear analysis is enough to help us to find the Lyapunov exponents both in the classical case and quantum mechanical case. This can be understood by recalling the standard method in computing the Lyapunov exponent of classical chaotic systems. Thus, a linearized theory in a black hole background is enough to check the "chaos bound" (it may be more suitable to call it "bound on Lyapunov exponent").
In order to check if the models (15), (16) and (18) really contain chaotic behaviors when the coupling constants are given by a black hole metric, we study the statistics of "nearest-neighbor level spacing", which is an other characteristic quantity of chaotic system. We denote the energy levels of outside Hamiltonian to be E i with E i < E i+1 , which are obtained by directly diagonalizing the Hamiltonian numerically (the cut-off in high energy levels is needed as high energy levels have low accuracy and are not trustworthy in physics). Assume that ∆ to be the mean value of E i+1 − E i and N to be the total number of energy levels. Then we define N P (s)δs to be number of energy levels The function P (s) is called 'nearest-neighbor level spacing" function. It has been shown that if the system is integrable, P (s) satisfies Poisson statistics P (s) = e −s [46]. If the system is chaotic, P (s) will deviate from the Poisson statistics. For Gaussian orthogonal ensemble or Gaussian unitrary ensemble, P (s) is given by Wigner distribution. For other general cases, the P (s) may be given by general Brody distribution approximately [47]. In Fig (3) we show the numerical results about P (s). For the case that x h = 0, the effective spacetime has no black hole and we find that P (s) is given by a Poisson statistics approximately. However, once x h = 0, we find that a dramatic change happens and P (s) is no longer a Poisson distribution, which suggests that the system is not integrable.
From the viewpoint of quantum many-body theory, Hamiltonians shown in Eqs. (14), (15) and (18) contain only nearest hopping and quadratic interactions no matter how to design coupling constants. When we design the coupling constant by setting x h = 0, it is not easy to understand why these models can exhibit exponential growth of OTOC and why the systems have "integrable-nonintegrale phase transition". However, if we use our framework to convert them into an effective black hole, we immediately understand these properties. When x h = 0, a black hole metric is encoded into the sit-dependent coupling κ n and so the system has exponential growth of OTOC and signal of chaos. This offers an example about how to use our models to get some insights about the quantum many-body systems from the back hole physics.

V. SUMMARY
In summary, we have shown that a massless scalar/Dirac field in the static 1+1 dimensional curved spacetime can be simulated by some basic models in condensed matter physics: the bosonic hopping model, free Hubbard model and XY model. We suggested a possible experimental realization in trapped ions system for this analogue gravity model. As some applications of the analogue gravity model, we have numerically shown that this model can be used to simulate Hawking radiation. We have also checked the quantum chaos behave of black hole and verified that black hole is one of the most chaotic systems and has the fastest scrambling in nature. These are predictions of AdS/CFT correspondence. In this sense our model provides the possibility to test experimentally the correspondence. In addition, our results show the site-dependent hopping is one-toone related to spacetime point of curved background. By a concrete example, we showed how this framework can help us get some insights about the quantum many-body systems from the back hole physics. This not only provides new motivation to study the site-dependent hopping model, but also indicates a large number of applications in the analogue gravity model. These will bring us new viewpoints and phenomena of quantum many-body systems, and also will enlighten us to deeply understand the nature of gravity.
Appendix A: Tunneling rate and Hawking temperature In this appendix, let us show how to use the picture "quantum tunneling" to obtain the tunneling rate and the Hawking temperature of black hole.
To obtain the tunneling rate, we need to find the solution Eq. (8) with the energy E. By a variable transformation we can find that The positive energy (measured by v) solution is As f (x) = 0 at the horizon, this solution is not continuous at the horizon. Let us separate the integration in the above equation as follows The function F (x) is continuous at the horizon. The divergence has been absorbed into the logarithm function. The solution (A3) can be separated into two pieces, for x > x h . The tunneling rate then reads, Following the argument in Ref. [43], the two pieces of the solution in Eqs. (A5) and (A6) should be connected continuously under the bottom half of complex plane. Treat the piece of x < x h as the starting point and analytically continue it into the region of x > x h , the logarithm function in Eq. (A5) will obtain an additional phase factor and so we can obtain the following relationship Take it into Eq. (A7) and we then obtain As expected, the tunneling rate and energy satisfy the blackbody spectrum and the temperature is just given by T H = g h /(2π).
In physics, Eq. (A3) implies an infinite momentum at horizon and so will break down our condition for discretization. This belongs to the question named "trans-Planckian problem", which widely exists in all discussions of Hawking radiation. A particle emitted from a black hole with a finite frequency (measured at infinity), if traced back to the horizon, must have an infinite momentum, and therefore a trans-Planckian wavelength. The trans-Planckian problem is a mathematical artifact of horizon calculations. In all analogue models, when the emitted particle is near the "horizon", the smooth approximation is invalid and so a truncation is needed. However, it has been showed that the details of truncation will not change the behavior of Hawking radiation in "low energy" region (the "low energy" means the energy is low at infinity), see Ref. [4], for example.

Appendix B: Details of numerical simulations on
Hawking radiation Let us first explain how to make the numerical simulation on Hawking radiation. We take f (x) = α tanh x Then we can see that the hopping amplitude reads There is a horizon at x = x h = 0 with the Hawking temperature T H = α/(4π). It is worth noting that κ n = 0 at horizon though f (x) is zero at horizon. The numerical computation needs a finite cut-off n = −L, −L + 1, · · · , L − 1, L. To match this cut-off, we have to set hopping amplitude κ n such that Without lose of generality, we set µ = 0 as the total particle number is conserved. The Hamiltonians for inner region and outer region are and Note that the total Hamiltonian is not the sum of inner part and outer part. In fact, we have Here H 0 is the contribution at the horizon which mixes the inner region and outer region. Assume N to be the total particle number. It is difficult to simulate the dynamics for large N and L. For example, in the case 2L + 1 = N = 13, the dimension of total Hilbert space is D ≈ 5 × 10 6 . To simplify the issue in numerical algorithm, let us choose N = 1 and so the dimension of Hilbert space is D = 2L + 1. In this case, we can choose the eigenvectors ofâ † nân as the basic vectors of Hilbert space |e −L = (1, 0, · · · , 0) T , |e −L+1 = (0, 1, 0, · · · , 0) T , · · · , |e L = (0, 0, 0, · · · , 1) T , which satisfy e l |â † nân |e k = δ nl δ nk (B6) and e l |â † nân−1 |e k = δ n−1,k δ l,n .
Then we can write down the matrix elements of Hamiltonian.

Appendix C: Details of numerical simulations on OTOC
The simulation on OTOC is similar. We take f (x) = x 2 (1 − x h /x). Then we can see that the hopping reads There is a horizon at x = x h with the Hawking temperature T H = x h /(4π). In this case we make the cut-off in the following way n = 1, 2, · · · , 2L + 1, with Ld = x h .