# Machine Learning Statistical Gravity from Multi-Region Entanglement Entropy

Jonathan Lam  Yi-Zhuang You  Department of Physics, University of California,
San Diego, CA 92093, USA

###### Abstract

The Ryu-Takayanagi formula directly connects quantum entanglement and
geometry. Yet the assumption of static geometry lead to an exponentially small
mutual information between far-separated disjoint regions, which does not hold
in many systems such as free fermion conformal field theories. In this work,
we proposed a microscopic model by superimposing entanglement features of an
ensemble of random tensor networks of different bond dimensions, which can be
mapped to a statistical gravity model consisting of a massive scalar field on
a fluctuating background geometry. We propose a machine-learning algorithm
that recovers the underlying geometry fluctuation from multi-region
entanglement entropy data by modeling the bulk geometry distribution via a
generative neural network. To demonstrate its effectiveness, we tested the
model on a free fermion system and showed mutual information can be mediated
effectively by geometric fluctuation. Remarkably, locality emerged from the
learned distribution of bulk geometries, pointing to a local statistical
gravity theory in the holographic bulk.

##  I Introduction

The holographic duality [1, 2, 3, 4, 5] is a duality between boundary
$d$-dimensional quantum field theories and bulk $(d+1)$-dimensional
gravitational theories in asymptotically anti-de Sitter (AdS) space. It
provides an appealing explanation for the emergence of spacetime geometry from
quantum entanglement[6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]. The connection
is manifested in the Ryu-Takayanagi (RT) formula [17, 18]
$S(A)=\frac{1}{4G_{N}}\min_{\gamma_{A}}|\gamma_{A}|$ that relates the
entanglement entropy $S(A)$ of a boundary region $A$ to the area of the
extremal surface $\gamma_{A}$ in the bulk that is homologous to the same
region $A$. Progress has been made to reconstruct the bulk geometry from the
boundary data in terms of geodesic lengths[19, 20, 21, 22], extremal areas[23,
24, 25] or entanglement entropies[26, 27]. A majority of the effort has been
focused on reconstructing a classical geometry from _single-region_
entanglement entropies (or independent extremal surfaces). However, _multi-
region_ entanglement entropies further encode the correlation among multiple
extremal surfaces, which could reveal how the bulk geometry fluctuations
around its classical background (assuming a semiclassical description of the
bulk gravity). In this work, we will explore the possibility to extract
information about fluctuating holographic bulk geometries from multi-region
entanglement entropies of a quantum system using generative models in machine
learning.

A feature of the holography entanglement entropy based on the RT formula is
that the mutual information $I_{A:B}=S_{A}+S_{B}-S_{AB}$ vanishes between two
disjoint boundary regions $A$ and $B$ that are far separated from each
other[28, 29], because the minimum surface enclosing the combined region $AB$
will be a disjoint union of $\gamma_{A}$ and $\gamma_{B}$ such that the
entropies simply add up as $S_{AB}=S_{A}+S_{B}$, leaving no room for mutual
information. While the vanishing mutual information is a correct feature of
holographic conformal field theories (CFT), it is not generally the case for
many other quantum systems (e.g. free-fermion CFT). One idea to remedy the
problem is to introduce bulk matter fields to mediate the mutual information
[30, 31, 32]. Another possibility is to consider statistical fluctuations of
bulk geometries such that $\gamma_{A}$ and $\gamma_{B}$ are correlated to
produce the finite mutual information. The statistical gravitational
fluctuation may be viewed as an effective description arising from tracing out
bulk matter fields. We will further explore the second possibility of
fluctuating geometry using a concrete model of random tensor network (RTN)[33,
34] with fluctuating bond dimensions. The bond dimension fluctuation
translates to the bulk geometry fluctuation in the context of tensor network
holography[35, 36, 37], which is presumably governed by some statistical
gravity model.

However, it is unclear what should be the appropriate bulk statistical gravity
model that best reproduces the entanglement feature of a given quantum system
on the boundary. To address this challenge, we propose to apply data-driven
and machine-learning approaches to uncover the statistical gravity model
behind the observational data of quantum many-body entanglement. What needs to
be learned is the joint probability distribution of bond dimensions (or bulk
geometries). Generative models[38, 39, 40, 41, 42, 43, 44] in machine learning
provide us precisely the tool to learn unknown distributions from data. In
particular, we apply a deep generative model[39, 41] to describe the bulk
geometry fluctuation. We train the model by matching the model predictions of
multi-region entanglement entropies with their actual values evaluated in the
given quantum many-body state. After training, the generative model should
tell us the statistical gravity model that emerges from learning.

The approach developed in this work extends the general idea of entanglement
feature learning[26], which aims to reconstruct the bulk geometry by learning
from the entanglement data on the boundary. Compare to the previous work, this
study makes significant progress in including the gravitational fluctuation in
the model, which will enable us to learn an emergent gravity theory rather
than a static classical geometry. We will focus on (1+1)D quantum systems, and
assume that the system admits an approximate semiclassical geometry
description in the holographic bulk. Based on a random tensor network model
with fluctuating bond dimensions, we first establish a holographic model for
quantum entanglement involving a scalar matter field on a statistically
fluctuating spatial geometry. Applying our approach to a free-fermion CFT
state with a large central charge, we uncover a statistical gravity model
governed by Weyl field fluctuations propagating on the hyperbolic background
geometry. We show that the Weyl field fluctuation has the emergent bulk
locality by studying its bulk correlation. We further analyze the spectrum and
the leading collective modes of the emergent gravity theory. We also show that
the matter field mass gets renormalized by the gravitational fluctuation as
expected.

##  II Holographic Models of Entanglement

###  II.1 Random Tensor Network Model

The random tensor network (RTN) model is an intuitive toy model for
holographic duality, which directly connects quantum states and emergent
geometries. The original proposal[33] of RTN assumes a fixed bond dimension on
every link of the tensor network. It can be generalized to include bond
dimension fluctuations (or more precisely, bond entanglement fluctuations)[34,
45]. The generalized RTN model in consideration is defined as follows: (i) A
planar graph $G=(V,E)$ is given to describe the background network geometry,
where $V$ denotes the vertex set and $E$ denotes the edge set.
$V=V_{\text{blk}}\cup V_{\text{bdy}}$ is divided into two subsets: the bulk
$V_{\text{blk}}$ and the boundary $V_{\text{bdy}}$ sets, see Fig. 1(b). (ii) A
local Hilbert space $\mathcal{H}_{v}^{e}$ is associated with each pair $(v,e)$
of vertex $v\in V$ and its adjacent edge $e\in E$ (for $e$ not adjacent to
$v$, the associated Hilbert space is considered trivial
$\mathcal{H}_{v}^{e}\cong\mathbb{C}$), see Fig. 1(a). (iii) A random pure
state
$|\psi_{v}\rangle\in\mathcal{H}_{v}\equiv\bigotimes_{e\in\mathrm{d}v}\mathcal{H}_{v}^{e}$
is defined on every bulk vertex $v\in V_{\text{blk}}$. (iv) A random entangled
state $|\phi_{e}\rangle\in\mathcal{H}^{e}\equiv\bigotimes_{v\in\partial
e}\mathcal{H}_{v}^{e}$ is defined across every edge $e\in E$. (v) RTN defines
an ensemble $\mathcal{E}_{\text{RTN}}=\\{|\Psi\rangle\\}$ of pure states in
the boundary Hilbert space $\mathcal{H}_{\text{bdy}}\equiv\bigotimes_{v\in
V_{\text{bdy}}}\mathcal{H}_{v}$ by taking a (partial) projection in the bulk
Hilbert space $\mathcal{H}_{\text{blk}}\equiv\bigotimes_{v\in
V_{\text{blk}}}\mathcal{H}_{v}$ as

$|\Psi\rangle=\langle\psi|\phi\rangle:|\psi\rangle=\bigotimes_{v\in
V_{\text{blk}}}|\psi_{v}\rangle,|\phi\rangle=\bigotimes_{e\in
E}|\phi_{e}\rangle.$ (1)

The probability measure of $|\Psi\rangle$ in the RTN ensemble
$\mathcal{E}_{\text{RTN}}$ is given by
$P(|\Psi\rangle)=P(|\psi\rangle)P(|\phi\rangle)$. The vertex state
distribution $P(|\psi\rangle)=\prod_{v\in V_{\text{blk}}}P(|\psi_{v}\rangle)$
is assumed to be factorized, and on each vertex, the distribution
$P(|\psi_{v}\rangle)$ is taken to be the Haar measure (i.e. uniform random
states in $\mathcal{H}_{v}$). The edge (link) state distribution
$P(|\phi\rangle)$ is generally a nontrivial joint distribution depending on
all $|\phi_{e}\rangle$ on all edges, which allows the quantum entanglement
across different edges to fluctuate collectively.

Figure 1: (a) Details of the RTN near an edge. (b) The planar graph $G$ on
which the RTN is defined. The vertices are classified into the bulk vertices
(in red) and the boundary vertices (in green). (c) The dual graph $\tilde{G}$
of the RTN graph $G$.

For any operator $O^{(k)}$ defined in $k$ copies of the boundary Hilbert space
$\mathcal{H}_{\text{bdy}}^{\otimes k}$, its expectation value in the product
state $|\Psi\rangle^{\otimes k}$ is defined to be

$\langle
O^{(k)}\rangle=\mathop{\mathbb{E}}_{|\Psi\rangle\in\mathcal{E}_{\text{RTN}}}\frac{\operatorname{Tr}\big{(}(|\Psi\rangle\langle\Psi|)^{\otimes
k}O^{(k)}\big{)}}{\langle\Psi|\Psi\rangle^{k}}.$ (2)

We assume that the correlation between denominator and numerator is not
important (which is generally valid in the semiclassical regime when
fluctuations are weak), so that we can approximate the ensemble average of the
ratio by the ratio of separate averages,

$\langle
O^{(k)}\rangle\simeq\frac{1}{\mathcal{N}_{k}}\mathop{\mathbb{E}}_{|\Psi\rangle\in\mathcal{E}_{\text{RTN}}}\operatorname{Tr}\big{(}(|\Psi\rangle\langle\Psi|)^{\otimes
k}O^{(k)}\big{)},$ (3)

where
$\mathcal{N}_{k}=\mathop{\mathbb{E}}_{|\Psi\rangle\in\mathcal{E}_{\text{RTN}}}\langle\Psi|\Psi\rangle^{k}$
is the $k$th moment of the state norm squared. For example, the 2nd Rényi
entanglement entropy $S_{A}$ (or more precisely, the purity $e^{-S_{A}}$) of
RTN states in a boundary region $A$ can be calculated by taking $k=2$ and
$O^{(k)}=X_{A}$ (the swap operator supported in region $A$),

$e^{-S_{A}}\propto\mathop{\mathbb{E}}_{|\Psi\rangle\in\mathcal{E}_{\text{RTN}}}\operatorname{Tr}\big{(}(|\Psi\rangle\langle\Psi|)^{\otimes
2}X_{A}\big{)}.$ (4)

We will suppress the Rényi index throughout this work, and use $S_{A}$ to
denote the 2nd Rényi entropy. The RTN model provides an effective description
of entanglement entropies of typical quantum states on the holographic
boundary, given the background geometry $G$ together with fluctuations of
states $|\psi_{v}\rangle,|\phi_{e}\rangle$ in the holographic bulk.

It worth mention that in modeling the 2nd Rényi entanglement entropy by Eq.
(4), the average over the RTN ensemble $\mathcal{E}_{\text{RTN}}$ is taken
neither on the state vector level (i.e. not a pure state superposition
$\mathop{\mathbb{E}}_{|\Psi\rangle\in\mathcal{E}_{\text{RTN}}}|\Psi\rangle$),
nor on the density matrix level (i.e. not a mixed state superposition
$\mathop{\mathbb{E}}_{|\Psi\rangle\in\mathcal{E}_{\text{RTN}}}|\Psi\rangle\langle\Psi|$),
but on the double density matrix level (as
$\mathop{\mathbb{E}}_{|\Psi\rangle\in\mathcal{E}_{\text{RTN}}}(|\Psi\rangle\langle\Psi|)^{\otimes
2}$). The same average strategy commonly appeared in random tensor
network/quantum circuit literatures[33, 34, 26, 46, 47, 48]. Such average may
not have direct physical realization, nevertheless it defines a RTN model for
entanglement entropy which can produce (i) positive mutual information
$I_{A:B}$ that does not vanish between distant regions and (ii) possibly
negative tripartite information $I_{A:B:C}$ (see Appendix A for a perturbative
proof). These features indicate that the generalized RTN model is expressive
enough to describe quantum chaotic states with information scrambling[49, 50]
and to capture mutual information between distant entanglement regions, which
goes beyond holographic CFT states.

###  II.2 Ising and Dual Ising Models

Evaluating the ensemble average in Eq. (4) following the approach developed in
Ref. 33, the RTN purity $e^{-S_{A}}$ can be map to the partition function of
an Ising model on the graph $G$ with fluctuating coupling constants

$e^{-S_{A}}=\sum_{[\sigma,J]}P[\sigma|J]P[J]\,\delta[\sigma_{\text{bdy}}\Leftrightarrow
A],$ (5)

with $P[\sigma|J]$ given by

$\begin{split}P[\sigma|J]&=\frac{e^{-E[\sigma|J]}}{Z[J]},\\\
E[\sigma|J]&=-\sum_{e\in E}\bigg{(}\frac{J_{e}}{2}\prod_{v\in\partial
e}\sigma_{v}\bigg{)},\\\
Z[J]&=\sum_{[\sigma]}e^{-E[\sigma|J]}\delta[\sigma_{\text{bdy}}\Leftrightarrow\emptyset].\end{split}$
(6)

and $P[J]$ given by

$P[J]=\int_{|\phi\rangle}\prod_{e\in
E}\delta\big{(}J_{e}-S(|\phi_{e}\rangle)\big{)}P(|\phi\rangle).$ (7)

Here $\sigma_{v}=\pm 1$ is the Ising variable defined on every vertex $v\in
V$, $J_{e}\geq 0$ is the ferromagnetic coupling strength on every edge $e\in
E$. $J_{e}$ is determined by $S(|\phi_{e}\rangle)$, the 2nd Rényi entropy of
the state $|\phi_{e}\rangle$ (entangled between the Hilbert spaces
$\mathcal{H}_{v_{+}e}$ and $\mathcal{H}_{v_{-}e}$ where $v_{\pm}$ are the two
vertices on the boundary of $e$). $J_{e}$ characterizes how much the tensors
are entangled with each other across the edge $e$ in the tensor network. It
corresponds to the notion of bond dimension when $|\phi_{e}\rangle$ is
maximally entangled. The distribution $P[J]$ describes the how the effective
bond dimension (bond entanglement) fluctuates in the RTN ensemble. Finally,
the partition function is subject to the boundary condition that is set by the
boundary region $A$ of $S_{A}$,

$\forall v\in V_{\text{bdy}}:\sigma_{v}=\left\\{\begin{array}[]{cc}+1&v\notin
A,\\\ -1&v\in A,\end{array}\right.$ (8)

which is denoted as $\delta[\sigma_{\text{bdy}}\Leftrightarrow A]$ in Eq. (5).
The partition function $Z[J]$ properly normalizes the Boltzmann weight of the
Ising model, such that $S_{A}=0$ when the entanglement region $A=\emptyset$ is
empty.

Given that $G$ is a planar graph111The model can be more expressive if the
planar graph assumption is lifted, however it will be challenging to make
connection to the RTN model on non-planar graphs., we can use the Kramers-
Wannier duality to rewrite the Ising model Eq. (5) on the dual lattice
$\tilde{G}=(\tilde{V},\tilde{E})$, as shown in Fig. 1(c), where $\tilde{V}$
corresponds to the set of faces in $G$ and $\tilde{E}\cong E$. The dual Ising
model takes the similar form

$\begin{split}e^{-S_{A}}=\sum_{[\tilde{\sigma},\tilde{J}]}\bigg{(}\prod_{\tilde{v}\in\partial
A}\tilde{\sigma}_{\tilde{v}}\bigg{)}P[\tilde{\sigma}|\tilde{J}]P[\tilde{J}],\end{split}$
(9)

with $P[\tilde{\sigma}|\tilde{J}]$ given by

$\begin{split}P[\tilde{\sigma}|\tilde{J}]&=\frac{e^{-E[\tilde{\sigma}|\tilde{J}]}}{Z[\tilde{J}]},\\\
E[\tilde{\sigma}|\tilde{J}]&=-\sum_{\tilde{e}\in\tilde{E}}\bigg{(}\frac{\tilde{J}_{\tilde{e}}}{2}\prod_{\tilde{v}\in\partial\tilde{e}}\tilde{\sigma}_{\tilde{v}}\bigg{)},\\\
Z[\tilde{J}]&=\sum_{[\tilde{\sigma}]}e^{-E[\tilde{\sigma}|\tilde{J}]},\end{split}$
(10)

and $P[\tilde{J}]$ related to $P[J]$ by

$P[\tilde{J}]=\bigg{(}\prod_{e}\frac{\partial
J_{e}}{\partial\tilde{J}_{\tilde{e}}}\bigg{)}P[J].$ (11)

Here $\tilde{\sigma}_{\tilde{v}}=\pm 1$ is the dual Ising variable and
$\tilde{J}_{\tilde{e}}=-\ln\tanh(J_{e}/2)$ is the dual coupling. The boundary
condition in the original Ising model translates to the insertion of the dual
Ising variable at every boundary point of entanglement region $A$ (i.e. at
every entanglement cut). The partition function $Z[\tilde{J}]$ on the
denominator ensures $S_{A}=0$ when the entanglement region $A=\emptyset$ is
empty, i.e. when there is no insertion of dual Ising variables. Both
$P[\tilde{\sigma}|\tilde{J}]$ and $P[\tilde{J}]$ are normalized probability
distributions, which defines the joint distribution
$P[\tilde{\sigma},\tilde{J}]=P[\tilde{\sigma}|\tilde{J}]P[\tilde{J}]$ for dual
Ising variables and their couplings. Therefore the purity of the RTN state can
be interpreted as the boundary correlation of dual Ising variables
$e^{-S_{A}}=\langle\prod_{\tilde{v}\in\partial
A}\tilde{\sigma}_{\tilde{v}}\rangle$ in an Ising model with fluctuating
couplings.

Figure 2: (a) Two-point correlation and (b) four-point correlation of dual
Ising spins.

The RT formula can be recovered in the classical limit when the RTN bond
dimensions are large and fixed, which corresponds to the deep ferromagnetic
phase of the original Ising model ($J_{e}\gg 1$) or equivalently the deep
paramagnetic phase of the dual Ising model ($\tilde{J}_{\tilde{e}}\ll 1$). In
such limit, the dual Ising correlation decays exponentially with the geodesic
distance $\langle\tilde{\sigma}_{1}\tilde{\sigma}_{2}\rangle\propto
e^{-|\gamma_{12}|/\xi}$, as illustrated in Fig. 2(a), which reproduces the RT
formula $S_{A}=|\gamma_{12}|/\xi$ with some appropriate choice of the
correlation length $\xi$. Multi-region entanglement entropies will correspond
to higher-point correlations functions, such as
$e^{-S_{AB}}\sim\langle\tilde{\sigma}_{1}\tilde{\sigma}_{2}\tilde{\sigma}_{3}\tilde{\sigma}_{4}\rangle\sim
e^{-|\gamma_{12}|/\xi}e^{-|\gamma_{34}|/\xi}$ in Fig. 2(b). Allowing the dual
Ising coupling $\tilde{J}$ to fluctuate collectively will introduce
perturbations to the geodesic distance
$|\gamma_{12}|\to|\gamma_{12}|+\delta|\gamma_{12}|$ in a correlated manner,
such that

$\begin{split}e^{-S_{AB}}&\sim\mathop{\mathbb{E}}\big{(}e^{-(|\gamma_{12}|+\delta|\gamma_{12}|)/\xi}e^{-(|\gamma_{34}|+\delta|\gamma_{34}|)/\xi}\big{)}\\\
&\sim
e^{-|\gamma_{12}|/\xi}e^{-|\gamma_{34}|/\xi}\;e^{\frac{1}{2\xi^{2}}\mathop{\mathbb{E}}\delta|\gamma_{12}|\delta|\gamma_{34}|}\\\
&\sim e^{-S_{A}}e^{-S_{B}}e^{I_{A:B}}.\end{split}$ (12)

Thus the correlated geometric fluctuation provides an effective mechanism to
generate the mutual information between far-separated regions $A$ and $B$
(beyond the classical RT formula). Therefore we anticipate the fluctuating RTN
model to be a more expressive holographic model for entanglement entropies.
However, it is not clear how the dual Ising coupling $\tilde{J}$ (or the
effective bond dimension $J$) should fluctuate precisely in order to
quantitatively reproduce all multi-region entanglement entropies of a given
quantum many-body state. The remaining task is learn the distribution
$P[\tilde{J}]$ (or other equivalent distributions) from data.

###  II.3 Effective Statistical Gravity Model

Suppose the fluctuation of $\tilde{J}$ is small around its static background
configuration, such that there is a meaningful notion of background geometry
in the bulk. The dual Ising model can be described by an effective field
theory in the continuum limit

$S[\phi|g]=\frac{1}{2}\int\mathrm{d}^{2}x\sqrt{g}(g^{ij}\partial_{i}\phi\partial_{j}\phi+m^{2}\phi^{2}),$
(13)

where the dual Ising variable $\tilde{\sigma}_{\tilde{v}}$ is coarse-grained
to a massive real scalar field $\phi(x)$, as the Ising model universally flows
to this massive Gaussian fixed point in the paramagnetic phase. The theory is
defined in the holographic space (without time dimension). The fluctuating
Ising coupling $\tilde{J}$ can be translated to a fluctuating bulk metric
tensor $g$ around a reference background geometry $\bar{g}$ 222Another
interpretation is to translate the fluctuating Ising coupling to the
fluctuating mass term, as the scalar field mass is the relevant perturbation
that drives the order-disorder transition, which plays the same role as the
Ising coupling. This alternative view turns out to be equivalent to the
fluctuating metric interpretation in two-dimension, as we will see soon.,
since a stronger local coupling creates a larger local correlation, which
effectively reduces the local distance measure
$\mathrm{d}s^{2}=g_{ij}\mathrm{d}x^{i}\mathrm{d}x^{j}$ between the correlated
Ising variables. Therefore the purity of RTN state Eq. (9) can be effectively
described by a statistical gravity model

$\begin{split}e^{-S_{A}}&=\int_{[\phi,g]}\bigg{(}\prod_{x\in\partial
A}\phi(x)\bigg{)}P[\phi|g]P[g],\\\
P[\phi|g]&=\frac{e^{-S[\phi|g]}}{Z[g]},Z[g]=\int_{[\phi]}e^{-S[\phi|g]},\end{split}$
(14)

where the gravity is “quenched” in the sense that the metric configuration is
generated with a probability distribution $P[g]$ independent of the scalar
field $\phi$ configuration.

In two-dimensional space, the metric tensor has three independent components.
However two of them can be removed by gauge transformation $g_{ij}\to
g_{ij}+\nabla_{i}\xi_{j}+\nabla_{j}\xi_{i}$. We can choose the conformal gauge
where the metric tensor $g_{ij}(x)$ is parametrized by a Weyl field
$\omega(x)$ that rescales a fixed background $\bar{g}_{ij}(x)$

$g_{ij}(x)=e^{2\omega(x)}\bar{g}_{ij}(x),$ (15)

such that each Weyl field configuration represents a physically distinct
geometry. As a result, the integration $\int_{[g]}P[g]$ can be replaced by
$\int_{[\omega]}P[\omega]$ in Eq. (14). The unknown joint distribution
$P[\omega]$ will be what we aim to learn from the entanglement entropy data.

To numerically evaluate the multi-point scalar field correlation, we can place
the bulk field theory back on a lattice, say on the dual graph
$\tilde{G}=(\tilde{V},\tilde{E})$. Using Regge calculus[53] to discretize the
action,

$S[\phi|\omega]=\sum_{\langle
xy\rangle\in\tilde{E}}\frac{A_{xy}}{2}\Big{(}\frac{\phi_{x}-\phi_{y}}{\ell_{xy}}\Big{)}^{2}+\sum_{x\in\tilde{V}}\frac{m^{2}A_{x}}{2}e^{2\omega_{x}}\phi_{x}^{2},$
(16)

where $\ell_{xy}$ can be interpreted as the geodesic distance between two
vertices $x$ and $y$ on the background geometry. $A_{x}$ and $A_{xy}$ are the
areas associated to the vertex $x$ and the edge $\langle xy\rangle$
respectively. $\ell_{xy},A_{x},A_{xy}$ are all fixed according to the choice
of background metric, which will be specified later. The statistical variables
in the model are the scalar field $\phi_{x}$ and the Weyl field $\omega_{x}$
in the holographic bulk. The model predicts the entanglement entropy on the
holographic boundary by

$\begin{split}e^{-S_{A}}&=\int_{[\phi,\omega]}\bigg{(}\prod_{x\in\partial
A}\phi_{x}\bigg{)}P[\phi|\omega]P[\omega],\\\
P[\phi|\omega]&=\frac{e^{-S[\phi|\omega]}}{Z[\omega]},Z[\omega]=\int_{[\phi]}e^{-S[\phi|\omega]}\end{split}$
(17)

which is the underlying lattice model that will be used in the machine
learning algorithm. The unknown distribution $P[\omega]$ will be parameterized
by a generative model. By matching the model prediction with the actual data
of entanglement entropies calculated from a quantum state, the algorithm can
reconstruct the distribution $P[\omega]$ and infer the statistical gravity
model behind the entanglement structure.

##  III Machine Learning Method

###  III.1 Generative Modeling

Generative modeling is about learning probability distributions[54]. We will
apply the simplest latent-variable generative model [55] in this work. The
basic idea is to start with a easy-to-sample prior distribution, such as a
Gaussian distribution. Draw a random vector $z\in\mathbb{R}^{n}$ (as a
collection of latent variables) from the prior distribution $P(z)$. Then
transform the latent variables $z$ by a deep neural network $G_{\vartheta}$
(parametrized by some variational parameter $\vartheta$) to the designated
random variable $\omega$, i.e. $z\to\omega=G_{\vartheta}(z)$. The mapping
$G_{\vartheta}$ is called the generator, which defines the distribution of
generated samples

$P_{\vartheta}[\omega]=\int\mathrm{d}z\;\delta(\omega-G_{\vartheta}(z))P(z).$
(18)

A large batch of $\omega$ can be sampled efficiently in parallel, when
hardware accelerators (e.g GPU or TPU) are available. If the neural network
$G_{\vartheta}$ is expressive enough, Eq. (18) will provide a sufficiently
expressive probability model $P_{\vartheta}[\omega]$ for the Weyl field
$\omega$ configuration.

The distribution $P_{\vartheta}[\omega]$ defines the model prediction of the
purity based on Eq. (17)

$e^{-S_{A}|_{\vartheta}}=\mathop{\mathbb{E}}_{[\omega]\sim
P_{\vartheta}}\langle\textstyle\prod_{\partial A}\phi\rangle_{\omega}.$ (19)

We will use $S_{A}|_{\vartheta}$ to denote the Rényi entropy predicted by the
machine learning model as it depends on the model parameters $\vartheta$. In
Eq. (19), we introduced the short-hand notation

$\langle{\textstyle\prod_{\partial
A}}\phi\rangle_{\omega}\equiv\int_{[\phi]}\bigg{(}\prod_{x\in\partial
A}\phi_{x}\bigg{)}P[\phi|\omega]$ (20)

to denote the scalar field correlation on a background Weyl field
configuration. The conditional distribution $P[\phi|\omega]$ is defined in Eq.
(17) with $S[\phi|\omega]$ given in Eq. (16). The scalar field correlation
$\langle{\textstyle\prod_{\partial A}}\phi\rangle_{\omega}$ can be efficiently
evaluated when $S[\phi|\omega]$ is a Gaussian action (which is the case here).

Figure 3: Flow diagram of the machine learning algorithm. Black arrows denotes
the forward evaluation of the loss function. Red dashed arrows denotes the
gradient back propagation to train the parameter.

The task is to learn the optimal Weyl field distribution
$P_{\vartheta}[\omega]$ that gives the best prediction of the purity data
based on Eq. (19). The dataset will contain the purity $\\{e^{-S_{A}}\\}$ of a
quantum state in different regions $A$. The distribution
$P_{\vartheta}[\omega]$ can be learned by optimizing model parameters
$\vartheta$ to minimize the following loss function (to be explained later)

$\mathcal{L}_{\vartheta}=\mathop{\operatorname{avg}}_{A}\big{(}1-e^{S_{A}|_{\vartheta}-S_{A}}\big{)}^{2}.$
(21)

As illustrated in Fig. 3, the training initiates from randomly choosing a
batch of entanglement regions $A$. On one hand, we query the dataset to get
the ground truth of $S_{A}$. On the other hand, a collection of Weyl field
configurations are sampled from the generative model, based on which the model
prediction $S_{A}|_{\vartheta}$ is estimated. Then the loss function is
calculated by comparing $S_{A}|_{\vartheta}$ with $S_{A}$, and the gradient
signal propagates back to train the parameters via gradient descent
$\vartheta\to\vartheta-r\partial_{\vartheta}\mathcal{L}_{\vartheta}$. After
some iterations, the parameters are expected to converge. In the following, we
will explain different modules in Fig. 3 in detail.

###  III.2 Entanglement Dataset

While efficient experimental approaches[56, 57] have been developed to
estimate Rényi entropies from randomized measurements, which enables the
acquisition of a large amount of entanglement data to drive the entanglement
feature learning, preparing an entanglement dataset by numerically computing
entanglement entropies from a given quantum many-body state remains difficult
in general. As a proof of concept, we choose to use the ground state of a free
fermion system for demonstration, on which entanglement entropies can be
efficiently calculated.

Consider $N$ copies of the (1+1)D massless Majorana fermion chain, described
by the Hamiltonian

$H=\sum_{a=1}^{N}\sum_{j}\mathrm{i}\chi_{j,a}\chi_{j+1,a},$ (22)

where $\\{\chi_{i,a},\chi_{j,b}\\}=\delta_{ij}\delta_{ab}$. Let $|\Psi\rangle$
be the ground state of $H$. The 2nd Rényi entropy can be efficiently computed
from the fermion correlation function,

$S_{A}=-\frac{1}{2}\operatorname{Tr}\ln(C_{A}^{2}+(1-C_{A})^{2}),$ (23)

where $C_{A,ij}=\langle\Psi|\chi_{i,a}\chi_{j,a}|\Psi\rangle$ (for $i,j\in A$)
is the two-point correlation function (matrix) of Majorana fermions restricted
inside the entanglement region $A$. The quantum system is critical and is
described by the free-fermion CFT at low energy.

To construct the dataset, we will take the Majorana fermion chain of 32 sites,
and randomly sample a large collection of single-region, two-region, and
three-region subsets. We then compute the entanglement entropy using Eq. (23)
for every region and record the results in the entanglement dataset.

###  III.3 Bulk Model Solver

The bulk model solver is expected to calculate the scalar field correlation
given the Weyl field background $\omega$ and the entanglement region $A$ that
specifies the scalar field inserting position on the boundary. We will use the
lattice model specified by the action in Eq. (16), which describe a free
scalar field $\phi$. The action can be written as the bilinear form

$S[\phi|\omega]=\frac{1}{2}\sum_{x,y\in\tilde{V}}\phi_{x}K^{(\phi)}_{xy}[\omega]\phi_{y}$
(24)

where $x,y$ label the vertices on the dual graph
$\tilde{G}=(\tilde{V},\tilde{E})$ on which the holographic model is defined.
The kernel matrix takes the form of $K^{(\phi)}[\omega]=\nabla^{2}+M[\omega]$,
with

$M_{xy}[\omega]=m^{2}A_{x}e^{2\omega_{x}}\delta_{xy}$ (25)

being the mass term, and

$\nabla^{2}_{xy}=\sum_{\langle
x^{\prime}y^{\prime}\rangle\in\tilde{E}}\frac{A_{x^{\prime}y^{\prime}}}{\ell_{x^{\prime}y^{\prime}}^{2}}(\delta_{xx^{\prime}}-\delta_{xy^{\prime}})(\delta_{yx^{\prime}}-\delta_{yy^{\prime}})$
(26)

being the discrete Laplace operator on the dual graph. The length $\ell_{xy}$
and area $A_{x},A_{xy}$ constants are fixed and are set by the background
geometry, as to be specified soon. The two-point correlation is given by the
inverse of the kernel matrix,

$\langle\phi_{x}\phi_{y}\rangle_{\omega}=\frac{1}{\beta}((K^{(\phi)}[\omega])^{-1})_{xy},$
(27)

where a trainable constant $\beta$ is introduced to take care of the field
renormalization. Higher-point correlations follow from Wick’s theorem. For
example,

$\begin{split}\langle\phi_{x}\phi_{y}\phi_{z}\phi_{w}\rangle_{\omega}&=\langle\phi_{x}\phi_{y}\rangle_{\omega}\langle\phi_{z}\phi_{w}\rangle_{\omega}\\\
&+\langle\phi_{x}\phi_{z}\rangle_{\omega}\langle\phi_{y}\phi_{w}\rangle_{\omega}\\\
&+\langle\phi_{x}\phi_{w}\rangle_{\omega}\langle\phi_{y}\phi_{z}\rangle_{\omega}.\end{split}$
(28)

We will treat $\beta$ and $m^{2}$ as trainable parameters, which will be
optimized (together with other model parameters for $P[\omega]$) to fit the
entanglement data.

Since we intend to apply our approach to entanglement data collected from
CFTs, following the idea of AdS/CFT correspondence, it is natural to choose
the two-dimension hyperbolic geometry (the spatial slice of AdS3) as the
background geometry. We use the following background metric

$\mathrm{d}s^{2}=\mathrm{d}\rho^{2}+\sinh^{2}\rho\;\mathrm{d}\theta^{2},$ (29)

where $0\leq\theta<2\pi$ and $\rho\leq\rho_{\text{bdy}}$ (the UV cutoff scale
is set by $\rho_{\text{bdy}}$, which is another parameter to learn). The
geodesic distance any two points on the boundary $\rho=\rho_{\text{bdy}}$
separated by $\theta$ is given by

$\begin{split}|\gamma|(\theta)&=\text{arccosh}(1+2\sinh^{2}\rho_{\text{bdy}}\sin^{2}(\theta/2))\\\
&\xrightarrow{e^{\rho_{\text{bdy}}}\gg
1}2\ln(\sin(\theta/2)+\rho_{\text{bdy}}).\end{split}$ (30)

Figure 4: Triangular lattice discretization of hyperbolic space in
$(\rho,\theta)$ coordinate. The lattice is divided into different layers along
the the radius direction. The $i$th layer corresponds to the radius
$\rho_{i}$.

Without loss of generality, we chose to discretize space using a triangular
lattice with periodic boundary condition along the $\theta$-direction. All
vertices in the same layer are of the same $\rho$-coordinate and their
$\theta$-coordinates are uniformly spaced, see Fig. 4. The geodesic distance
$\ell_{xy}$ between two vertices $x$ and $y$ is given by

$\begin{split}\cosh\ell_{xy}&=\cosh\rho_{x}\cosh\rho_{y}\\\
&-\sinh\rho_{x}\sinh\rho_{y}\cos(\theta_{x}-\theta_{y}).\end{split}$ (31)

The area of an elementary triangle in the $i$th layer reads

$\begin{split}\tan\frac{A_{\bigtriangleup,i}}{4}&=\tanh\Big{(}\frac{\rho_{i+1}-\rho_{i}}{2}\Big{)}\tanh\Big{(}\frac{b_{i}}{4}\Big{)},\\\
\tan\frac{A_{\bigtriangledown,i}}{4}&=\tanh\Big{(}\frac{\rho_{i}-\rho_{i-1}}{2}\Big{)}\tanh\Big{(}\frac{b_{i}}{4}\Big{)},\\\
\cosh b_{i}&=\cosh^{2}\rho_{i}-\sinh^{2}\rho_{i}\cos\Delta\theta,\end{split}$
(32)

which defines the vertex and edge areas in a barycentric scheme. Specifically,
the vertex area $A_{x}$ is given by

$\begin{split}A_{x}=\tfrac{1}{3}(2A_{\bigtriangleup,i}+2A_{\bigtriangledown,i}+A_{\bigtriangleup,i-1}+A_{\bigtriangledown,i+1}),\end{split}$
(33)

for $\rho_{x}=\rho_{i}$. The edge area $A_{xy}$ is given by

$A_{xy}=\left\\{\begin{array}[]{ll}\tfrac{1}{3}(A_{\bigtriangleup,i}+A_{\bigtriangledown,i})&\rho_{x}=\rho_{y}=\rho_{i};\\\
\tfrac{1}{3}(A_{\bigtriangleup,i}+A_{\bigtriangledown,i+1})&\rho_{x}=\rho_{i},\rho_{y}=\rho_{i+1};\\\
\tfrac{1}{3}(A_{\bigtriangleup,i-1}+A_{\bigtriangledown,i})&\rho_{x}=\rho_{i},\rho_{y}=\rho_{i-1}.\\\
\end{array}\right.$ (34)

These equations defines $\ell_{xy}$, $A_{x}$ and $A_{xy}$ used in the lattice
model Eq. (16), which all rely on the values of $\rho_{i}$ for different
layers. The discretization scheme in the radial dimension is specified by how
$\rho_{i}$ is spaced from $0$ to $\rho_{\text{bdy}}$. A bad choice of the
discretization scheme may cause some triangle elements to have high aspect
ratios, reducing the quality of the triangulation in approximating the
continuous background geometry. We will take a data-driven approach to learn
the optimal discretization scheme by treating $\\{\rho_{i}\\}$ as trainable
parameters.

To summarize, the bulk model solver contains the following parameters: the
normalization $\beta$ and the squared mass $m^{2}$ associated with the scalar
field dynamics, and the radial coordinates $\\{\rho_{i}\\}$ associated with
the discretization of background geometry. These parameters will be trained
together with other neural network parameters (see Sec. III.4) to optimize the
model prediction of the entanglement entropy data.

###  III.4 Neural Network Design

The central goal is to learn the Weyl field distribution $P[\omega]$ using a
latent-variable generative model $P_{\vartheta}[\omega]$, recall Eq. (18). The
key component of the generative model is a generator $G_{\vartheta}$ that maps
the latent variable $z$ to a Weyl field configuration
$\omega=G_{\vartheta}(z)$. The generator is realized as a deep neural network
consists of consecutive layers of simpler maps

$G_{\vartheta}(z)=g_{N}\circ\cdots g_{2}\circ g_{1}(z),$ (35)

where each layer $g_{n}(z)$ is an affine transformation followed by some non-
linearity such as ReLU [58]. The weight and bias parameters are introduced to
parametrize the affine transformations, which constitute part of the training
parameters $\vartheta$.

It is both practical and theoretically motivated to enforce the neural
network’s architecture such that the learned distribution
$P_{\vartheta}[\omega]$ will respects certain symmetries, i.e. to construct an
equivariant neural network [59]. Let $Q$ be a symmetry transformation that we
wish to impose. The sufficient condition for the generated distribution to be
symmetric (i.e. $P_{\vartheta}[Q\omega]=P_{\vartheta}[\omega]$) is to require
(i) $Qg_{n}(z)=g_{n}(Qz)$ and (ii) $P(Qz)=P(z)$. The symmetries in
consideration are

  1. 1.

$\omega_{(\rho,\theta)}\to\omega_{(\rho,\theta+a)}$ [translation],

  2. 2.

$\omega_{(\rho,\theta)}\to\omega_{(\rho,-\theta)}$ [reflection].

The translation symmetry can be imposed by parameter sharing between relation-
related weights and biases, making the affine transformation in each layer
effectively a convolution along the translation direction. The reflection
symmetry can be imposed by using a reflection symmetric convolution kernel.
The prior distribution $P(z)=\prod_{x}P(z_{x})$ automatically satisfies the
symmetry condition as it factorizes to identical independent Gaussian
distributions on every site.

We would like to emphasize that although each layer looks like a convolutional
layer under the symmetry constraint, we do not restrict the convolution kernel
to be local (the kernel size extends to the whole lattice), because we do not
want to impose locality by hand. As we will see, a sense of locality could
emerge in the neural network as the holographic model gets trained, which
corresponds to the emergent locality in the bulk gravity theory.

###  III.5 Loss Function Design

The loss function is designed to evaluate the average difference between the
purity $e^{-S_{A}|_{\vartheta}}$ predicted by the holographic model and the
purity $e^{-S_{A}}$ given by the entanglement data. A straightforward option
would be the mean squared error (MSE) loss

$\mathcal{L}_{\vartheta,\text{MSE}}=\mathop{\operatorname{avg}}_{A}(e^{-S_{A}|_{\vartheta}}-e^{-S_{A}})^{2}.$
(36)

The model prediction $e^{-S_{A}|_{\vartheta}}$ should be evaluated according
to Eq. (19), which involves the ensemble expectation
$\mathop{\mathbb{E}}_{[\omega]\sim P_{\vartheta}}$. In practice, the
expectation value can only be estimated by sampling a finite number of Weyl
field configurations from the generative model $P_{\vartheta}$ and take the
average

$e^{-S_{A}|_{\vartheta}}=\frac{1}{N_{\omega}}\sum_{[\omega]\sim
P_{\vartheta}}\langle{\textstyle\prod}_{\partial A}\phi\rangle_{\omega},$ (37)

where $N_{\omega}$ denotes the number of Weyl field samples. With the help of
modern GPU, $\langle{\textstyle\prod}_{\partial A}\phi\rangle_{\omega}$ can be
computed in parallel efficiently. The sample size $N_{\omega}$ is thus
ultimately limited by the GPU memory. In our case, $N_{\omega}$ ranges from
512 to 2048.

For any finite sample size $N_{\omega}$, the finite average
$e^{-S_{A}|_{\vartheta}}$ will have a finite variance, which bias the MSE loss

$\mathcal{L}_{\vartheta,\text{MSE}}=\mathop{\operatorname{avg}}_{A}\bigg{(}(\mathop{\mathbb{E}}e^{-S_{A}|_{\vartheta}}-e^{-S_{A}})^{2}+\frac{\mathop{\operatorname{var}}e^{-S_{A}|_{\vartheta}}}{N_{\omega}}\bigg{)},$
(38)

causing the parameter to converge to a wrong saddle point. The bias can be
corrected by assigning a larger weight to the prediction with a higher
precision, i.e.

$\mathop{\operatorname{avg}}_{A}\frac{(e^{-S_{A}|_{\vartheta}}-e^{-S_{A}})^{2}}{\mathop{\operatorname{var}}e^{-S_{A}|_{\vartheta}}},$
(39)

which can also be argued from the maximum-likelihood estimation. The variance
is generally proportional to the square of the purity
$\mathop{\operatorname{var}}e^{-S_{A}|_{\vartheta}}\propto(e^{-S_{A}|_{\vartheta}})^{2}$,
which leads to the mean squared relative error (MSRE) loss

$\begin{split}\mathcal{L}_{\vartheta,\text{MSRE}}&=\mathop{\operatorname{avg}}_{A}\Big{(}\frac{e^{-S_{A}|_{\vartheta}}-e^{-S_{A}}}{e^{-S_{A}|_{\vartheta}}}\Big{)}^{2}\\\
&=\mathop{\operatorname{avg}}_{A}\big{(}1-e^{S_{A}|_{\vartheta}-S_{A}}\big{)}^{2}.\end{split}$
(40)

We numerically test the loss function by generating some data using a model
with known parameters, and train new models with different loss function on
the generated data to see if the parameter converges to the known result. Our
test shows that Eq. (40) indeed converges better compare to Eq. (36).
Therefore, we will use the MSRE loss function to train the model, as mentioned
in Eq. (21).

##  IV Numerical Results

###  IV.1 Fitting Entanglement Data with Static and Fluctuating Geometry

We apply the proposed machine learning approach to learn the entanglement
feature of a Majorana fermion chain of 32 sites (16 unit cells) with a
relatively large central charge $c=8$. The entanglement data is partitioned
into the training set and the test set that does not overlap. Within the
training/test set, the data can be further classified by the number of
subregions of the entanglement region, including the single-region, double-
region, and triple-region entanglement. To demonstrate the effect of
introducing gravitational fluctuations, we will compare two holographic
models: (i) the _fluctuating_ model, i.e. the model
$e^{-S_{A}|_{\vartheta}}=\mathop{\mathbb{E}}_{[\omega]\sim
P_{\vartheta}}\langle\textstyle\prod_{\partial A}\phi\rangle_{\omega}$
proposed in Eq. (19) with fluctuating geometries, (ii) the _static_ model,
i.e. the model $e^{-S_{A}|_{\vartheta}}=\langle\textstyle\prod_{\partial
A}\phi\rangle_{\omega\equiv 0}$ with a fixed static geometry. We train both
models using the MSRE loss in Eq. (21). The algorithm is implemented in the
TensorFlow[61] framework using the ADAM[62] optimizer. Upon convergence, the
MSRE loss is evaluated on the test set to characterize the performance of the
model. The result is summarized in Tab. 1

Model | static | static | fluctuating
---|---|---|---
Training set | single | single+double | single+double
Test set | single | $8.7\times 10^{-6}$ | $2.1\times 10^{-2}$ | $1.5\times 10^{-3}$
double | $1.1\times 10^{-1}$ | $3.9\times 10^{-2}$ | $5.7\times 10^{-3}$
triple | $7.5\times 10^{-1}$ | $6.0\times 10^{-1}$ | $3.1\times 10^{-1}$
Table 1: MSRE loss on test sets for different models and training sets. The
model can be static geometry or fluctuating geometry. The training set can
include on single-regions or both single- and double-regions. The test set can
be either single-, double-, or triple-regions separately.

If we train the static geometry model with single-region data only, the model
can easily achieve high accuracy ($\text{MSRE}\sim 10^{-5}$) in predicting
single-region entanglement, as also shown in Fig. 5(a). But the prediction of
multi-region entanglement is rather inaccurate ($\text{MSRE}\sim 10^{-1}$),
meaning that the static geometry model overfits the single-region data and can
not be generalized to multi-region data. If we include the double-region data
in the training set, and train the static geometry model with both single- and
double-region entanglement, the model will learn to predict double-region
entanglement better at the price of losing the accuracy in predicting single-
region entanglement, with the MSRE saturates at the $\sim 10^{-2}$ level. This
implies an intrinsic conflict for the static geometry model in modeling the
single- and multi-region entanglement simultaneously.

However, by introducing gravitational fluctuations to the model, the
fluctuating geometry model achieves one order of magnitude improvement in the
prediction accuracy of both single- and double-region entanglement, as the
MSRE drops to the $\sim 10^{-3}$ level, which is also manifest in Fig. 5(b).
This indicates that the gravitational fluctuation indeed helps to reconcile
the conflict between single- and multi-region entanglements in the classical
gravity model (RT formula) (see Appendix B for an analytic analysis of how the
conflict can be reconciled in principle). Moreover, the prediction accuracy in
triple-regions is also improved significantly, even if the model is never
trained on the triple-region data. This speaks for the better generalizability
of the fluctuating geometry model.

Figure 5: Model predicted entanglement entropy v.s. ground truth for (a)
static model trained on single-region entanglement data, (b) static and
fluctuating models trained on both single-region and double-region
entanglement data. Figure 6: Mutual information $I_{A:B}$ between two equal-
sized regions $A$ and $B$ for different inter-region separations. Data points
show predictions by the static model (red) and the fluctuating model (blue).
Both models are trained with _both_ single- and double-region entanglement
data.

As argued previously, the static geometry model suffers from the problem of
vanishing mutual information between far separated regions. One motivation to
introduce gravitational fluctuations is to mediate the mutual information
between distant regions through the holographic bulk. Indeed, as shown in Fig.
6, by allowing the geometry to fluctuate, the model can better capture the
behavior of mutual information. In particular, the static model fails to
produce the non-vanishing mutual information between distant regions, which is
most obviously seen in Fig. 6(a), where the regions are most far separated
(compare to their sizes). However, the fluctuating model fixes this problem,
demonstrating the importance of introducing the geometric fluctuation in
modeling the multi-region entanglement.

###  IV.2 Weyl Field Correlation and Effective Bulk Gravity Theory

After training, we want to open up the model and see what bulk gravity theory
has been learned. With the trained generative model $P_{\vartheta}[\omega]$
that describes the statistical fluctuation of the Weyl field $\omega$, we can
explore various statistical properties of the distribution
$P_{\vartheta}[\omega]$ to gain a deeper understanding of the optimal bulk
gravity theory that emerges from learning the boundary entanglement data.

Figure 7: (a) Weyl field local covariance $\Sigma^{(\omega)}_{xx}$ v.s. the
radius coordinate $\rho_{x}$. (b) Intra-layer correlation $C^{(\omega)}_{xy}$
of the Weyl field (in logarithmic scale) v.s. the bulk geodesic distance
$d_{xy}$, where $x$, $y$ points are taken from the same layer. The distance
between points on a larger radius (or a higher layer) appears farther due to
the hyperbolic background geometry, but the inverse correlation length
$\Delta$ (the slope) remains roughly on the same order of magnitude across
layers.  Figure 8: Weyl field correlation in the bulk, between a marked
reference point and the remaining points. Results are shown for the reference
point placed in different layers. Figure 9: Spectrum of the Weyl field
fluctuation. (a) Leading eigenvalues $\lambda^{(\alpha)}$ of the covariance
function (in logarithmic scale). (b) Selected eigenmodes $\omega^{(\alpha)}$
labeled by the principal and angular quantum numbers.

We first study the covariance function $\Sigma^{(\omega)}_{xy}$ of the Weyl
field $\omega$, defined as

$\Sigma^{(\omega)}_{xy}=\mathop{\mathbb{E}}_{[\omega]\sim
P_{\vartheta}}\omega_{x}\omega_{y}=\int_{[\omega]}P_{\vartheta}[\omega]\omega_{x}\omega_{y}.$
(41)

We observe that its diagonal elements $\Sigma^{(\omega)}_{xx}$ (i.e. the local
covariance) grows with the radius $\rho_{x}$ coordinate (as $x$ approaches the
boundary), see Fig. 7(a). This is because the discretization scale is changing
along the radius direction. In our discretization scheme as shown in Fig. 4,
the hyperbolic space is finer discretized towards the center of the bulk,
therefore the field $\omega$ will appear to be stiffer near the bulk center,
and hence its covariance is smaller. To eliminate this influence of the
discretization scheme, we normalize (standardize) the covariance and define
the correlation function

$C^{(\omega)}_{xy}=\frac{\Sigma^{(\omega)}_{xy}}{\sqrt{\Sigma^{(\omega)}_{xx}\Sigma^{(\omega)}_{yy}}}.$
(42)

We found that the Weyl field correlation $C^{(\omega)}_{xy}$ decays
exponentially with respect to the geodesic distance $d_{xy}$ in the
holographic bulk

$C^{(\omega)}_{xy}\sim\exp(-\Delta d_{xy}),$ (43)

where the inverse correlation length $\Delta$ remains almost the same across
different layers in the bulk, as shown in Fig. 7(b). The short-ranged nature
of the Weyl field correlation is more obviously shown in Fig. 8, which is an
unequivocal sign of locality. In other words, the machine has learned from the
entanglement data that the Weyl field fluctuation can be described by a local
model (as the correlation is short-ranged) in the bulk. This emergent locality
is remarkable since locality was never explicitly given to the generative
model at the architecture level: the neural network in the generator
$G_{\vartheta}$ was fully connected, which in principle allows non-local /
long-ranged correlation of $\omega$ across the bulk, yet a short-ranged
correlation emerges from learning the entanglement data.

Further more, we can learn about the leading modes of gravitational
fluctuations in the machine-learned distribution $P_{\vartheta}[\omega]$ by
computing the spectral decomposition of the covariance function

$\Sigma^{(\omega)}_{xy}=\sum_{\alpha}\lambda^{(\alpha)}\omega^{(\alpha)}_{x}\omega^{(\alpha)}_{y},$
(44)

where $\lambda^{(\alpha)}$ is the $\alpha$th eigenvalue and
$\omega^{(\alpha)}$ is the corresponding eigenmode. The result is shown in
Fig. 9. The long wave-length collective fluctuations emerges as the leading
(low-energy) modes of gravitational fluctuation automatically. Using the
covariant function $\Sigma^{(\omega)}_{xy}$, one can reconstruct the effective
gravitational action to the quadratic order (at Gaussian level)

$S[\omega]=\frac{1}{2}\sum_{x,y}\omega_{x}\Sigma^{(\omega)}_{xy}\omega_{y}+\cdots,$
(45)

such that $P[\omega]\propto e^{-S[\omega]}$ approximately. In this way, the
machine-learning model helps us to extract a statistical gravity theory (in
terms of the Weyl field theory) from the entanglement data, demonstrating a
data-driven approach to establish the holographic duality.

###  IV.3 Matter Field Mass Renormalization Effect

As we have seen, geometric fluctuation effectively introduces interactions
between the bulk scalar field $\phi$, which generates the desired behavior for
mutual information. As a consequence, the bare mass $m$ of the scalar field
should also be renormalized by the gravitational interaction. Remarkably, we
can observe such a renormalization effect in our holographic model, by
comparing the static model (without geometric fluctuation) and the fluctuating
model (with geometric fluctuation). The mass parameter $m$ is trainable in
both models, but their optimal values are different due to the renormalization
effect. We train the static model on the single-region entanglement data, and
the fluctuating model on both the single- and double-region entanglement data.
For a range of total central charge $11/2\leq c\leq 8$ studied, we observe
that the trained value of the (bare) mass $m$ in the fluctuating model is
systematically larger than that $m_{0}$ in the static model, as shown in Fig.
10.

Figure 10: Trained values of the scalar field mass $m$ for different central
charge $c$ of the free fermion CFT, based on the static model (red) and the
fluctuating model (blue).

The mass renormalization effect can be understood heuristically by consider a
single-region entanglement. In the static model, the entanglement entropy is
modeled by $e^{-S_{A}}\sim e^{-m_{0}|\gamma_{A}|}$ where $\gamma_{A}$ is the
geodesic connecting the entanglement cuts of $A$ through the static bulk. With
geometric fluctuation $|\gamma_{A}|\to|\gamma_{A}|+\delta|\gamma_{A}|$ (where
$\delta|\gamma_{A}|$ is the additional geodesic length due to the Weyl field
$\omega$), the entanglement entropy will be modeled by

$e^{-S_{A}}\sim\mathop{\mathbb{E}}_{\omega}e^{-m(|\gamma_{A}|+\delta|\gamma_{A}|)}\simeq
e^{-m(|\gamma_{A}|-\frac{m}{2}\mathop{\mathbb{E}}_{\omega}(\delta|\gamma_{A}|)^{2})}.$
(46)

For these two models to match, we must have $m>m_{0}$, which qualitatively
explains our observation.

##  V Summary and Discussion

We present a machine-learning approach to extract the holography statistical
gravity theory from the data of multi-region entanglement entropy in a quantum
many-body system. Our work advances both the field of tensor network
holography and the field of machine learning holography. (i) On the tensor
network holography side, we generalize the random tensor network (RTN) model
to incorporate the bond dimension fluctuation, which makes the model more
expressive in capturing features of multi-region entanglement. We derive the
holographic bulk theory for the RTN with bond dimension fluctuation and show
that the dual gravity theory consists of a massive scalar field on a
fluctuating background geometry. The idea of using Ising duality in the
derivation is also quite original, which provides an alternative view of the
bulk theory that has not been presented in literature, as we are aware of.
(ii) On the machine learning holography side, our work goes beyond the
previous approaches[63, 26, 64, 65, 66, 67, 68] of inferring only a static
background geometry from the boundary quantum data. By modeling the bulk
geometric fluctuation with a generative model, our approach can extract a
statistical gravity theory from the quantum entanglement data. Remarkably, we
found that the machine-constructed gravity theory exhibit an emergent
locality, which reveals the hidden bulk locality behind the non-local quantum
entanglement on the boundary.

Our work provides a novel data-driven approach to explore the emergent gravity
from quantum entanglement. Combining with the recent development of efficient
numerical methods to simulate entanglement dynamics in quantum many-body
systems[47, 69, 48], we can further explore the corresponding gravity dynamics
in the holographic bulk, which will deepen our understanding of emergent
gravity from quantum entanglement. On the practical side, our algorithm will
boost the efficiency to model the entanglement structure of quantum many-body
systems, which will find applications in quantum algorithm optimization and
quantum circuit design.

###### Acknowledgements.

We acknowledge the helpful discussions with John McGreevy, Xiao-Liang Qi,
Zhenbin Yang. The authors are supported by a startup fund from UCSD and a UC
Hellman Fellowship.

## References

  * Brown and Henneaux [1986] J. D. Brown and M. Henneaux, Central charges in the canonical realization of asymptotic symmetries: an example from three-dimensional gravity, Comm. Math. Phys. 104, 207 (1986).
  * Witten [1998a] E. Witten, Anti-de Sitter space and holography, Advances in Theoretical and Mathematical Physics 2, 253 (1998a), arXiv:hep-th/9802150 [hep-th] .
  * Witten [1998b] E. Witten, Anti-de Sitter Space, Thermal Phase Transition, And Confinement In Gauge Theories, Adv. Theor. Math. Phys. 2, 505 (1998b), hep-th/9803131 .
  * Gubser  _et al._ [1998] S. S. Gubser, I. R. Klebanov, and A. M. Polyakov, Gauge theory correlators from non-critical string theory, Physics Letters B 428, 105 (1998), hep-th/9802109 .
  * Maldacena [1999] J. Maldacena, The Large-N Limit of Superconformal Field Theories and Supergravity, International Journal of Theoretical Physics 38, 1113 (1999), hep-th/9711200 .
  * Van Raamsdonk [2009] M. Van Raamsdonk, Comments on quantum gravity and entanglement, arXiv e-prints , arXiv:0907.2939 (2009), arXiv:0907.2939 [hep-th] .
  * van Raamsdonk [2010] M. van Raamsdonk, Building up spacetime with quantum entanglement, General Relativity and Gravitation 42, 2323 (2010), arXiv:1005.3035 [hep-th] .
  * Maldacena and Susskind [2013] J. Maldacena and L. Susskind, Cool horizons for entangled black holes, Fortschritte der Physik 61, 781 (2013), arXiv:1306.0533 [hep-th] .
  * Jensen and Karch [2013] K. Jensen and A. Karch, Holographic dual of an einstein-podolsky-rosen pair has a wormhole, Phys. Rev. Lett. 111, 211602 (2013).
  * Balasubramanian  _et al._ [2013] V. Balasubramanian, B. D. Chowdhury, B. Czech, and J. de Boer, The entropy of a hole in spacetime, Journal of High Energy Physics 2013, 220 (2013), arXiv:1305.0856 [hep-th] .
  * Qi [2013] X.-L. Qi, Exact holographic mapping and emergent space-time geometry, ArXiv e-prints (2013), arXiv:1309.6282 [hep-th] .
  * Balasubramanian  _et al._ [2014] V. Balasubramanian, P. Hayden, A. Maloney, D. Marolf, and S. F. Ross, Multiboundary Wormholes and Holographic Entanglement, arXiv e-prints , arXiv:1406.2663 (2014), arXiv:1406.2663 [hep-th] .
  * Susskind [2014] L. Susskind, ER=EPR, GHZ, and the Consistency of Quantum Measurements, arXiv e-prints , arXiv:1412.8483 (2014), arXiv:1412.8483 [hep-th] .
  * Balasubramanian  _et al._ [2014] V. Balasubramanian, B. D. Chowdhury, B. Czech, J. de Boer, and M. P. Heller, Bulk curves from boundary data in holography, Phys. Rev. D 89, 086004 (2014).
  * Czech and Lamprou [2014] B. Czech and L. Lamprou, Holographic definition of points and distances, Phys. Rev. D 90, 106005 (2014).
  * Cao  _et al._ [2017] C. Cao, S. M. Carroll, and S. Michalakis, Space from Hilbert space: Recovering geometry from bulk entanglement, Phys. Rev. D 95, 024031 (2017), arXiv:1606.08444 [hep-th] .
  * Ryu and Takayanagi [2006a] S. Ryu and T. Takayanagi, Holographic Derivation of Entanglement Entropy from the anti de Sitter Space/Conformal Field Theory Correspondence, Phys. Rev. Lett. 96, 181602 (2006a), arXiv:hep-th/0603001 [hep-th] .
  * Ryu and Takayanagi [2006b] S. Ryu and T. Takayanagi, Aspects of holographic entanglement entropy, Journal of High Energy Physics 2006, 045 (2006b), arXiv:hep-th/0605073 [hep-th] .
  * Porrati and Rabadan [2004] M. Porrati and R. Rabadan, Boundary Rigidity and Holography, Journal of High Energy Physics 2004, 034 (2004), arXiv:hep-th/0312039 [hep-th] .
  * Hammersley [2006] J. Hammersley, Extracting the bulk metric from boundary information in asymptotically AdS spacetimes, Journal of High Energy Physics 2006, 047 (2006), arXiv:hep-th/0609202 [hep-th] .
  * Bilson [2008] S. Bilson, Extracting spacetimes using the AdS/CFT conjecture, Journal of High Energy Physics 2008, 073 (2008), arXiv:0807.3695 [hep-th] .
  * Cao  _et al._ [2020] C. Cao, X.-L. Qi, B. Swingle, and E. Tang, Building bulk geometry from the tensor Radon transform, Journal of High Energy Physics 2020, 33 (2020), arXiv:2007.00004 [hep-th] .
  * Bilson [2011] S. Bilson, Extracting spacetimes using the AdS/CFT conjecture. Part II, Journal of High Energy Physics 2011, 50 (2011), arXiv:1012.1812 [hep-th] .
  * Alexakis  _et al._ [2017] S. Alexakis, T. Balehowsky, and A. Nachman, Determining a Riemannian Metric from Minimal Areas, arXiv e-prints , arXiv:1711.09379 (2017), arXiv:1711.09379 [math.DG] .
  * Bao  _et al._ [2019a] N. Bao, C. Cao, S. Fischetti, and C. Keeler, Towards bulk metric reconstruction from extremal area variations, Classical and Quantum Gravity 36, 185002 (2019a), arXiv:1904.04834 [hep-th] .
  * You  _et al._ [2018] Y.-Z. You, Z. Yang, and X.-L. Qi, Machine learning spatial geometry from entanglement features, Phys. Rev. B 97, 045153 (2018).
  * Roy and Sarkar [2018] S. R. Roy and D. Sarkar, Bulk metric reconstruction from boundary entanglement, Phys. Rev. D 98, 066017 (2018), arXiv:1801.07280 [hep-th] .
  * Headrick [2010] M. Headrick, Entanglement Rényi entropies in holographic theories, Phys. Rev. D 82, 126010 (2010), arXiv:1006.0047 [hep-th] .
  * Headrick [2019] M. Headrick, Lectures on entanglement entropy in field theory and holography, arXiv e-prints , arXiv:1907.08126 (2019), arXiv:1907.08126 [hep-th] .
  * Faulkner  _et al._ [2013] T. Faulkner, A. Lewkowycz, and J. Maldacena, Quantum corrections to holographic entanglement entropy, Journal of High Energy Physics 2013, 74 (2013), arXiv:1307.2892 [hep-th] .
  * Engelhardt and Wall [2015] N. Engelhardt and A. C. Wall, Quantum extremal surfaces: holographic entanglement entropy beyond the classical regime, Journal of High Energy Physics 2015, 73 (2015), arXiv:1408.3203 [hep-th] .
  * Dong  _et al._ [2020] X. Dong, X.-L. Qi, Z. Shangnan, and Z. Yang, Effective entropy of quantum fields coupled with gravity, Journal of High Energy Physics 2020, 52 (2020), arXiv:2007.02987 [hep-th] .
  * Hayden  _et al._ [2016] P. Hayden, S. Nezami, X.-L. Qi, N. Thomas, M. Walter, and Z. Yang, Holographic duality from random tensor networks, Journal of High Energy Physics 2016, 9 (2016), arXiv:1601.01694 [hep-th] .
  * Qi  _et al._ [2017] X.-L. Qi, Z. Yang, and Y.-Z. You, Holographic coherent states from random tensor networks, Journal of High Energy Physics 8, 60 (2017), arXiv:1703.06533 [hep-th] .
  * Swingle [2012a] B. Swingle, Constructing holographic spacetimes using entanglement renormalization, ArXiv e-prints (2012a), arXiv:1209.3304 [hep-th] .
  * Swingle [2012b] B. Swingle, Entanglement renormalization and holography, Phys. Rev. D 86, 065007 (2012b), arXiv:0905.1317 [cond-mat.str-el] .
  * Pastawski  _et al._ [2015] F. Pastawski, B. Yoshida, D. Harlow, and J. Preskill, Holographic quantum error-correcting codes: toy models for the bulk/boundary correspondence, Journal of High Energy Physics 2015, 149 (2015), arXiv:1503.06237 [hep-th] .
  * Salakhutdinov [2015] R. Salakhutdinov, Learning deep generative models, Annual Review of Statistics and Its Application 2, 361 (2015), https://doi.org/10.1146/annurev-statistics-010814-020120 .
  * Jimenez Rezende and Mohamed [2015] D. Jimenez Rezende and S. Mohamed, Variational Inference with Normalizing Flows, arXiv e-prints , arXiv:1505.05770 (2015), arXiv:1505.05770 [stat.ML] .
  * van den Oord  _et al._ [2016a] A. van den Oord, N. Kalchbrenner, and K. Kavukcuoglu, Pixel Recurrent Neural Networks, arXiv e-prints , arXiv:1601.06759 (2016a), arXiv:1601.06759 [cs.CV] .
  * Dinh  _et al._ [2016] L. Dinh, J. Sohl-Dickstein, and S. Bengio, Density estimation using Real NVP, arXiv e-prints , arXiv:1605.08803 (2016), arXiv:1605.08803 [cs.LG] .
  * Kingma  _et al._ [2016] D. P. Kingma, T. Salimans, R. Jozefowicz, X. Chen, I. Sutskever, and M. Welling, Improving Variational Inference with Inverse Autoregressive Flow, arXiv e-prints , arXiv:1606.04934 (2016), arXiv:1606.04934 [cs.LG] .
  * van den Oord  _et al._ [2016b] A. van den Oord, S. Dieleman, H. Zen, K. Simonyan, O. Vinyals, A. Graves, N. Kalchbrenner, A. Senior, and K. Kavukcuoglu, WaveNet: A Generative Model for Raw Audio, arXiv e-prints , arXiv:1609.03499 (2016b), arXiv:1609.03499 [cs.SD] .
  * Papamakarios  _et al._ [2017] G. Papamakarios, T. Pavlakou, and I. Murray, Masked Autoregressive Flow for Density Estimation, arXiv e-prints , arXiv:1705.07057 (2017), arXiv:1705.07057 [stat.ML] .
  * Vasseur  _et al._ [2018] R. Vasseur, A. C. Potter, Y.-Z. You, and A. W. W. Ludwig, Entanglement Transitions from Holographic Random Tensor Networks, arXiv e-prints , arXiv:1807.07082 (2018), arXiv:1807.07082 [cond-mat.stat-mech] .
  * Bao  _et al._ [2019b] N. Bao, G. Penington, J. Sorce, and A. C. Wall, Holographic Tensor Networks in Full AdS/CFT, arXiv e-prints , arXiv:1902.10157 (2019b), arXiv:1902.10157 [hep-th] .
  * Kuo  _et al._ [2019] W.-T. Kuo, A. A. Akhtar, D. P. Arovas, and Y.-Z. You, Markovian Entanglement Dynamics under Locally Scrambled Quantum Evolution, arXiv e-prints , arXiv:1910.11351 (2019), arXiv:1910.11351 [cond-mat.dis-nn] .
  * Fan  _et al._ [2021] R. Fan, S. Vijay, A. Vishwanath, and Y.-Z. You, Self-organized error correction in random unitary circuits with measurement, Phys. Rev. B 103, 174309 (2021), arXiv:2002.12385 [cond-mat.stat-mech] .
  * Hosur  _et al._ [2016] P. Hosur, X.-L. Qi, D. A. Roberts, and B. Yoshida, Chaos in quantum channels, Journal of High Energy Physics 2016, 4 (2016).
  * Seshadri  _et al._ [2018] A. Seshadri, V. Madhok, and A. Lakshminarayan, Tripartite mutual information, entanglement, and scrambling in permutation symmetric systems with an application to quantum chaos, arXiv e-prints , arXiv:1806.00113 (2018), arXiv:1806.00113 [quant-ph] .
  * Note [1] The model can be more expressive if the planar graph assumption is lifted, however it will be challenging to make connection to the RTN model on non-planar graphs.
  * Note [2] Another interpretation is to translate the fluctuating Ising coupling to the fluctuating mass term, as the scalar field mass is the relevant perturbation that drives the order-disorder transition, which plays the same role as the Ising coupling. This alternative view turns out to be equivalent to the fluctuating metric interpretation in two-dimension, as we will see soon.
  * Regge [1961] T. Regge, General relativity without coordinates, Nuovo Cim. 19, 558 (1961).
  * Goodfellow  _et al._ [2016] I. Goodfellow, Y. Bengio, and A. Courville,  _Deep Learning_ (MIT Press, 2016).
  * Goodfellow  _et al._ [2014] I. J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, Generative Adversarial Networks, arXiv e-prints , arXiv:1406.2661 (2014), arXiv:1406.2661 [stat.ML] .
  * Brydges  _et al._ [2019] T. Brydges, A. Elben, P. Jurcevic, B. Vermersch, C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F. Roos, Probing Rényi entanglement entropy via randomized measurements, Science 364, 260 (2019), arXiv:1806.05747 [quant-ph] .
  * Huang  _et al._ [2020] H.-Y. Huang, R. Kueng, and J. Preskill, Predicting many properties of a quantum system from very few measurements, Nature Physics 16, 1050 (2020), arXiv:2002.08953 [quant-ph] .
  * Agarap [2018] A. F. Agarap, Deep Learning using Rectified Linear Units (ReLU), arXiv e-prints , arXiv:1803.08375 (2018), arXiv:1803.08375 [cs.NE] .
  * Cohen and Welling [2016] T. Cohen and M. Welling, Group equivariant convolutional networks, in  _Proceedings of The 33rd International Conference on Machine Learning_, Proceedings of Machine Learning Research, Vol. 48, edited by M. F. Balcan and K. Q. Weinberger (PMLR, New York, New York, USA, 2016) pp. 2990–2999.
  * Note [3] We test the loss function by generating the data using a model with known parameters, and train new models with different loss function to see if the parameter converges to the known result.
  * Abadi  _et al._ [2015] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, TensorFlow: Large-scale machine learning on heterogeneous systems (2015), software available from tensorflow.org.
  * Kingma and Ba [2014] D. P. Kingma and J. Ba, Adam: A Method for Stochastic Optimization, ArXiv e-prints (2014), arXiv:1412.6980 [cs.LG] .
  * Gan and Shu [2017] W.-C. Gan and F.-W. Shu, Holography as deep learning, International Journal of Modern Physics D 26, 1743020 (2017), arXiv:1705.05750 [gr-qc] .
  * Dong and Zhou [2018] X. Dong and L. Zhou, Spacetime as the optimal generative network of quantum states: a roadmap to QM=GR?, arXiv e-prints , arXiv:1804.07908 (2018), arXiv:1804.07908 [quant-ph] .
  * Hashimoto  _et al._ [2018] K. Hashimoto, S. Sugishita, A. Tanaka, and A. Tomiya, Deep learning and the ads/cft correspondence, Phys. Rev. D 98, 046019 (2018).
  * Hashimoto  _et al._ [2018] K. Hashimoto, S. Sugishita, A. Tanaka, and A. Tomiya, Deep learning and holographic QCD, Phys. Rev. D 98, 106014 (2018), arXiv:1809.10536 [hep-th] .
  * Hashimoto [2019] K. Hashimoto, AdS/CFT as a deep Boltzmann machine, arXiv e-prints , arXiv:1903.04951 (2019), arXiv:1903.04951 [hep-th] .
  * Hashimoto  _et al._ [2020] K. Hashimoto, H.-Y. Hu, and Y.-Z. You, Neural ODE and Holographic QCD, arXiv e-prints , arXiv:2006.00712 (2020), arXiv:2006.00712 [hep-th] .
  * Akhtar and You [2020] A. A. Akhtar and Y.-Z. You, Multiregion entanglement in locally scrambled quantum dynamics, Phys. Rev. B 102, 134203 (2020), arXiv:2006.08797 [cond-mat.dis-nn] .

##  Appendix A Perturbative Analysis of Mutual and Tripartite Information

The random tensor network model points to a bulk theory described by the
following action

$S[\phi,\omega]=S[\phi|\omega]+S[\omega],$ (47)

where
$S[\phi|\omega]=\frac{1}{2}\sum_{x,y}\phi_{x}K^{(\phi)}_{xy}[\omega]\phi_{y}$
follows from Eq. (24), and we take a quadratic action
$S[\omega]=\frac{1}{2}\sum_{x,y}\omega_{x}K^{(\omega)}_{xy}\omega_{y}$ for
simplicity. In the perturbative limit, we assume that the fluctuation of the
field $\omega$ is small, such that we can expand
$K^{(\phi)}_{xy}[\omega]=K^{(\phi)}_{xy}+g\omega_{x}\delta_{xy}$, where
$K^{(\phi)}_{xy}$ denotes the bare kernel of $\phi$ on the $\omega=0$
background. Therefore, the bulk theory becomes

$S[\phi,\omega]=\frac{1}{2}\sum_{x,y}(\phi_{x}K^{(\phi)}_{xy}\phi_{y}+\omega_{x}K^{(\omega)}_{xy}\omega_{y})+\frac{g}{2}\sum_{x}\omega_{x}\phi_{x}^{2}.$
(48)

Define the field theory average as

$\langle\cdots\rangle=\frac{1}{Z}\int_{[\phi,\omega]}\cdots
e^{-S[\phi,\omega]},$ (49)

then the entanglement entropy $S_{A}$ of a signle-region $A$ that ends at the
dual sites $(x_{1},x_{2})$ can be written as

$e^{-S_{A}}=\langle\phi_{x_{1}}\phi_{x_{2}}\rangle.$ (50)

The entanglement entropy for multi-regions are modeled similarly as multi-
point covariance of the $\phi$ field among all boundary points.

Now we consider three regions $A$, $B$ and $C$ boundaried by $(x_{1},x_{2})$,
$(x_{3},x_{4})$ and $(x_{5},x_{6})$ respectively. The mutual information
$I_{A:B}=S_{A}+S_{B}-S_{AB}$ and the tripartite information
$I_{A:B:C}=S_{A}+S_{B}+S_{C}-S_{AB}-S_{BC}-S_{AC}+S_{ABC}$ can be evaluated by
the following ratios of covariance function

$\begin{split}e^{-I_{A:B}}&=\frac{\langle\phi_{x_{1}}\phi_{x_{2}}\rangle\langle\phi_{x_{3}}\phi_{x_{4}}\rangle}{\langle\phi_{x_{1}}\phi_{x_{2}}\phi_{x_{3}}\phi_{x_{4}}\rangle},\\\
e^{-I_{A:B:C}}&=\frac{\langle\phi_{x_{1}}\phi_{x_{2}}\rangle\langle\phi_{x_{3}}\phi_{x_{4}}\rangle\langle\phi_{x_{5}}\phi_{x_{6}}\rangle\langle\phi_{x_{1}}\phi_{x_{2}}\phi_{x_{3}}\phi_{x_{4}}\phi_{x_{5}}\phi_{x_{6}}\rangle}{\langle\phi_{x_{1}}\phi_{x_{2}}\phi_{x_{3}}\phi_{x_{4}}\rangle\langle\phi_{x_{1}}\phi_{x_{2}}\phi_{x_{3}}\phi_{x_{4}}\rangle\langle\phi_{x_{1}}\phi_{x_{2}}\phi_{x_{3}}\phi_{x_{4}}\rangle}.\end{split}$
(51)

Figure 11: A few useful graphs.

To simplify the notation in the following discussion, we introduce a few
graphs in Fig. 11. Let $\mathcal{G}_{4}$ be the complete graph over
$x_{1},x_{2},x_{3},x_{4}$, and $\mathcal{G}_{6}$ be the complete graph over
$x_{1},x_{2},x_{3},x_{4},x_{5}$. Further denote $\mathcal{G}^{\prime}_{2n}$
graph (with a prime) to be the graph with edges $(x_{2k-1}x_{2k})$ (for
$k=1,\cdots,n$) removed from $\mathcal{G}_{2n}$. Define the set of perfect
matchings on a graph $\mathcal{G}$ by $\mathcal{M}[\mathcal{G}]$ (where each
perfect matching is a subset of edges such that every vertex is covered and
only covered by one edge). We define the bare propagators (the covariance
functions) $\Sigma^{(\phi)}$ and $\Sigma^{(\omega)}$ from the inverses of the
bare kernels $K^{(\phi)}$ and $K^{(\omega)}$ for both $\phi$ and $\omega$
fields respectively,

$\Sigma^{(\phi)}_{xy}=((K^{(\phi)})^{-1})_{xy},\quad\Sigma^{(\omega)}_{xy}=((K^{(\omega)})^{-1})_{xy}.$
(52)

Using perturbative field theory (treating $g$ in Eq. (47) as perturbation), to
the 2nd order in $g$ (and keeping only the tree level diagrams), we can
calculate the covariance functions

$\begin{split}\langle\phi_{x_{1}}\phi_{x_{2}}\rangle&=\Sigma^{(\phi)}_{x_{1}x_{2}},\end{split}$
(53)
$\begin{split}\langle\phi_{x_{1}}\phi_{x_{2}}\phi_{x_{3}}\phi_{x_{4}}\rangle&=\sum_{(x_{i}x_{j})(x_{k}x_{l})\in\mathcal{M}[\mathcal{G}_{4}]}\bigg{(}\Sigma^{(\phi)}_{x_{i}x_{j}}\Sigma^{(\phi)}_{x_{k}x_{l}}+g^{2}\sum_{y_{1},y_{2}}\Sigma^{(\omega)}_{y_{1}y_{2}}\Sigma^{(\phi)}_{x_{i}y_{1}}\Sigma^{(\phi)}_{y_{1}x_{j}}\Sigma^{(\phi)}_{x_{k}y_{2}}\Sigma^{(\phi)}_{y_{2}x_{l}}\bigg{)},\end{split}$
(54)
$\begin{split}\langle\phi_{x_{1}}\phi_{x_{2}}\phi_{x_{3}}\phi_{x_{4}}\phi_{x_{5}}\phi_{x_{6}}\rangle&=\sum_{(x_{i}x_{j})(x_{k}x_{l})(x_{m}x_{n})\in\mathcal{M}[\mathcal{G}_{6}]}\bigg{(}\Sigma^{(\phi)}_{x_{i}x_{j}}\Sigma^{(\phi)}_{x_{k}x_{l}}\Sigma^{(\phi)}_{x_{m}x_{n}}\\\
&\hskip
120.0pt+g^{2}\sum_{y_{1},y_{2}}\Sigma^{(\omega)}_{y_{1}y_{2}}(\Sigma^{(\phi)}_{x_{i}y_{1}}\Sigma^{(\phi)}_{y_{1}x_{j}}\Sigma^{(\phi)}_{x_{k}y_{2}}\Sigma^{(\phi)}_{y_{2}x_{l}}\Sigma^{(\phi)}_{x_{m}x_{n}}\\\
&\hskip
180.0pt+\Sigma^{(\phi)}_{x_{i}y_{1}}\Sigma^{(\phi)}_{y_{1}x_{j}}\Sigma^{(\phi)}_{x_{k}x_{l}}\Sigma^{(\phi)}_{x_{m}y_{2}}\Sigma^{(\phi)}_{y_{2}x_{n}}\\\
&\hskip
180.0pt+\Sigma^{(\phi)}_{x_{i}x_{j}}\Sigma^{(\phi)}_{x_{k}y_{1}}\Sigma^{(\phi)}_{y_{1}x_{l}}\Sigma^{(\phi)}_{x_{m}y_{2}}\Sigma^{(\phi)}_{y_{2}x_{n}})\bigg{)}.\end{split}$
(55)

Substitute the correlation functions Eq. (53)-Eq. (55) to Eq. (51), we find
(to the $g^{2}$ order)

$\begin{split}I_{A:B}=&\sum_{(x_{i}x_{j})(x_{k}x_{l})\in\mathcal{M}[\mathcal{G}^{\prime}_{4}]}\frac{\Sigma^{(\phi)}_{x_{i}x_{j}}\Sigma^{(\phi)}_{x_{k}x_{l}}}{\Sigma^{(\phi)}_{x_{1}x_{2}}\Sigma^{(\phi)}_{x_{3}x_{4}}}\\\
&+g^{2}\sum_{(x_{i}x_{j})(x_{k}x_{l})\in\mathcal{M}[\mathcal{G}_{4}]}\sum_{y_{1},y_{2}}\Sigma^{(\omega)}_{y_{1}y_{2}}\frac{\Sigma^{(\phi)}_{x_{i}y_{1}}\Sigma^{(\phi)}_{y_{1}x_{j}}\Sigma^{(\phi)}_{x_{k}y_{2}}\Sigma^{(\phi)}_{y_{2}x_{l}}}{\Sigma^{(\phi)}_{x_{1}x_{2}}\Sigma^{(\phi)}_{x_{3}x_{4}}},\end{split}$
(56)
$\begin{split}I_{A:B:C}=&-\sum_{(x_{i}x_{j})(x_{k}x_{l})(x_{m}x_{n})\in\mathcal{M}[\mathcal{G}^{\prime}_{6}]}\frac{\Sigma^{(\phi)}_{x_{i}x_{j}}\Sigma^{(\phi)}_{x_{k}x_{l}}\Sigma^{(\phi)}_{x_{m}x_{n}}}{\Sigma^{(\phi)}_{x_{1}x_{2}}\Sigma^{(\phi)}_{x_{3}x_{4}}\Sigma^{(\phi)}_{x_{5}x_{6}}}\\\
&-g^{2}\sum_{(x_{i}x_{j})(x_{k}x_{l})(x_{m}x_{n})\in\mathcal{M}[\mathcal{G}^{\prime}_{6}]}\sum_{y_{1},y_{2}}\Sigma^{(\omega)}_{y_{1}y_{2}}\Big{(}\frac{\Sigma^{(\phi)}_{x_{i}y_{1}}\Sigma^{(\phi)}_{y_{1}x_{j}}\Sigma^{(\phi)}_{x_{k}y_{2}}\Sigma^{(\phi)}_{y_{2}x_{l}}\Sigma^{(\phi)}_{x_{m}x_{n}}}{\Sigma^{(\phi)}_{x_{1}x_{2}}\Sigma^{(\phi)}_{x_{3}x_{4}}\Sigma^{(\phi)}_{x_{5}x_{6}}}+\frac{\Sigma^{(\phi)}_{x_{i}y_{1}}\Sigma^{(\phi)}_{y_{1}x_{k}}\Sigma^{(\phi)}_{x_{j}y_{2}}\Sigma^{(\phi)}_{y_{2}x_{l}}\Sigma^{(\phi)}_{x_{m}x_{n}}}{2\Sigma^{(\phi)}_{x_{1}x_{2}}\Sigma^{(\phi)}_{x_{3}x_{4}}\Sigma^{(\phi)}_{x_{5}x_{6}}}\\\
&\hskip
170.0pt+\frac{\Sigma^{(\phi)}_{x_{i}y_{1}}\Sigma^{(\phi)}_{y_{1}x_{j}}\Sigma^{(\phi)}_{x_{k}x_{l}}\Sigma^{(\phi)}_{x_{m}y_{2}}\Sigma^{(\phi)}_{y_{2}x_{n}}}{\Sigma^{(\phi)}_{x_{1}x_{2}}\Sigma^{(\phi)}_{x_{3}x_{4}}\Sigma^{(\phi)}_{x_{5}x_{6}}}+\frac{\Sigma^{(\phi)}_{x_{i}y_{1}}\Sigma^{(\phi)}_{y_{1}x_{m}}\Sigma^{(\phi)}_{x_{k}x_{l}}\Sigma^{(\phi)}_{x_{j}y_{2}}\Sigma^{(\phi)}_{y_{2}x_{n}}}{2\Sigma^{(\phi)}_{x_{1}x_{2}}\Sigma^{(\phi)}_{x_{3}x_{4}}\Sigma^{(\phi)}_{x_{5}x_{6}}}\\\
&\hskip
170.0pt+\frac{\Sigma^{(\phi)}_{x_{i}x_{j}}\Sigma^{(\phi)}_{x_{k}y_{1}}\Sigma^{(\phi)}_{y_{1}x_{l}}\Sigma^{(\phi)}_{x_{m}y_{2}}\Sigma^{(\phi)}_{y_{2}x_{n}}}{\Sigma^{(\phi)}_{x_{1}x_{2}}\Sigma^{(\phi)}_{x_{3}x_{4}}\Sigma^{(\phi)}_{x_{5}x_{6}}}+\frac{\Sigma^{(\phi)}_{x_{i}x_{j}}\Sigma^{(\phi)}_{x_{k}y_{1}}\Sigma^{(\phi)}_{y_{1}x_{m}}\Sigma^{(\phi)}_{x_{l}y_{2}}\Sigma^{(\phi)}_{y_{2}x_{n}}}{2\Sigma^{(\phi)}_{x_{1}x_{2}}\Sigma^{(\phi)}_{x_{3}x_{4}}\Sigma^{(\phi)}_{x_{5}x_{6}}}\Big{)}.\end{split}$
(57)

In the case that $\Sigma^{(\phi)}_{xy}\geq 0$ and $\Sigma^{(\omega)}_{xy}\geq
0$ (which is typically the case), we can ensure $I_{A:B}\geq 0$ and
$I_{A:B:C}\leq 0$. The result proves that random tensor network model can
produce a negative tripartite information $I_{A:B:C}$, which is a unique
feature of quantum many-body entanglement that can not be achieved in
classical systems. A negative tripartite information is an indication of
quantum information scrambling and chaotic quantum dynamics in the quantum
system. Although the bulk theory is a classical statistical gravity model, it
can still model the quantum chaotic entanglement features on the holographic
boundary, which speak for the strong expression power of the random tensor
network model.

##  Appendix B Necessity and Expected Behavior of Weyl Field Fluctuation

We would like to take a closer look at the mutual information. It would more
intuitive to present the diagrammatic representation of Eq. (56)

$I_{A:B}=\raisebox{-34.0pt}{\includegraphics[height=72.0pt]{fig_MIdiag}}.$
(58)

where points on the boundary correspond to $x_{1},x_{2},x_{3},x_{4}$
(following the arrangement of vertices in the $\mathcal{G}_{4}$ graph shown in
Fig. 11) and the small circles in the bulk correspond to $y_{1},y_{2}$ that
should be summed over. The black lines represent $\Sigma^{(\phi)}_{xy}$ and
the gray lines represent $\Sigma^{(\omega)}_{xy}$. The perturbation $g$
parameterizes the coupling strength of the bulk scalar field $\phi$ to the
background gravitational fluctuation (the Weyl field $\omega$). Setting $g=0$
will decouple the gravitational fluctuation, which effectively corresponds to
a static bulk model (because gravitational fluctuation will have no effect in
the decoupled limit). Let us consider the case when regions $A$ and $B$ are
far separated, meaning that the spacings $|x_{1}-x_{2}|$ and $|x_{3}-x_{4}|$
are small. In this case, we expect the mutual information to decay with the
inter-region spacing in a power-law manner with the power set by the smallest
scaling dimension of the critical field in the quantum system on the
holographic boundary, because the mutual information upper bounds all
correlation function between regions $A$ and $B$, which can not decay faster
than the lightest critical field. As we will see, this behavior can only be
reproduce via the bulk model gravitational fluctuation is included.

To argue the necessity of including gravitational fluctuation, we first
consider the decoupled limit (i.e. $g=0$) to demonstrate why it fails to
capture the correct behavior of mutual information. In the $g=0$ limit, the
first two terms can still contribute to a finite mutual information that
decays with the inter-region separation in a power-law manner, but the power
will be set by the central charge of the quantum system on the holographic
boundary. Because the power-law comes from the $\phi$-field correlation
$\Sigma^{(\phi)}_{xy}$, whose scaling is determined by the single-region
entanglement entropy, as $-\ln\Sigma^{(\phi)}_{x_{1}x_{2}}\sim S_{A}\sim
c\ln|A|$ (such that $x_{1},x_{2}$ are boundary points of the region $A$).
However, the total central charge $c=N/2$ can be as large as we wish in the
large $N$ limit. Therefore, although the first two terms (the $g^{0}$ terms)
in Eq. (58) can produce a power-law decay mutual information, but the power
will typically be too large (i.e. the mutual information will decay too fast).
This reflects the internal inconsistency in describing both single- and
double-region entanglements using a holographic bulk model without
gravitational fluctuation.

An obvious solution is to introduce a different field from $\phi$ to mediate
the mutual information across the holographic bulk. Then we will have an
independent freedom, such that we can tune its scaling dimension to match that
of the lightest critical field. This is one major motivation to introduce the
gravitational fluctuation (or to couple the scaler field $\phi$ to a
fluctuating Weyl field $\omega$). As we turn on the coupling $g$, the mutual
information will be dominated by

$I_{A:B}=g^{2}\sum_{y_{1},y_{2}}\Sigma^{(\omega)}_{y_{1}y_{2}}\frac{\Sigma^{(\phi)}_{x_{1}y_{1}}\Sigma^{(\phi)}_{y_{1}x_{2}}\Sigma^{(\phi)}_{x_{3}y_{2}}\Sigma^{(\phi)}_{y_{2}x_{4}}}{\Sigma^{(\phi)}_{x_{1}x_{2}}\Sigma^{(\phi)}_{x_{3}x_{4}}},$
(59)

whose long range behavior scales with $\sim
g^{2}\Sigma^{(\omega)}_{y_{1}y_{2}}$. The Weyl field $\omega$ has a different
propagator $\Sigma^{(\omega)}_{xy}$, which can be independently tuned to make
$-\ln\Sigma^{(\omega)}_{xy}\sim 2\Delta_{\text{min}}\ln|x-y|$ with
$\Delta_{\text{min}}$ being the smallest scaling dimension in the quantum
critical theory. Here $|x-y|$ denotes the distance between two boundary points
$x$ and $y$ measured using the bound metric. Translate $|x-y|$ into the bulk
distance $d_{xy}$ assuming the bulk has a hyperbolic background geometry, we
have $d_{xy}\sim\ln|x-y|$, which implies $\Sigma^{(\omega)}_{xy}\sim
e^{-2\Delta_{\text{min}}d_{xy}}$. This indicates that Weyl field must be heavy
in the bulk to produce the exponential decay of its correlation function with
the bulk distance. Indeed, such a massive Weyl field fluctuation does emerge
in the machine-learnt bulk gravity theory.