Holographic duality from random tensor networks

Tensor networks provide a natural framework for exploring holographic duality because they obey entanglement area laws. They have been used to construct explicit toy models realizing many interesting structural features of the AdS/CFT correspondence, including the non-uniqueness of bulk operator reconstruction in the boundary theory. In this article, we explore the holographic properties of networks of random tensors. We find that our models naturally incorporate many features that are analogous to those of the AdS/CFT correspondence. When the bond dimension of the tensors is large, we show that the entanglement entropy of boundary regions, whether connected or not, obey the Ryu-Takayanagi entropy formula, a fact closely related to known properties of the multipartite entanglement of assistance. Moreover, we find that each boundary region faithfully encodes the physics of the entire bulk entanglement wedge. Our method is to interpret the average over random tensors as the partition function of a classical ferromagnetic Ising model, so that the minimal surfaces of Ryu-Takayanagi appear as domain walls. Upon including the analog of a bulk field, we find that our model reproduces the expected corrections to the Ryu-Takayanagi formula: the minimal surface is displaced and the entropy is augmented by the entanglement of the bulk field. Increasing the entanglement of the bulk field ultimately changes the minimal surface topologically in a way similar to creation of a black hole. Extrapolating bulk correlation functions to the boundary permits the calculation of the scaling dimensions of boundary operators, which exhibit a large gap between a small number of low-dimension operators and the rest. While we are primarily motivated by AdS/CFT duality, our main results define a more general form of bulk-boundary correspondence which could be useful for extending holography to other spacetimes.


Introduction
Tensor networks have been proposed [1] as a helpful tool for understanding holographic duality [2][3][4] due to the intuition that the entropy of a tensor network is bounded by an area law that agrees with the Ryu-Takayanagi (RT) formula [5]. In general, the area law only gives an upper bound to the entropy [1], which for particular tensor networks and choices of regions has been shown to be saturated [6]. Tensor networks can also be used to build holographic mappings or holographic codes [6][7][8], which are isometries between the Hilbert space of the bulk and that of the boundary. In particular, some of us have recently proposed bidirectional holographic codes built from tensors with particular properties, socalled pluperfect tensors [8]. These codes simultaneously satisfies several desired properties, including the RT formula for a subset of boundary states, error correction properties of bulk local operators [9], a kind of bulk gauge invariance, and the possibility of sub-AdS locality. The perfect and pluperfect tensors defined in Refs. [6] and [8], respectively, have entanglement properties that are idealized version of large-dimensional random tensors, which is part of the motivation why it is natural to study these tensor networks. In this work, we will show that by directly studying networks of large dimensional random tensors, instead of their "idealized" counterpart, their properties can be computed more systematically. Specifically, we will assume that each tensor in the network is chosen independently at random. We find that the computation of typical Rényi entropies and other quantities of interest in the corresponding tensor network states can be mapped to the evaluation of partition functions of classical statistical models, namely generalized Ising models with boundary pinning fields. When each leg of each tensor in the network has dimension D, these statistical models have inverse temperature β ∝ log D. For large enough D, they are in the long-range ordered phase, and we find that the entropies of a boundary region is directly related to the energy of a domain wall between different domains of the order parameter. The minimal energy condition for this domain wall naturally leads to the RT formula. 1 Besides yielding the RT formula for general boundary subsystems, the technique of random state averaging allows us to study many further properties of a random tensor network: 1. Effects of bulk entanglement. Using the random tensor network as a holographic mapping rather than a state on the boundary, we derive a formula for the entropy of a boundary region in the presence of an entangled state in the bulk. As a special example of the effect of bulk entanglement, we show how the behavior of minimal surfaces (which are minimal energy domain walls in the statistical model) is changed qualitatively by introducing a highly entangled state in the bulk. When the state is sufficiently highly entangled, no minimal surface penetrates into this region, so that the topology of the space has effectively changed. This phenomenon is analogous to the change of spatial geometry in the Hawking-Page transition [11,12], where the bulk geometry changes from perturbed AdS to a black hole upon increasing temperature.
2. Bidirectional holographic code and code subspace. By calculating the entanglement entropy between a bulk region and the boundary in a given tensor network, we can verify that the random tensor network defines a bidirectional holographic code (BHC). When the bulk Hilbert space has a higher dimension than the boundary, we obtain an approximate isometry from the boundary to the bulk. When we restrict the bulk degrees of freedom to a smaller subspace ("code subspace", or "low energy subspace") which has dimension lower than the boundary Hilbert space dimension, we also obtain an approximate isometry from this bulk subspace to the boundary. This bulk-to-boundary isometry satisfies the error correction properties defined in Ref. [9]. To be more precise, all bulk local operators in the entanglement wedge of a boundary region can be recovered from that boundary region. 2 3. Correlation spectrum. In addition to entanglement entropies, we can also study properties of boundary multi-point functions. In particular, we show that the boundary two-point functions are determined by the bulk two-point functions and the properties of the statistical model. When the bulk geometry is a pure hyperbolic space, the boundary two-point correlations have power-law decay, which defines the scaling dimension spectrum. We show that in large-dimensional random tensor networks there are two kinds of scaling dimensions, those from the bulk "low energy" theory which do not grow with the bond dimension D, and those from the tensor network itself which grow ∝ log D. This confirms that the holographic mapping defined by a random tensor network maps a weakly-interacting bulk state to a boundary state with a scaling dimension gap, consistent with the expectations of AdS/CFT.
The use of random matrix techniques has a long and rich history in quantum information theory (see, e.g., the recent review [13] and references therein). Previous work on random tensor network states has originated from a diverse set of motivations, including the construction of novel random ensembles that satisfy a generalized area law [14,15], the relationship between entropy and the decay of correlations [16], and the maximum entropy principle [17]. The relation between the Schmidt ranks of tensor network states and minimal cuts through the network has been investigated in [18]. While the primary motivation for this work is to better understand holographic duality, its methods and even the nature of many of its conclusions place it squarely in this earlier tradition. In the holographic context, it was in fact previously shown that using a class of pseudo-random tensors known as quantum expanders in a MERA tensor network would reproduce the qualitative scaling of the Ryu-Takayanagi formula [19].
The remainder of the paper is organized as follows. In Section 2 we define the random tensor networks. We show how the calculation of the second Rényi entropy is mapped to the partition function of a classical Ising model. In Section 3 we investigate the RT formula in the large dimension limit of the random tensors, and discuss the effect of bulk entanglement.
As an explicit example we study the minimal surfaces for a highly entangled (volume-law) bulk state and discuss the transition of the effective bulk geometry as a function of bulk entropy density. In Section 4 we study the properties of the holographic mapping defined by random tensor networks, including boundary-to-bulk isometries and bulk-to-boundary isometries for the code subspace, and we discuss the recovery of bulk operators from boundary regions. In Section 5 we generalize the calculation of the second Rényi entropy to higher Rényi entropies. We show that the n-th Rényi entropy calculation is mapped to the partition function of a statistical model with a Sym n permutation group element at each vertex. The same technique also enables us to compute other averaged quantities involving higher powers of the density operator. In Section 6 we use this technique to study the boundary two-point correlation functions. We show that the boundary correlation functions are determined by the bulk correlations and the tensor network, and that a gap in the scaling dimensions opens at large D in the case of AdS geometry. In Section 7 we bound the fluctuations around the typical values calculated previously and discuss the effect of finite bond dimensions. Section 8 explains the close relationship between the random tensors networks of this paper and optimal multipartite entanglement distillation protocols previously studied in the quantum information theory literature. In Section 9 we consider other ensembles of random states. We find that the RT formula can be exactly satisfied in tensor networks built from random stabilizer states, which allows for the construction of exact holographic codes. Finally, Section 10 is devoted to conclusion and discussion.

Definition of random tensor networks
We start by defining the most general tensor network states in a language that is suitable for our later discussion. A rank-n tensor has components T µ 1 µ 2 ...µn with µ k = 1, 2, . . . , D k . We can define a Hilbert space H k with dimension D k for each leg of the tensor, and consider the index µ k as labeling a complete basis of states |µ k in this Hilbert space. In this language, T µ 1 µ 2 ...µn (with proper normalization) corresponds to the wavefunction of a quantum state |T = {µ k } T µ 1 µ 2 ...µn |µ 1 ⊗ |µ 2 ⊗ · · · ⊗ |µ n defined in the product Hilbert space n k=1 H k . A tensor network is obtained by connecting tensors, i.e., by contracting a common index. For purposes of illustration, a small tensor network is shown in Fig. 1 (a). Before connecting the tensors, each tensor corresponds to a quantum state, so that the collection of all tensors can be considered as a tensor product state x |V x . Here, x denotes all vertices in the network, and |V x is the state corresponding to the tensor at vertex x. Each leg of a tensor corresponds to a Hilbert space. We will denote the Hilbert space corresponding to a leg from x to another vertex y by H xy , and its dimension by D xy . If a leg is dangling, i.e., not connected to any other vertex, we will denote the corresponding Hilbert space by H x∂ and its dimension by D x∂ . (Without loss of generality we can assume there is at most A tensor network that defines a mapping from bulk legs (red) to boundary legs (blue). An arbitrary bulk state (orange triangle) is mapped to a boundary state. (For simplicity, we have drawn a pure state in the bulk. For a mixed state the map needs to be applied to both indices of the bulk density operator.) (c) The internal lines of the tensor network can always be combined with the bulk state and viewed as a state in an enlarged Hilbert space (enclosed by the dashed hexegon). In this view, each tensor acts independently on this generalized bulk state and maps it to the boundary state.
one dangling leg at each vertex.) Connecting two tensors at x, y by an internal line then corresponds to a projection in the Hilbert space H xy ⊗ H yx onto a maximally entangled state |xy = 1 √ Dxy Dxy µ=1 |µ xy ⊗ |µ yx . Here |µ xy denotes a state in the Hilbert space H xy and similarly for |µ yx . By connecting the tensors according to the internal lines of the tensor network, we thus obtain the state (2.1) in the Hilbert space corresponding to the dangling legs, x∈∂ H x∂ . We note that |Ψ is in general not normalized. Tensor network states defined in this way are often referred to as projected entangled pair states (PEPS) [20]. As has been discussed in previous works [6][7][8], tensor networks can be used to define not only quantum states but also holographic mappings, or holographic codes, which map between the Hilbert space of the bulk and that of the boundary. Fig. 1 (b) shows a very simple "holographic mapping" which maps the bulk indices (red lines) to boundary indices (blue lines), with internal lines (black lines) contracted. A bulk state (orange triangle in the figure) is mapped to a boundary state by this mapping. Such a boundary state can also be written in a form similar to Eq. (2.1). Instead of viewing the tensor network as defining a mapping, we can equivalently consider it as a quantum state in the Hilbert space H b ⊗ H ∂ , which is a direct product of the bulk Hilbert space H b and the boundary Hilbert space H ∂ .
Denoting the bulk state as |Φ b , the corresponding boundary state is From this expression one can see that the internal lines of the tensor network can actually be viewed as part of the bulk state. As is illustrated in Fig. 1 (c), one can view the maximally entangled states on internal lines together with the bulk state |Φ b as a state in the enlarged "bulk Hilbert space". This point of view will be helpful for our discussion. More generally, one can also have a mixed bulk state with density operator ρ b , instead of the pure state |Φ b . The most generic form of the boundary state is given by the density operator Here the partial trace tr P is carried over the bulk and internal legs of all tensors (i.e., over all but the dangling legs). In this compact form, one can see that the state ρ is a linear function of the independent pure states of each tensor |V x V x |.
In this work, we study tensor network states of the form (2.3), where the tensors |V x are unit vectors chosen independently at random from their respective Hilbert spaces. We will mostly use the "uniform" probability measure that is invariant under arbitrary unitary transformations. Equivalently, we can take an arbitrary reference state |0 x and define |V x = U |0 x with U a unitary operator. The random average of an arbitrary function f (|V x ) of the state |V x is then equivalent to an integration over U according to the Haar probability measure dU f (U |0 x ), with normalization dU = 1.
All nontrivial entanglement properties of such a tensor network state are induced by the projection, i.e., the partial trace with ρ P . However, the average over random tensors can be carried out before taking the partial trace, since the latter is a linear operation. This is the key insight that enables the computation of entanglement properties of random tensor networks.

Calculation of the second Rényi entropy
We will now apply this technique to study the second Rényi entropies of the random tensor network state ρ defined in Eq. (2.3). For a boundary region A with reduced density matrix ρ A , the second Rényi entropy S 2 (A) is given by e −S 2 (A) = tr ρ 2 A /(tr ρ) 2 . 3 It is helpful to write this expression in a different form by using the "swap trick", (2.5) 3 In the quantum information theory literature, the Rényi entropy is usually defined with logarithm in base 2, Sn Here we use base e to keep the notation consistent with the condensed matter and high energy literature.
Here we have defined a direct product ρ ⊗ ρ of two copies of the original system, and the operator F A is defined on this two-copy system and swaps the states of the two copies in the region A. To be more precise, its action on a basis state of the two-copy Hilbert space is given by F A (|n A 1 ⊗ |mĀ 1 ⊗ |n A 2 ⊗ |m Ā 2 ) = |n A 1 ⊗ |mĀ 1 ⊗ |n A 2 ⊗ |m Ā 2 , whereĀ denotes the complement of A on the boundary.
We are now interested in the typical values of the entropy. Denote the numerator and denominator resp. of Eq. (2.5) by These are both functions of the random states |V x at each vertex. We would like to average over all states in the single-vertex Hilbert space. The variables Z 1 and Z 0 are easier to average than the entropy, since they are quadratic functions of the single-site density matrix |V x V x |. The entropy average can be expanded in powers of the fluctuations δZ 1 = Z 1 − Z 1 and δZ 0 = Z 0 − Z 0 : (2.8) We will later show in Section 7 that for large enough bond dimensions D xy the fluctuations are suppressed. Thus we can approximate the entropy with high probability by the separate averages of Z 1 and Z 0 : Throughout this article we use for asymptotic equality as the bond dimensions go to infinity. In the following we will compute Z 1 and Z 0 separately and use (2.9) to determine the typical entropy. To compute Z 1 , we insert Eq. (2.3) into Eq. (2.6) and obtain (2.10) In this expression we have combined the partial trace over bulk indices in the definition of the boundary state ρ and the trace over the boundary indices in Eq. (2.6) into a single trace over all indices. In the expression it is now transparent that the average over states, one at each vertex, can be carried out independently before couplings between different sites are introduced by the projection. The average over states can be done by taking an arbitrary reference state |0 x and setting |V x = U x |0 x . Then the average is equivalent to an integration over U x ∈ SU (D x ) with respect to the Haar measure. The result of this integration can be obtained using Schur's lemma (see, e.g., Ref. [21]): Here, I x denotes the identity operator and F x the swap operator defined in the same way as F A described above, swapping the two copies of Hilbert space of the vertex x (which means all legs connecting to x). The Hilbert space dimension is D x = y n.n. x D xy , the product of the dimensions corresponding to all legs adjacent to x, including the boundary dangling legs. It is helpful to represent Eq. (2.11) graphically as in Fig. 2 (a) and (b). Carrying out the average over states at each vertex x, Z 1 then consists of 2 N terms if there are N vertices, with an identity operator or swap operator at each vertex. We can then introduce an Ising spin variable s x = ±1, and use s x = 1 (s x = −1) to denote the choice of I x and F x , respectively. In this representation, Z 1 becomes a partition function of the spins {s x }: For each value of the Ising variables {s x }, the operator being traced is now completely factorized into a product of terms, since F x acts on each leg of the tensor independently. This fact is illustrated in Fig. 2 Then the trace at a boundary leg x∂ gives D . Taking a product of these two kinds of terms in the trace, we obtain the Ising action The form of the action can be further simplified by recalling that ρ P has the direct product form in Eq. (2.4). Therefore the second Rényi entropy factorizes into that of the bulk state ρ b and that of the maximally entangled states at each internal line xy. The latter is a standard Ising interaction term, since the entropy of either site is log D xy while the entropy of the two sites together vanishes. Therefore (2.13) Here we have omitted the details of the constant term since it plays no role in later discussions. Eq. (2.13) is the foundation of our later discussion. The same derivation applies to the average of the denominator Z 0 = tr [ρ ⊗ ρ] in Eq. (2.5), which leads to the same Ising partition function with a different boundary condition h x = 1 for all boundary sites, since  there is no swap operator F A applied. One can define F 1 = − log Z 1 , F 0 = − log Z 0 , such that F 1 and F 0 are the free energy of the Ising model with different boundary conditions. 4 Then Eq. (2.9) reads That is, the typical second Rényi entropy is given by the difference of the two free energies, i.e., the "energy cost" induced by flipping the boundary pinning field to down (−1) in region A, while keeping the remainder of the system with a pinning field up (+1).
In summary, what we have achieved is that the second Rényi entropy is related to the partition function of a classical Ising model defined on the same graph as the tensor network. Besides the standard two-spin interaction term, the Ising model also has an additional term in its energy contributed by the second Rényi entropy of the bulk state ρ b , and the Ising spins at the boundary vertices are coupled to a boundary "pinning field" h x determined by the boundary region A. If the bulk contribution from ρ b is small (which means major part of quantum entanglement of the boundary states is contributed by the tensor network itself), one can see that the parameters log D xy and log D x∂ determine the effective temperature of the Ising model. For simplicity, in the following we assume D xy = D x∂ = D for all internal legs and boundary dangling legs. In this case we can take β = 1 2 log D as the inverse temperature of the classical Ising model.

Ryu-Takayanagi formula
Once the mapping to the classical Ising model is established, it is easy to see how the Ryu-Takayanagi formula emerges. In the large D limit, the Ising model is in the low-temperature long-range ordered phase (as long as the bulk has spatial dimension ≥ 2), so that the Ising action can be estimated by the lowest energy configuration. The boundary pinning field h x leads to the existence of an Ising domain wall bounding the boundary region A, and in the absence of a bulk contribution the minimal energy condition of the domain wall is exactly the RT formula. In this section we will analyze this emergence of the Ryu-Takayanagi formula and corrections due to bulk entanglement in more detail.

Ryu-Takayanagi formula for a bulk direct-product state
We first consider the simplest situation with the bulk state a pure direct-product state In this case one can contract the bulk state at each site with the tensor of that site, which leads to a new tensor with one fewer legs. Since each tensor is a random tensor, the new tensor obtained from contraction with the bulk state is also a random tensor. Therefore the holographic mapping with a pure direct-product state in the bulk is equivalent to a purely in-plane random tensor network, similar to a MERA, or a "holographic state" defined in Ref. [6]. The second Rényi entropy of such a tensor network state is given by the partition function of Ising model in Eq. (2.13) without the ρ b term. Omitting the constant terms that appears in both Z 0 and Z 1 , the Ising action can be written as In the large D limit, the Ising model is in the low temperature limit, and the partition function is dominated by the lowest energy configuration. As illustrated in Fig. 3 (a), the "energy" of an Ising configuration is determined by the number of links crossed by the domain wall between spin-up and spin-down domains, with the boundary condition of the domain wall set by the boundary field h x . For the calculation of denominator Z 0 , h x = +1 everywhere, so that the lowest energy configuration is obviously s x = +1 for all x, with energy F 0 = 0. For F 1 , the nontrivial boundary field h x = −1 for x ∈ A requires the existence of a spin-down domain. Each link xy with spins anti-parallel leads to an energy cost of log D. Therefore the Rényi entropy in large D limit is The minimization is over surfaces Σ such that Σ ∪ A form the boundary of a spin-down domain, and |Σ| denotes the area of Σ, i.e., the number of edges that cross the surface. Therefore the minimal area surface, denoted by γ A , is the geodesic surface bounding A region. Here we have assumed that the geodesic surface is unique. More generally, if there are k degenerate minimal surfaces (as will be the case for a regular lattice in flat space), F 1 is modified by − log k. With this discussion, we have proved that Ryu-Takayanagi formula applies to the second Rényi entropy of a large dimensional random tensor network, with the area of geodesic surface given by the graph metric of the network. As will be discussed later in Section 5, the higher Rényi entropies take the same value in the large D limit, and it can also be extended to the von Neumann entropy, at least if the minimal geodesics are unique (see Section 7). However, the triumph that the second Rényi entropy is equal to the area of the minimal surface in the graph metric is in fact a signature that the random tensor construction deviates from the holographic theory. The holographic calculation of the second Rényi entropy amounts to evaluating the Euclidean action of the two-fold replica geometry, which satisfies the Einstein equation everywhere in the bulk. Thus, in general, the second Rényi entropy does not exactly correspond to the area of the minimal surface in the original geometry. Due to the back-reaction of the gravity theory, the n-fold replica geometry is in general different from the geometry constructed by simply gluing n copies of the original geometry around the minimal surfaces, the discrepancy between which can be seen manifestly from the n-dependence of the holographic Rényi-n entropy. We will see in Section 5 that our random tensor model can reproduce the correct Rényi entropies for a single boundary region if we replace the bond states |xy by appropriate short-range entangled states with non-trivial entanglement spectrum. However, this does not resolve the problem for multiple boundary regions, for which we will have a more detailed discussion in Section 5.
To compare with the RT formula defined on a continuous manifold, one can consider a triangulation of a given spatial manifold and define a random tensor network on the graph of the triangulation. (See [22, Appendix A] for further discussion of the construction of the triangulation graph.) Denoting by l g the length scale of the triangulation (the average distance between neighboring triangles), the area |γ A | in our formula is dimensionless and the area |γ c A | defined on the continuous Riemann manifold is given by |γ c A | = l d−1 g |γ A | (when the spatial dimension of bulk is d). Therefore S(A) = l 1−d g log D |γ c A |, and we see that l 1−d g log D corresponds to the gravitational coupling constant 1 4G N . Compared to previous results about the RT formula in tensor networks [6,8], our proof of RT formula has the following advantages: Firstly, our result does not require the boundary region A to be a single connected region on the boundary. Since the entropy in the large D limit is always given by the Ising spin configuration with minimal energy, the result applies to multiple boundary regions. Secondly, our result does not rely on any property of the graph structure, except for the uniqueness of the geodesic surface (if this is not satisfied then the entropy formula acquires corrections as discussed above; cf. Section 9). If we obtain a graph by triangulation of a manifold, our formula applies to manifolds with zero or positive curvature, even when the standard AdS/CFT correspondence does not apply. In addition to these two points, we will also see in later discussions that our approach allows us to study corrections to the RT formula systematically. Notice that we are not limited to two-dimensional manifolds. One can consider a higher dimensional manifold and construct a graph approximating its geometry. It follows from our results that the entropy of a subregion of the boundary state is given by the size of the minimum cut on the graph, i.e., the area of the minimal surface in the bulk homologous to the boundary region.

Ryu-Takayanagi formula with bulk state correction
If we do not assume the bulk state to be a pure direct-product state, the bulk entropy term in Eq. (2.13) is nonzero. If we still take the D → ∞ limit, the Ising model free energy is still determined by the minimal energy spin configuration, which is now determined by a balance between the area law energy log D |Σ| for a domain wall Σ, and the energy cost from bulk entropy. We can define the spin-down region in such a minimal energy configuration as E A , which bounds the boundary region A, and corresponds to the region known as the entanglement wedge in the literature [23,24]. The second Rényi entropy is then given by The bulk contribution has two effects. First it modifies the position of the minimal energy domain wall |γ A |, and thus modified the area law (RT formula) term of the entropy. Second it gives an additional contribution to the entanglement entropy of the boundary region. This is similar to how bulk quantum fields contribute corrections to the RT formula in AdS/CFT [25].
To understand the consequence of this bulk correction, we consider an example shown in Fig. 3 (b) and (c), where A and B are two distant disjoint regions on the boundary. If the bulk entanglement entropy vanishes, the RT formula applies and the entanglement wedges E A and E B are disjoint. Therefore we find that S 2 (A)+S 2 (B) = S 2 (AB) and so the "mutual information" between the two intervals I 2 (A : B) = S 2 (A) + S 2 (B) − S 2 (AB) vanishes in the large D limit. 5 When the bulk state is entangled, if we assume the entanglement is not too strong, so that the entanglement wedges remain disjoint, the minimal energy domain walls γ A and γ B may change position, but remain disconnected. Therefore: From this equation, we see that even if a small bulk entanglement entropy may only lead to a minor correction to the minimal surface location, it is the only source of mutual information between two far-away regions in the large D limit. (If we consider a large but finite D, and include spin fluctuations of the Ising model, we obtain another source of mutual information between far-away regions, which vanishes exponentially with log D.) The suppression of mutual information between two far-away regions implies that the correlation functions between boundary regions A and B are suppressed, even if each region has a large entanglement entropy in the large D limit. In the particular case when the bulk geometry is a hyperbolic space, the suppression of two-point correlations discussed here translates into the scaling dimension gap of boundary operators, which is known to be a required property for CFTs with gravity duals [26][27][28]. A more quantitative analysis of the behavior of two-point correlation functions and scaling dimension gap will be postponed to Section 6.

Phase transition of the effective bulk geometry induced by bulk entanglement
We have shown that a bulk state with nonzero entanglement entropy gives rise to corrections to the Ryu-Takayanagi formula. In the discussion in Section 3.2, we assumed that the bulk entanglement was small enough that the topology of the minimal surfaces remained the same as those in the absence of bulk entanglement. Alternatively, one can also consider the opposite situation when the bulk entanglement entropy is not a small perturbation compared to the area law term log D |γ A |, in which case the behavior of the minimal surfaces may change qualitatively. In this subsection, we will study a simple example of this phenomenon, with the bulk state being a random pure state in the Hilbert space of a subregion in the bulk. As is well-known, a random pure state is nearly maximally entangled [29], which we will use as a toy model of a thermal state (i.e., of a pure state that satisfies the eigenstate thermalization hypothesis [30,31]). The amount of bulk entanglement can be controlled by the dimension of the Hilbert space D b of each site. We will show that the topologies of minimal surfaces experience phase transitions upon increasing D b which qualitatively reproduces the transition of the bulk geometry in the Hawking-Page phase transition [11,12]. To be more precise, the entropy of the boundary region receives two contributions: the area of the minimal surfaces in the AdS background and the bulk matter field correction. However, above a critical value of D b , the minimal surface tends to avoid the highly entangled region in the bulk, such that there is a region which no minimal surface ever penetrates into, and the minimal surface jumps discontinuously from one side of the region to the other side as the boundary region size increases to half of the system. This is qualitatively similar to how a black hole horizon emerges from bulk entanglement.
(A black hole cannot be identified conclusively in the absence of causal structure, however, so our conclusions in this section are necessarily tentative.) We consider a tensor network which is defined on a uniform triangulation of a hyperbolic disk. Each vertex is connected to a bulk leg with dimension D b in addition to internal legs between different vertices. Then we take a disk-shaped region, as shown in Fig. 4 (a). We define the bulk state to be a random state in the disk region, and a direct-product state outside: The second Rényi entropy of a boundary region is determined by the Ising model partition function with the action (2.13). 6 The bulk contribution S 2 ({s x = −1} ; ρ b ) for a random state with large dimension only depends on the volume of the spin-down domain in the disk region, since all sites a play symmetric role. After an average over random states, the entropy of a bulk region with N sites is given by [32] in which N T is the total number of sites in the disk region. Therefore the Ising action contains two terms, an area law term and the bulk term which is a function of the volume of spin-down domain. For simplicity, we can consider a fine-grained triangulation and approximate the area and volume by that in the continuum limit. If we denote the average distance between two neighboring vertices as l g , as in previous subsections, we obtain Here M ↓ is a spin-down region bounding a boundary region A, and ∂M ↓ is the boundary of this region in the bulk (which does not include A). V T = N T l 2 g is the total volume of the disk region in the bulk.
Consider the Poincaré disk model of hyperbolic space, with the metric ds 2 = 4(dr 2 + r 2 dθ 2 )/(1 − r 2 ) 2 . The boundary is placed at r = 1 − with > 0 a small cutoff parameter. The disk region is defined by r ≤ b. Choose a boundary region θ ∈ [−ϕ, ϕ], with ϕ ≤ π/2 so that the boundary region is smaller than half the system size. (Boundary regions that exceed half the system size have the same entropy as their complement, since the whole system is in a pure state.) If we assume the minimal surface ∂M ↓ to be a curve described by r = r(θ) (i.e., for each θ there is only one r value), the volume and area of this curve can be written explicitly as For fixed l −1 g log D, when we gradually increase l −2 g log D b , there are three distinct phases: the perturbed AdS phase, the small black hole phase, and the maximal black hole phase. The phase diagram can be obtained numerically, as shown in Fig. 4 (b). In the calculation, we fix b = tanh(1/2), which means that the radius of the disk in proper distance is 1 (i.e., the AdS radius). In the perturbed AdS phase, although the minimal surfaces are deformed due to the existence of the bulk random state, there is no topological change in the behavior of minimal surfaces. As the size of the boundary region increases, the minimal surface swipes through the whole bulk continuously ( Fig. 5 (a)). In the small black hole phase, the minimal surface experiences a discontinuous jump as the boundary region size increases. There exists a region with radius 0 < r c < b that cannot be accessed by the minimal surfaces of any boundary regions ( Fig. 5 (b)). Qualitatively, the minimal surfaces therefore behave like those in a black hole geometry, which always stay outside the black hole horizon. As log D b increases, r c increases until it fills the whole disk (r c = b). Further increase of log D b does not change the entanglement property of the boundary anymore, since the entropy in the bulk disk region has saturated at its maximum. This is the maximal black hole phase (Fig. 5 (c)).
More quantitatively, the two phase boundaries in Fig. 4 The square root behavior of the blue line can be understood by taking the maximal boundary region of half the system size ϕ = π 2 . At the critical l −2 g log D b , the diameter of the Poincaré disk goes from the minimal surface bounding the half system to a local maximum. For more detailed discussion, see Appendix A. The second transition at the red line is roughly where the entanglement entropy of the bulk region reaches its maximum. However, more work is required to obtain the correct coefficient (1 + b 2 )/2b, as we show in Appendix A. In Fig. 6 (a), we present the evolution of the black hole size r c /b when l −2 g log D b increases and l −1 g log D = 10 is fixed. Fig. 6 (b) provides another diagnostic to differentiate the geometry with and without the black hole. The entanglement entropy S 2 (ϕ) is plotted as a function of the boundary region size. In the perturbed AdS phase (blue curve), S 2 (ϕ) is a smooth function of ϕ, just like in the pure AdS space. In the small black hole phase (black curve) and the maximal black hole phase (red curve), there is a cusp in the function S 2 (ϕ) at ϕ = π 2 , as a consequence of the discontinuity of the minimal surface. For ϕ ≤ π 2 , S 2 (ϕ) shows a crossover from the AdS space behavior (which corresponds to the entanglement entropy of a CFT ground state) to a volume law. Such behavior of S 2 (ϕ) is qualitatively consistent with the behavior of a thermal state (more precisely a pure state with finite energy density) on the boundary.
In summary, we see that a random state in the bulk region is mapped by the random tensor network to qualitatively different boundary states depending on the entropy density of the bulk. This is a toy model of the transition between a thermal gas state in AdS space and a black hole. In a more realistic model of the bulk thermal gas, the thermal entropy is mainly at the IR region (around the center of the Poincaré disk), but there is no hard cutoff. Therefore there is no sharp transition between small black hole phase and maximal black hole phase. The size of black hole will keep increase as a function of temperature. In contrast, the lower phase transition between perturbed AdS phase and the small black hole phase remains a generic feature, since the minimal surface will eventually skip some region in the bulk when the volume law entanglement entropy of the bulk states is sufficiently high. From this simple example we see how the bulk geometry defined by a random tensor network has nontrivial response to the variation of the bulk quantum state. Finding a more systematic and quantitative relation between the bulk geometry and bulk entanglement properties will be postponed to future works.
At last, we comment on the case of two-sided black holes. As is well-known, an eternal black-hole in AdS space is the holographic dual of a thermofield double state [33], which is an entangled state between two copies of CFTs, such that the reduced density matrix of each copy is thermal. As a toy model of the eternal black hole we consider a mixed bulk state with density matrix Here ρ pure x is a pure state density matrix while ρ mix x is a mixed state with finite entropy. This density matrix described a bulk state in which all qudits in the disk region |x| < b are entangled with some thermal bath. The behavior of the geometry can be tuned by the entanglement entropy of ρ mix 2. Similar to the single-sided case, there is a second phase transition where further increase of bulk entropy density does not change the boundary entanglement feature any more. The transition point for two-sided case occurs at a slightly different value When the bulk entropy exceeds this value, the boundary state is a mixed state with entropy l −1 g log D 4πb 1−b 2 , which is given by the boundary area of the disk region in the bulk. The boundary of the disk plays the role of the black hole horizon.
While the behavior observed here is consistent with black hole formation, it is important to stress that the conclusion is actually ambiguous. Geodesics can be excluded from regions of space even in the absence of a black hole. 7 The presence of a black hole is ultimately a feature of the causal structure, so resolving the ambiguity would require introducing time into our model.

Random tensor networks as bidirectional holographic codes
In the previous section we discussed the entanglement properties of the boundary quantum state obtained from random tensor networks. In this section we will investigate the properties of random tensor networks interpreted as holographic mappings (or holographic codes).
In Ref. [8], the concept of a bidirectional holographic code (BHC) was introduced, which is a holographic mapping with two different kinds of isometry properties. A BHC is a tensor network with boundary legs and bulk legs. We denote the number of boundary legs as L and the number of bulk legs (i.e., the number of bulk vertices) as V , and denote the dimension of each boundary leg as D and that of each bulk leg as D b . The first isometry is defined from the boundary Hilbert space with dimension D L to the bulk Hilbert space with dimension D V b . The physical Hilbert space is identified with the image of this isometry from the boundary to the bulk, so that the full bulk Hilbert space is redundant in the sense that it contains many non-physical states. The condition identifying these physical states can be formulated as a gauge symmetry. The second isometry is defined from a subspace of the bulk Hilbert space to the boundary. The physical interpretation of this subspace is as the low energy subspace of the bulk theory. The bulk theory is intrinsically nonlocal in the space of all physical states, but locality emerges in the low energy subspace. More precisely, the degrees of freedom at different locations of the low energy subspace are all independent, and a local operator acting in the low energy subspace can be recovered from certain boundary regions, satisfying the so-called "error-correction property" [6,9]. For this reason, the low energy subspace is also referred to as the code subspace.
In this section, we will investigate the properties of random tensor networks and show that they satisfy the BHC conditions in the large D limit and moreover have properties that are even better than the BHC constructed using pluperfect tensors in Ref. [8].

Code subspace
We start from the holographic mapping in the low energy subspace, or "code subspace" in the language of quantum error correction [9]. Physically, the code subspace is a subspace of the Hilbert space which corresponds to small fluctuations around a classical geometry in the bulk. More precisely, the criterion of "small fluctuations" states that these states are described well by a bulk quantum field theory with the given geometrical background. In other words, in the code subspace the bulk fields (operators) at different spatial locations are independent and the Hilbert space seems to factorize with respect to the bulk position. The fact that one cannot take the code subspace to be the entire Hilbert space, i.e. that locality in the bulk fails if we consider the entire Hilbert space, is the essential feature of a theory of quantum gravity (defined as the holographic dual of a boundary theory), as compared to an ordinary quantum field theory in the bulk.
In general, the choice of code subspace is not unique. However, the random tensor network approach allows for a simple and explicit choice. We define the code subspace to be the tensor product of lower-dimensional subspaces at each vertex of the graph: The holographic mapping restricted to this subspace is simply a tensor network with a smaller bond dimension D b for each bulk leg. In the following, we investigate the condition for the bulk-to-boundary map to be an isometry, which thus determines the value of D b that makes such a subspace an eligible code subspace.
When we view the tensor network as a linear map M from the bulk to the boundary, the isometry condition means M † M = I is the identity operator. To apply the results we obtained for the second Rényi entropy, it is more convenient to view the tensor network as a pure state. Choose an orthonormal basis {|α } of the bulk and a basis {|a } for the boundary. The linear map M with matrix element M αa = α|M |a can then be identified with the pure quantum state In terms of the state, the requirement that M † M = I is equivalent to the statement that the bulk reduced density matrix I is maximally mixed. Therefore, the isometry condition can be verified by an entropy calculation.
For that purpose we calculate the second Rényi entropy of the whole bulk. In the large D limit, this is mapped to an Ising model partition function in the same way as in the RT formula discussion, except that there is now a pinning field everywhere in the bulk, in addition to the boundary: For computation of the bulk-boundary entanglement entropy, we should take b x = −1 for all x, and h x = +1 for all boundary sites. (We have written Eq. (4.2) in this general form because other configurations of h x , b x will be used in our later discussion.) In this action, the effect of the bulk pinning field b x competes with the boundary pinning field h x . The relative strength of these two pinning fields is determined by the log D, the lowest energy configuration will be the one with all spins pointing up. In the opposite limit log D b log D, all spins point down. For the purpose of defining a code subspace with isometry to the boundary, we consider the limit log D b log D. In that case all spins are pointing up, and the only energy cost in the Ising action (4.2) comes from the last term, leading to the entropy which is the maximum possible for a state on the bulk Hilbert space since its dimension is D V b . In the limit log D b log D, D → ∞, the bulk is therefore in a maximally mixed state, so the corresponding holographic mapping from the bulk to the boundary is isometric. The isometry condition is equivalent to the condition that the lowest energy configuration of the Ising model has all spins pointing up.
Instead of requiring log D b log D, we can write down more precisely the isometry condition by requiring that the all-up configuration has the lowest energy. Consider a generic spin configuration with a spin-down domain Ω. The energy of this configuration is Here |Ω| and |∂Ω| are the volume and the surface area of Ω, respectively. In order for the all-up configuration to be stable, we need A(Ω) > V log D b for all nontrivial Ω, which requires For example, if the bulk is a (triangulation of) hyperbolic space (with curvature radius R = 1), a disk with boundary area |∂Ω| = 2πR/l g has volume |Ω| = 2π Here we have measured both area and volume by the triangulation scale l g . Therefore the isometry condition requires There is a finite range of D b which satisfies the isometry condition, which is a consequence of the fact that the area/volume ratio is finite in hyperbolic space. For comparison, the same discussion for a disk in flat space with boundary area 2πR will require log D b log D < 2 R l g . Therefore the ratio log D b / log D must scale inversely with the size of the whole system R max . 8 A useful remark is that the isometry condition (4.4) (or more precisely, a slightly weaker condition with < replaced by ≤) is obviously necessary by a counting argument: In order for an operator defined in region Ω to be mapped to the boundary isometrically, it needs to be first mapped to the boundary of Ω, so that the dimension of the Hilbert space at the boundary D |∂Ω| must be at least as large as the dimension of the bulk Hilbert space D |Ω| b . With this observation, what we see from the Ising model representation is that the large-D random tensor network is an optimal holographic code, in the sense that an isometry is defined as long as the counting argument does not exclude it. Of course one should keep in mind that this optimal property is only true asymptotically in the large D limit.

Entanglement wedges and error correction properties
Having shown that the holographic mapping M defines an isometry from the bulk to the boundary degrees of freedom for suitable ratios log D b / log D, it is natural to ask whether this isometry has the error correction properties proposed in Ref. [9], i.e., whether operators in the bulk can be recovered from parts of the boundary instead of from the whole boundary. Specifically, consider an operator φ C in the bulk which only acts nontrivially in a region C. Denote the complement of C in the bulk by C. We say that φ C can be recovered from a boundary region A if there exists a boundary operator O A such that [6] We note that condition (4.6) is composable: For example, if φ C and φ C can be recovered from A and A , respectively, then for any bulk state ρ b and the corresponding boundary state ρ = M ρ b M † . In the same way, an arbitrary n-point function in the bulk can be obtained from a corresponding correlation function on the boundary. In the language of quantum error correction, Eq. (4.6) states that the logical operator φ C acting on the degrees of freedom in C can be realized by an equivalent physical operator acting on the degrees of freedom in A only. We are now interested in understanding when all operators φ C in the region C can be recovered from A. That is, we would like the quantum information stored in subsystem C to be protected against erasure of the degrees of freedom in B, the complement of A on the boundary. This amounts to another entropic condition, namely, that in the pure state |Ψ M defined in Eq. (4.1) there is no mutual information between C and the region BC [34], which ensures that the mutual information between A and C is maximal: For the reader's convenience, we recount a short proof of this fact in Appendix B.
In general it is important that Eq. (4.7) is evaluated in terms of von Neumann entropies rather than Rényi entropies. In the limit of large D, however, both entropies are closely approximated by the Ryu-Takayanagi formula as long as the minimal surfaces are unique (see Section 7). What is more, we may even arrange for the Ryu-Takayanagi formula to be satisfied exactly, without any assumption on the uniqueness of minimal surfaces, by using ensembles of random stabilizer states instead of Haar random states (see Section 9). In the following we shall therefore evaluate the quantum error correction condition (4.7) in terms of second Rényi entropies and assume (for simplicity) that the RT formula holds exactly.
To understand when the error correction condition holds, we consider the configuration shown in Fig. 7. The calculation of S 2 (C) is straightforward. Given the isometry condition (4.4), the whole bulk is in a maximally mixed state after tracing over the boundary, so that S 2 (C) also takes the maximal value |C| log D b . In the calculation of S 2 (BC), the pinning field is set to b x = −1 for x ∈ B and h x = −1 for x ∈ C. The boundary spin-down field in B will pin a spin-down domain (orange region in Fig. 7). We consider the case when C is in the spin-up (blue) domain, in which case the energy cost gives the entropy Here γ A is the domain wall bounding region A, and E A is the spin-up domain, which is the entanglement wedge of A. The first term is the area law energy cost of the domain wall, and the second term is the volume law energy cost. S 2 (BCC) can be computed similarly by flipping the pinning field in C to downwards. Due to the isometry condition (4.4), flipping the field in C does not create new spin-down domains, so that the only difference between S 2 (BCC) and S 2 (BC) is an additional energy cost in the C region that is exactly S 2 (C). Therefore condition (4.7) holds, and the operators in C can be recovered from A. As a final note, observe that the domain wall γ A is generally not the minimal surface, due to the presence of the bulk pinning field, but our conclusion holds as long as C is in the spin-up domain and is disconnected from γ A .
For comparison, we can consider the same configuration in Fig. 7 and ask whether operators in C can be recovered from B. This requires the calculation of S 2 (C) + S 2 (AC) − S 2 (ACC). Following an analysis similar to the previous paragraph, one can obtain Here E B is the complement of E A in the bulk, which is the entanglement wedge of B. Therefore the mutual information I 2 (C : AC) = 2 S 2 (C) > 0, so that C cannot be recovered from B.
From the two cases studied above, we can see that operators in a bulk region C can be recovered from a boundary region A if and only if C is included in the entanglement wedge E A of A. It should be noted that this statement only applies to small bulk D b , or for sufficiently small regions C if D b is larger, when the entanglement wedge E A (spin-down domain in the Ising model) is independent of the direction of the pinning field in C.

Gauge invariance and absence of local operators
In the two subsections above, we showed how a large D and small D b random tensor network defines bulk-to-boundary isometries with error correction properties. In this subsection we would like to investigate the other direction of the BHC, i.e., the boundary-to-bulk isometry. To define this isometry, we need to require that the boundary-bulk entanglement entropy be equal to |A| log D, which is the maximum possible entropy for the boundary. This requires the opposite condition from Eq. (4.4): To satisfy this condition, we can take Ω = {x} as a single site in the bulk, for which the condition is reduced to D b > D nx , with n x the number of links connected to x. If this condition is satisfied for each site, Eq. (4.8) also applies to other regions, since |∂B| ≤ x∈B n x always holds. Therefore the condition ensuring a boundary-to-bulk isometry is This is similar to the condition proposed in Ref. [8], with the difference that Ref. [8] has D b = D nx because each tensor is required to be rigorously a unitary mapping from the in-plane legs to the bulk leg. When this isometry condition is satisfied, the boundary-to-bulk isometry maps each boundary state isometrically to a bulk state in a larger Hilbert space with dimension D V b . It should be clarified that the physical Hilbert space is always that of the boundary, and that the D V b -dimensional Hilbert space, which is factorizable into a direct product of each bulk site, is just an auxiliary tool. The situation is very similar to a gauge theory, in which one can embed gauge invariant states into a larger auxiliary Hilbert space by treating the gauge vector potential as a physical field. In fact, it was shown in Ref. [8] that the physical Hilbert space -the image of the boundary Hilbert space under the holographic mappingcan be defined by a gauge invariance condition. The discussion also applies to the random tensor network satisfying condition (4.9).
The main property of the boundary-to-bulk isometry is that the bulk theory is intrinsically nonlocal. To be more precise, consider an arbitrary region C that disconnected from the boundary, as shown in Fig. 8. We would like to show that any operator φ C supported in C is mapped to the boundary trivially, i.e., Here, we have denoted the whole boundary as region A, while I A is the identity operator on the boundary, and c is a constant. This statement is equivalent to the statement I(A : C) = S(A) + S(C) − S(AC) = 0, which means there is no mutual information between C and the whole boundary. Following an argument similar to that of the previous subsection, and using condition (4.8) one can easily conclude that as is illustrated in Fig. 8. Therefore all purely bulk operators are trivial, and only those in regions adjacent to the boundary contain nontrivial information about boundary physical operators. As was discussed in Ref. [8], this property is a consequence of the gauge symmetry of the tensor network. For all tensor networks, there is a gauge symmetry induced by acting unitarily on each internal leg while preserving the physical state after contraction. However, for tensor networks with the boundary-to-bulk isometry property, this gauge symmetry is isometrically mapped to constraints on the bulk legs.
In summary, we have shown that a BHC can be built from a large D random tensor network with bulk leg dimension D b satisfying condition (4.9). The boundary theory is mapped isometrically to a nonlocal theory in the bulk, with the physical (boundary) Hilbert space defined by gauge constraints. A code subspace is defined by a local projection at every bulk site to a smaller subspace with dimension D b which satisfies condition (4.4). A bulk-to-boundary isometry is defined in the code subspace, and a bulk local operator in the code subspace can be recovered from a boundary region as long as the entanglement wedge of this region encloses the support of this bulk operator. In this way, random tensor networks can be used to define a bulk theory with intrinsic nonlocality and emergent locality in a subspace, as is desired for a theory of quantum gravity.

Higher Rényi entropies
In this section, we will generalize the second Rényi entropy calculation to higher Rényi entropies, and show that the higher Rényi entropies of a random tensor network are also mapped to partition functions of classical spin models, with the spin now living in a different target space, the permutation group Sym n of {1, . . . , n}. For n = 2, the permutation group Sym 2 = Z 2 reduces to the target space of the Ising model.
The derivation is in exact parallel with that for the second Rényi entropy in Section 2.2 For the random tensor network state given by Eq. (2.3), the n-th Rényi entropy is: Again we use the natural logarithm to define higher Rényi entropies. We now define: Here ρ ⊗n denotes the direct product of n-copies of ρ, and C

(n)
A is the permutation operator that permutes the n copies cyclically in A region. For a basis |m A of region A, a basis of the direct product space is given by |m 1A ⊗ |m 2A ⊗ · · · ⊗ |m nA , and the action of C As in Section 2.2, we approximate the typical Rényi entropy by and compute Z The average of |V x V x | ⊗n results in a projector onto the symmetric subspace of the n-fold tensor power Hilbert space (e.g., [21]): Here g x runs over all permutation group elements and we identify g x with its action on the n-copy single site Hilbert space. This action is defined by permuting the n copies of systems, similar to the definition of C becomes a partition function of classical spin model: The statistical weight of a configuration {g x } is determined by the expectation value of this permutation (multiplied by the additional cyclic permutation C

(n)
A on the boundary) in the state ρ ⊗n P . In general, this expectation value is not related to the Rényi entropy of any bulk region, which is a key difference from the case of second Rényi entropy. One can view such expectation values of permutation operators as generalized multi-partite entanglement measures that contains more information than Rényi entropies. 9 In our case, ρ P = ρ b ⊗ xy |xy xy| by Eq. (2.4). Thus the action A (n) [{g x }] becomes a sum of bond contributions and contributions of the bulk state ρ b , similar to (2.13) in the case of the second Rényi entropy: Here, χ(g) denotes the number of cycles in a permutation g (including cycles of length one). The boundary pinning field h x takes the value x the cyclic permutation acting on site x. For n = 2, it is straightforward to check that (5.6) reduces to the Ising action (2.13). The EPR pairs on the internal legs contribute a two-spin interaction energy which vanishes only if g x = g y . In this sense, the interaction is "ferromagnetic", which prefers all g x to align. We first consider the case when the bulk is a pure direct-product state, so that the contribution of ρ b to Eq. (5.6) vanishes. We also take D xy = D x∂ = D for simplicity, as in previous sections. The action (5.6) describes a Sym n -spin model with ferromagnetic interaction and a boundary pinning field, at inverse temperature β = log D. Some examples of {g x } configurations are shown in Fig. 9. Each domain wall between two different values of g x has an energy cost which is proportional to the area of the domain wall and χ(g −1 x g y ). Up to this point, the derivation applies to arbitrary values of D. In the large D limit, the partition function is dominated by the lowest energy contribution. If the entanglement wedge E A is unique (i.e., if the Ising model used to evaluate the second Rényi entropy has a unique minimal energy configuration) then the spin model with action (5.6) likewise has a unique minimal energy configuration. It is given by setting g x equal to the cyclic permutation throughout region E A and to the identity elsewhere (see Fig. 9 (b) for an illustration in the case n = 3). We give a detailed proof of this fact in Appendix C. Since the boundary of E A in the bulk is the geodesic surface γ A , we obtain 9 This is because they are invariant under local unitary transformations that only act on a domain with the same permutation value gx = g. Such quantities are known as LU invariants in the quantum information literature (see, e.g., Ref. [35] and references therein).
where the constant prefactor is independent from the choice of region A and will be canceled by the same factor in the denominator Z (n) 0 . The factor (1 − n) comes from the fact that the cyclic permutation contains one loop. Therefore we conclude that the typical n-th Rényi entropy of the random tensor network state is given by If the bulk is in an entangled state ρ b then we need to consider the corresponding contribution to the action (5.6), which is given by As long as the bulk dimension is not too large, we may think of (5.8) as a perturbation to the statistical model for the direct product case. In the minimal energy configuration of the unperturbed model, the contribution (5.8) is precisely equal to (n−1)S n (E A ; ρ b ), i.e., (n−1) times the n-th Rényi entropy of the reduced density matrix of the entanglement wedge in the bulk state. For a general configuration {g x }, however, (5.8) cannot be interpreted as an entropy. In fact, the corresponding statistical weight tr ρ ⊗n b x g x can even be a complex number, so that the interpretation of the action (5.6) requires suitable care. However, we note that the partition function (5.4) is by definition an average of the positive quantities (5.1) and therefore always positive. The choice of branch for the logarithm in (5.8) is also irrelevant for the resulting statistical weight and so does not concern us further.
The key observation now is that |tr ρ ⊗n b x g x | ≤ 1 by the Cauchy-Schwarz inequality. Thus the real part of (5.8) is always non-negative: The bulk correction only ever increases the real part of the energy levels. In particular, the only way that the (real part of the) energy gap can decrease in the perturbed model is due to the bulk corrections in the minimal energy configuration of the unperturbed model. Since S n (E A ; ρ b ) ≤ log D b |E A |, the energy gap can therefore be lower bounded by log D − (n − 1) log D b |E A |. As long as this gap diverges for large D, the minimal energy configuration remains unchanged and dominates the partition sum, so that We conclude that, for an entangled bulk state of sufficiently low dimension and large D, the typical n-th Rényi entropy of the random tensor network state is given by It should be noted that these conclusions only hold when there is a unique minimal geodesic surface, as discussed above. When there are multiple degenerate minimal surfaces, the entropy is reduced by a factor log N with N the number of minimal energy configurations. An example of degenerate minimal energy configurations are shown in Fig. 9 (c) for a square lattice in flat space. The fact that in leading order all Rényi entropies approach the same value in the large D limit tells us an important difference between a large D random tensor network state if the minimal surface is degenerate. and a large central charge CFT ground state. 10 When A is a length-l interval, the Rényi entropy of A in a (1 + 1)-dimensional CFT ground state is given by [36,37] S n (A) = 1 + 1 n c 6 log l, (5.10) which shows that in the large central charge limit, the n dependence remains nontrivial. In term of the eigenvalue spectrum of the reduced density matrix ρ A , this difference tells us that in a random tensor network state, the eigenvalues of ρ A are more strongly concentrated than that in a CFT ground state, although the density of states is also highly peaked in the latter case. From the point of view of the dual gravity theory, the nontrivial n-dependence is enforced by the requirement that the dual geometry of tr[ρ n ] ought to satisfy the equations of motion. To be more specific, if we naively constructed the dual geometry of tr[ρ n ] by gluing n copies of the original bulk geometry around the minimal surfaces (in view of Eq. (5.1), this is quite literally what the calculation of the Rényi entropy of the tensor network state amounts to) then there will be no n-dependence of the Rényi entropy [38]. The problem here is that this naively replicated geometry does not satisfy the equations of motion. In other words, the geometry does not backreact and converge to the saddle point of some gravitational action.
If we are interested in modifying the random tensor network to realize the same Rényi entropy behavior (5.10) as a CFT ground state, the simplest way is by replacing the maximally entangled EPR pair state |xy at each internal leg by a more generic state. For a more generic link state |L xy , the calculation of Z (n) 1 still applies, with the link terms (5.7) in the action (5.6) replaced by − log tr |L xy L xy | ⊗n g x ⊗ g y . (5.11) Since such terms are always non-negative, and only vanish if g x = g y , the qualitative behavior of the Sym n -spin model remains ferromagnetic, and the lowest energy configuration in large D limit still only contains a single domain wall with minimal area bounding A, as shown in Fig. 9 (b). In this case one obtains the following formula for the Rényi entropies, with S n (|L xy ) the Rényi entropy of the bond state |L xy with the partition between x and y. As a reminder, l g is the area occupied by each leg of the tensor, in units of the AdS radius R. Therefore the Rényi entropy behavior is identical to that of a CFT ground state if the tensor network is a triangulation of the hyperbolic space, and the bond state |L xy satisfies S n (|L xy ) = l g c 12 1 + 1 n . For example, one can take l g = 2 c , and define the state |L xy as a thermofield double state of the free boson CFT. In this case, the reduced density matrix ρ x = tr y |L xy L xy | is a thermal density matrix of the free boson CFT on the torus with aspect ratio β L of order one. The aspect ratio can be chosen to fit the entropy S n = 1 6 1 + 1 n . The random tensors on each site impose random projections acting on these simple free CFT states, each of which is defined on a small circle. This random projection defines a state on the boundary which for single intervals has the Rényi entropy behavior of a strongly correlated CFT in 1 + 1 dimensions. However, for more complicated subsystems, such as a disjoint union of two distant intervals [39], the Rényi entropy S n exhibits a dependence on n that cannot be accommodated by a suitable choice of link state (which only affects the Rényi entropy per bond but not the minimal surface or the replica geometry itself). 11 Besides, it is not quite clear whether the above modification will reproduce the correct Rényi entropies even for a single region if we go to the higher dimensions. More systematic investigation of the comparison between random tensor networks and CFT ground states will be reserved for future works.
The mapping we derived from the calculation of n-th Rényi entropy to classical Sym nspin models applies to more general situations. For example, studying the second-order correction terms in the calculation of average second Rényi entropy, Eq. (2.8), involves the computation of δZ 2 1 and δZ 2 0 . These quantities are quartic in x |V x V x |, so that one can apply formula (5.3) and translate δZ 2 1 into a partition function of Sym 4 -spin model. The only difference between δZ 2 1 and Z 1 in the 4-th Rényi entropy calculation is the value of boundary field. In the calculation of δZ 2 1 the value of boundary field should be chosen as permutation (12) (34) in region A, and identity elsewhere. We will use this strategy in Section 7 to bound the fluctuations of the Rényi entropies around their semiclassical value (5.9). Another important application of this mapping with n > 2 is the calculation of boundary two-point functions, which will be discussed in next section.

Boundary two-point correlation functions
As we discussed in Section 3, the mutual information between two distant regions ( Fig. 3 (b)) does not grow with log D (if the bulk state ρ b remains D-independent), although the entropy 11 We would like to thank Xi Dong for explaining this to us. of each region is proportional to log D. This observation indicates that two-point correlation functions between A and B are suppressed, as a consequence of strong multi-partite correlation in the random tensor network state. In this section we will investigate the behavior of two-point correlation functions more systematically, making use of the state averaging techniques.
We consider two small regions A and B, whose entanglement wedges are disconnected, as was shown in Fig. 3 (b). Here the bulk has a given state ρ b and the entanglement wedges are defined with respect to the Ising action (2.13) is determined by the correlation matrix. Therefore M αβ contains complete information about correlation functions between A and B.
To define a basis-independent measure of correlation, one natural choice is the singular value spectrum of M αβ . Denote the singular value decomposition of M as M αβ = s U αs λ s V βs , with λ s ≥ 0 real singular values, and U, V unitary matrices. This decomposition tells us that there is a particular set of operators which satisfies This set of operators can be considered as the analogs of the quasi-primary fields in a conformal field theory, and the singular values λ s are basis-independent measures of two-point correlations between A and B.
Instead of directly carrying out singular value decomposition of M and studying λ s , it is more convenient to consider the following quantity: Knowing C 2n for all integers n determines the singular values λ s , in the same way that the eigenvalue spectrum of a density matrix is determined by all Rényi entropies. On the other hand, using the orthonormality condition (6.1), C 2n can be reexpressed into the following form: Here X A and Y B are two permutation operators which means X A permutes each copy with an odd label 2k − 1 with copy 2k, and Y B permutes each copy 2k + 1 with copy 2k. The details of this derivation are presented in Appendix D. In this way, we have expressed C 2n in a form similar to the 2n-th Rényi entropy, with a different permutation operator. Once we obtain Eq. (6.6) it is straightforward to perform the state average, which maps C 2n in large D limit to the same classical Sym 2n -spin model as in the 2n-th Rényi entropy calculation: is the same partition function defined in Eq. (5.6), with a different boundary field h x . h x takes the value of the two permutations in Eq. (6.7) for x in A and B respectively, and identity elsewhere. The denominator is the same as that of the 2n-th Rényi entropy. In the large D limit, the minimal energy configuration that dominates Fig. 10. The minimal energy domain walls are the same as in the 2n-th Rényi entropy, which are minimal sufaces bounding A and B. However, the prefactor of the area law term is different, since tr X A = tr Y B = D n . Therefore we obtain Here X E A and Y E B are the same permutations as X A and Y B in Eq. (6.7), respectively, but acting in the bulk regions E A and E B . Interestingly, the bulk state contribution is exactly the same expression as C 2n in Eq. (6.6), but for the bulk state ρ b . In other words, we can define an orthonormal basis φ α E A and φ β E B in the bulk regions E A and E B , and define the bulk correlation matrix Following the same derivation as Eq. (6.5) and Eq. (6.6) we obtain the bulk correlation moments where λ s,bulk are the singular values of the bulk correlation matrix M b . (Note that tr (ρ b ) = 1 so that the denominator for C bulk 2n is trivial.) Therefore Eq. (6.9) can be interpreted as the following relation between the boundary correlation matrix and the boundary one: between the two entanglement wedges E A and E B , multiplied by a constant factor that is independent of the distance between the two regions (as long as their joint entanglement wedge E AB stays disconnected).
Eq. (6.12) has several important consequences. Firstly, it tells us that in a proper basis choice, there is a one-to-one correspondence between bulk two-point correlators and boundary ones. When we take the limit that both A and B are small (compared with the extrinsic curvature radius of the boundary), the entanglement wedges E A and E B become narrow regions near the boundary. In this limit the two-point correlation functions between E A and E B can be viewed as the boundary limit of bulk two-point functions. Therefore Eq. (6.12) shows that the boundary two-point functions between local operators are, up to a constant prefactor, equal to the bulk two-point functions with both points approaching the boundary. In other words, the holographic mapping defined by a random tensor network gives a bulk-boundary correspondence consistent with the usual "dictionary" of holographic duality. Secondly, if ρ b is taken to be independent from D, the bulk correlation spectrum {λ s,bulk } is D-independent. Therefore the boundary correlation spectrum is also D-independent (except for the prefactor), although the total number of operators in A and B are both increasing with D.
To understand the consequences of this observation more explicitly, we consider the special case that the bulk is the Poincaré patch of hyperbolic space, and the state ρ b is also invariant with respect to the isometry group of the bulk geometry. If we take the limit of small A and B (much smaller than the AdS radius), the bulk entanglement wedges approach the boundary, and the bulk two-point functions all decay as a power law of the boundary distance due to scale invariance. Therefore, in this limit λ s,bulk = C s |x − y| ∆s , (6.13) with {∆ s } defining the spectrum of scaling dimensions. According to Eq. (6.12), the boundary two-point functions also decay as a power law, with the same set of scaling dimensions. Compared with the situation in the AdS/CFT corresondence, we see that the boundary operators with scaling dimension ∆ s are analogs of low-dimensional operators with scaling dimensions independent of N . The number of low-dimensional operators is determined by the bulk theory. It is natural to ask whether there are also high-dimensional operators in the random tensor network state, which are the analog of "stringy" operators in AdS/CFT with scaling dimensions growing with N . To address this question, one needs to consider the finite D fluctuations. In the following we will provide some arguments about finite D corrections to the correlation spectrum which are not rigorous but may be helpful for physical understanding. At finite D the partition function Z (2n) 1 receives a contribution from other spin configurations with higher energy. Many low energy spin configurations are separate deformations of the domain walls bounding A and that bounding B. Although such fluctuations will renormalize the correlation functions, they do not change the correlation length since there is no distance dependence. The lowest energy configuration which contributes nontrivially to the distance dependence of correlation function is the one shown in Fig. 10 (b). This configuration contains two domains E 1 and E 2 with permutations the same as X A and Y B in Eq. (6.7). The boundary of E 1 ∪ E 2 consists of the connected geodesics bounding A and B, and the interface between E 1 and E 2 is chosen as the narrowest "throat" between the two geodesics. There are many configurations with similar energy, so that it is difficult to give a quantitative estimate of the finite D correction to C 2n . However as a rough estimate if we only consider the contribution of this configuration, we obtain The constant term is of order 1, given by , with |W | the width of the throat. Since E 1 and E 2 are adjacent to each other, the bulk correlation term tr[ρ ⊗2n b (X E 1 ⊗ Y E 2 )] will be dominated by short-range correlations, and thus does not decay with the distance d AB . For hyperbolic space, at long distance d AB ∝ 1 lg log |x − y|, with l g the discretization scale. Therefore the finite D correction due to this domain configuration contributes new power laws with the scaling dimension ∆ = 1 lg log D. With this new contribution to C 2n , the scaling dimension spectrum of the boundary state now contains ∆ s , 1 lg log D , which consists of the low-lying scaling dimensions ∆ s that are D-independent, and the high scaling dimension that grows linearly with log D. Such a separation in scaling dimensions is consistent with the requirement in AdS/CFT for CFT's with a gravitational dual, known as the scaling dimension gap [26][27][28]. Although the analysis here is clearly incomplete, it is reasonable to believe that the separation of two types of operators remains valid in a more detailed analysis, since there are two different origins of power law correlations, those from the bulk state and those from the spin fluctuations in the classical statistical model.

Fluctuations and corrections for finite bond dimension
In preceding sections, we have shown that the unnormalized state averages Z (n) 1 and Z (n) 0 are mapped to Ising partition functions with inverse temperature β = log D and different boundary conditions. In the large D limit, these Ising partition function are dominated by the contribution of the lowest energy spin configuration, which gave rise to the Ryu-Takayanagi formula for the Rényi entropies assuming that S n (A) = log Z 0 . In this section we will make this step precise and quantify how well the Rényi entropies S n (A) are approximated by the Ryu-Takayanagi formula. Before going into the details, we first present our conclusion: 1. For a system with volume (i.e., number of bulk vertices) V , for an arbitrary small deviation δ > 0, one can define a critical bond dimension with α and c 2n constants independent from the volume. The meaning of the exponent c 2n will be explained below. In the limit D D c the deviation satisfies is the RT formula for the n-th Rényi entropy, including the bulk correction. We will always assume that the bulk dimension D b is finite, so that in the large D limit the minimal surface γ A is determined by minimizing the area.
2. We subsequently show that under a plausible physical assumption on the free energy of the statistical models, the bound given in Eq. (7.1) can be improved by reducing the critical bond dimension to with α a non-universal constant. The meaning of the exponent ∆ 2n will be explained below.

The general bound on fluctuations
To start, we denote the Ising action (5.6) of the minimal energy spin configuration with a boundary field h x by A We shall assume throughout this section that the minimal energy configuration is unique (otherwise see Section 9 and Appendix F). In the large D limit, Z 1 , respectively. We note that A (n) 3) and the text below it. Thus the RT formula (7.2) can also be written as To bound the fluctuations of Z , we consider where we have used that Z since at finite temperature the partition function receives contributions from all spin configurations, not just the minimal energy configuration.
The key insight now is that the second moment of Z . Thus the ground state energy of this Sym 2n -spin model is essentially two times that of the Sym n -Ising model. More precisely, it follows from (7.5) that To bound the right-hand side term, we use that by assumption the minimal energy configuration is unique; all other configurations incur an additional energy cost of at least log D − (2n − 1) log D b V (cf. the discussion before Eq. (5.9)). Since there are (2n)! configurations at each bulk site this leads to the conservative upper bound where c 2n ≡ log (2n)! + (2n − 1) log D b . By combining Eqs. (7.4), (7.6) and (7.7) we obtain that The same conclusion holds for Z (n) 0 = (tr ρ) n (corresponding to a boundary field with h x = I everywhere). By Markov's inequality, it follows that and likewise for Z where we have used δ ≤ 2 such that log(1 ± δ/4) ≤ δ/2. We have thus proved that the desired bound (7.1) holds with probability at least 1 − Dc D , where D c = 32δ −2 e c 2n V . Interestingly, the above results for the Rényi entropies can be used to show corresponding statements for the von Neumann entropy. For a bulk direct product state, this is easy to see: Here, the Ryu-Takayanagi formula amounts to S RT n (A) ≡ log D |γ A |. Since S(A) ≥ S n (A) for any quantum state and S(A) ≤ log rank ρ A ≤ log D|γ A | in any tensor network state, we have essentially matching upper and lower bounds for the von Neumann entropy, and hence S(A) log D|γ A | with high probability. This result can be established more generally even in the presence of an entangled bulk state as long as D b D by adapting the techniques of [40] (cf. Section 8).

Improvement of the bound under a physical assumption
The above results establish rigorously that the entropies approximate the Ryu-Takayanagi formula in the limit of large D. However, the technique lead to a rather conservative estimate of the finite D correction, since it only proves that the entropy is close to the RT value for exponentially large bond dimension D e c 2n V . In this subsection we would like to argue based on a plausible physical assumption that actually the RT formula applies to a much larger range of D, as long as D is bigger than some power law function of V .
To start, let us reinvestigate Eq. (7.7), which was the basis of the general bound (7.1). In obtaining Eq. (7.7) we replaced the energy of all higher energy spin configurations by their minimum log D . This leads to a very conservative bound since most configurations certainly have an energy much higher than that. Since the statistical model has a local action, the number of excitations with lowest energy is actually proportional to V rather than exponential of V . Although the number of slightly higher energy excitations are super-extensive, it is still true that the free energy of the spin model is extensive at finite temperature. Furthermore, the free energy approaches the ground state energy in the lo temperature (large D) limit exponentially, since the probability of lowest energy excitation with energy E g is suppressed by the Boltzman weight e −βEg = D −Eg/2 . Using these plausible physical observations we can write the asymptotic form of the free energy Note that there is generically a power law term (log D) a multiplying the exponential factor in the free energy density. However, this power law correction is not important for our bound, since we can choose an energy ∆ 2n slightly smaller than E g , such that (7.10) In the limit D V 1/∆ 2n this implies that by choosing a constant C slightly larger than C. If we substitute this estimate for the conservative bound (7.7) then Eq. (7.8) becomes and likewise for Z (n) 0 . We may now proceed as above and conclude that, under the assumption (7.10) on the Ising models, the Rényi entropies satisfy the RT formula to arbitrary precision and with arbitrarily high probability if D V 2/∆ 2n . This improves the dependency of the bond dimension on the system size from an exponential function of V to a power law.
To illustrate the behavior of the free energy (7.10) in an explicit example, we consider an Ising model on the square lattice. (It should be noted that the Ising model case does not directly apply to the discussion above since Sym 2n -spin models are used there. However the behavior of free energy is generic for gapped spin models.) Specifically, we shall consider a cylindrical geometry given by an M × N square lattice with periodic boundary conditions along the first direction and open boundary conditions along the second one. In this setup, the minimal surface bounding a boundary region is unique. As the boundary region we choose a single interval of length L < M . Z 1 /Z ∞ 1 and Z 0 /Z ∞ 0 can be computed exactly using Onsager's solution [41,42]. The asymptotic behavior for large D is given by Therefore Eq. (7.10) holds with exponent ∆ 2n = 2. More details about this calculation are presented in Appendix E.

Possible effects of even smaller bond dimension
When D does not satisfy the condition D (CV ) 1/∆ , the deviation of the entropy from the RT value can be large. An interesting question is whether the correction to the RT formula is simply a renormalization of the coefficient of the area law, or if there is a qualitative change. For the second Rényi entropy, the quantity − log Z 1 /Z 0 is the free energy cost induced by the boundary pinning field h x in the Ising model. The behavior of this free energy cost depends on the strength of the fluctuations of the domain wall configuration. If the domain wall only fluctuates mildly around the minimal energy configuration, one can naturally expect the energy cost of the domain wall is still proportional to its area, although the coefficient may be renormalized to be different from the bare value given by the lowest energy configuration. In contrast, if the domain wall is strongly fluctuating, the energy of the domain wall may have a qualitatively different dependence in the minimal area |γ A |. Interestingly, the behavior of the domain wall in the Ising model was studied a long time ago. If the bulk spatial dimension d ≥ 3, it was found that there is a finite critical temperature T r (which is lower than the phase transition temperature T c of the Ising model), below which the fluctuations of a domain wall configuration have a finite range. The transition at T r is known as the roughening transition [43,44]. 12 It is natural to expect that the RT formula for second Rényi entropy applies for any log D > T −1 r , which remains finite even if the system size V goes to infinity. However, it is not clear how to bound the deviation of S 2 (A) from − log Z 1 Z 0 , given by the fluctuation terms in Eq. (2.8). When the bulk spatial dimension is d = 2, the domain wall is one-dimensional, and thus the fluctuation of its position is always strong. Consequently, the RT formula does not apply to any finite D if we take V → ∞ first.

Relation to random measurements and the entanglement of assistance
The average over random tensors that has played a central role in this work has appeared previously in the quantum information literature, but with a very different motivation. The definition of the boundary state |Ψ in Eq. (2.2) involves contracting the random vertex states x |V x at the bulk vertices with a bulk state |Φ b as well as a collection of Bell pairs xy |xy for the internal edges and ⊗ x |x∂ x connecting boundary vertices to their boundary connecting points ∂ x . To obtain a new physical interpretation for the state |Ψ , one can start with the state and perform a random measurement at every bulk vertex x. The post-measurement state on the unmeasured boundary vertices will then have the same distribution as |Ψ . Note Figure 11. (a) A simple tensor network consisting of two vertices V 1 and V 2 of degree 3, a pure bulk state and two boundary dangling legs ∂ 1 and ∂ 2 . (b) The construction of the Section 2. 2). The reason lies in the change of perspective; in Section 2 the random vertex states were being projected to Bell pairs and the bulk state, but here the Bell pairs and the bulk state are being projected to the random states, and therefore we need a larger Hilbert space to get a non-empty Hilbert space after projection. See Fig. (11). These two perspectives are mathematically equivalent in our examples. The post-measurement state on the unmeasured boundary vertices will then have the same distribution as |Ψ . From this point of view, boundary entanglement is being induced by performing a suitable measurement on a joint bulk-boundary state.
One of the basic problems of quantum information theory is how to establish as much high-quality entanglement as possible between spatially separated parties. One scenario that had been considered was to start with a pure state |Φ ABC of three systems and to ask how much entanglement could be induced on average between A and B upon measuring C, optimized over all possible C measurements. Because the party in possession of C is helping A and B establish entanglement, this quantity is known as the entanglement of assistance E A (A; B) Φ [45]. Concavity of the entropy implies a trivial upper bound: , and likewise for Y . In a remarkable paper, Smolin et al. showed that this upper bound was asymptotically achievable [46]: Going further, one can imagine partitioning C into subsystems C 1 , C 2 , . . . , C m and allowing only local measurements of each C j instead of joint measurements of the entire C system. From an engineering perspective, such a scenario could arise naturally if A and B are distant and the C j represent intermediate "repeater" stations in a network [47]. The additional locality restrictions will reduce the amount of entanglement that can be induced between A and B. While the concavity upper bound still applies, it can be applied here with a bit more finesse. If we choose any subset S ⊆ {C 1 , . . . , C m }, then the bound implies that this multipartite version of the entanglement of assistance, E multi A (A; B) Φ , will be bounded above by S(Φ AS c ) since the total entanglement generated between A and B will be no more than the entanglement between AS c and B after measuring S but prior to measuring S c . Likewise, E multi Consider now the special case in which |Φ has the form of Eq. (8.1) used in this paper. Set A to be any boundary region, B the complement A c of A in the boundary and identify the different subsystems C j with the bulk vertices x. The righthand side of Eq. (8.4) is then nothing other than the Ryu-Takayanagi formula with corrections due to the bulk state |Φ b , since minimizing over subsets S amounts to minimizing over cuts in the tensor network: (8.5) where E A is the bulk region corresponding to the minimizing set S and |γ A | is the size of the cut separating E A from its complement. This matches Eq. (3.3) up to the substitution of the von Neumann entropy for the second Rényi entropy. (The reason for taking the k → ∞ limit in (8.4) is essentially to make all Rényi entropies equal after suitable small perturbations to the state. For reasonable physical choices of |Φ b such as quantum field theory ground states, it should be sufficient to take k = 1 and include a small correction on the righthand side of (8.4). This has been shown, for example, in the case that A is an interval in a 1+1 dimensional CFT [48].) While the original proof of the multipartite entanglement of assistance formula used classically-inspired random coding information theory techniques, subsequent proofs proceeded by performing appropriate isotropic measurements of the C subsystems [40,49]. Because of the equivalence between contracting random tensors and performing random measurements, the analyses in the quantum information theory literature are mathematically very similar to the calculations in this article. The analog of the calculations justifying reconstruction of a bulk operator contained in the entanglement wedge of a boundary region A has even appeared, again with a different motivation, as the "split-transfer protocol" [40,50].
One could go as far as to rename the one-shot multipartite entanglement of assistance formula of [50] the "fully-quantum Ryu-Takayanagi" formula, in that it captures the essence of Ryu-Takayanagi without making any prior assumptions about the geometrical interpretation of the bulk state. Aside from connecting to pre-existing literature, one virtue of this change of perspective is that it suggests a possible physical justification for the random tensor networks in our model. One could imagine taking the state in a quantum theory of gravity and measuring the Planckian degrees of freedom of a large "bulk" subset, leaving some "boundary" degrees of freedom and bulk fields unmeasured. If the Hilbert spaces are large and the measurements generic, then the measurements should reveal almost no information about the bulk, inducing a nontrivial mapping between non-Planckian bulk degrees of freedom and the boundary. In this way, fixing the bulk Planckian degrees of freedom in the bulk-boundary state through measurement generically produces a holographic correspondence. Expanding around a particular background geometry in this picture amounts to choosing a bulk-boundary state with the correct area law entropy and randomly fixing the Planckian degrees of freedom through projection.

Random tensor networks from 2-designs
In the construction of our random tensor network state (2.1), the tensors |V x were chosen to be Haar-random, i.e., drawn from to the unitarily invariant ensemble of pure states. However, our calculations for the second Rényi entropy in Section 2.2 made use only of the second moments of the Haar measure. This calculation led to the emergence of a classical Ising model and thereby to the Ryu-Takayanagi formula. It is therefore natural to consider other ensembles of pure states whose first two moments agree with those of the Haar measure, known collectively as complex projective 2-designs [51].
It follows from the discussion in Section 7 that for a tensor network state with Haarrandom tensors and bulk direct product state ρ b , the Ryu-Takayanagi formula S(A) log D|γ A | will be satisfied with high probability in the limit of large D if the minimal geodesic is unique. This conclusion was obtained from considering higher moments of the Haar measure and therefore does not apply for a general 2-design. Another complication arises from the fact that the tensor network state can be zero (i.e., ρ = 0) with nonzero probability, in which case its entropies are not well-defined. In Appendix F we show that for any 2-design the boundary state is nonzero with high probability and that, moreover, where k denotes the number of minimal geodesics, where S RT (A) = log D|γ A | since we consider the case of a direct product bulk state, and where we write S 2 (A) =0 for the average second Rényi entropy conditioned on the boundary state being nonzero. We note that, since the lower bound in (9.1) matches the deterministic upper bound up to a constant, it follows that S(A) is at most constantly away from S RT (A) with high probability.
One random ensemble of particular interest is given by stabilizer states. Stabilizer states, defined as common eigenvectors of generalized Pauli operators, are quantum states that can be highly entangled, but whose particular algebraic structure allows for efficient simulation and effective reasoning [52]. It has been shown in [53] that pure stabilizer states in prime power dimension D = p n form a 2-design when drawn uniformly at random. Thus (9.1) applies to the entropies of the corresponding tensor network state (2.1) constructed from random stabilizer states. Such a state is again a stabilizer state, as we argue in Appendix G. The particular algebraic structure of stabilizer states implies that their reduced density matrices not only have flat spectrum (so that all Rényi entropies agree with the von Neumann entropy) but in fact that all their entropies are quantized in units of log p. It follows that for large p the Ryu-Takayanagi formula will hold exactly with high probability (even in the presence of multiple minimal geodesics).
In particular, we may use this construction to obtain random holographic codes and evaluate their error correcting properties by using (4.7) and the exact Ryu-Takayanagi formula (9.2) purely from the structure of the tensor network. In [6], holographic codes were constructed from perfect tensors, i.e., states that are maximally entangled across any bipartition, and it was shown that under certain circumstances this already implies the Ryu-Takayanagi formula (such as for single intervals in nonpositively curved space). Random stabilizer states are perfect tensors with high probability, 13 and so the analysis and results of [6] can likewise be applied to our random tensor networks constructed from stabilizer states with high bond dimension. However, our tensors are not only perfect or pluperfect [8] but also generically so and therefore can achieve the Ryu-Takayanagi formula for arbitrary subsystems.
Another consequence of (9.2) is that any entropy inequality that is valid for arbitrary quantum states, or even just for stabilizer states [54,55], is also valid for the Ryu-Takayanagi entropy formula, thereby establishing a conjecture from [22]. This can be understood as consistency check of the Ryu-Takayanagi formula, generalizing [56], where the validity of strong subadditivity was verified for the Ryu-Takayanagi formula. We refer to [57] for a detailed analysis of the entanglement properties of tensor networks built from random stabilizer states.

Conclusion and discussion
In this work we have studied the quantum information theoretic properties of random tensor networks with large bond dimension. In the following we will revisit our method from a more general perspective and summarize our findings. Viewing each tensor as a quantum state |V x , the tensor network state ρ = ρ (|V x V x |) obtained by contracting these tensors is a linear function of each tensor. Denote by f n (ρ) an arbitrary function that is a monomial function of the state ρ with degree n. Then the state average of f n over all possible choices of |V x is exactly mapped to the partition function of an classical spin model, with degrees of freedom in the permutation group Sym n , with the spins defined on the vertices of the same graph that underlies the tensor network. Different physical quantities can be translated to different functions f n (ρ). When the tensor network is used as a quantum state of the boundary, one can consider tr tr A ρ n for an arbitrary region A, which corresponds to the n-th Renyi entropy of A. When the tensor network is used as a linear map, it can be viewed as a "holographic mapping" between two parts of the degrees of freedom (boundary and bulk, respectively). In this case, in addition to the Renyi entropies one can study the entanglement entropy of a given region while another region is projected to a certain quantum state. For example, one can project the bulk into a given quantum state and study the entanglement properties of the resulting boundary state. We can also define basis-independent measures of correlation functions and relate that to a calculation of monomial functions, which allows us to study the behavior of two-point functions in the boundary state.
The mapping between the random state average and the spin-model partition function has rich consequences. For a random tensor network state, in the large D limit the Ryu-Takayanagi formula can be proven for all Renyi entropies, where the minimal surface area condition comes naturally from minimizing the energy of the spin model with given boundary conditions. The Ryu-Takayanagi formula also generalizes naturally to include bulk state corrections when there is nontrivial quantum entanglement in the bulk. As a particular example, we study the behavior of minimal surfaces in the presence of a bulk random state, and show how the minimal surface behavior can change topologically upon increase of the bulk entanglement entropy, in a way that is qualitatively consistent with black hole formation. In addition to entanglement entropy, we also studied the behavior of two-point correlation functions. The boundary correlation functions between two regions are directly determined by bulk correlation functions between two corresponding regions known as the entanglement wedges of the boundary regions. In the special case of hyperbolic space, our results on correlation functions imply that the boundary theory has power law correlations with a large scaling dimension gap. In the large D limit there are two types of scaling dimensions, those which does not scale with D coming from the bulk quantum state, and those which scale with D coming from the tensor network contribution. Such behavior of the scaling dimension gap is consistent with those of CFT ground states with a gravity dual, although the condition is necessary but not sufficient.
Random tensor networks provide a new framework for understanding holographic duality. Besides the properties studied in this paper, many other physical properties can be evaluated by the mapping to classical spin models. Compared to other tensor network models, properties of the random tensor networks can be studied much more systematically. The large dimension D limit is an analog of the large N limit in gauge theories. The fact that a random tensor network with large dimension automatically satisfies many desired properties for holographic duality further supports the point of view that semi-classical gravity is deeply related to scrambling and chaos.
There are a several open questions that shall be studied in future works. One question is whether it is possible to use a random tensor network to describe the ground state of a conformal field theory. The underlying graph of random tensor networks on hyperbolic space is invariant under a subgroup of discrete isometries of the bulk which do not involve transformation in time. Therefore we expect the distribution of tensor network states on the boundary to remain invariant under the subgroup of boundary conformal transformations that correspond to the bulk discrete isometries, modulo complications arising from the cutoff.
It is an open question whether we can modify the tensor network state to preserve the whole conformal symmetry. Related to the discussion of Rényi entropies, this may require modification of the state on links between vertices. It would also be interesting to consider random tensor network models where the same tensor is placed at each vertex. 14 Another question is how to generalize this formalism to include dynamics. What Hamiltonians of the boundary theory can be mapped to local Hamiltonians in the bulk "low energy" subspace? How to see that conserved currents on the boundary correspond to massless fields in the bulk? The answers to these questions will also be essential for understanding how the bulk gravity equation emerges.
the bulk random state. However, when this solution becomes a local maximum instead of a minimum, it means that the minimal surfaces of all the boundary regions would avoid the center of the Poincaré disk. In other words, there exists a region in the bulk inaccessible to any measurements from the boundary smaller than half system size. For convenience, we use (x, y) coordinates instead of (r, θ) in this problem. Thus what we care about is where Θ(x) = 1 when x > 0 and 0 otherwise. It is obvious that, if we treat the above expression as a matrix, the first term is always a positive definite matrix after integrating by parts of the derivative term, and the second term is a negative definite matrix, which corresponds to the fact that y = 0 minimizes the area contribution from the domain wall and maximizes the volume contribution from the bulk random pure state. Although it is hard to analytically diagonalize δ 2 S 2 (π/2) , it is straightforward to observe that the instability happens at l −1 2 . Now we turn to the second question, the transition between the small black hole phase and the maximal black hole phase. In order to understand the formation of the maximal black hole, we need a more detailed investigation of Eq. (3.4). We first focus on the random pure state region r ≤ b, and assume the minimal surface enters this region at angle ϕ and −ϕ. The minimization problem in Eq. (3.4) can be solved by asking (r (θ)) 2 + r 2 (θ) The above variational equation contains both the derivative and the integration (contained in V r(θ) ) of r(θ). But in the large D limit, which indicates that the transition happens when << D V T b near transition point. Thus in this limit, the above equation can be simplified with only r(θ) and its derivatives left.
The trick we use to solve this equation is to transform it back to a minimization problem, I[r(θ)] is the objective function to be minimized with respect to r(θ). Figure 12. Construction of the boundary operator O A corresponding to a bulk operator φ C and illustration of the recovery equation (4.6).
In the following it will be crucial to distinguish the input systems C and C of the bulk-to-boundary isometry M from the corresponding subsystems of the pure state |Ψ M defined in (4.1). We will thus denote the latter by C and C , so that where |φ + CC and |φ + C C denote maximally entangled states between C and C and between C and C , respectively. Eq. (4.7), which becomes I(C : BC ) = 0, implies at once that where τ C = tr ABC (Ψ M ) is a maximally mixed state (since M is an isometry). By definition, |Ψ M is a purification of (B.1), but we can also find a purification that respects the product structure |φ + CC ⊗ |γ BC E , obtained by purifying τ C to a maximally entangled state and tr AC (Ψ M ) to an arbitrary pure state |γ BC E . If we choose the dimension of E to be sufficiently large then the two purifications can be related by an isometry V from A to CE: It can now be readily verified that any bulk operator φ C can be recovered from A by using the boundary operator O A = V † φ C V . Indeed, Eq. (4.6), which states that O A M = M φ C , is a direct consequence of the following calculation: where we have used (B.2) and that φ C |φ + CC = φ T C |φ + CC (twice). We refer to Fig. 12 for an illustration.

C Uniqueness of minimal energy configuration for higher Rényi models
In this appendix we give a formal proof of the assertion made in Section 5 that the spin model with action (5.6) has a unique minimal energy configuration, given by setting the entanglement wedge E A to the cyclic permutation C (n) and its complement to the identity, provided the entanglement wedge E A is unique. For simplicity, we assume that D xy = D x∂ = D (but it is easy to see that the same conclusions hold true if all bond dimensions are powers of a fixed integer), and we consider the equivalent spin model with energy where the g x are variables in Sym n , with x and y ranging over both bulk and boundary vertices, subject to the boundary conditions g x = C (n) in A and g x = 1 inĀ (cf. Section 8).
The first observation is that n − χ(g) is equal to the minimal number of transpositions (i.e., permutations that exchange only two indices) required to write a permutation g. This implies that d(g x , g y ) := n − χ(g −1 x g y ) defines a metric. In particular, it satisfies the triangle inequality. The second ingredient is that, by the integral flow theorem, we can decompose a maximal flow between A andĀ into edge-disjoint paths. Each path starts in A, ends inĀ, and by the max-flow/min-cut theorem there are |γ A | many such paths P 1 , . . . , P |γ A | . Now consider an arbitrary configuration {g x } that satisfies the boundary conditions. We can bound its energy by looking only at those edges that occur in one of the paths, resulting in the lower bound (C.1) Along each path P k , the first spin is assigned the cyclic permutation C (n) and the last spin the identity permutation 1. Therefore, the triangle inequality (invoked once for each path) implies that Note that the right-hand side is just the energy cost of the configuration where we assign C (n) to the spins in E A and 1 to all other spins. We claim that this is the unique minimal energy configuration. To see this, suppose that {g x } is an arbitrary configuration that achieves this energy cost. Case 1: The only permutations that appear in {g x } are C (n) and 1. Then the domain where g x = C (n) is a minimal cut between A andĀ, i.e., an entanglement wedge for A. Since we have assumed that the entanglement wedge is unique, it must be equal to E A . Thus {g x } is the configuration described above.
Case 2: The configuration {g x } contains some other permutations. Since it is a minimal energy configuration, both inequalities (C.1) and (C.2) above must be tight. The fact that the first inequality is tight means that if an edge is not contained in any of the paths P k then the configuration {g x } necessarily assigns the same permutation in Sym n to its endpoints. It follows that the first inequality remains tight if we modify the configuration {g x } by changing an entire domain from one permutation to another. For the second inequality, we can use the triangle inequality to see that the sequence of permutations in any path P k must always be of the form C (n) , . . . , C (n) , * * * , 1, . . . , 1, where * * * denotes a sequence of permutations that are neither C (n) nor 1. Indeed, if this were not the case then the energy cost of the corresponding path would be higher than (n − 1). But this implies that by either changing all other permutations to C (n) , or by changing all of them to 1, we obtain two distinct minimal energy configurations that only contain C (n) and 1. By case 1, this is a contradiction.

D Calculation of C 2n in Section 6
In this appendix, we will present the derivation from Eq. (6.5) to Eq. (6.6) in Section 6. We first calculate M † M using the orthonormality condition (6.1). in which a, b, c, d are indices in the Hilbert space of A, and m, n, k, l are those in B. The best way of visualizing this calculation is by introducing a diagrammatic representation, as is shown in Fig. 13. In the trace of (M † M ) n , there are 2n copies of the density matrix ρ. The n contractions of A indices lead to pairwise permutations between pairs of density matrices 1 ↔ 2, 3 ↔ 4,. . . ,2n − 1 ↔ 2n. Similarly, the contractions of B indices lead to pairwise permutations between 2 ↔ 3, 4 ↔ 5, . . . , 2n ↔ 1. This concludes the proof of Eq. (6.6).

E Partition function of Ising model on the square lattice
In this appendix, we calculate the partition function of the Ising model in the large D (low temperature) limit on a 2D rectangular lattice of size M × N , with periodic boundary conditions along the first direction and open boundary conditions along the second one. We will use several results that can be found in [42]. As in the main text, we let 2β = log D.
Thus the leading order correction in the large D limit is M N D −4 . Now we turn to Z 1 (β), the partition function of the Ising model at temperature 1/β, with boundary pinning field down everywhere except for in a single interval of length L. Similarly as above, we denote the corresponding zero-temperature limit by Z ∞ 1 = Z 1 (β → ∞). Using the duality of Ising model, we know that Z 1 (β)/Z 0 (β) = S 0,0 S 0,L (β ), where e −2β = tanh β. Here, S 0,0 S 0,L (β ) denotes the two-point correlation function on the boundary of the dual lattice at temperature 1/β , whose analytical form is also provided in [42]. When L 1, we can expand Z 1 (β)/Z ∞ 1 to leading order in D and L, Thus the leading order correction is 2LD −1 .

F Average second Rényi entropy for 2-designs
In the following, we will show that (9.1) holds for an arbitrary 2-design in the limit of large bond dimension. We recall from Section 7 that the inequality S 2 (A) ≤ S(A) holds for arbitrary quantum states, while S(A) ≤ log rank ρ A ≤ log D|γ A | in any tensor network state. Therefore it remains to prove the lower bound on the average of the second Rényi entropy.
The first moments of the Haar measure are given by |V x V x | = I/D x , and so T = 1/ x D x . Together with our calculation in Sections 2 and 3 it follows that in the large D limit, where have introduced T = tr ρ and recall that Z 1 = tr ρ 2 A . We now bound the fluctuations in the trace T . Noting that T 2 = Z 0 , we can bound the variance of T /T as follows: where A 0 refers to the Ising action in its original form (2.13) and A to the simplified form (3.1) with constants removed. But any nontrivial spin configuration incurs an energy cost of at least log D, so that we obtain the upper bound By Chebyshev's inequality, it follows that, for any ε > 0, We now condition on the event that T /T ≥ 1 − ε. Writing X good for corresponding averages, we obtain the following bound using concavity of the logarithm, This implies both that K is a subgroup of G × H and that ϕ K is a group homomorphism. Thus L ≡ ϕ(K) is a (commutative) subgroup of the Pauli group; it follows that Thus |ψ B is indeed a subnormalized stabilizer state, as we set out to show. Since maximally entangled states are stabilizer states, it follows at once that a tensor network state (2.1) constructed by contracting stabilizer states |V x is again a stabilizer state.