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

TL;DR

This paper introduces a framework linking quantum field theories in curved spacetime to quantum many-body systems, enabling simulations of phenomena like Hawking radiation and black hole chaos using experimental setups such as trapped ions.

Contribution

It establishes a one-to-one correspondence between quantum fields in curved spacetime and specific quantum many-body models, with potential experimental realizations and applications in simulating black hole physics.

Findings

Massless scalar fields map to site-dependent bosonic hopping models.

Massless Dirac fields map to Hubbard or XY models.

Black holes are identified as the fastest scramblers and most chaotic systems.

Abstract

This paper proposes a new general framework to build a one-to-one correspondence between quantum field theories in static 1+1 dimensional curved spacetime and quantum many-body systems. We show that a massless scalar field in an arbitrary 2-dimensional static spacetime is always equivalent to a site-dependent bosonic hopping model, while a massless Dirac field is equivalent to a site-dependent free Hubbard model or a site-dependent isotropic XY model. A possible experimental realization for such a correspondence in trapped ions system is suggested. As applications of the analogue gravity model, we show that they can be used to simulate Hawking radiation of black hole and to study its entanglement. We also show in the analogue model that black holes are most chaotic systems and the fastest scramblers in nature. We also offer a concrete example about how to get some insights about quantum many-body systems from back hole physics.

Equations

$$\text{d}{s}^2=f(x)\text{d}{t}^2-f(x{)}^{-1}\text{d}{x}^2.$$

$$g_h=\frac{1}{2}{f}^{\prime }(x_h)>0.$$

$$t\to v,s.t.,v=t+\int f(x{)}^{-1}\text{d}x.$$

$$\text{d}{s}^2=f\text{d}{v}^2-2\text{d}v\text{d}x.$$

$$m^2\varphi -2{\partial }_v{\partial }_x\varphi -f^{\prime }{\partial }_x\varphi -f{\partial }_x^2\varphi =0.$$

$$m\phi =2{\partial }_v\varphi +f{\partial }_x\varphi ,$$

$${\partial }_v\varphi =-\frac{f}{2}{\partial }_x\varphi +\frac{m\phi }{2},{\partial }_x\phi =m\varphi .$$

$${\partial }_{v}w=-\frac{f}{2}{\partial }_{x}w-\frac{f^{\prime }}{4}w+\frac{m\phi }{2\sqrt{f}},{\partial }_x\phi =mw\sqrt{f}.$$

$${\partial }_{v}w=-\frac{1}{4}[{\partial }_x(fw)+f{\partial }_{x}w].$$

$$i{\gamma }^a{e^{\mu }}_a{\partial }_{\mu }\psi +\frac{i}{2}{\gamma }^a\frac{1}{\sqrt{-g}}{\partial }_{\mu }(\sqrt{-g}{e^{\mu }}_a)\psi -m\psi =0.$$

$${e^{\mu }}_a=\left[\begin{array}{cc}-1,& 1\\ -\frac{f}{2}+\frac{1}{2},& \frac{f}{2}+\frac{1}{2}\end{array}\right]$$

$$\psi =\frac{1}{\sqrt{2}}\left[\begin{array}{c}u+w\\ u-w\end{array}\right]$$

$${\partial }_{v}w=-\frac{f}{2}{\partial }_{x}w-\frac{f^{\prime }}{4}w+\frac{i}{2}mu,{\partial }_{x}u=-imw.$$

$$f_n=f(nd),w_n(v)=w(v,x_n).$$

$$i\frac{\text{d}}{\text{d}v}{\tilde{w}}_n=-{\kappa }_n{\tilde{w}}_{n-1}-{\kappa }_{n+1}{\tilde{w}}_{n+1}-\mu {\tilde{w}}_n.$$

$${\kappa }_n=\frac{f_n+f_{n-1}}{8d}\approx \frac{f[(n-1/2)d]}{4d}.$$

$$[{\hat{a}}_n,{\hat{a}}_m^{†}]={\delta }_{nm},[{\hat{a}}_n,{\hat{a}}_m]=[{\hat{a}}_n^{†},{\hat{a}}_m^{†}]=0.$$

$$i\frac{\text{d}}{\text{d}v}{\hat{a}}_n=-{\kappa }_n{\hat{a}}_{n-1}-{\kappa }_{n+1}{\hat{a}}_{n+1}-\mu {\hat{a}}_n.$$

$$\mathcal{H}=\sum _n\left[-{\kappa }_n\left({\hat{a}}_n^{†}{\hat{a}}_{n-1}+{\hat{a}}_{n-1}^{†}{\hat{a}}_n\right)-\mu {\hat{a}}_n^{†}{\hat{a}}_n\right].$$

$$\{{\hat{c}}_n^{†},{\hat{c}}_m\}={\delta }_{nm},\{{\hat{c}}_n,{\hat{c}}_m\}=\{{\hat{c}}_n^{†},{\hat{c}}_m^{†}\}=0,$$

$$\mathcal{H}=\sum _n\left[-{\kappa }_n\left({\hat{c}}_n^{†}{\hat{c}}_{n-1}+{\hat{c}}_{n-1}^{†}{\hat{c}}_n\right)-\mu {\hat{c}}_n^{†}{\hat{c}}_n\right].$$

$${\sigma }_n^{+}=\exp \left[i\pi \sum _{j=1}^{n-1}{\hat{c}}_j^{†}{\hat{c}}_j\right]{\hat{c}}_n^{†},{\sigma }_n^{-}=\exp \left[-i\pi \sum _{j=1}^{n-1}{\hat{c}}_j^{†}{\hat{c}}_j\right]{\hat{c}}_n$$

$$\mathcal{H}=\sum _n\left[-{\kappa }_n\left({\sigma }_n^{+}{\sigma }_{n-1}^{-}+{\sigma }_{n-1}^{+}{\sigma }_n\right)+\frac{1}{2}\mu {\sigma }_n^z\right].$$

$${\sigma }_n^x={\sigma }_n^{+}+{\sigma }_n^{-},{\sigma }_n^y=-i({\sigma }_n^{+}-{\sigma }_n^{-})$$

$$\mathcal{H}=\frac{1}{2}\sum _n\left[-{\kappa }_n\left({\sigma }_n^x{\sigma }_{n-1}^x+{\sigma }_n^y{\sigma }_{n-1}^y\right)+\mu {\sigma }_n^z\right].$$

$$V_0=\frac{1}{2}m\sum _{i=1}^N({\omega }_x^2{x}_i^2+{\omega }_y^2{y}_i^2+{\omega }_z^2{z}_i^2)$$

$$V_C=\sum _{i>j}^N\frac{e^2}{\sqrt{(z_i-z_j{)}^2+(x_i-x_j{)}^2+(y_i-y_j{)}^2}}$$

$$V_L=\sum _{j>i}^N{t}_{i,j}(a_i^{†}{a}_j+a_i{a}_j^{†})$$

$$H=\sum _{i=1}^N{\omega }_x{a}_i^{†}{a}_i+\sum _{j>i}^N{t}_{i,j}(a_i^{†}{a}_j+a_i{a}_j^{†}),$$

$$P(E)\propto e^{-E/T_H}.$$

Full text

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 Jacobson (2003); Hollands and Wald (2015) 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 Unruh (1981, 1995), 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 Weinfurtner et al. (2011), Bose-Einstein condensate (BEC) Garay et al. (2000); Steinhauer (2016); de Nova et al. (2019), optic systems Sheng et al. (2013); Bekenstein et al. (2017); Drori et al. (2019); Sheng et al. (2018); Zhong et al. (2018), ultracold atoms in optical lattices Pedernales et al. (2018) and so on. See Refs. Barcelo et al. (2005); Faccio et al. (2013); Barceló (2018) 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 Tasaki (1998), 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. Hu et al. (2019); Trautmann and Hauke (2016); Schutzhold and Unruh (2013)). Let us notice that over the past decades, a remarkable progress in gravity and relevant fields is the proposal of AdS/CFT correspondence Maldacena (1999); Witten (1998); Gubser et al. (1998), 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 Ryu and Takayanagi (2006); Nishioka et al. (2009). Recently based on the AdS/CFT correspondence, quantum scrambling has been suggested as a powerful tool for characterizing chaos in black holes Hayden and Preskill (2007); Shenker and Stanford (2014), and Refs. Sekino and Susskind (2008); Shenker and Stanford (2014); Maldacena et al. (2016) 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-to-one 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,

[TABLE]

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 $x=x_{h}$ such that $f(x_{h})=0$ but

[TABLE]

$g_{h}$ is the surface gravity of the horizon, which gives the Hawking temperature $T_{H}=g_{h}/(2\pi)$ 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,

[TABLE]

The metric (1) in the infalling Eddington-Finkelstein coordinates $\{v,x\}$ becomes

[TABLE]

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

[TABLE]

By introducing the variable $\varphi$

[TABLE]

Eq. (4) can be rewritten into two coupled 1st order equations

[TABLE]

Now make variable transformation $\phi=w\sqrt{f}$ and we can rewrite the above equations into the following forms111This equation is singular at horizon due to $f(x_{h})=0$. A discussion about this point can be found in our appendix A.

[TABLE]

In the massless limit $m\rightarrow 0$, the above two equations decouple and there is only one independent evolutional equation,

[TABLE]

A similar result can also be obtained for Dirac field. The Dirac equation with the general vielbein ${e^{\mu}}_{a}$ and metric $g_{\mu\nu}$ can be written as Mann et al. (1991); Pedernales et al. (2018)

[TABLE]

Here $g$ is the determinate of metric $g_{\mu\nu}$. The $\gamma$-matrices in the two-dimensional case are chosen such that $\gamma^{a}=(\sigma_{z},i\sigma_{y})$. Choose the vielbein to be

[TABLE]

and take the decomposition

[TABLE]

into account, we find that there are two independent equations

[TABLE]

In the massless limit $m\rightarrow 0$, there is only one evolutional equation remained, which is the same as Eq. (8).

III Map into quantum many-body systems

III.1 Theory model

Now let us discretize the system. The spatial position is discretized as $x=x_{n}=nd$ with $n\in\mathbb{N}$ and $d\ll\lambda_{0}$, where $\lambda_{0}$ is the effective average wavelength in the system. The functions in the fixed spacetime are then transformed into discrete forms as follows,

[TABLE]

The spatial derivatives in Eq. (8) are approximated by central differences. Upon a variable transformation $w_{n}=(-i)^{n}e^{-i\mu v}\tilde{w}_{n}$, Eq. (8) can be rewritten into the following form

[TABLE]

with

[TABLE]

Here $\mu$ 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 field $\tilde{w}_{n}$ into operator. For bosonic field, we use the replacement $\tilde{w}_{n}\rightarrow\hat{a}_{n}/\sqrt{d}$ and introduce bosonic commutators such that

[TABLE]

The evolutional equation for the field operator then reads,

[TABLE]

Considering the evolutional equation in Heisenberg picture $i\partial_{v}\hat{a}_{n}=[\hat{a}_{n},\mathcal{H}]$, Eq. (13) implies a following Hamiltonian

[TABLE]

This Hamiltonian describes a bosonic hopping model and can be treated as a limit case of a certain of different well-studied quantum systems. For example, in condensed matter systems, it is the Bose-Hubbard model Gersch and Knollman (1963); Ma et al. (1986); Giamarchi and Schulz (1988); Fisher et al. (1989) with site-dependent hopping amplitude and zero on-site self-interaction.

For the Dirac field, we can do the similar thing. Take the replacement $\tilde{w}_{n}\rightarrow\hat{c}_{n}/\sqrt{d}$ and introduce anti-commutators such that

[TABLE]

we can obtain a following Hamiltonian form

[TABLE]

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 Hensgens et al. (2017); Tarruell and Sanchez-Palencia (2018); Salfi et al. (2016) , for instance.

The Hamiltonian (15) can also be rewritten into other well-studied model in condensed matter physics: the isotropic XY model Barouch et al. (1970); Barouch and McCoy (1971). To do that, let us introduce the following operators according to Jordan-Wigner transformation

[TABLE]

and $\sigma_{n}^{z}=1-2\hat{c}_{n}^{\dagger}\hat{c}_{n}$ with the periodic/anti-periodic boundary condition. Upon a constant, Hamiltonian (15) can be rewritten as

[TABLE]

Now introduce the Pauli matrices

[TABLE]

then the above Hamiltonian reads,

[TABLE]

This is nothing, but the isotropic XY model with site-dependent hopping.

III.2 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 $H=V_{0}+V_{C}+V_{L}$, where

[TABLE]

[TABLE]

[TABLE]

where $\omega_{\alpha},\alpha=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 $\omega_{x,y}\gg\omega_{z}$, the ions form a chain along the $z$ axis and occupy equilibrium positions. Phonons in the $z$ direction can be described approximately by D. and I. (2004)

[TABLE]

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.

IV Applications in black hole physics

IV.1 Hawking radiation and its entanglements

In this section, we will use the above quantum many-body 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)=\alpha\tanh x$ and $d=0.1$. In this case, there is a horizon at $x=x_{h}=0$ with the Hawking temperature $T_{H}=\alpha/(4\pi)$. It is worth noting that $\kappa_{n}\neq 0$ at the horizon though $f(x)$ is zero at the horizon. Without lose of generality, we set $\mu=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 $|\Psi(0)\rangle=|e_{n_{0}}\rangle$ and choose $n_{0}=-2/d$ as an example. It describes an initial particle which is localized at the $n_{0}$-th 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. Damour and Ruffini (1976); Parikh and Wilczek (2000); Arzano et al. (2005), 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,

[TABLE]

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\lesssim\mathcal{O}(2\pi\alpha/(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 $S(v)=-\text{Tr}[\rho(v)\ln\rho(v)]$ and the reduced density matrix for outer region is given by $\rho(v)=\text{Tr}_{\text{inner}}(|\Psi(v)\rangle\langle|\Psi(v)|)$. 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.

IV.2 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 AdS2 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 $\hat{W}(t)$ and $\hat{V}(t)$ in the Heisenberg picture, their OTOC is typically defined as

[TABLE]

Here $\hat{W}(0)$ and $\hat{V}(0)$ can be same or different, $\langle\cdot\rangle$ stands for average in an initial state. Ref. Maldacena et al. (2016) 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

[TABLE]

where $\lambda_{L}$ is the Lyapunov exponent and satisfies following “chaos bound”,

[TABLE]

Here $T$ is the temperature of the system. The exponential growth (25) will be broken after the “scrambling time”

[TABLE]

Here $N_{f}$ stands for the degrees of freedom of the system. It is conjectured in Sekino and Susskind (2008); Shenker and Stanford (2014); Maldacena et al. (2016) 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,

[TABLE]

Here $\hat{N}_{n_{0}}$ is a local operator associated to particle number operator at the $n_{0}$-th site,

[TABLE]

Here $l_{0}$ has the length scale and stands for the width of distribution of $\hat{N}_{n_{0}}$. The time evolutional operator $\hat{N}_{n}$ is given by Heisenberg picture $\hat{N}_{n}(v)=\exp(-i\mathcal{H}v)\hat{N}_{n}\exp(i\mathcal{H}v)\,.$ The reason we use the Eq. (29) to define the local operator $\hat{N}_{n_{0}}$ rather than $\hat{N}_{n_{0}}=\hat{a}_{n_{0}}^{\dagger}\hat{a}_{n_{0}}$ is that Eq. (29) is a well-defined smooth local operator in the continuous limit $d\rightarrow 0$. Instead, $\hat{N}_{n_{0}}=\hat{a}_{n_{0}}^{\dagger}\hat{a}_{n_{0}}$ will become a $\delta$-like function in continuous, which is singular. The initial state is a thermal state with the temperature same as the temperature of black hole

[TABLE]

Here $Z$ is the normalized factor which insures Tr$(\rho)=1$ and the summation contains all the positive energy modes of outside Hamiltonian $\mathcal{H}_{\text{out}}$ (as the negative modes are assumed to fall into black hole). $\mathcal{H}_{\text{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 define $\tilde{C}(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\pi 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}$ Maldacena et al. (2016), where $G$ is the Newton’s constant. Here we neglect the backreaction of matters on geometry, which means the limit of $G\rightarrow 0$. Thus, in principle, the OTOC will increase forever, i.e., $t_{*}\rightarrow\infty$. 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}$ of $\hat{N}_{n_{0}}$. This suggests that the time scale of exponential growth will become infinity in the continuous limit $d\rightarrow 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 $\Delta$ to be the mean value of $E_{i+1}-E_{i}$ and $\mathcal{N}$ to be the total number of energy levels. Then we define $\mathcal{N}P(s)\delta s$ to be number of energy levels $E_{i}$ which satisfy $s\leq(E_{i+1}-E_{i})/\Delta\leq s+\delta s$. 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}$ Guhr et al. (1998). 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 Jafarizadeh et al. (2012). 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}\neq 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}\neq 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}\neq 0$, a black hole metric is encoded into the sit-dependent coupling $\kappa_{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-to-one 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

[TABLE]

we can find that

[TABLE]

The positive energy (measured by $v$) solution is

[TABLE]

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

[TABLE]

The function $F(x)$ is continuous at the horizon. The divergence has been absorbed into the logarithm function. The solution (33) can be separated into two pieces,

[TABLE]

for $x<x_{h}$ and

[TABLE]

for $x>x_{h}$. The tunneling rate then reads,

[TABLE]

Following the argument in Ref. Damour and Ruffini (1976), the two pieces of the solution in Eqs. (35) and (36) 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. (35) will obtain an additional phase factor and so we can obtain the following relationship

[TABLE]

Take it into Eq. (37) and we then obtain

[TABLE]

As expected, the tunneling rate and energy satisfy the blackbody spectrum and the temperature is just given by $T_{H}=g_{h}/(2\pi)$.

In physics, Eq. (33) 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. Unruh (1995), 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)=\alpha\tanh x$ Then we can see that the hopping amplitude reads

[TABLE]

There is a horizon at $x=x_{h}=0$ with the Hawking temperature $T_{H}=\alpha/(4\pi)$. It is worth noting that $\kappa_{n}\neq 0$ at horizon though $f(x)$ is zero at horizon. The numerical computation needs a finite cut-off $n=-L,-L+1,\cdots,L-1,L$. To match this cut-off, we have to set hopping amplitude $\kappa_{n}$ such that

[TABLE]

Without lose of generality, we set $\mu=0$ as the total particle number is conserved.

The Hamiltonians for inner region and outer region are

[TABLE]

and

[TABLE]

Note that the total Hamiltonian is not the sum of inner part and outer part. In fact, we have

[TABLE]

Here $\mathcal{H}_{0}$ is the contribution at the horizon

[TABLE]

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\approx 5\times 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 $\hat{a}_{n}^{\dagger}\hat{a}_{n}$ as the basic vectors of Hilbert space

[TABLE]

which satisfy

[TABLE]

and

[TABLE]

Then we can write down the matrix elements of Hamiltonian.

[TABLE]

[TABLE]

and

[TABLE]

For a given initial state $|\Psi(0)\rangle$, the time evolutional state is given by $|\Psi(v)\rangle=e^{-i\mathcal{H}v}|\Psi(0)\rangle$. In Fig. 1, we take parameters $d=0.1,L=300$ and $\alpha=10$. The results are similar if we choose $N=2$. Due to the technical difficulties addressed above, we cannot explore larger $N$. Roughly speaking, as we consider the free theory in fixed background, many particles can be understood as a collection of single particle. Thus, the simplification here does not loss the essential physics.

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

[TABLE]

There is a horizon at $x=x_{h}$ with the Hawking temperature $T_{H}=x_{h}/(4\pi)$. In this case we make the cut-off in the following way

[TABLE]

Similar to the case in Hawking radiation, if we choose $N=1$, we still have Eqs. (44) and (45). Then we can write down the matrix elements of Hamiltonian

[TABLE]

and

[TABLE]

For a given initial state $|\Psi(0)\rangle$, the time evolutional state is given by $|\Psi(v)\rangle=e^{-i\mathcal{H}v}|\Psi(0)\rangle$ and $\hat{N}_{n}(v)=\exp(-i\mathcal{H}v)\hat{N}_{n}\exp(i\mathcal{H}v)$. Then we can obtain the OTOC in Eq. (27).

Acknowledgements.

R.-G. Cai was supported in part by the National Natural Science Foundation of China Grants Nos. 11690022, 11821505, 11851302, 11947302, 11745006, and by the Strategic Priority Research Program of CAS Grant No.XDB23030100, and by the Key Research Program of Frontier Sciences of CAS. H. Liu was supported by the National Natural Science Foundation of China Grants No.11690033. L. Luo is a member of the Indiana University Center for Spacetime Symmetries (IUCSS). L. Luo received supports from National Natural Science Foundation of China under Grant No. 11774436, Sun Yat-sen University Discipline Construction Fund, Sun Yat-sen University Three Major Construction Fund, and Guangdong Province Pearl River Youth Talent Fund.