# Emergent gravity from the correlation of spin-$\displaystyle\tfrac{1}{2}$
systems coupled with a scalar field

Quentin Ansel<EMAIL_ADDRESS>Institut UTINAM, CNRS UMR 6213,
Université Bourgogne Franche-Comté, Observatoire des Sciences de l’Univers
THETA, 41 bis avenue de l’Observatoire, F-25010 Besançon, France

(August 28, 2024)

###### Abstract

This paper introduces several ideas of emergent gravity, which come from a
system similar to an ensemble of quantum spin-$\displaystyle\tfrac{1}{2}$
particles. To derive a physically relevant theory, the model is constructed by
quantizing a scalar field in curved space-time. The quantization is based on a
classical discretization of the system, but contrary to famous approaches,
like loop quantum gravity or causal triangulation, a Monte-Carlo based
approach is used instead of a simplicial approximation of the space-time
manifold. This avoids conceptual issues related to the choice of the lattice.
Moreover, this allows us to easily encode the geometric structures of space,
given by the geodesic length between points, into the mean value of a
correlation operator between two spin-like systems. Numerical investigations
show the relevance of the approach, and the presence of two regimes: a
classical and a quantum regime. The latter is obtained when the density of
points reaches a given threshold. Finally, a multi-scale analysis is given,
where the classical model is recovered from the full quantum one. Each step of
the classical limit is illustrated with numerical computations, showing the
very good convergence towards the classical limit and the computational
efficiency of the theory.

###### Contents

  1. 1 Introduction
  2. 2 From continuous to discontinuous field theory with a Monte-Carlo approximation
  3. 3 Discrete quantum scalar field theory in curved space-time
    1. 3.1 Hamiltonian of a scalar field
    2. 3.2 Discretization procedure
  4. 4 Quantum space-time
  5. 5 Multi-scales analysis
    1. 5.1 Classical and non-classical gravitational states
    2. 5.2 Classical limit of the quantum scalar field
  6. 6 Conclusion
  7. A Normal law in a curved space-time
  8. B Details on the numerical optimizations
  9. C Covariant Monte-Carlo discretization
  10. D Hamiltonian and equation of motion of discretized systems

##  1 Introduction

Emergent gravity is a field of growing interest. The idea that the curved
space-time of general relativity is not fundamental but comes from an
effective theory has been explored in plenty of different approaches, such as
in string theory [1, 2, 3, 4], quantum graphity [5, 6, 7], the work of Carroll
and a.l. [8, 9, 10] and many other [11, 12, 13, 14]. Each approach is usually
well motivated, and it shields light on some key elements of what could be a
complete theory of quantum gravity. Unfortunately, these approaches are not
fully satisfactory yet [11]. If the lack of experimental evidence is left
aside,the current issues come partially from computational problems: from the
main equation it becomes hardly possible to make physical predictions that can
be tested [15], or it is arduous to connect altogether the main ingredients
expected in a complete theory of quantum gravity [16, 17, 18]. In this paper,
I propose yet another step forward in the direction of a complete theory of
emergent quantum gravity.

The ideas presented here do not cast exactly in the main theories, but it can
be seen at the crossing between several of them. One of the main sources of
inspiration is the work of Carroll and a.l. [8, 9, 10] In their series of
papers, they develop the idea that the Einstein field equation is encoded in
some way by the correlation between parts of the wave function describing the
entire universe. They postulate that the correlation, quantified by the
quantum mutual information, gives us areas between adjacent region, and they
are able to recover, at least partially, Einstein field equations. The
relation between the mutual information and space-time areas is strongly
motivated by different arguments, the main ones being that the entropy of the
horizon is proportional to its area, and the mutual information follows an
approximated area law [19]. Unfortunately, the inverse radon transform used to
recover the metric from the data of all areas is far from being simple and not
well understood in dimensions larger than 2. Moreover, the coupling with
matter fields is not straightforward, and it seems difficult to recover
quantum field theory in curved space-time with an adiabatic elimination [3].
Nonetheless, the idea that geometric quantities are encoded in the correlation
between subparts of a quantum state is interesting, and the fundamental
discrete nature of space-time follows other approaches, such as causal-
networks [20], loop quantum gravity [21, 22, 16], and quantum graphity [5, 6,
7]. In these approaches, space-time is given by "atoms of space" connected
together in a lattice. One drawback of these approaches is that the
connectivity between these atoms of space is not imposed by fundamental
principles, and one usually assumes 4 valent nodes (discretization with
tetrahedrons) or 6-valent nodes (discretization with cuboids), but this is
only motivated by a computational convenience. This discretization can be
thought of as a quantum gravity analog of the nearest neighbors approximation
of solid-state physique, but the ultimate theory must be free of the choice of
a lattice.

Based on these observations, I propose a model of emergent gravity for which
the main ingredients are the following: (i) the structure of space-time is
encoded in the correlation of subparts of a quantum state, each subpart being
tough as an "atom of space", a quantum region of space (ii) Correlations are
related to the geodesic distance between two quantum regions, not the area
(iii) Quantum field theory in curved space-time is recovered in the limit when
gravitational degrees of freedom are eliminated adiabatically [3, 4] and when
the size of the wave-packet is large compared to the Planck scale, and finally
(iv) the quantum state describing the quantum space-time is similar to an
ensemble of spin-$\displaystyle\tfrac{1}{2}$, a quantum region being described
by a spin, and geometric quantities are given by a spin-spin correlation. The
use of an ensemble of two-level systems is interesting because it limits the
size of the Hilbert space, but we can imagine rewriting the theory with
arbitrary spin-$\displaystyle j$ or even with harmonic oscillators.

This paper is organized as follows. In sec. 2, the Monte-Carlo discretization
procedure is introduced. In Sec. 3, the scalar field is introduced and it is
quantized using the material developed in Sec. 2. In Sec. 4, a quantum
description of the gravitational field is developed so that quantum field
theory in curved space-time is recovered in the appropriate limit. This part
of the theory is still significantly incomplete, but the minimum is made to
provide a well-defined theory for which explicit calculations en physical
predictions can be made. In Sec. 5, a multi-scale analysis is provided, where
the classical scalar field theory is recovered from the full quantum theory by
means of successive approximations. Each approximation step is illustrated
with numerical simulations. A conclusion and prospective views are given in
Sec. 6.

##  2 From continuous to discontinuous field theory with a Monte-Carlo
approximation

As outlined in the introduction, space-time is assumed fundamentally discrete,
and the continuous theory is obtained at a large scale, without assumptions on
a lattice structure at the Planck scale. Here, quantum regions are not
connected with the nearest neighbors, they are all connected together with
more or less weight. This is quite different from the typical situation
encountered in lattice quantum field theory (such as in Wilson’s approach of
path integrals [23, 24, 25, 26, 27], loop quantum gravity [22],…). The
underlying idea is illustrated in Fig. 1. Here, a probabilistic approach is
used, and integrals are discretized with a Monte-Carlo method [28].

a)  b)

Figure 1: a) A triangulated space. b) A discretized space with the vertices
connected all together (only the edges between $\displaystyle x_{2}$ and the
other vertices are shown) the weight between the vertices is a function of the
distance between the points. In the continuum limit, the weights define
density probability functions, depicted with colored gradients.

In the following, the points defining space are assumed uniformly distributed
on a manifold $\displaystyle\Sigma$, such that for any function $\displaystyle
f:\Sigma\rightarrow\mathbb{C}$, and any $\displaystyle R\subset\Sigma$ we have
[29, 30]

$\int_{R\in\Sigma}d^{3}x\sqrt{q(x)}f(x)=\lim_{N\rightarrow\infty}\frac{V_{R}}{N}\sum_{n=1}^{N}f(x_{n})$
(1)

where $\displaystyle V_{R}$ is the volume of $\displaystyle R$, and
$\displaystyle q$ is the metric on $\displaystyle\Sigma$. In this paper, the
Lebesgue measure is replaced by a probability distribution $\displaystyle
p_{x}(y)$. Using the fact that
$\displaystyle\int_{\Sigma}d^{3}y\sqrt{q(y)}p_{x}(y)=1$, the normalization
factor for an integration around the point $\displaystyle x$ is deduced to be

$\frac{V_{x}}{N}=\left(\sum_{n=1}^{N}p_{x}(y_{n})\right)^{-1}.$ (2)

This result is written for a distribution with support on
$\displaystyle\Sigma$, but it can be easily restricted to $\displaystyle
R\subset\Sigma$. This leads to a simple formula for the calculation of
integral without reference to a coordinate system:

$\begin{split}\mathbb{E}[f]_{x}&=\int_{\Sigma}d^{3}y\sqrt{q(y)}p_{x}(y)f(y)\\\
&\approx\frac{\sum_{n=1}^{N}f(y_{n})p_{x}(y_{n})}{\sum_{n=1}^{N}p_{x}(y_{n})}.\end{split}$
(3)

To be true, $\displaystyle p_{x}(y)$ must be gauge invariant. A natural
definition for such a quantity is based on the geodesic length between
$\displaystyle x$ and $\displaystyle y$, which can be computed without too
much difficulty as soon as the points are not too far from each other. In the
following, it is assumed that

$p_{x}(y)=\frac{\mathcal{N}(R)}{(\sqrt{2\pi}\varepsilon)^{3}}\exp\left(-\frac{\sigma(x,y)}{\varepsilon^{2}}\right)$
(4)

where $\displaystyle\varepsilon$ has the dimension of a length, and
$\displaystyle\sigma$ is Synge’s world function [31], which is equal to half
the square of the geodesic length, and $\displaystyle\mathcal{N}(R)$ is a
normalization coefficient that depends on the space-time curvature. It is
equal to 1 in a flat space. In this setting, $\displaystyle\varepsilon$
provides a natural notion of distance. It is not fixed by first principles,
but it is assumed of the order of the Planck length, i.e.,
$\displaystyle\varepsilon^{2}\sim\hbar G/c^{3}$. At a first sight, it is not
obvious that Eq. (4) gives a valid probability distribution in a curved space-
time. This aspect is explored in appendix A with a taylor expansion of Synge’s
world function in Riemann normal coordinates.

The discretization procedure described here above can be easily generalized to
the case when the distribution of points is not homogeneously distributed on
$\displaystyle\Sigma$, but is distributed homogeneously with respect to a
given coordinate system. In such a case, the previous formulas can be adapted
using the replacement rule $\displaystyle
p_{x}(y)\rightarrow\sqrt{q(y)}p_{x}(y)$.

##  3 Discrete quantum scalar field theory in curved space-time

###  3.1 Hamiltonian of a scalar field

Now that the general setting for the discretization of a field theory is made,
a specific example can be studied. As a toy model, a real scalar field
$\displaystyle\varphi$ is considered. The starting point of the presentation
is the Lagrangian density [32, 33]:

$\mathcal{L}=\frac{\sqrt{-g}}{2}\left[g^{\mu\nu}(\partial_{\mu}\varphi)(\partial_{\nu}\varphi)+m^{2}\varphi^{2}\right]$
(5)

where $\displaystyle g_{\mu\nu}$ is the metric tensor of the space-time
manifold $\displaystyle\mathcal{M}$, $\displaystyle\sqrt{-g}\equiv\sqrt{-\det
g}$, and $\displaystyle m$ is the mass. Both $\displaystyle g$ and
$\displaystyle\varphi$ depend on the space-time position $\displaystyle x$. In
the following, a 3+1 splitting is performed, to define a Hamiltonian [34, 21].
With this splitting, the metric is rewritten as follows: $\displaystyle
g_{00}=-\mathsf{N}+\mathsf{N}^{a}\mathsf{N}_{a}$, $\displaystyle
g_{0a}=\mathsf{N}_{a}$, and $\displaystyle g_{ab}=q_{ab}$, where coordinate
indices with Latin characters run over space coordinates only (i.e.,
$\displaystyle a,b=1,2,3$), and $\displaystyle q$ is the metric of a space
sub-manifold $\displaystyle\Sigma\subset\mathcal{M}$. With this notation, the
volume density is $\displaystyle\sqrt{-g}=\sqrt{q}\mathsf{N}$. The next step
is to determine the conjugated momentum of $\displaystyle\varphi$. Using
$\displaystyle\Pi=\tfrac{\partial\mathcal{L}}{\partial(\partial_{0}\varphi)}$,
one can determine

$\Pi=\sqrt{q}\mathsf{N}g^{0\nu}\leavevmode\nobreak\ \partial_{\nu}\varphi,$
(6)

and the Hamiltonian reads

$H=\frac{1}{2}\int_{\Sigma}d^{3}x\leavevmode\nobreak\
\left(\frac{\Pi^{2}}{\sqrt{q}\mathsf{N}}+\sqrt{q}\mathsf{N}\left[q^{ab}\partial_{a}\varphi\partial_{b}\varphi+m^{2}\varphi^{2}\right]\right).$
(7)

The Klein-Gordon equation in curved space-time is then recovered using
Hamilton’s equations $\displaystyle\tfrac{\delta
H}{\delta\varphi}=-\partial_{0}\Pi$ and $\displaystyle\tfrac{\delta
H}{\delta\Pi}=\partial^{0}\varphi$. To simplify the analysis, the coordinate
gauge choice $\displaystyle\mathsf{N}=1$ and $\displaystyle\mathsf{N}_{a}=0$
is made. This avoids the use of many subtleties of algebraic QFT in curved
space-time [32, 33], and thus, the derivation of a quantum theory is
significantly simplified. A first look at a covariant approach to the theory
is presented in Appendix C. The natural next step is to introduce complexifed
field variables. They are defined by

$\displaystyle\displaystyle
a(x)=\sqrt{\frac{m}{2}}\left(\varphi+\frac{i}{m\sqrt{q}}\Pi\right)$ (8)
$\displaystyle\displaystyle
a^{\dagger}(x)=\sqrt{\frac{m}{2}}\left(\varphi-\frac{i}{m\sqrt{q}}\Pi\right).$
(9)

With these new variables, the Hamiltonian takes the form

$\begin{split}H=&\int_{\Sigma}d^{3}x\sqrt{q}\leavevmode\nobreak\
\left(m\left(a^{\dagger}a+\frac{1}{2}\right)\right.\\\
&\left.+\frac{1}{4m^{2}}q^{ab}\partial_{a}(a^{\dagger}+a)\partial_{b}(a^{\dagger}+a)\right).\end{split}$
(10)

Contrary to the standard approach, the Hamiltonian is not diagonalized by the
complexified variables [35]. This enables the development of a quantum theory
where the gravitational field plays a role only in the coupling terms between
different space-time positions (the term with derivatives).

###  3.2 Discretization procedure

So far, the theory is given by a continuous space-time background. This
section is devoted to its discretization. The most subtle point concerns the
discretization of the coupling Hamiltonian
$\displaystyle\partial^{a}\varphi\partial_{a}\varphi$, which must return a
term of the form $\displaystyle\nabla_{a}\nabla^{a}\varphi$ in the equation of
motion, $\displaystyle\nabla$ being the covariant derivative 111In a curved
space-time, the covariant derivative $\displaystyle\nabla$ replaces the usual
derivative $\displaystyle\partial$, to take into account the modification of
the tangent space $\displaystyle T_{x}\mathcal{M}$ when a vector is moved in
the manifold. For a scalar field, $\displaystyle\nabla$ reduces to
$\displaystyle\partial$, but
$\displaystyle\partial_{a}\varphi=\nabla_{a}\varphi$ is a vector and its
derivative must be computed with $\displaystyle\nabla$ [34].. To apply the
procedure of Sec. 2, Stokes’s theorem [34] must be used and a field vanishing
on the boundary of the domain must be assumed. Therefore,

$\int d^{3}x\leavevmode\nobreak\ \sqrt{q}\mathsf{N}\leavevmode\nobreak\
\partial_{a}\varphi\partial^{a}\varphi\rightarrow-\int
d^{3}x\leavevmode\nobreak\ \sqrt{q}\mathsf{N}\leavevmode\nobreak\
\varphi\leavevmode\nobreak\ \nabla_{a}\nabla^{a}\varphi.$ (11)

with this setting, the Laplacian is related to the average of the field around
a given position. More precisely, one has

$\begin{split}\mathbb{E}[\varphi]_{x}&=\int_{\Sigma}d^{3}y\sqrt{q}\leavevmode\nobreak\
p_{x}(y)\varphi(y)\\\
&=\varphi(x)+\int_{\Sigma}d^{3}y\sqrt{q}\leavevmode\nobreak\
p_{x}(y)\leavevmode\nobreak\ y^{a}\partial_{a}\varphi(x)\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
+\frac{1}{2}\int_{\Sigma}d^{3}y\sqrt{q}\leavevmode\nobreak\
p_{x}(y)\leavevmode\nobreak\ y^{a}y^{b}\nabla_{a}\nabla_{b}\varphi(x)+...\\\
&=\varphi(x)+\frac{\varepsilon^{2}}{2}\nabla_{a}\nabla^{a}\varphi(x)+O(\varepsilon^{3}),\end{split}$
(12)

and as a consequence,

$\nabla_{a}\nabla^{a}\varphi(x)=\frac{2}{\varepsilon^{2}}\left[\mathbb{E}[\varphi]_{x}-\varphi(x)\right]+O(\varepsilon^{3}).$
(13)

This result is easily obtained in a flat space, but a few more cautions must
be taken in curved space. The validity of the $\displaystyle
O(\varepsilon^{3})$ approximation is explored in appendix A.

Next, the last equation is plugged into the integral over
$\displaystyle\Sigma$, and a Monte-Carlo approximation is made. The result is

$\int_{\Sigma}d^{3}y\sqrt{q}\mathsf{N}\varphi\nabla_{a}\nabla^{a}\varphi\approx\frac{2V}{\varepsilon^{2}N}\sum_{n=1}^{N}\mathsf{N}_{n}\left[\varphi_{n}\mathbb{E}[\varphi]_{n}-\varphi^{2}_{n}\right].$
(14)

The subscript $\displaystyle n$ is used to specify that quantities are
evaluated at the point $\displaystyle x_{n}$.

Complexified variables can be restored and Eq. (14) can be inserted into the
Hamiltonian (10). The theory is quantized straightforwardly by promoting the
complexified field variables into creation and annihilation operators acting
on a Fock space. The quantum Hamiltonian is

$\begin{split}\hat{H}=&\frac{V}{N}\sum_{n=1}^{N}\left(m\left(\hat{a}^{\dagger}_{n}\hat{a}_{n}+\frac{1}{2}\right)\right.\\\
&\left.-\frac{1}{2\varepsilon^{2}m}\left[(\hat{a}^{\dagger}_{n}+\hat{a}_{n})\mathbb{E}[(\hat{a}^{\dagger}+\hat{a})]_{x_{n}}-(\hat{a}^{\dagger}_{n}+\hat{a}_{n})^{2}\right]\right),\end{split}$
(15)

with $\displaystyle\hat{a}^{\dagger}_{n}$ and $\displaystyle\hat{a}_{n}$ the
creation/annihilation operators at the position $\displaystyle x_{n}$. Despite
the relatively simple expression of the Hamiltonian, the average
$\displaystyle\mathbb{E}$ hides a factor $\displaystyle V/N$, which plays the
role of an unknown coupling constant between oscillators. To circumvent this
issue, the Hamiltonian is factorized with a global factor $\displaystyle V/N$,
and at each position $\displaystyle x_{n}$ the volume is evaluated with the
relation (2). The global factor is only responsible for a rescale of the time
coordinate in the Schrödinger equation. The final result is:

$\begin{split}\hat{H}=&\frac{V^{2}}{N^{2}}\sum_{n=1}^{N}\left(\sum_{k=1}^{n}p_{x_{n}}(x_{k})\right)\left(m\left(\hat{a}^{\dagger}_{n}\hat{a}_{n}+\frac{1}{2}\right)\right.\\\
&\left.+\frac{1}{2\varepsilon^{2}m}(\hat{a}^{\dagger}_{n}+\hat{a}_{n})^{2}\right)\\\
&-\frac{1}{2\varepsilon^{2}m}\left((\hat{a}^{\dagger}_{n}+\hat{a}_{n})\sum_{k=1}^{N}p_{x_{n}}(x_{k})(\hat{a}^{\dagger}_{k}+\hat{a}_{k})\right).\end{split}$
(16)

This Hamiltonian shares very strong similarities with the Bose-Hubbard model
[36, 37, 38, 39, 40]. In particular, this model was used to simulate quantum
fields in curved space-times [41]. While the standard approach is to start
from a given lattice configuration of the Bose-Hubbard model and to compute
the resulting effective space-time, here the reverse approach is made.
However, Eq. (16) is free from the nearest neighbor approximation, and it
provides a simple mapping between the space-time geometry and the parameters
of the Hamiltonian.

##  4 Quantum space-time

In the previous sections, a formalism that allows us to discretize a quantum
scalar field in space was developed, where gravitational degrees of freedom
have an influence only through an interaction Hamiltonian that couples the
fields at different locations. With these materials at hand, a (partial)
quantum theory for the gravitational field consistent with QFT in curved
space-time can be constructed.

A closer look at Eq. (16) highlights that in the interaction Hamiltonian, the
geometry of $\displaystyle\mathcal{M}$ is hidden in $\displaystyle
P_{x_{n}}(x_{k})$. These weights are the quantities that must be quantized.

The quantization relies on several assumptions. In particular, the geometry of
the gravitational state is encoded in the correlation between different sites
of an ensemble of spin-$\displaystyle\tfrac{1}{2}$, a spin being associated
with a point $\displaystyle x_{n}$. Then, a general quantum state
$\displaystyle|\Psi_{G}\rangle$ is a linear combination of vectors of
$\displaystyle\mathcal{H}_{G}=\mathbb{C}^{2N}$.

Many different operators can be used to describe a two-point correlation
function. Here, the following choice is made:

$\hat{C}_{nk}=\hat{\sigma}_{+}(x_{n})\hat{\sigma}_{-}(x_{k})+\hat{\sigma}_{+}(x_{k})\hat{\sigma}_{-}(x_{n})$
(17)

with $\displaystyle\hat{\sigma}_{-}(x_{n})$ and
$\displaystyle\hat{\sigma}_{+}(x_{n})$ the spin creation/annihilation
operators associated with the point $\displaystyle x_{n}$. The mean value of
this operator, $\displaystyle
C_{nk}=\langle\Psi_{G}|\hat{C}_{nk}|\Psi_{G}\rangle$, gives us information on
how much two regions are correlated, and more specifically how much two
regions are entangled since a superposition of different states is necessary
to get $\displaystyle C_{nk}\neq 0$.

To relate correlations and geometric data of the manifold, I propose to
postulate the following relation between the correlation and the geodesic
distance:

$C_{nk}^{2}=\alpha\leavevmode\nobreak\
\exp\left(-\frac{\sigma(x_{n},x_{k})}{\varepsilon^{2}}\right),$ (18)

with $\displaystyle\alpha$ a proportionality coefficient. This definition,
motivated by a numerical investigations, requires the exclusion of the terms
$\displaystyle n=k$ in the sum. This is because
$\displaystyle\hat{\sigma}_{+}(x_{n})\hat{\sigma}_{-}(x_{n})+\hat{\sigma}_{+}(x_{n})\hat{\sigma}_{-}(x_{n})=\textrm{1}\\!\\!\\!\\!1_{2}$,
and thus, by taking the expectation value, it is not possible to recover the
coefficient $\displaystyle\alpha$. As a consequence, the condition
$\displaystyle p_{x_{n}}(x_{n})=0$ must be imposed in Eq. (16). This is not a
problem for the convergence of the Monte-Carlo approximation scheme since it
concerns a modification of the probability distribution whose Lebesgue measure
is equal to zero.

Eq. (18) has a very simple physical interpretation. The larger the distance
between $\displaystyle x_{n}$ and $\displaystyle x_{k}$, the smaller the
correlation. In order to suppress the quantum fuzziness at large scales, the
correlation must decay sufficiently fast. A normal law (see Eq. (4)) is
assumed for convenience, but it is not excluded that another distribution can
be more adapted to the problem. For example, $\displaystyle
C_{nk}^{2}\propto\exp(-\sqrt{2\sigma_{nk}}/\varepsilon)$ could be also a
viable choice. The validity of this ansazt is explored numerically in Fig. 2.
The approach of the numerical investigation is the following. A cost function
$\displaystyle\mathcal{F}=\frac{1}{N}\sum_{n=1}^{N}\sum_{n=1}^{N}|\alpha
e^{-\sigma_{nk}/\varepsilon^{2}}-C_{nk}^{2}(\Psi_{G})|^{2}$ is minimized with
respect to the vector components of $\displaystyle|\Psi_{G}\rangle$. For a
different number of points (i.e., a different number of spins), a quantum
state is optimized numerically such that $\displaystyle C_{nk}^{2}$ is as
close as possible to $\displaystyle\alpha e^{-\sigma_{nk}/\varepsilon^{2}}$.
The value of $\displaystyle\alpha$ is fixed by hand, to simplify the
optimization process. Specific details on these numerical computations are
gathered in Appendix B. Different kinds of configurations and spaces are
investigated, with flat and curved spaces, and with points organized randomly
or regularly. In all situations, very low values of $\displaystyle\mathcal{F}$
are obtained, except when the density of points is too large. For example, in
the case of Fig. 2 b), the value of $\displaystyle\mathcal{F}$ drops
significantly when $\displaystyle L$ reaches a given threshold ($\displaystyle
L/\varepsilon\approx 1.5$ for the cube and $\displaystyle L/\varepsilon\approx
1.2$ for the 3-simplex). The threshold depends on the number of spins and it
seems to correspond to $\displaystyle C_{nk}\approx 1/\sqrt{N}$. Further
investigations may clarify this point, but this numerical result has an
important physical consequence: there is a natural minimum scale for which we
can obtain a classical geometry. The value of the threshold depends on
$\displaystyle\alpha$, and thus it is not clear at this point in which exact
circumstances it is possible to encounter a quantum space-time. This issue
will be investigated elsewhere.

Figure 2: Test of the validity of the gravitational wave function ansazt,
using the numerical minimization of the cost function
$\displaystyle\mathcal{F}=\frac{1}{N}\sum_{n=1}^{N}\sum_{n=1}^{N}|\alpha\exp(-\sigma(x_{n},x_{k})/\varepsilon^{2})-C_{nk}^{2}(\Psi_{G})|^{2}$.
The parameters to optimize are vector components of
$\displaystyle|\Psi_{G}\rangle$. Different case studies are considered. In the
first one (panel a), the locations of $\displaystyle N$ points are generated
randomly in a cube of edge length $\displaystyle L/\varepsilon$. In the second
case (panel b), the number of points is fixed, and they are located at the
vertices of a regular polyhedron (a cube or a 3-simplex), the edge length
being varied. In the third case (panel c), the situation is similar to the
first one, but the cube is replaced by a 3-sphere of radius $\displaystyle r$.
In the last case (panel d), points are randomly sampled in a Hayward space-
time [42]. Further details concerning the numerical calculations are given in
Appendix B.

The Hamiltonian of the scalar field (16) can be reconstructed with the
information given by $\displaystyle C_{nk}$. However, the square of
$\displaystyle C_{nk}$ is required, to ensure a positive quantity that can be
assimilated with a probability. This cannot be realized as the mean value of a
linear operator acting on $\displaystyle|\Psi_{G}\rangle$, the quadraticity
can be recovered with a duplicate of the state and a tensor product. We
introduce
$\displaystyle|\tilde{\Psi}_{G}\rangle=|\Psi_{G}\rangle\otimes|\Psi_{G}\rangle$,
and with this state, $\displaystyle C_{nk}^{2}$ is obtained with the relation

$\langle\tilde{\Psi}_{G}|\hat{C}_{nk}\otimes\hat{C}_{nk}|\tilde{\Psi}_{G}\rangle=C_{nk}^{2}.$
(19)

The full quantum Hamiltonian is obtained by replacing $\displaystyle
P_{x_{n}}(x_{k})$ by $\displaystyle\hat{C}_{nk}\otimes\hat{C}_{nk}$ in Eq.
(16). The parameter $\displaystyle\alpha$ is omitted since it can be
factorized with the term $\displaystyle V^{2}/N^{2}$, which is already a
global scaling factor of the Hamiltonian.

The dynamics of the full quantum system are given by the Schrödinger equation

$\frac{d}{dt}|\Psi\rangle=-\frac{iN}{V}\left(\hat{H}+\hat{H}_{G}\right)|\Psi\rangle,$
(20)

where $\displaystyle\hat{H}$ is the Hamiltonian [34, 21] given in Eq. (16),
and the origin of factor $\displaystyle N/V$ is explained in appendix D.
$\displaystyle\hat{H}_{G}$ is a gravitational Hamiltonian. Deriving the
explicit expression of $\displaystyle\hat{H}_{G}$ is a nontrivial key point of
the theory. It is left unresolved in this article. Both $\displaystyle\hat{H}$
and $\displaystyle\hat{H}_{G}$ act on gravitational degrees of freedoms, but
in the low coupling limit, $\displaystyle|\Psi\rangle$ approximately take the
form of a tensor product at any time, i.e.
$\displaystyle|\Psi\rangle\approx|\tilde{\Psi}_{G}\rangle\otimes|\Psi_{S.F.}\rangle$
(Born approximation). This limit is not obtained by $\displaystyle
m\rightarrow 0$, but when
$\displaystyle\|\hat{H}|\Psi\rangle\|\ll\|\hat{H}_{G}|\Psi\rangle\|$. By
taking the trace over gravitational degrees of freedom, one gets:

$\frac{d}{dt}|\Psi_{S.F.}\rangle=-\frac{iN}{V}\left(\langle\hat{H}\rangle_{\tilde{\Psi}_{G}}+\langle\hat{H}_{G}\rangle_{\tilde{\Psi}_{G}}\textrm{1}\\!\\!\\!\\!1\right)|\Psi_{S.F.}\rangle$
(21)

The mean value of $\displaystyle\hat{H}_{G}$ returns only an irrelevant phase
from the point of view of the scalar field, and thus, the Schrödinger equation
can be redefined as

$\frac{d}{dt}|\Psi_{S.F.}\rangle=-\frac{iN}{V}\langle\hat{H}\rangle_{\tilde{\Psi}_{G}}|\Psi_{S.F.}\rangle.$
(22)

It gives us the dynamics of the scalar field in the limit of quantum fields in
curved space-time. This is a mean-field Schrödinger equation similar to the
one encountered in quantum optics [43, 44, 45]. From the point of view of open
quantum systems, this corresponds to a first-order adiabatic elimination [46].
The second-order expansion introduces non-unitary dynamics and an effective
Lindblad equation must replace Eq. (22). The second-order correction goes
beyond the standard quantum field theory in curved space-time and it
incorporates the notion of gravitational decoherence [47, 48]. The rest of
this paper is focused on the first-order approximation and only Eq. (22) is
considered.

##  5 Multi-scales analysis

In the previous sections, a quantum theory for the scalar field and the
gravitational field was developed, such that quantum field theory in curved
space-time is recovered when suitable limits are taken. Several aspects have
been omitted and the goal of this section is to go deeper in the analysis,
using numerical calculations.

###  5.1 Classical and non-classical gravitational states

Figure 3: Panel a) shows the IPR of $\displaystyle|\Psi_{G}\rangle$, given by
$\displaystyle
IPR=\left(\sum_{m=1}^{dim{\mathcal{H}_{G}}}|\Psi_{G,m}|^{4}\right)^{-1}$ as a
function of the edge length $\displaystyle L$ for the optimized quantum states
whose final cost function are given in Fig. 2 b). In this situation, the
optimized quantum state aims to describe a classical geometry with 8 points
located at the vertices of a cube. We recall that when $\displaystyle
L/\varepsilon\leq 1.5$, the state fails to describe the target geometry, but
when $\displaystyle L/\varepsilon>1.5$ the geometry is well described by the
$\displaystyle|\Psi_{G}\rangle$. Panel b) is the same as panel a), but it
shows the mean, max. and min. values of the purity
$\displaystyle\text{Tr}[\hat{\rho}(x_{n})^{2}]$ of each spin reduced states.
Panel c) and d) show the geodesic distance between each point, as given
respectively by the optimized state at $\displaystyle L/\varepsilon=2$ and
$\displaystyle L/\varepsilon=0.1$. The geodesic distance is deduced from Eq.
(18), with value $\displaystyle\alpha=(\sqrt{2\pi}\varepsilon)^{3}$ (fixed for
the numerical optimization), and it reads: $\displaystyle
l_{nk}=\sqrt{-2\varepsilon\log(C_{nk}^{2})}$.

First, a few properties of the gravitational quantum state are investigated by
computing a few relevant quantities. A first one is the inverse participation
ratio (IPR) [49], defined by $\displaystyle
IPR=\left(\sum_{m=1}^{dim{\mathcal{H}_{G}}}|\Psi_{G,m}|^{4}\right)^{-1}$, with
$\displaystyle\Psi_{G,m}$ the vector components of
$\displaystyle|\Psi_{G}\rangle$ in the canonical basis of $\displaystyle
H_{G}$. Roughly, the IPR gives us the number of basis states on which the
state is decomposed with equal weights. This provides us some information on
the complexity of the quantum state, and how many basis states are used to
encode the data of a classical space-time. A second interesting quantity is
the purity of the reduced state associated with each spin.

The purity associated with a point $\displaystyle x_{n}$ is defined by
$\displaystyle\text{Tr}[\hat{\rho}(x_{n})^{2}]$ [44, 45] , with
$\displaystyle\hat{\rho}(x_{n})$ the reduced density matrix of the
corresponding spin. Since the system is a spin-$\displaystyle\tfrac{1}{2}$,
the reduced density matrix can be easily computed using
$\displaystyle\hat{\rho}(x_{n})=(\textrm{1}\\!\\!\\!\\!1_{2}+\sum_{i=1}^{3}\hat{\sigma}_{i}\langle\Psi_{G}|\hat{\sigma}_{i}(x_{n})|\Psi_{G}\rangle)/2$,
with $\displaystyle\sigma_{i}$ the Pauli matrices. The purity is a quantity
describing how the density matrix is close to a pure state. The purity of a
pure state is equal to 1, and the minimum value is 0.5 for the state
$\displaystyle\hat{\rho}=\textrm{1}\\!\\!\\!\\!1_{2}/2$.

The IPR and the averaged purity of each $\displaystyle x_{n}$ are given in
Fig. 3 a) and b), for the quantum states optimized to describe a cube, whose
cost functions $\displaystyle\mathcal{F}$ are given in Fig. 2 b). With these
states, there are the ones that describe a classical geometry, and the ones
that do not. The physical relevance of this second kind of state is less
important than the first one, because they fail in the description of a
classical object, and there is no guarantee that they describe a true physical
situation. However, they give us elements of comparison with the states
describing a classical geometry. We observe that the IPR is usually quite high
(between 30% and 50% of the total dimension of the Hilbert space) and the
purity is low, close to the minimum value of 0.5. The local states are
therefore highly statistical mixtures. The relatively high IPR and low purity
suggest that the states are quite highly entangled. Interestingly, the data of
Fig. 3 a) and b) clearly indicates the change of behavior near $\displaystyle
L/\varepsilon=1.5$. This observation, already commented in the previous
section, requires a deeper analysis to explain its origin.

Having a large-scale entanglement and locally mixed states seem to be
necessary conditions to obtain classical geometries, but we see that there
exist non-classical states that have similar behaviors. Only the correlation
function allows us to distinguish a classical to a non-classical geometry. In
figure Fig. 3 c) and d) are plotted the geodesic distance between the
$\displaystyle x_{n}$ for the optimized state computed with $\displaystyle
L/\varepsilon=2$ and $\displaystyle L/\varepsilon=0.1$. They are respectively
classical and non-classical states. The geodesic distance is computed from
$\displaystyle C_{nk}$ by inverting Eq. (18), and using
$\displaystyle\alpha=(\sqrt{2\pi}\varepsilon)^{3}$ (fixed for the numerical
optimization). It reads $\displaystyle
l_{nk}=\sqrt{-2\varepsilon\log(C_{nk}^{2})}$. In the classical case, the
distances correspond to the ones between the vertices of a cube of edge length
2 (they are equal to $\displaystyle 2$, $\displaystyle 2\sqrt{2}\approx 2.8$
and $\displaystyle 2\sqrt{3}\approx 3.5$). In the nonclassical case, all the
distances are $\displaystyle\approx 1.5$. This is not possible unless the
geometric figure is a hypercube of dimension equal to 8, or if the
interpretation of well-located points in a 3D space is lost. This ensemble
could also be seen as a kind of fuzzy object. This object is divided into
well-identifiable subsystems but do not have well defined positions in a 3D
space. They are identically close to one another, violating the basis axiom of
normed 3D spaces.

###  5.2 Classical limit of the quantum scalar field

Figure 4: Panels a) and b) show the time evolution of the field at different
space positions, computed using the Schrödinger equation and the Hamiltonian
Eq. (16) for the quantum case, or the classical Hamiltonian and Hamilton
equations for the classical case given in Appendix D. The system is defined
with $\displaystyle 5$ points in Minkowski space whose locations are given in
panel d). The initial quantum state condition is a coherent state
$\displaystyle|\Psi_{S.F.}\rangle=|0.1\rangle\otimes|0\rangle\otimes|0\rangle\otimes|0\rangle\otimes|0\rangle$,
and the initial classical state is $\displaystyle\varphi_{n}=(\langle
0.1|\hat{\varphi}|0.1\rangle,0,0,0,0)$ and
$\displaystyle\tilde{\Pi}_{n}=(0,0,0,0,0)$. The mass is fixed such that
$\displaystyle m\varepsilon=2$. Panel c) shows the error made by the semi-
classical approximation. It consists of the 2-norm in
$\displaystyle\mathbb{R}^{N}$ of the difference between the vector
$\displaystyle\langle\hat{\varphi}_{n}\rangle$ and $\displaystyle\varphi_{n}$.
The accumulation of error is due to the small nonclassical effects, but also
to some extent numerical errors due to the truncation of the Hilbert space
(only the first four states of the Fock space are kept, leading to a total
Hilbert space dimension of 1024).

The next step of the multi-scales analysis is to consider a regime where the
gravitational field is in its semi-classical limit such that the dynamics of
the scalar field are governed by the Schrödinger equation Eq. (22). In this
limit, the mean value of the Hamiltonian with respect to the gravitational
state gives us the coefficients $\displaystyle p_{x_{k}}(x_{n})$ in Eq. (16).
We can thus forget the quantum nature of gravity and the coefficients
$\displaystyle p_{x_{k}}(x_{n})$ become inputs of the model. As previously
outlined in the text, the Hamiltonian is in the form of a Bose-Hubbard
Hamiltonian, and it is well defined and already studied in the literature.
However, it is not so easy to recover standard quantum field theory, and even
harder to recover the classical field. Of course, we can follow the reverse
procedure of the one followed in Sec. 3, and come back to the initial point,
but it is quite difficult to visualize the transition from the quantum world
to the semi-classical one. To this end, different semi-classical
approximations can be performed. Here, we are interested in an approach
particularly well suited for numerical computations. First, we consider the
regime where $\displaystyle\langle\hat{\varphi}\rangle\approx\varphi$ and
$\displaystyle\langle\hat{\Pi}\rangle\approx\Pi$ [44], where the dynamics of
$\displaystyle\varphi$ and $\displaystyle\Pi$ are computed using Hamilton’s
equations, with the classical Hamiltonian as a starting point (see appendix
D). This allows us to replace the Schrödinger equation with a set of classical
equations of motion, with a small number of degrees of freedom (compared to
the size of the truncated Hilbert space used in a numerical simulation) [50,
51, 52]. Physically, the classical system corresponds to a set of coupled
oscillators. Equations of motions are non-linear, but their numerical
integration is not particularly difficult. Next, we take the limit when the
number of discretization points tends to infinity, to recover the continuum
limit.

The regime with $\displaystyle\langle\hat{\varphi}\rangle\approx\varphi$ and
$\displaystyle\langle\hat{\Pi}\rangle\approx\Pi$ is obtained when the quantum
state is of the form

$|\Psi_{S.F.}\rangle=\bigotimes_{n=1}^{N}|a_{n}\rangle$ (23)

with $\displaystyle|a_{n}\rangle$ a coherent state defined by [44]

$|a_{n}\rangle=e^{-|a_{n}|^{2}/2}\sum_{k=0}^{\infty}\frac{(a_{n})^{k}}{\sqrt{k!}}|k_{n}\rangle,$
(24)

$\displaystyle|k_{n}\rangle$ being the Fock basis state of the oscillator
$\displaystyle n$. A coherent state has the particularity that
$\displaystyle\langle a_{n}|\hat{a}_{n}|a_{n}\rangle=a_{n}$ and
$\displaystyle\langle a_{n}|\hat{a}_{n}^{\dagger}|a_{n}\rangle=a_{n}^{*}$, and
thus $\displaystyle\langle\hat{\varphi}\rangle=\varphi$ and
$\displaystyle\langle\hat{\Pi}\rangle=\Pi$. The Hamiltonian has the
particularity to preserve quite well the coherent state structure of the
state, and thus, a coherent state remains coherent as a function of time. As a
consequence, if $\displaystyle|\Psi_{S.F.}\rangle$ is a coherent state, it
behaves like a classical ensemble of coupled oscillators. This point is
illustrated in Fig. 4. We see that quantum and classical dynamics are very
close to each other, even after very long times (i.e., after a large number of
oscillations). Physically, this system is quite far from a classical
continuous scalar field because the number of points taken in this example is
quite small. The number of quantum oscillators that can be considered in a
numerical simulation is limited by the size of the Hilbert space.
Nevertheless, a feature of the continuous limit can already be observed: an
initially localized excitation will tend to propagate towards the nearby
points. This is clearly visible in panels a) and b), where the field amplitude
in $\displaystyle x_{1}$ decreases while the amplitude in the other points
increases.

Next, the state of the scalar field is assumed to stay coherent at any time,
such that dynamics can be computed using only the classical equation of motion
(see appendix D). In this regime, the number of $\displaystyle x_{n}$ in the
numerical simulation can be increased considerably, but the Monte Carlo method
requires quite a lot of points to obtain a good convergence, and thus, the
illustrative example is limited to a 1D flat space. In practice, the
expression of the Hamiltonian remains the same as in 3D, but the coefficients
$\displaystyle p_{x_{n}}(x_{k})$ are now defined by. $\displaystyle
p_{x_{n}}(x_{k})=e^{-(x_{n}-x_{k})^{2}/(2\varepsilon^{2})}/(\sqrt{2\pi}\varepsilon)$.
In a flat 1D space, the Klein-Gordon equation is given by [35]

$\frac{\partial^{2}\varphi}{\partial
x^{2}}-\frac{\partial^{2}\varphi}{\partial t^{2}}-m^{2}\varphi=0.$ (25)

It can be derived readily from the Lagrangian Eq. (5), with $\displaystyle
g^{\mu\nu}=\eta^{\mu\nu}$ (the Minkowski metric), and using Euler Lagrange
equations with respect to the field variable $\displaystyle\varphi$. Solutions
of Eq. (25) are given by [35]

$\varphi(x,t)=\int\frac{dk}{4\pi\omega(k)}\left(a(k)e^{i(-\omega(k)t+kx)}+C.C.\right),$
(26)

with $\displaystyle\omega(k)=\sqrt{k^{2}+m^{2}}$. In order to satisfy the
condition that $\displaystyle\varphi$ vanish on the boundary of the domain
(remember that this is a necessary condition to derive the quantum Hamiltonian
in the form of Eq. (16)), the density of modes is chosen to be

$a(k)=\frac{\omega(k)}{\kappa}\exp\left(-\frac{k^{2}}{2\kappa^{2}}\right).$
(27)

Plugging Eq. (27) into Eq. (26) leads us to

$\varphi(x,t)=\int\frac{dk}{2\pi\kappa}\exp\left(-\frac{k^{2}}{2\kappa^{2}}\right)\cos(kx-\omega(k)t).$
(28)

This integral can be easily evaluated numerically with very good precision.
Even if a short numerical computation is necessary for its evaluation, Eq.
(28) is called below "the analytic solution" of the Klein-Gordon equation. In
the limit when $\displaystyle m\rightarrow\infty$, a simple approximated
expression can be derived, $\displaystyle\varphi(x,t)\approx
e^{-x^{2}/2\kappa^{2}}\cos(mt)/\sqrt{2\pi}$. The shape of the wave packet is a
Gaussian oscillating at a frequency $\displaystyle m$.

Figure 5: Panel a) shows the analytic solution of the Klein Gordon equation,
given by Eq. (28) with $\displaystyle m\varepsilon=0.8$ and
$\displaystyle\kappa\varepsilon=0.15$. Panel b) shows the numerical solution,
computed with the same parameters and an initial state given by the analytic
solution at $\displaystyle t=0$. In the numerical solution, time must be
rescaled due to the global factorization in the Hamiltonian Eq. (16). Panel c)
gives the deviation (also called the error) of the numerical solution from the
Klein-Gordon equation. It is given by the sup norm of
$\displaystyle(\nabla^{2}-m^{2})\varphi=\partial_{x}^{2}\varphi-\partial_{t}^{2}\varphi-m^{2}\varphi$,
the norm being computed over the variable $\displaystyle x$.

In Fig. 5, the analytic solution of the Klein-Gordon equation is compared with
the discrete model, defined by the classical version of Eq. (16). The equation
of motion of the second model being entirely integrated numerically, the
resulting solution is called below the "numerical solution". The two dynamics
are very close to each other. In both cases, two observations can be made, the
field oscillates at the frequency $\displaystyle m$ and there is a time
deformation of the wave packet because $\displaystyle m\varepsilon$ is of the
order of unity. The validity of the numerical solution can be characterized
more precisely by computing its deviation from the Klein-Gordon equation. The
number of points $\displaystyle\varphi(x,t)$ in the numerical solution being
sufficiently high, it is possible to compute an analytic interpolation of the
data points. The interpolated function is then inserted in the Klein-Gordon
equation, and the error at a given time $\displaystyle t$ is defined by the
sup norm, computed over the variable $\displaystyle x$. The error of the
numerical solution is given in Fig. 5 c). It is divided by the sup norm of
$\displaystyle\varphi(t=0)$ to fix an irrelevant scaling amplitude. The error
remains small and bounded ($\displaystyle\lesssim$ 2 percent). Better accuracy
of the method can be obtained by increasing the number of points, increasing
the region of integration, and increasing the size of the wave packet (the
size of the wave packet must be large to ensure a good convergence of the
approximation scheme of the Laplacian, as given in Eq. (13)).

To finalize the illustration of the method, a short example of propagation in
a curved space-time is provided in Fig. 6. The chosen one is a Hayward space-
time [42] where space-slices are restricted to a single dimension. In Lemaitre
type coordinates, the line element is given by:

$ds^{2}=-dt^{2}+\frac{r_{s}r^{2}}{r^{3}+l^{2}r_{s}}d\rho^{2}$ (29)

with $\displaystyle t$ the time coordinate, and $\displaystyle\rho$ a space
coordinate. $\displaystyle r_{s}$ correspond to the Schwarzschild radius, and
$\displaystyle l$ is a length scale of the central body. The function
$\displaystyle r(\rho,t)$ is a non-analytic, and defined through the
differential equation $\displaystyle d\rho-
dt=dr\left(\tfrac{r_{s}r^{2}}{r^{3}+l^{2}r_{s}}\right)^{-1/2}$. It can be
easily evaluated numerically or approximated analytically using a Taylor
expansion. Since the coordinates are related by a differential equation, an
arbitrary choice of origin is required. Here, $\displaystyle\rho-t=0$ is
chosen to correspond with $\displaystyle r=2^{1/3}l^{2/3}r_{s}^{1/3}$. In this
case study, the space-time has two horizons [42]. The wave packet is initially
located near the inner horizon. Due to the non-zero energy density at the
origin of this space-time, the wave packet propagates towards the horizon
(free fall towards the center of the body, located at $\displaystyle r=0$),
and a very strong distorsion happens near this latter. Here, the goal is just
to illustrate that the method allows us modelize the non-trivial properties of
a curved space-time222The second order expansion of $\displaystyle\sigma$ used
for its estimation allows us to show the existance of curvature effects, but
it may also leads to numerical errors that could be corrected with more
accurate estimations of $\displaystyle\sigma$.. No further details are
investigated, this could be the topic of another independent paper.

The very good results given by the numerical solution show the relevance of
the discretization method. It is, therefore, an interesting for quantum
gravity purposes, but it can also be an alternative discretization scheme in
quantum field theory in curved space-time for which most of the current
techniques are based on a simplicial discretization of the space-time manifold
[23, 24, 25, 26, 27]. The simplicial discretization requires a nontrivial
discretization scheme of the differential operators which may be difficult to
compute in large dimensions. With the Monte-Carlo based approach the
nontrivial task is to compute the geodesic length, but once this task is
achieved, potentially with some approximations, then, the implementation of
the scheme is straightforward.

Figure 6: Propagation of a scalar field near the inner horizon of a Hayward
space-time. The dashed line corresponds to this horizon, and the solid line
corresponds to the radial coordinate $\displaystyle r=(2rs)^{1/3}l^{2/3}$. The
parameters used in the simulation are $\displaystyle l/\varepsilon=10$,
$\displaystyle r_{s}/\varepsilon=60$, $\displaystyle m\varepsilon=0.8$, and
$\displaystyle\kappa\varepsilon=0.4$.

##  6 Conclusion

In this paper, several ideas of quantum gravity have been explored, such as
the possibility to encode space-time geometric distances in the correlation of
a spin-like system, and to use a discretization scheme based on Monte-Carlo
principles instead of a lattice-based discretization. These two ingredients
offer interesting advantages compared to other common quantum gravity theories
which are based on elementary areas (instead of distances) and/or simplicial
lattices for the discretization scheme. The use of distances leads to a
simpler connection with general relativity and our everyday world, for which
elementary quantities are distances. The Monte Carlo method simplifies the
discretization scheme of differential operators, and it avoids conceptual
problems induced by the arbitrary choice of a lattice.

To construct a physically relevant theory, the starting point was to quantize
a scalar field in curved space-time. Working with such a test system has led
to a simple link between general ideas and a well-defined mathematical model.
In this model, the gravitational field is described by a quantum system which
is similar to an ensemble of spin-$\displaystyle\tfrac{1}{2}$ particles. At
each spin corresponds a point of space. The geodesic length between two of
these points is given by the quantum state through the relation $\displaystyle
C_{nk}^{2}=\langle\hat{\sigma}_{+}(x_{n})\hat{\sigma}_{-}(x_{k})+\hat{\sigma}_{+}(x_{k})\rangle_{\Psi_{G}}^{2}\propto
p_{x_{n}}(x_{k})$. This relation is a postulate of the theory, and its
physical relevance has been investigated with numerical tests. Quantum states
that describe various classical geometries have been found with a very high
numerical precision. This is a first positive result supporting this very
strong postulate. The numerical simulations have also shown the existence of
two different regimes associated with two different length scales. Bellow a
specific threshold, the classical description does not hold anymore and one
obtains something that can be called "a quantum space-time". This is a
nontrivial consequence of the model, which remotely resembles the quantization
of areas and volumes in loop quantum gravity. In loop quantum gravity, the
minimal scale is given by the lowest eigenvalues of area and volume operators
[21, 22, 16]. This can be interpreted as the existence of minimal building
blocks that can be used to construct the space-time. However, in the theory of
this paper, the interpretation is slightly different. The classical space-time
is observed only if the distance between two points is sufficiently far from
each other. If the density of points becomes too short, there is a threshold
where a classical picture is not possible anymore and the resulting object
cannot be interpreted as a space-like sub-manifold. This effect may have
interesting consequences in black hole physics and cosmology. The relevance of
the model has also been investigated from another point of view, which has
consisted of recovering, by means of successive approximation, the classical
theory of a scalar field. Each step of the procedure has been illustrated by
numerical examples. These examples show clearly the different regimes of the
theory, but they also demonstrate its computational efficiency.

In these few pages, many issues have been glossed over. Here is a (non-
exhaustive) list of the points that shall be investigated:

  * $\displaystyle\bullet$

One of the most important ones concerns the explicit construction of the
gravitational Hamiltonian $\displaystyle H_{G}$, which is compatible with the
Einstein-field equation (at least in some limit cases, compatible with current
experimental observations). Another issue concerns the choice of the
probability density $\displaystyle p_{x_{n}}(x_{k})$, which has been assumed
to be a normal law. In fact, the form of the density may be fixed by the
theory itself, by $\displaystyle\hat{H}_{G}$ and its eigenvectors. The use of
a normal Law is interesting because at small distances we have $\displaystyle
C_{nk}^{2}\propto 1-q_{ab}dx^{a}dx^{b}/(2\varepsilon^{2})$, with
$\displaystyle dx$ a vector connecting $\displaystyle x_{n}$ to $\displaystyle
x_{k}$. The relation with the metric tensor is therefore very simple and quite
natural, but many different distributions can have locally this form and very
different behaviors at large scales.

  * $\displaystyle\bullet$

In this paper, it was sufficient to provide a quantum mechanical description
of the 3-metric $\displaystyle q_{ab}$. However, this is not sufficient to
describe the entire space-time manifold $\displaystyle\mathcal{M}$, the
extrinsic curvature $\displaystyle k_{ab}$ is also necessary, and it must be
encoded somewhere in $\displaystyle|\tilde{\Psi}_{G}\rangle$.

  * $\displaystyle\bullet$

A similar study with a Dirac field instead of a scalar field could also be
interesting to test the coupling operator $\displaystyle\hat{C}_{nk}$. This
operator (see Eq. (17)) has been chosen because it provides interesting
numerical results, but this choice is not unique. For a Dirac field, it could
be more interesting to choose an operator of the form
$\displaystyle\hat{C}_{nk}=\sum_{i,j\in\\{x,y,z\\}}w_{ij}\hat{\sigma}_{i}(x_{n})\hat{\sigma}(x_{k})$,
with $\displaystyle w_{ij}\in\mathbb{R}$, with
$\displaystyle\hat{\sigma}_{i}(x_{n})$ the Pauli matrices.

  * $\displaystyle\bullet$

Another point to clarify concerns the possibility of interpreting
$\displaystyle\hat{\sigma}_{i}(x_{n})$ as a quantized 3D reference frame. If
such a link can be consistently made with the ideas of the paper, this may
provide a very interesting connection with current theories of quantum
reference frame.

  * $\displaystyle\bullet$

Finally, a last point to investigate with deeper details, together with the
construction of $\displaystyle H_{G}$, is the physical consequence of the use
of $\displaystyle|\tilde{\Psi}_{G}\rangle=|\Psi_{G}\rangle^{\otimes 2}$. This
condition has a pure mathematical origin, it could be relaxed, such that
$\displaystyle|\tilde{\Psi}_{G}\rangle=|\Psi_{1}\rangle\otimes|\Psi_{2}\rangle$.
The two states would correspond only in the classical limit.

## Acknowledgement

The author acknowledges David Viennot for useful discussions on emergent
gravity.

## Data Availability

Numerical codes and data can be provided upon reasonable request.

##  Appendix A Normal law in a curved space-time

######  Theorem 1.

Let Synge’s world function be expressed in terms of Riemann normal coordinates
[31, 53]:

$\displaystyle\sigma(x_{n},x_{k})=\frac{1}{2}\left(\delta_{ab}y^{a}y^{b}-\frac{\lambda^{2}}{3}R_{acbd}y^{a}y^{b}y^{c}y^{d}+O(\lambda^{3})\right)$

with the coordinate system defined such that $\displaystyle x_{n}$ is located
at the origin, and $\displaystyle x_{k}$ located at $\displaystyle y^{a}$.
$\displaystyle\lambda$ is a scale parameter that must be sufficiently small,
but also sufficiently large in front of $\displaystyle\varepsilon$, and
$\displaystyle R_{abcd}$ is Riemann tensor evaluated at $\displaystyle y=0$.
Then, we have,

$\displaystyle\int d^{3}y\sqrt{q}\leavevmode\nobreak\
\frac{\mathcal{N}(R)}{(\sqrt{2\pi}\varepsilon)^{3}}e^{-\sigma/\varepsilon^{2}}=1+O(\lambda^{3})$

with
$\displaystyle\mathcal{N}(R)=1+R\varepsilon^{2}\lambda^{2}/6+O(\lambda^{3})$,
$\displaystyle R$ being the scalar curvature, and

$\displaystyle\int d^{3}y\sqrt{q}\leavevmode\nobreak\
y^{a}y^{b}\frac{\mathcal{N}(R)}{(\sqrt{2\pi}\varepsilon)^{3}}e^{-\sigma/\varepsilon^{2}}=\varepsilon^{2}+O(\varepsilon^{3})+O(\lambda^{3}).$

Proof  The proof is based in both cases on a series expansion in powers of
$\displaystyle\lambda$, followed by the integration. To perform the
calculation $\displaystyle\sqrt{q}$ must be also Taylor expanded. Following
anew [53], one get

$\sqrt{q}=1-\frac{\lambda^{2}}{6}R_{ab}y^{a}y^{b}+O(\lambda^{3}),$ (30)

with $\displaystyle R_{ab}$ the Ricci tensor. Plugging all together the
expansions of $\displaystyle\mathcal{N}(R),\sigma$ and
$\displaystyle\sqrt{q}$, the integral becomes:

$\begin{split}&\int d^{3}y\sqrt{q}\leavevmode\nobreak\
\frac{\mathcal{N}(R)}{(\sqrt{2\pi}\varepsilon)^{3}}e^{-\sigma/\varepsilon^{2}}\\\
&=\int
d^{3}y\frac{e^{-\delta_{ab}y^{a}y^{b}/\varepsilon^{2}}}{(\sqrt{2\pi}\varepsilon)^{3}}\left(1+\frac{\lambda^{2}}{6}(R\varepsilon^{2}-R_{ab}y^{a}y^{b}\right.\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\
\left.+\frac{R_{abcd}}{\varepsilon^{2}}y^{a}y^{b}y^{c}y^{d})\right)+O(\lambda^{3})\end{split}$

The integral has non zero contributions only if the powers in $\displaystyle
y^{a}$ are all even (otherwise positive and negative positions cancel each
other). This leads us to the following replacement rules:

$\displaystyle\displaystyle\int
d^{3}y\frac{e^{-\delta_{ab}y^{a}y^{b}/\varepsilon^{2}}}{(\sqrt{2\pi}\varepsilon)^{3}}y^{a}y^{b}$
$\displaystyle\displaystyle=\varepsilon^{2}\delta^{ab}$ (31)
$\displaystyle\displaystyle\int
d^{3}y\frac{e^{-\delta_{ab}y^{a}y^{b}/\varepsilon^{2}}}{(\sqrt{2\pi}\varepsilon)^{3}}(y^{a})^{2}(y^{b})^{2}$
$\displaystyle\displaystyle=\varepsilon^{4}(3-2\delta^{ab}).$ (32)

Moreover, the Ricci tensor is symmetric and the Riemann tensor is skew-
symmetric on some of its indices. As a consequence, one get $\displaystyle
R_{ab}y^{a}y^{b}\rightarrow R_{ab}\delta^{ab}=R$ (the scalar curvature being
evaluated at $\displaystyle y^{a}=0$), and $\displaystyle
R_{abcd}y^{a}y^{b}y^{c}y^{d}\rightarrow 0$. As a result, the second order term
in $\displaystyle\lambda$ in (A) is canceled, and the result is $\displaystyle
1+O(\lambda^{3})$.

For the second relation of the theorem, the starting point is the same, but
now, the minimum number of powers in $\displaystyle y^{a}$ is at least one.
Once again, the term with the Riemann tensor is canceled, and the only
surviving terms are the ones with $\displaystyle y^{a}y^{b}$, $\displaystyle
R\varepsilon^{2}\lambda^{2}y^{a}y^{b}/6$ , and $\displaystyle
R_{cb}\lambda^{2}y^{a}y^{b}y^{c}y^{d}/6$. The result of the integral, is thus
$\displaystyle\varepsilon^{2}-\varepsilon^{4}\lambda^{2}(R+2R^{ab})/6$. Then,
this leads to a final expression of
$\displaystyle\varepsilon^{2}+O(\varepsilon^{3})+O(\lambda^{3})$.

##  Appendix B Details on the numerical optimizations

In Sec. 4, the assumptions at the core of the quantum description of space-
time is tested numerically by searching numerically quantum states that can
describe a classical space-time.

As earlier commented in the main text, the goal is to determine the vector
component of a quantum state that minimize a cost function
$\displaystyle\mathcal{F}$. Several cases study are considered, and results
are plotted in Fig. 2. In the first one (panel a), the locations of
$\displaystyle N$ points are generated randomly in a cube of edge length
$\displaystyle L/\varepsilon$. In the second case (panel b), the number of
points is fixed, and they are located at the vertices of a regular polyhedron
(a cube or a 3-simplex), the edge length being varied. In the third case
(panel c), the situation is similar to the first one, but the cube is replaced
by a 3-sphere of radius $\displaystyle r$. In the last case (panel d), points
a randomly sampled in a Hayward space-time [42]. With these examples, flat and
two kind of curved spaces are investigated. The first curved space of interest
is a 3-sphere, for which it is easy to compute geodesic lengths. The second
curved space has, contrary to the 3-sphere, a position dependent scalar
curvature. Hence, the normalization coefficient of the probability
distribution (see Eq. (4)) depends on the position.

The Hayward333A Hayward space-time describes a Black hole without singularity.
space-time is defined by the line element (written in Lemaitre-type
coordinates)

$ds^{2}=-dt^{2}+\frac{r_{s}r^{2}}{r^{3}+l^{2}r_{s}}d\rho^{2}+r^{2}d\Omega^{2},$
(33)

with $\displaystyle t$ the time coordinate, $\displaystyle\rho$ the radial
coordinate, and $\displaystyle\Omega$ is a solid angle coordinate,
$\displaystyle r_{s}$ correspond to the Schwarzschild radius, and
$\displaystyle l$ is a length scale of the central body. $\displaystyle
r(\rho,t)$ is a non-analytic function of $\displaystyle\rho$ and
$\displaystyle t$, defined by $\displaystyle d\rho-
dt=dr\left(\tfrac{r_{s}r^{2}}{r^{3}+l^{2}r_{s}}\right)^{-1/2}$. It can be
easily evaluated numerically or approximated analytically using a Taylor
expansion. Since the coordinates are related by a differential equation, an
arbitrary choice of origin is required. Here, $\displaystyle\rho-t=0$ is
chosen to correspond with $\displaystyle r=2^{1/3}l^{2/3}r_{s}^{1/3}$. In this
paper only a limited portion of the space-time is considered. In the case of
Fig. 2, $\displaystyle t=0$ and $\displaystyle\rho/\varepsilon\in[-1,2]$,
$\displaystyle l/\varepsilon=0.5$, and $\displaystyle r_{s}/\varepsilon=3$.

To perform the calculations, several key points must be considered. First of
all, the minimization of $\displaystyle\mathcal{F}$ is performed with a two-
step procedure. The first step is a global minimization with the algorithm
JAYA [54]. With this algorithm, convergence towards a minimum may not be
achieved (mostly when the number of spins is large), but the resulting quantum
state is used as an input for a second optimization using a gradient descent
algorithm [55].

The second point is that $\displaystyle C_{nk}^{2}\propto p_{x_{n}}(x_{k})$,
and the proportionality coefficient is unknown. The value of this coefficient,
namely $\displaystyle\alpha$, is fixed by hand in the numerical optimization.
In principle it can be optimized simultaneously with the quantum state, but at
the price of an important computational cost. A third important point concerns
the computation of geodesic lengths. This can be done exactly in the simplest
cases, but otherwise they are only estimated approximately.

  * $\displaystyle\bullet$

In the first two cases $\displaystyle\alpha$ is simply fixed to
$\displaystyle\alpha=1$, and the geodesic length is given by the Euclidean
distance. We thus have $\displaystyle\sigma=\|x_{n}-x_{k}\|^{2}/2$.

  * $\displaystyle\bullet$

In the case of the 3-sphere, the space is curved,
$\displaystyle\alpha=(1+R\varepsilon^{2}\lambda^{2}/6)(\sqrt{2\pi}\varepsilon)^{-3}$,
with $\displaystyle\lambda=3r/5$ (fixed by hand, such that
$\displaystyle\varepsilon<\lambda<r$). Recall that $\displaystyle r$ is the
radius of the sphere, and the scalar curvature is given by $\displaystyle
R=6/r^{2}$. The geodesic length is given by the arc length, so that
$\displaystyle\sigma=r^{2}\arccos(u_{n}.u_{k}/r^{2})^{2}/2$, with
$\displaystyle u_{n}$ the point position in $\displaystyle\mathbb{R}^{4}$.

  * $\displaystyle\bullet$

In the case of the Hayward space,
$\displaystyle\alpha=(\sqrt{2\pi}\varepsilon)^{-3}$. The geodesic length is
estimated using a second order expansion of $\displaystyle\sigma$ [31, 53]

##  Appendix C Covariant Monte-Carlo discretization

In this appendix, the ideas introduced in Sec. 3.2 are reused to define a
covariant quantum theory of the scalar field. In such a framework, the action
is discretized, and it can be inserted in a second step into a path-integral.
The precise definition of the path integral is not discussed here, the
following few lines focus on the discretization of the action. The usual
action for a scalar field is given by:

$S=\int
d^{4}x\frac{\sqrt{-g}}{2}\left[g^{\mu\nu}(\partial_{\mu}\varphi)(\partial_{\nu}\varphi)+m^{2})\varphi^{2}\right]$
(34)

It can be discretized with a Monte-Carlo approach, as given in Sec. 2, but
using the space-time metric $\displaystyle g$ instead of the metric
$\displaystyle q$ of a 3D hypersurface. The nontrivial point would be to
estimate the term $\displaystyle
g^{\mu\nu}(\partial_{\mu}\varphi)(\partial_{\nu}\varphi)$.

One can notice that:

$\begin{split}\mathbb{E}[\varphi^{2}]_{x}&=\int
d^{4}y\sqrt{-g}\leavevmode\nobreak\ p_{x}(y)\varphi(x)^{2}\\\
&\approx\varphi(x)^{2}+\frac{\varepsilon^{2}}{2}\varphi(x)\nabla_{\mu}\nabla^{\mu}\varphi(x)+\varepsilon^{2}(\partial^{\mu}\varphi)(\partial_{\mu}\varphi)\end{split}$
(35)

The computation steps are basically the same as in Eq. (12), but using the
Taylor expansion of $\displaystyle\varphi(x)^{2}$ over space-time. Using the
relation of Eq. (13), generalized to 4D, it is possible to rewrite Eq. (35)
into

$\mathbb{E}[\varphi^{2}]_{x}\approx\varepsilon^{2}\left(\frac{1}{2}\varphi(x)\mathbb{E}[\varphi]_{x}+(\partial^{\mu}\varphi)(\partial_{\mu}\varphi)\right),$
(36)

and thus,

$(\partial^{\mu}\varphi)(\partial_{\mu}\varphi)\approx\frac{1}{\varepsilon^{2}}\mathbb{E}[\varphi^{2}]_{x}-\frac{1}{2}\varphi(x)\mathbb{E}[\varphi]_{x}.$
(37)

As a consequence, with the Monte-Carlo approximation, the action is
approximately

$S\approx\frac{V}{2N}\sum_{n=1}^{N}m^{2}\varphi_{n}^{2}+\frac{V}{N}\sum_{k=1}^{N}p_{x_{n}}(x_{k})\left(\frac{1}{\varepsilon^{2}}\varphi_{k}^{2}-\frac{1}{2}\varphi_{n}\varphi_{k}\right)$
(38)

One can also factorize by $\displaystyle V/N$, as in Eq. (16), and the final
result is

$\begin{split}S\approx&\frac{V^{2}}{2N^{2}}\sum_{n=1}^{N}\left(\sum_{k=1}^{N}p_{x_{n}}(x_{k})\right)m^{2}\varphi_{n}^{2}\\\
&+\sum_{k=1}^{N}p_{x_{n}}(x_{k})\left(\frac{1}{\varepsilon^{2}}\varphi_{k}^{2}-\frac{1}{2}\varphi_{n}\varphi_{k}\right).\end{split}$
(39)

Once again, this formula is evaluated in 4D, which means that $\displaystyle
p_{x_{n}}(x_{k})$ is a function of the geodesic length in
$\displaystyle\mathcal{M}$, not in $\displaystyle\Sigma$.

##  Appendix D Hamiltonian and equation of motion of discretized systems

In the main text of this paper, repeated references to the classical
Hamiltonian and the classical equation of motion are made, without detailed
explanations on the subject. A few details on the subject are reported in this
appendix.

The classical Hamiltonian and the classical equation of motions are already
given in Sec. 3, in the case of the continuous field. However, the formula
must be slightly adapted in the discretized versions.

In the continuous case, Hamilton’s equations are given by functional
derivatives [56]. In the case of $\displaystyle\delta$-peacked variation of a
field, a functional $\displaystyle I=\int d\mu(x)f(\varphi(x))$ has for
variation:

$\frac{\delta I}{\delta\varphi(y)}=\int d\mu(x)\frac{\partial
f}{\partial\varphi}\delta(x-y)=\frac{\partial f}{\partial\varphi}(y).$ (40)

In the Monte-Carlo approach, the Lebesgue measure 444for simplicity, a flat
space is assumed $\displaystyle\mu$ is replaced by a sum of Dirac
distribution, i.e.

$I\approx I^{\prime}=\frac{V}{N}\sum_{n}\int
d\mu(x)\delta(x-x_{n})f(\varphi(x)).$ (41)

The functional derivative is, therefore, 555to get this result, a product of
Dirac distributions must be carefully computed using the Colombeau algebra
[57].

$\frac{\delta I^{\prime}}{\delta\varphi(y)}=\frac{V}{N}\sum_{n}\frac{\partial
f}{\partial\varphi}(x_{n})\delta(x_{n}-y).$ (42)

This is a distribution. To identify the variations of $\displaystyle
I^{\prime}$ to the variations of $\displaystyle I$,
$\displaystyle\tfrac{\delta I^{\prime}}{\delta\varphi(y)}$ must be
regularized. The regularization can be performed by different means. A simple
one is to consider:

$\frac{\delta
I}{\delta\varphi(y)}=\frac{N}{V}\int_{\mathcal{B}_{y}}d\mu(x)\frac{\delta
I^{\prime}}{\delta\varphi(y)}$ (43)

with $\displaystyle\mathcal{B}_{y}$ a ball centered on $\displaystyle y$ with
a radius sufficiently small, so that a single Dirac distribution is
integrated. The two functional derivatives are thus proportional, with a
factor $\displaystyle N/V$.

As a consequence, the proportionality coefficient is recovered in all the
formulas adapted from the continuous case. Then, Hamilton’s equations for the
discretized system are given by:

$\displaystyle\displaystyle\frac{d\varphi_{n}}{dt}=\frac{N}{V}\frac{\partial
H}{\partial\Pi_{n}}$ (44)
$\displaystyle\displaystyle\frac{d\Pi_{n}}{dt}=-\frac{N}{V}\frac{\partial
H}{\partial\varphi_{n}}.$ (45)

with $\displaystyle H$ the discretized Hamiltonian, and the Schrödinger
equation is

$\frac{d}{dt}|\psi\rangle=-i\frac{N}{V}\hat{H}|\psi\rangle.$ (46)

Now, a few details on the classical version of Eq. (16) are given. The
classical Hamiltonian is obtained by reversing the canonical quantization
procedure, which means that $\displaystyle\hat{a}_{n}$ and
$\displaystyle\hat{a}_{n}^{\dagger}$ are replaced by $\displaystyle\mathbb{C}$
numbers $\displaystyle a_{n}$ and $\displaystyle a_{n}^{\dagger}$. The
corresponding Hamiltonian $\displaystyle H$ is therefore a function on the
complex phase space $\displaystyle\mathbb{C}^{N}$. Equations of motion are
then obtained by complexified Hamilton equations $\displaystyle
id_{t}a_{n}=\tfrac{N}{V}\partial_{a_{n}^{\dagger}}H$. We can also come back to
the real representation of the system, using the change of variable introduced
in Eq. (8) and (9). After the change of variable, the classical Hamiltonian
reads:

$\begin{split}H=&\frac{V^{2}}{2N^{2}}\sum_{n=1}^{N}\left(\sum_{k=1}^{n}p_{x_{n}}(x_{k})\right)\left(\Pi_{n}^{2}+\left[m^{2}+\frac{1}{\varepsilon^{2}}\right]\varphi_{n}^{2}\right)\\\
&-\frac{1}{\varepsilon^{2}}\left(\varphi_{n}\sum_{k=1}^{N}p_{x_{n}}(x_{k})\varphi_{k}\right).\end{split}$
(47)

## References

  * [1] Daniela Klammer and Harold Steinacker.  “Fermions and emergent noncommutative gravity”.  Journal of High Energy Physics 2008, 074 (2008).
  * [2] Harold Steinacker.  “Emergent geometry and gravity from matrix models: an introduction”.  Classical and Quantum Gravity 27, 133001 (2010).
  * [3] David Viennot.  “Emergent gravity and d-brane adiabatic dynamics: emergent lorentz connection”.  Classical and Quantum Gravity 38, 245004 (2021).
  * [4] David Viennot.  “Fuzzy schwarzschild (2 + 1)-spacetime”.  Journal of Mathematical Physics 63, 082302 (2022).
  * [5] Tomasz Konopka, Fotini Markopoulou, and Simone Severini.  “Quantum graphity: A model of emergent locality”.  Phys. Rev. D 77, 104029 (2008).
  * [6] Francesco Caravelli and Fotini Markopoulou.  “Properties of quantum graphity at low temperature”.  Phys. Rev. D 84, 024002 (2011).
  * [7] James Q. Quach, Chun-Hsu Su, Andrew M. Martin, and Andrew D. Greentree.  “Domain structures in quantum graphity”.  Phys. Rev. D 86, 044001 (2012).
  * [8] ChunJun Cao, Sean M. Carroll, and Spyridon Michalakis.  “Space from Hilbert space: Recovering geometry from bulk entanglement”.  Phys. Rev. D 95, 024031 (2017).
  * [9] Ning Bao, Sean M. Carroll, and Ashmeet Singh.  “The Hilbert space of quantum gravity is locally finite-dimensional”.  Int. J. Mod. Phys. D 26, 1743013 (2017).
  * [10] ChunJun Cao and Sean M. Carroll.  “Bulk entanglement gravity without a boundary: Towards finding einstein’s equation in hilbert space”.  Phys. Rev. D 97, 086003 (2018).
  * [11] Thinakkal Padmanabhan.  “Emergent gravity paradigm: recent progress”.  Modern Physics Letters A 30, 1540007 (2015).
  * [12] Andrea Mondino and Stefan Suhr.  “An optimal transport formulation of the einstein equations of general relativity”.  Journal of the European Mathematical Society 25, 933–994 (2022).
  * [13] Jonathan Gorard.  “Some quantum mechanical properties of the wolfram model”.  Complex Syst.29 (2020).
  * [14] Jonathan Gorard.  “Some relativistic and gravitational properties of the wolfram model”.  Complex Systems 29, 599–654 (2020).
  * [15] Lisa Glaser and Sebastian Steinhaus.  “Quantum gravity on the computer: Impressions of a workshop”.  Universe5 (2019).
  * [16] Abhay Ashtekar and Eugenio Bianchi.  “A short review of loop quantum gravity”.  Reports on Progress in Physics 84, 042001 (2021).
  * [17] Alexander Y Yosifov.  “Aspects of semiclassical black holes: Development and open problems”.  Advances in High Energy Physics 2021, 1–13 (2021).
  * [18] Suddhasattwa Brahma, Robert Brandenberger, and Samuel Laliberte.  “Bfss matrix model cosmology: Progress and challenges”.  Physics 5, 1–10 (2023).
  * [19] Michael M. Wolf, Frank Verstraete, Matthew B. Hastings, and J. Ignacio Cirac.  “Area laws in quantum systems: Mutual information and correlations”.  Phys. Rev. Lett. 100, 070502 (2008).
  * [20] Jan Ambjorn, Andrzej Görlich, Jerzy Jurkiewicz, and Renate Loll.  “Causal dynamical triangulations and the search for a theory of quantum gravity”.  Int. J. Mod. Phys. D 22, 1330019 (2013).
  * [21] Carlo Rovelli.  “Quantum Gravity”.  Cambridge Monographs on Mathematical Physics. Cambridge University Press. (2004).
  * [22] Carlo Rovelli and Francesca Vidotto.  “Covariant loop quantum gravity: An elementary introduction to quantum gravity and spinfoam theory”.  Cambridge University Press. (2014).
  * [23] N. H. Christ, R. Friedberg, and T. D. Lee.  “Weights of links and plaquettes in a random lattice”.  Nuclear Physics B 210, 337–346 (1982).
  * [24] Hai-cang Ren.  “Matter fields in lattice gravity”.  Nuclear Physics B 301, 661–684 (1988).
  * [25] F. L. Teixeira.  “Random Lattice Gauge Theories and Differential Forms” (2013).  url: http://arxiv.org/abs/1304.3485.
  * [26] Arata Yamamoto.  “Lattice QCD in curved spacetimes”.  Phys. Rev. D 90, 054510 (2014).
  * [27] Richard C. Brower, George Fleming, Andrew Gasbarro, Timothy Raben, Chung-I Tan, and Evan Weinberg.  “Quantum Finite Elements for Lattice Field Theory”.  PoS LATTICE2015, 296 (2016).  arXiv:1601.01367.
  * [28] Russel E. Caflisch.  “Monte carlo and quasi-monte carlo methods”.  Acta Numerica 7, 1–49 (1998).
  * [29] Emilio Zappa, Miranda Holmes-Cerfon, and Jonathan Goodman.  “Monte carlo on manifolds: sampling densities and integrating functions”.  Communications on Pure and Applied Mathematics 71, 2609–2647 (2018).
  * [30] Stefano De Marchi and Giacomo Elefante.  “Quasi-monte carlo integration on manifolds with mapped low-discrepancy points and greedy minimal riesz s-energy points”.  Applied Numerical Mathematics 127, 110–124 (2018).
  * [31] Eric Poisson, Adam Pound, and Ian Vega.  “The motion of point particles in curved spacetime”.  Living Reviews in Relativity 14, 7 (2011).
  * [32] Igor Khavkine and Valter Moretti.  “Algebraic QFT in curved spacetime and quasifree Hadamard states: an introduction”.  In Advances in algebraic quantum field theory.  Pages 191–251.  Springer (2015).
  * [33] Christian Gérard.  “An introduction to quantum field theory on curved spacetimes”.  Page 171–218.  London Mathematical Society Lecture Note Series. Cambridge University Press. (2018).
  * [34] Charles W. Misner, Kip S. Thorne, and John Archibald Wheeler.  “Gravitation”.  W. H. Freeman. (1973).  1st edition.
  * [35] Claude Itzykson and Jean-Bernard Zuber.  “Quantum field theory”.  Courier Corporation. (2012).
  * [36] Matthew P. A. Fisher, Peter B. Weichman, G. Grinstein, and Daniel S. Fisher.  “Boson localization and the superfluid-insulator transition”.  Phys. Rev. B 40, 546–570 (1989).
  * [37] Omjyoti Dutta, Mariusz Gajda, Philipp Hauke, Maciej Lewenstein, Dirk-Sören Lühmann, Boris A Malomed, Tomasz Sowiński, and Jakub Zakrzewski.  “Non-standard hubbard models in optical lattices: a review”.  Reports on Progress in Physics 78, 066001 (2015).
  * [38] Farokh Mivehvar, Francesco Piazza, Tobias Donner, and Helmut Ritsch.  “Cavity qed with quantum gases: new paradigms in many-body physics”.  Advances in Physics 70, 1–153 (2021).
  * [39] Daniel P Arovas, Erez Berg, Steven A Kivelson, and Srinivas Raghu.  “The hubbard model”.  Annual review of condensed matter physics 13, 239–274 (2022).
  * [40] Monika Aidelsburger, Luca Barbiero, Alejandro Bermudez, Titas Chanda, Alexandre Dauphin, Daniel González-Cuadra, Przemysław R Grzybowski, Simon Hands, Fred Jendrzejewski, Johannes Jünemann, et al.  “Cold atoms meet lattice gauge theory”.  Philosophical Transactions of the Royal Society A 380, 20210064 (2022).
  * [41] F. Bemani, R. Roknizadeh, and M.H. Naderi.  “Quantum simulation of discrete curved spacetime by the bose–hubbard model: From analog acoustic black hole to quantum phase transition”.  Annals of Physics 388, 186–196 (2018).
  * [42] Sean A. Hayward.  “Formation and evaporation of nonsingular black holes”.  Phys. Rev. Lett. 96, 031103 (2006).
  * [43] Jin-Sheng Peng.  “Introduction To Modern Quantum Optics”.  Wspc. Singapore (1998).  1st edition.
  * [44] Crispin Gardiner, Peter Zoller, and Peter Zoller.  “Quantum noise: a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics”.  Springer Science & Business Media. (2004).
  * [45] Heinz-Peter Breuer and Francesco Petruccione.  “The Theory of Open Quantum Systems”.  Oxford University Press. Oxford (2007).
  * [46] Remi Azouit, Alain Sarlette, and Pierre Rouchon.  “Adiabatic elimination for open quantum systems with effective Lindblad master equations” (2016).  url: http://arxiv.org/abs/1603.04630.
  * [47] Angelo Bassi, André Großardt, and Hendrik Ulbricht.  “Gravitational decoherence”.  Classical and Quantum Gravity 34, 193002 (2017).
  * [48] David Viennot and Olivia Moro.  “Adiabatic transport of qubits around a black hole”.  Classical and Quantum Gravity 34, 055005 (2017).
  * [49] B. Kramer and A. MacKinnon.  “Localization: theory and experiment”.  Rep. Prog. Phys. 56, 1469–1564 (1993).
  * [50] Subir Sachdev.  “Quantum phase transitions”.  Physics world 12, 33 (1999).
  * [51] E. M. Graefe and H. J. Korsch.  “Semiclassical quantization of an $\displaystyle n$-particle bose-hubbard model”.  Phys. Rev. A 76, 032116 (2007).
  * [52] Hagar Veksler and Shmuel Fishman.  “Semiclassical analysis of bose–hubbard dynamics”.  New Journal of Physics 17, 053030 (2015).
  * [53] Leo Brewin.  “Riemann normal coordinate expansions using cadabra”.  Classical and Quantum Gravity 26, 175017 (2009).
  * [54] Ravipudi Venkata Rao.  “Jaya: A simple and new optimization algorithm for solving constrained and unconstrained optimization problems”.  International Journal of Industrial Engineering Computations 7, 19–34 (2016).
  * [55] Joseph-Frédéric Bonnans, Jean Charles Gilbert, Claude Lemaréchal, and Claudia A Sagastizábal.  “Numerical optimization: theoretical and practical aspects”.  Springer Science & Business Media. (2006).
  * [56] Izrail Moiseevitch Gelfand, Richard A Silverman, and others.  “Calculus of variations”.  Courier Corporation. (2000).
  * [57] J. F. Colombeau.  “Elementary Introduction to New Generalized Functions”.  Elsevier. (2011).