# Random insights into the complexity of two-dimensional tensor network
calculations

Sofía González-García  Perimeter Institute for Theoretical Physics, Waterloo,
Ontario N2L 2Y5, Canada  University of Waterloo, Waterloo, Ontario, N2L 3G1,
Canada  Google Quantum AI, Santa Barbara, CA 93111, USA  Shengqi Sang
Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada
University of Waterloo, Waterloo, Ontario, N2L 3G1, Canada  Timothy H. Hsieh
Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada
Sergio Boixo  Google Quantum AI, Santa Barbara, CA 93111, USA  Guifré Vidal
Google Quantum AI, Santa Barbara, CA 93111, USA  Andrew C. Potter  Department
of Physics and Astronomy, and Quantum Matter Institute, University of British
Columbia, Vancouver, BC, Canada V6T 1Z1  Romain Vasseur  Department of
Physics, University of Massachusetts, Amherst, MA 01003, USA

###### Abstract

Projected entangled pair states (PEPS) offer memory-efficient representations
of some quantum many-body states that obey an entanglement area law, and are
the basis for classical simulations of ground states in two-dimensional (2d)
condensed matter systems. However, rigorous results show that exactly
computing observables from a 2d PEPS state is generically a computationally
hard problem. Yet approximation schemes for computing properties of 2d PEPS
are regularly used, and empirically seen to succeed, for a large subclass of
(‘not too entangled’) condensed matter ground states. Adopting the philosophy
of random matrix theory, in this work we analyze the complexity of
approximately contracting a 2d random PEPS by exploiting an analytic mapping
to an effective replicated statistical mechanics model that permits a
controlled analysis at large bond dimension. Through this statistical-
mechanics lens, we argue that: $i)$ although approximately sampling wave-
function amplitudes of random PEPS faces a computational-complexity phase
transition above a critical bond dimension, $ii)$ one can generically
efficiently estimate the norm and correlation functions for any finite bond
dimension. These results are supported numerically for various bond-dimension
regimes. It is an important open question whether the above results for random
PEPS apply more generally also to PEPS representing physically relevant ground
states.

Tensor network states (TNS) provide a compact representation of quantum states
with spatially-local entanglement, such as ground states of local
Hamiltonians. One-dimensional (1d) TNS, matrix product states (MPS), can be
efficiently contracted, enabling dramatic progress in studying 1d many-body
ground states White (1992); Fannes  _et al._ (1992); Verstraete  _et al._
(2008); Schollwöck (2011). Projected entangled pair states (PEPS) Verstraete
and Cirac (2004) are a higher dimensional generalization of MPS tensor
networks. Contrary to their 1d counterpart, contracting PEPS is a
$\\#\mathsf{P}$-complete task Schuch  _et al._ (2007) (the complexity class of
hard counting problems) even for simple square lattice PEPS with constant bond
dimension, $D$. Moreover, approximately contracting a 2d PEPS with bond
dimension $D\sim{\rm poly}(L)$ that scales polynomially in the linear
dimension $L$ of the system was shown to be average-case hard even for
calculating simple physical quantities like expectation values of local
observables Haferkamp  _et al._ (2020). Yet, in contrast to these complexity
results, the practical experience of PEPS practitioners Corboz and Mila
(2014); Corboz  _et al._ (2014); Niesen and Corboz (2017); Zheng  _et al._
(2017); Ponsioen  _et al._ (2019); Chen  _et al._ (2020) is that the standard
algorithm Verstraete and Cirac (2004); Jordan  _et al._ (2008) (which we
review below) for approximating the physical properties of PEPS seems to work
efficiently for finitely-correlated ground states of 2d lattice models in many
condensed matter problems.

(a) (b) (c)

Figure 1:  Boundary MPS contraction of 2d PEPS: (a) Tensor network
representation of the norm $\mathcal{N}=\braket{\Psi}{\Psi}$, for a finite
$L_{x}\times L_{y}$ PEPS $\ket{\Psi}$ (top) and one iteration of the infinite
bMPS approximate contraction algorithm (bottom): a row of $TT^{*}$ tensors is
multiplied by the boundary MPS whose bond dimension is then truncated back to
$\chi$. (b) Schmidt coefficients at a bipartition of the infinite bMPS for
$\chi=72$ for a clean random iPEPS for bond dimension $D\in\\{2,3,4,5,6\\}$.
(c) Change $\Delta\varrho(\chi)$ in the physical single-site density matrix
$\varrho(\chi)$ as a function of $\chi$, where we take $\varrho(\chi_{\rm
max}=72)$ as a reference.

Systematically clarifying the behavior of PEPS and their approximate
calculation efficiency is crucial not only for computational physics, but also
for refining the targets for quantum computational advantage in materials and
chemistry simulation problems Lee  _et al._ (2023), which are one of the
leading prospective applications for quantum computers Ma  _et al._ (2020). A
key challenge is a lack of systematic analytical tools for analyzing the
complexity of approximate contraction schemes. To this end, we adopt the
philosophy of statistical mechanics and random matrix theory which show that,
in many instances Guhr  _et al._ (1998), the statistical properties of random
ensembles can be systematically understood, even when individual cases cannot
be directly analyzed.

With this motivation in mind, we study the “typical” contraction complexity of
2d PEPS with _random_ tensors. This enables us to systematically address the
statistical complexity of random PEPS through complementary analytical and
numerical approaches. Analytically, we exploit a replica-trick-based mapping
between entanglement features of random disordered PEPS and partition
functions of classical lattice-“spin” models. We present evidence that
standard PEPS contraction schemes can efficiently approximate correlations of
a 2d random PEPS, including both of local observables and many non-local
observables such as string order parameters. This contrasts the behavior of
computing global properties such as overlaps
$\langle\Psi^{\prime}|\Psi\rangle$ between two PEPS states, which includes
individual wave-function amplitudes $\langle s_{1},s_{2},\dots|\Psi\rangle$,
where one expects a complexity-phase-transition tuned by $D$ Vasseur  _et al._
(2019), which has been observed numerically Levy and Clark (2021); Yang  _et
al._ (2022). We emphasize that the arguments presented are not rigorous
proofs, but rather physical arguments about the replica statistical mechanics
that are grounded in standard paradigms of statistical mechanics and critical
phenomena.

To validate the predictions of this analytic derivation, we numerically
explore the $D$ dependence of contracting both “clean” (translation-invariant)
infinite PEPS (iPEPS) Jordan  _et al._ (2008), as well as “clean” and
“disordered” (see below) stabilizer (Clifford) PEPS, which allow for a much
larger bond dimension. The numerical results are consistent with the analytic
predictions, including quantitative agreement in the detailed large-$D$
asymptotic behavior. While the analytical findings are limited to disordered
PEPS, the numerical results suggest that the predictions also hold for “clean”
translation-invariant PEPS.

It remains an important open question to determine the relevance of these
random PEPS results to PEPS representing ground states of Hamiltonians
relevant to condensed matter, materials science, and chemistry. For example,
it is known that despite being highly entangled, at large $D$ the random
square-lattice PEPS exhibit only short-range correlations for local
observables Lancien and Pérez-García (2021) (see also Appendix C.5), whereas
the ground-states of physical systems can exhibit relatively long correlation
lengths, potentially making them harder to contract. On the other hand, the
statistical mechanics mapping applies more broadly to a large class of network
structures, including examples with arbitrarily-long-range correlations (see
Appendix D.2).

##### PEPS and boundary MPS method

We begin by briefly reviewing standard schemes to approximately contract 2d
PEPS. Consider a 2d square lattice made of $L_{x}\times L_{y}\equiv N$ sites,
$L_{x}\leq L_{y}$, where each site is represented by a $d$-dimensional vector
space $\mathbb{C}^{d}$, and a bulk state
$\ket{\Psi}\in\bigotimes_{i=1}^{N}\mathbb{C}^{d}$ of the lattice that accepts
a PEPS representation of the form
$|\Psi\rangle=\sum_{s_{1},\dots,s_{N}=1}^{d}\mathcal{C}\left[T_{[1]}^{s_{1}}\dots
T_{[N]}^{s_{N}}\right]|s_{1}\dots s_{N}\rangle$. Here, $T_{[r]}$ denotes the
5-index PEPS tensor on lattice site $r\in\\{1,2,\cdots,N\\}$, with complex
components $\left(T_{[r]}\right)^{s}_{ijkl}$, where $s=1,\cdots,d$ and
$i,j,k,l=1,\cdots,D$ are known as physical and bond indices, respectively, and
$\mathcal{C}$ denotes contraction of the bond indices connecting nearest
neighbor tensors.

An archetypal PEPS computation is the contraction of the 2d tensor network in
Figure 1a, of size $L_{x}\times L_{y}$, which corresponds to the overlap
$\mathcal{N}\equiv\braket{\Psi}{\Psi}$. Other closely related 2d networks,
such as those corresponding to $k$-point correlators of local operators
$\langle\Psi|O_{1}O_{2}\dots O_{k}|\Psi\rangle$, are contracted similarly. An
exact contraction incurs a computational cost $O(L_{y}\exp(L_{x}))$, but a
number of efficient algorithms exist, most notably the boundary MPS (bMPS)
method Verstraete and Cirac (2004); Jordan  _et al._ (2008), which scales
linearly both in $L_{x}$ and $L_{y}$ at the price of introducing
approximations111The corner transfer matrix (CTM) Baxter (2007) method is
another prominent approximate contraction scheme for the 2d tensor network
$\mathcal{N}$. CTM and bMPS are closely related and similar conclusions
regarding boundary entanglement, as discussed in this paper, apply to both..
The bMPS method considers the virtual 1d lattice made of $L_{x}$ sites, each
described by a $D^{2}$-dimensional vector space $\mathbb{C}^{D^{2}}$
(corresponding to two PEPS bond indices), obtained from a horizontal cut of
the 2d tensor network $\mathcal{N}$, and a so-called boundary state
$|\psi\rrangle\in\bigotimes_{i=1}^{L_{x}}\mathbb{C}^{D^{2}}$ that accepts an
MPS representation,
$|\psi\rrangle=\sum_{\alpha_{1},\cdots,\alpha_{L_{x}}}\mathcal{C}\left(B_{[1]}^{\alpha_{1}}\cdots
B_{[L_{x}]}^{\alpha_{L_{x}}}\right)|\alpha_{1}\cdots\alpha_{L_{x}}\rrangle$.
Here the double brackets $|\rrangle$ emphasize that the boundary state is
actually a vectorized density operator on the PEPS bond space, $B_{[x]}$ is a
3-index MPS tensor on site $x\in\\{1,2,\cdots,L_{x}\\}$ with complex
components $\left(B_{[x]}\right)^{\alpha}_{\beta\gamma}$, where
$\alpha=1,\cdots,D^{2}$ and $\beta,\gamma=1,\cdots,\chi$, and $\chi$ is the
bMPS bond dimension. In order to approximately contract the 2d tensor network
$\mathcal{N}$, which we now regard as made of $L_{y}$ rows of tensors where
each row is labelled by an integer $y=1,\cdots,L_{y}$, we first represent the
boundary state $|\psi(1)\rrangle$ corresponding to the top row ($y=1$) exactly
as a bMPS with bond dimension $D^{2}$, which we then truncate down to $\chi$
(if $D^{2}>\chi$). For increasing value of $y=1,2,\cdots$, given a (possibly
truncated, approximate) bMPS representation of the boundary state
$|\psi(y)\rrangle$, we produce an approximate bMPS for $|\psi(y+1)\rrangle$ by
contracting the row $y+1$ with the existing bMPS, see Figure 1a, then
truncating the resulting bond dimension $\chi^{\prime}$ down to $\chi$ (if
$\chi^{\prime}>\chi$). Finally, the value of $\mathcal{N}$ is obtained by
contracting the bMPS for $|\psi(L_{y}-1)\rrangle$ with the bottom row (i.e.
$y=L_{y}$). The leading computational cost of the bMPS method with respect to
index dimensions $d,D,\chi$ is $O(D^{4}\chi^{3}+dD^{6}\chi^{2})$ Lubasch  _et
al._ (2014) (or $O(D^{6}\chi^{3}+dD^{8}\chi^{2})$ in the infinite case Jordan
_et al._ (2008), $L_{x},L_{y}\to\infty$).

The above sequential contraction of the 2d tensor network $\mathcal{N}$ is
analogous to a dynamical evolution of the state of a 1d system (represented by
the bMPS) by a transfer matrix implementing a completely positive map made of
a row of $TT^{*}$ tensors. This setting is similar to 1d noisy quantum
circuits Noh  _et al._ (2020); Li  _et al._ (2023).

##### Random PEPS models

We consider applying the boundary MPS method to random PEPS, whose tensors $T$
are multi-dimensional arrays of complex entries, with each entry sampled
independently and identically distributed (i.i.d.) from the complex normal
(Gaussian) distribution with zero mean and unit variance. We consider two
different ensembles of random PEPS: i) “clean”, translation-invariant PEPS
with a single random instance of $T$ used for each site, and ii) “disordered”
PEPS with a distinct, random $T_{[r]}$ for each site, $r$. From the
perspective of spatial symmetries, clean and disordered random PEPS are akin
to crystalline and disordered materials respectively.

##### Boundary MPS entanglement and random PEPS contraction

The resulting value of $\mathcal{N}=\braket{\Psi}{\Psi}$ is approximate due to
the errors introduced at each truncation typically implemented by preserving
the $\chi$ largest Schmidt coefficients $\\{\sigma_{i}\\}$ (or coefficients
$\\{\sigma_{i}^{2}\\}$ in the entanglement spectrum) assigned to each bond of
the bMPS (see Appendix A). Intuitively, if most of the $O(D^{2y})$ Schmidt
coefficients (for $y\leq L_{x}/2$) are ‘large’, then an accurate estimate of
$\mathcal{N}$ requires using a bMPS with a bond dimension $\chi$ that grows
exponentially in $y$, and thus the method is not efficient. However, if only a
small number of Schmidt coefficients are significant, then retaining only a
constant value of $\chi$ may suffice for an accurate bMPS representation, thus
enabling efficiently approximate contraction of the network.

Figure 1b shows the entanglement spectrum $\\{\sigma_{i}^{2}\\}$ of the fixed-
point bMPS $\left|\tilde{\psi}\right\rrangle$ for clean random iPEPS, where we
exploit translation invariance to access the $L_{x},L_{y}\to\infty$ limit, and
$\left|\tilde{\psi}\right\rrangle$ represents the large-$y$ fixed-point of the
bMPS sequence $\\{\left|\psi(y)\right\rrangle\\}$ (see Appendix A for further
details). Notice the exponential decay of $\sigma_{i}^{2}$ as a function of
$i$, which implies that an accurate approximation of the boundary state can be
obtained with a small bond dimension $\chi$. In turn, this results in an
accurate approximation $\varrho(\chi)$ to the reduced density matrix
$\varrho=\text{tr}_{\bar{R}}|\Psi\rangle\langle\Psi|/\mathcal{N}$ on a local
region $R$ of the 2d lattice and thus also to the expectation value
$\braket{\Psi}{O}{\Psi}/\mathcal{N}=\text{tr}(\varrho
O)\approx\text{tr}(\varrho(\chi)O)$ for any local observable $O$ supported on
$R$. Indeed, consider the distance $\Delta\varrho(\chi)$ between the density
matrix $\varrho(\chi)$ and a reference density matrix $\varrho(\chi_{\rm
max})$ as given by the largest eigenvalue of $|\varrho(\chi)-\varrho(\chi_{\rm
max})|$, where we make the key assumption that $\varrho(\chi_{\rm max})$ is a
good approximation to the exact $\varrho=\varrho(\chi\to\infty)$. As shown in
Figure 1c for a single-site $\varrho$, $\Delta\varrho(\chi)$ decays
exponentially with $\chi$, with similar behavior observed when $R$ consists of
more than one site (see Figure 12b). We therefore conclude that an accurate
approximation to $\braket{\Psi}{O}{\Psi}/\mathcal{N}$ can be obtained
efficiently, for these instances of clean random iPEPS.

In practice, the bMPS entanglement is often characterized using the Renyi
entropy of order $n$,

$\displaystyle
S^{(n)}_{A}=\frac{1}{1-n}\log\left(\text{tr}\left[\left(\frac{\rho_{A}}{\text{tr}\rho}\right)^{n}\right]\right),~{}~{}~{}~{}$
(1)

where $\rho\equiv|\psi\rrangle\llangle\psi|$ and
$\rho_{A}=\text{tr}_{\bar{A}}|\psi\rrangle\llangle\psi|$ is the reduced
density matrix for a 1d sub-region, $A$, of the evolved bMPS, and the
denominator is included to explicitly normalize the state. Specifically, for
$n\rightarrow 1$, we expect an inefficient contraction when $|\psi(y)\rrangle$
obeys an entanglement volume law, where $S^{(n)}_{A}$ is proportional to the
size $|A|$ of region $A$, and an efficient one when it obeys an entanglement
area law, where $S^{(n)}_{A}$ is upper bounded by a constant. An advantage of
using a single statistic, such as $S^{(n)}_{A}$, rather than the full
entanglement spectrum is that its ensemble average for $n\geq 1$ can be
evaluated using a statistical mechanics model, as we will see below.

Figure 2:  Statistical Mechanics Model (schematic) – The bMPS entanglement of
an random PEPS maps onto the difference between free-energies
$\mathcal{F}_{A}-\mathcal{F}_{0}$ that differ in boundary conditions as shown
(note in the replica limit free-energy and partition functions coincide). Each
site of the random PEPS is replaced by a replica-permutation “spin” (shown as
colored circles; while there are $(2Q)!$ different spin values, we show only
the three corresponding to the bulk and boundary fields here). Solid lines
along bonds of the random PEPS represent ferromagnetic interactions. Lines
ending in a dot denote single-site fields. At large-$D$, the dominant
configuration is uniform polarization along $e$ (gray). Fluctuating domains of
spins have both a linear boundary tension (dashed lines) due to the
interactions, and a surface-area tension due to the bulk $e$ fields.

##### Statistical mechanics model mapping

The entanglement features of random tensor network contractions can be mapped
onto free-energy cost of boundary domains in a classical statistical mechanics
(stat-mech) model Hayden  _et al._ (2016); Vasseur  _et al._ (2019); Lopez-
Piqueres  _et al._ (2020); Nahum  _et al._ (2021). This mapping has been used
extensively to study both random holographic tensor networks (tensor networks
with physical legs only at the boundary) Hayden  _et al._ (2016); Vasseur  _et
al._ (2019), and to explore entanglement growth Nahum  _et al._ (2017, 2018);
Zhou and Nahum (2019), measurement-induced phase transitions Skinner  _et al._
(2019); Li  _et al._ (2018); Jian  _et al._ (2020); Bao  _et al._ (2020);
Potter and Vasseur (2022); Fisher  _et al._ (2023); Napp  _et al._ (2022) in
random quantum circuits and channels Li  _et al._ (2023). Here, we adapt it to
study the entanglement of the evolved boundary state involved in the
contraction of a random, disordered PEPS. In the main text, we merely
summarize the key elements of the stat-mech mapping, and refer readers to the
supplemental material in the Appendix for details. We focus on the calculation
of the PEPS norm, $\mathcal{N}=\langle\Psi|\Psi\rangle$, for a disordered,
random, square lattice PEPS $|\Psi\rangle$, and will later comment to how this
is modified for other observables. As above, we denote the density matrix of
the evolved bMPS in the PEPS contraction procedure described above as
$\rho=|\psi\rrangle\llangle\psi|$.

The mapping exploits a replica trick, $\log f=\lim_{m\rightarrow
0}\frac{1}{m}(f^{m}-1)$, to perform averages of $n^{\rm th}$-Renyi entropies
(Equation 1). The objects of interest are averages over powers of the reduced
density matrix of the boundary MPS for a region, $A$:
$\mathbb{E}\left[\left({\rm tr}\rho_{A}^{n}\right)^{m}\right]$, where
$\mathbb{E}[\dots]$ denotes averaging over the Gaussian-distribution of tensor
entries. For integer $m,n$, this average can be carried out using Wick’s
theorem. These quantities contain $Q=nm$ copies of $\rho$. Since each copy of
$\rho$ contains $2\times$ $T$ and $2\times$ $T^{*}$ tensors, this results in a
sum over all possible Wick pairings of the $2Q=2mn$ replicas of the tensor
network. A given Wick pairing is defined by a set of permutation elements
$g_{r}\in S_{2Q}$ where $S_{2Q}$ is the symmetric group (group of
permutations) on $2Q$ elements, for each site, $r$. The resulting entanglement
entropy can be written as a free-energy difference Vasseur  _et al._ (2019):

$\displaystyle S^{(n)}(A)=\lim_{m\rightarrow
0}\frac{1}{m(n-1)}(\mathcal{F}_{A}-\mathcal{F}_{0}),$ (2)

where $\mathcal{F}_{A,0}=-\log\mathcal{Z}_{A,0}$ are the free energies of a
classical statistical mechanics model with partition function
$\mathcal{Z}=\sum_{\\{g_{r}\in S_{2Q}\\}}e^{-H[g]}$, defined by Hamiltonian:

$\displaystyle\mathcal{H}[g]$ $\displaystyle=-\sum_{\langle
rr^{\prime}\rangle}JC(g_{r},g_{r^{\prime}})-\sum_{r,g\in
S_{2Q}}h_{g,r}C(g,g_{r}).$ (3)

The labels $A$ and $0$ in the free energies $\mathcal{F}_{A,0}$ refer to
different boundary fields that will be specified below. Here, $\langle
rr^{\prime}\rangle$ denotes neighboring pairs of sites, the function
$C(g,g^{\prime})=C(g^{-1}g^{\prime})$ counts the number of cycles in the
permutation $g^{-1}g^{\prime}$, and the coupling constant is related to bond
dimension: $J=\log D$. Intuitively, $C(g,g^{\prime})$ measures how similar the
permutations $g$ and $g^{\prime}$ are to each other, with $C(g,g)=2Q$, and
minimal value $C\geq 1$. The coupling between neighboring spins is
ferromagnetic, favoring aligned spins. This competes with entropic
fluctuations of the spins.

The second term can be thought of as local “fields” that point in a different
“directions” (in $S_{2Q}$) in the bulk and boundary regions. For
$\mathcal{F}_{A}$, the fields are:

$\displaystyle h_{g,r}=\begin{cases}\delta_{g,e}\log d&r\in{\rm bulk}\\\
\delta_{g,g_{A}}\log D&r\in\text{boundary region $A$}\\\ \delta_{g,g_{0}}\log
D&r\in\text{boundary region $\bar{A}$}\end{cases}$ (4)

where $e$ is the identity permutation, and $g_{0,A}$ are permutations
corresponding to the trace and cyclic permutation boundary conditions in the
entanglement region $A$ and its complement, $\bar{A}$, respectively (see
Appendix C for details). In contrast, the boundary fields in $\mathcal{F}_{0}$
are uniform at the boundary: $h_{g,{\rm boundary}}=\delta_{g,g_{0}}\log D$
[corresponding to Equation 4 with only boundary-region $\bar{A}$].

##### Global observables: Complexity transition

For computing global properties of a random PEPS, such as individual wave-
function components, $\langle s_{1}\dots s_{N}|\Psi\rangle$, or overlaps
between distinct random PEPS, $\langle\Psi^{\prime}|\Psi\rangle$, the tensors
are not completely positive, and tensor bulk $h_{e}$ fields vanish. Absent the
bulk $h\sim e$ terms, the model has $Q$ copies of each tensor $T$ on each
site, and a (bulk) left/right $S_{Q}\times S_{Q}$ symmetry. There are two
possible phases: a disordered phase with intact symmetry, and an ordered phase
where this symmetry is spontaneously broken and the spins polarize towards a
spontaneously chosen $g\in S_{Q}$. This model also describes the entanglement
of random holographic tensor networks Vasseur  _et al._ (2019), in which
physical degrees of freedom arise only at the boundary of a higher-dimensional
tensor network with tensors having only virtual bond legs. In the ordered
phase the boundary fields representing the entanglement cut force a $g_{A}$ to
$g_{0}$ domain wall at the boundary, which costs energy proportional to the
length of the domain, corresponding to volume law scaling of entanglement. By
contrast, in the disordered phase, domain walls have a vanishing line tension,
leading to area-law entanglement scaling. The universal properties of this
complexity transition remain unsolved, but are closely related to those of
(forced) measurement-induced entanglement transitions in random unitary
circuit dynamics Nahum  _et al._ (2021). Numerical analysis Levy and Clark
(2021) for square PEPS indicates that the critical bond dimension for this
complexity transition is close to $D_{c}\approx 2$, meaning that essentially
all practical calculations are expected to lie on the hard (volume-law) side
of the transitions.

##### Local observables: area-law bMPS

For local observables such as the wave-function norm,
$\langle\Psi|\Psi\rangle$, or $k$-point correlators of local operators:
$\langle\Psi|O_{1}O_{2}\dots O_{k}|\Psi\rangle$, the situation is rather
different. Here, the bulk identity-permutation ($h_{e}$) fields are present,
and explicitly break the $S_{2Q}$ symmetry. This erases the distinction
between the ordered and disordered phases of the model, and as we will now
argue leaves only the area-law phase behind. In fact, we will see that this
conclusion holds even for many non-local observables, such as string-order
parameters or Wilson loop observables, so long as they lack support on some
finite density of sites.

The boundary fields, $g_{A,0}$, associated with the entanglement region $A$,
and its complement region $\bar{A}$, compete with the bulk $e$ fields. The
resulting competition is straightforward to analyze at large $D$, where domain
wall configurations are strongly suppressed, and the partition functions are
dominated by the lowest energy configuration. We refer to this approximation
as “min-cut”. When the $L_{y}$ is large (to approach the thermodynamic limit,
we are interested in $L_{y}\rightarrow\infty$), the number of bulk fields
(which scales with system volume) overwhelms that of the boundary fields
(which scale with linear dimension), forcing the spins to polarize along $e$.
Consequently the linear free-energy per unit length of the boundary regions is
proportional to $C(g_{A,0},e)\log D$ in the $A,\bar{A}$ regions respectively.
Crucially, while the bulk fields break explicitly break the $S_{2Q}$, there is
a residual bulk $S_{Q}$ symmetry that guarantees that $C(g_{A},e)=C(g_{0},e)$,
i.e. that the leading contribution to the boundary-field energies cancels in
$\mathcal{F}_{A}-\mathcal{F}_{0}$, leading to area-law scaling.

In fact, in the $D\rightarrow\infty$ limit, this cancellation is exact, and
there is strictly zero operator entanglement. Generically, corrections to the
min-cut approximation for $S_{A}$ come from domains that span across the
entanglement cut between $A$ and $\bar{A}$. Since these fluctuations are
massive, (the probability of getting a large domain of size $\ell$ is expected
to be exponentially small in $\ell$), these fluctuations are only sensitive to
the local change in boundary fields between $A$ and $\bar{A}$, and hence can
give only an area-law contribution to $S_{A}$. In Appendix C.4, we estimate
that, beyond the infinite-$D$ limit, the leading corrections to the area-law
entanglement come at order $1/D^{2}$.

We remark that these results highlight a counterintuitive feature of random
tensor networks. While highly-correlated states require large bond dimension,
_typical_ large-$D$ PEPS have rather short-range correlations (despite being
quite entangled), and the set of highly-correlated large-$D$ PEPS is in fact
rare (by the measure of our Gaussian-random ensemble). See Appendix D for a
more detailed discussion of potential differences between random PEPS and
physical groundstates relevant to condensed matter. In Appendix D.2, we also
discuss different tensor network geometries for which _typical_ states exhibit
longer range correlations.

We emphasize that, while the calculations described were done at leading order
in the large-$D$ limit, the general symmetry principles described above
suggest that the prediction of area-law entanglement for the boundary-MPS hold
also at any $D$ and that there is no expected singular change (phase
transition) in the entanglement as a function of $D$.

Correlators and overlaps –  While we have focused on the computation of the
PEPS norm, this is closely related to the calculation of correlation functions
of local operators: $\langle\Psi|O_{1}O_{2}\dots O_{k}|\Psi\rangle$. The
insertion of local operators amounts to replacing the tensors on those sites
with non-positive tensors, which in the statistical mechanics model
corresponds to removing the local $e$-field on that site. Yet, so long as a
finite-density of sites are unaffected by the operator insertions, there
remains a net bulk $e$-field that explicitly breaks the replica symmetry
breaks the replica-permutation symmetry and removes the complexity phase
transition. Hence, we expect our “easiness” result to a very large class of
observables including even exotic non-local string correlation functions to
used detect topological phases and confinement in gauge theories.

Figure 3:  Entanglement barrier  – Mid-cut entanglement entropy, $S(y)$, of
the bMPS $|\psi(y)\rrangle$ for disordered (left) and clean (right) random
stabilizer PEPS as a function of layer number $y$, for various choices of
$(D,d)=(173^{k_{D}},173^{k_{d}})$. A schematic representation of the dominant
spin configurations at large-$D$ in the corresponding stat-mech model are
shown underneath the plot. $S(y)$ initially increases linearly up to
lengthscale $\ell_{*}\sim\log D/\log d$ (corresponding to a growing vertical
domain wall) before dropping to a small constant for $y\gtrsim\ell_{*}$
(corresponding to the bulk fields expelling the domain wall out to the
system’s boundary). In both cases, the simulated $S(y)$ agrees well with the
stat-mech model in the large-$D$ regime (solid gray lines), with very small
instance to instance variance.

##### Entanglement barrier and random stabilizer PEPS

We have so far considered the thermodynamic limit of the boundary-MPS
entanglement, where the statistical mechanics model predicts that the evolved
boundary MPS eventually reaches an area-law steady state. We can also use the
stat-mech mapping to examine the transient dynamics of the boundary MPS at
early stages in contraction of a PEPS with open boundaries at $y=0$. When $y$
is small, $y\ll\ell^{*}=\log D/\log d$, for $D>d$, the boundary fields
dominate and the spins $g_{r}$ polarize along the direction of the fields at
their closest boundary. For $\mathcal{F}_{A}$, this adds an extra
$g_{A}$-to-$g_{0}$ domain will along the y-direction with free energy cost
$\sim y$ (see Figure 3). As a result, the statistical mechanics model predicts
that entanglement grows as:

$\displaystyle S^{A}(y)\sim\begin{cases}y\log D^{2}&y\ll\ell^{*}\\\ {\rm
constant}&y\gg\ell^{*}\end{cases},~{}~{}~{}~{}~{}~{}\ell^{*}\equiv\log D/\log
d,$ (5)

achieving a maximum of $S_{\rm max}\sim\log^{2}D/\log d$ at
$y\approx\ell^{*}=\log D/\log d$, corresponding to a maximum bMPS bond
dimension of $\chi_{\rm max}\sim\exp\left[\log^{2}D/\log d\right]$. We refer
to this phenomenon as an “entanglement barrier”. We note that a similar
entanglement barrier arises in the classical simulation of random quantum
circuits in the presence of decoherence Noh  _et al._ (2020); Li  _et al._
(2023), which are relevant to “quantum supremacy” experiments Martinis et al’
(2019). For circuits, the entanglement barrier reflects the initial build up
of quantum correlations and entanglement before noise and errors overwhelm the
system making it essentially indistinguishable from a maximally-mixed state.
The evolution of the boundary MPS is similar to the dynamics of circuits with
decoherence (even though the bMPS evolution is not generally trace-
preserving), which provides a physical picture for the low (area law)
entanglement of the MPS steady state Noh  _et al._ (2020); Li  _et al._
(2023). In this correspondence, the PEPS physical dimension $d$ represents the
strength of noise, and the barrier height decreases with increasing noise
strength.

We remark that, in practice, the entanglement barrier may not be an obstacle
for extracting bulk properties of a PEPS. As illustrated by the iPEPS numerics
above, it may not be necessary to accurately track the transient evolution of
the boundary MPS through the entanglement barrier. Instead, one may be able to
directly find an accurate representation of the area-law steady state bMPS,
with $\chi$ reflecting the much smaller steady-state entropy deep in the bulk.
However, for select tasks, such as simulating edge states of topological
materials, the entanglement barrier may have practical consequences. We also
note that this prediction is consistent with the average-case hardness proof
for computing norms of finite-size random PEPS Haferkamp  _et al._ (2020) with
$D\sim{\rm poly(L)}$: in the stat-mech description this corresponds to an
entanglement barrier with $\chi_{\rm max}$ scaling in a (very-weakly) super-
polynomial fashion with $L$.

We can also use the entanglement barrier prediction simply as a test for the
validity of the min-cut approximation and replica limit for the statistical
mechanics model. To easily access large-$D$, we focus on random stabilizer
(Clifford) PEPS, which can be efficiently simulated via the stabilizer
formalism with $\sim{\rm poly}(\log D$) complexity (see the Appendix B for the
definition of stabilizer PEPS and the algorithm for contracting them). For
stabilizer PEPS, each of $T_{i}$’s bond dimensions needs to be an integer
power of some prime number $p$: $d=p^{k_{d}},~{}D=p^{k_{D}}$. In simulations,
we reach the large-$D$ regime by picking $p=173$ (and confirm that
qualitatively similar results also hold for $p=2$).

We performed simulations of both disordered and clean stabilizer PEPS with
periodic boundary conditions, and show results in Figure 3. The simulation
results agree quantitatively with the theoretical prediction in Equation 5,
suggesting the validity of the statistical mechanics model in capturing the
entanglement features of the boundary state. The agreement for the clean case
is noteworthy since the statistical mechanics model was derived for disordered
PEPS, and suggests that the conclusions of the stat-mech model also apply to
translation invariant PEPS.

##### Discussion

We have presented numerical and analytical evidence that the complexity of
(approximate) 2d random PEPS calculations greatly depends on the nature of the
object to be contracted. On the one hand, overlaps between a 2d random PEPS
and a product state $\braket{s_{1}\cdots s_{N}}{\Psi}$ or, more generally,
overlaps $\braket{\Psi^{\prime}}{\Psi}$ between two distinct 2d random PEPS,
face a complexity phase transition as a function of bond dimension. On the
other hand, the overlap $\braket{\Psi}{\Psi}$ of a single random PEPS as well
as (computationally related) physically relevant properties such as the
expectation value $\braket{\Psi}{O}{\Psi}$ or correlation functions
$\braket{\Psi}{O_{1}O_{2}\cdots O_{k}}{\Psi}$ of local operators, can be
efficiently computed. The statistical mechanics model provides an intuitive
picture for this result.

The chief objective hazard of our analytic analysis is the reliance on the
replica trick: away from the infinite-$D$ limit, we are able to analyze
properties only for positive integer number of replicas $n=2,3,\dots$. By
contrast (see Appendix A), errors in observables are best captured by the von-
Neumann entanglement entropy (which can be obtained by taking the limit of
Renyi index $n\rightarrow 1^{+}$). We note that stabilizer PEPS have a fine-
tuned, completely flat entanglement spectrum and do not constitute an
independent check of the replica limit convergence. Therefore, in Appendix A,
we numerically explore finite-$\chi$ scaling of the $n$ dependence of
bipartite Renyi entropies for clean random iPEPS. We observe a smooth and
well-converged Renyi-index dependence up to the largest accessible $D=6$,
implying no signs of trouble for the replica limit.

Our statistical mechanics mapping, when generalized to a 3d random PEPS
$\ket{\Phi_{3d}}$, predicts again an area-law entangled boundary state for the
overlap $\braket{\Phi_{3d}}{\Phi_{3d}}$ and expectation value and correlation
functions of local observables. This 2d area-law entangled boundary state
might again be efficiently approximated with a 2d PEPS $\ket{\Psi_{2d}}$. It
is thus plausible that approximate expectation values and correlators of local
observables can be efficiently evaluated also for 3d random PEPS.

An important question is to what extent the present results obtained for
random PEPS may also apply to the contraction of PEPS that represent ground
states of interest in condensed matter, materials science or quantum
chemistry. Elucidating such a question would be useful both for the classical
simulation of such systems, as well as to determine whether quantum computers
may provide an exponential speed-up for such problems Lee  _et al._ (2023). We
note that, even if classical TNS methods could efficiently calculate
properties of finitely-correlated states in 2d and 3d, many challenging tasks
remain ripe for quantum advantage, including simulating highly-entangled
states such as metals or phase-transitions, or calculating dynamical responses
such as transport coefficients or optical spectra. Progress in addressing this
question will likely require a combination of analytical work and numerical
simulations directly targeting “physically-relevant” systems of interest.

Acknowledgments.—We thank Sarang Gopalakrishnan, Andreas Ludwig, Yi-Zhuang You
and Yuri D. Lensky for insightful discussions. We acknowledge support from the
US Department of Energy, Office of Science, Basic Energy Sciences, under Early
Career Award No. DE-SC0019168 (R.V.) and DE-SC0022102 (A.C.P.), as well as the
Alfred P. Sloan Foundation through Sloan Research Fellowships (A.C.P. and
R.V.). R.V., A.C.P., and G.V. thank the Kavli Institute of Theoretical Physics
(KITP) for hospitality. KITP is supported in part by the National Science
Foundation under Grant No. NSF PHY-1748958. S.G.G., S.S., T.H. and G.V.
acknowledge support by the Perimeter Institute for Theoretical Physics (PI),
Natural Sciences and Engineering Research Council of Canada (NSERC), and
Compute Canada. Research at PI is supported in part by the Government of
Canada through the Department of Innovation, Science and Economic Development
and by the Province of Ontario through the Ministry of Colleges and
Universities. G.V. is a CIFAR associate fellow in the Quantum Information
Science Program, a Distinguished Invited Professor at the Institute of
Photonic Sciences (ICFO), and a Distinguished Visiting Research Chair at
Perimeter Institute.

## References

  * White (1992) Steven R. White, “Density matrix formulation for quantum renormalization groups,” Physical Review Letters 2863, 69 (1992).
  * Fannes  _et al._ (1992) M. Fannes, B. Nachtergaele, and R. F. Werner, “Finitely correlated states on quantum spin chains,” Communications in Mathematical Physics 144, 443 – 490 (1992).
  * Verstraete  _et al._ (2008) F. Verstraete, V. Murg, and J.I. Cirac, “Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems,” Advances in Physics 57, 143–224 (2008).
  * Schollwöck (2011) Ulrich Schollwöck, “The density-matrix renormalization group in the age of matrix product states,” Annals of Physics 326, 96–192 (2011).
  * Verstraete and Cirac (2004) F. Verstraete and J. I. Cirac, “Renormalization algorithms for quantum-many body systems in two and higher dimensions,”  (2004).
  * Schuch  _et al._ (2007) Norbert Schuch, Michael M. Wolf, Frank Verstraete, and J. Ignacio Cirac, “Computational complexity of projected entangled pair states,” Physical Review Letters 98 (2007), 10.1103/physrevlett.98.140506.
  * Haferkamp  _et al._ (2020) Jonas Haferkamp, Dominik Hangleiter, Jens Eisert, and Marek Gluza, “Contracting projected entangled pair states is average-case hard,” Phys. Rev. Res. 2, 013010 (2020).
  * Corboz and Mila (2014) Philippe Corboz and Frédéric Mila, “Crystals of bound states in the magnetization plateaus of the shastry-sutherland model,” Physical Review Letters 112 (2014), 10.1103/physrevlett.112.147203.
  * Corboz  _et al._ (2014) Philippe Corboz, T. M. Rice, and Matthias Troyer, “Competing states in the t-j model: uniform d-wave state versus stripe state,” Physical Review Letters 113 (2014), 10.1103/physrevlett.113.046402.
  * Niesen and Corboz (2017) Ido Niesen and Philippe Corboz, “A tensor network study of the complete ground state phase diagram of the spin-1 bilinear-biquadratic heisenberg model on the square lattice,” SciPost Physics 3 (2017), 10.21468/scipostphys.3.4.030.
  * Zheng  _et al._ (2017) Bo-Xiao Zheng, Chia-Min Chung, Philippe Corboz, Georg Ehlers, Ming-Pu Qin, Reinhard M. Noack, Hao Shi, Steven R. White, Shiwei Zhang, and Garnet Kin-Lic Chan, “Stripe order in the underdoped region of the two-dimensional hubbard model,” Science 358, 1155–1160 (2017).
  * Ponsioen  _et al._ (2019) Boris Ponsioen, Sangwoo S. Chung, and Philippe Corboz, “Period 4 stripe in the extended two-dimensional hubbard model,” Physical Review B 100 (2019), 10.1103/physrevb.100.195141.
  * Chen  _et al._ (2020) Ji-Yao Chen, Sylvain Capponi, Alexander Wietek, Matthieu Mambrini, Norbert Schuch, and Didier Poilblanc, “$su(3)_{1}$ chiral spin liquid on the square lattice: a view from symmetric peps,” Physical Review Letters 125 (2020), 10.1103/physrevlett.125.017201.
  * Jordan  _et al._ (2008) J. Jordan, R. Orús, G. Vidal, F. Verstraete, and J. I. Cirac, “Classical simulation of infinite-size quantum lattice systems in two spatial dimensions,” Physical Review Letters 101 (2008), 10.1103/physrevlett.101.250602.
  * Lee  _et al._ (2023) Seunghoon Lee, Joonho Lee, Huanchen Zhai, Yu Tong, Alexander M Dalzell, Ashutosh Kumar, Phillip Helms, Johnnie Gray, Zhi-Hao Cui, Wenyuan Liu, _et al._ , “Evaluating the evidence for exponential quantum advantage in ground-state quantum chemistry,” Nature Communications 14, 1952 (2023).
  * Ma  _et al._ (2020) He Ma, Marco Govoni, and Giulia Galli, “Quantum simulations of materials on near-term quantum computers,” npj Computational Materials 6 (2020), 10.1038/s41524-020-00353-z.
  * Guhr  _et al._ (1998) Thomas Guhr, Axel Müller-Groeling, and Hans A Weidenmüller, “Random-matrix theories in quantum physics: common concepts,” Physics Reports 299, 189–425 (1998).
  * Vasseur  _et al._ (2019) Romain Vasseur, Andrew C. Potter, Yi-Zhuang You, and Andreas W. W. Ludwig, “Entanglement transitions from holographic random tensor networks,” Phys. Rev. B 100, 134203 (2019).
  * Levy and Clark (2021) Ryan Levy and Bryan K. Clark, “Entanglement Entropy Transitions with Random Tensor Networks,” arXiv e-prints , arXiv:2108.02225 (2021), arXiv:2108.02225 [cond-mat.stat-mech] .
  * Yang  _et al._ (2022) Zhi-Cheng Yang, Yaodong Li, Matthew P. A. Fisher, and Xiao Chen, “Entanglement phase transitions in random stabilizer tensor networks,” Physical Review B 105 (2022), 10.1103/physrevb.105.104306.
  * Lancien and Pérez-García (2021) Cécilia Lancien and David Pérez-García, “Correlation length in random MPS and PEPS,” Annales Henri Poincaré 23, 141–222 (2021).
  * Baxter (2007) R J Baxter, “Corner transfer matrices in statistical mechanics,” Journal of Physics A: Mathematical and Theoretical 40, 12577–12588 (2007).
  * Lubasch  _et al._ (2014) Michael Lubasch, J Ignacio Cirac, and Mari-Carmen Bañuls, “Unifying projected entangled pair state contractions,” New Journal of Physics 16, 033014 (2014).
  * Noh  _et al._ (2020) Kyungjoo Noh, Liang Jiang, and Bill Fefferman, “Efficient classical simulation of noisy random quantum circuits in one dimension,” Quantum 4, 318 (2020).
  * Li  _et al._ (2023) Zhi Li, Shengqi Sang, and Timothy H. Hsieh, “Entanglement dynamics of noisy random circuits,” Phys. Rev. B 107, 014307 (2023).
  * Hayden  _et al._ (2016) Patrick Hayden, Sepehr Nezami, Xiao-Liang Qi, Nathaniel Thomas, Michael Walter, and Zhao Yang, “Holographic duality from random tensor networks,” Journal of High Energy Physics 2016, 9 (2016).
  * Lopez-Piqueres  _et al._ (2020) Javier Lopez-Piqueres, Brayden Ware, and Romain Vasseur, “Mean-field entanglement transitions in random tree tensor networks,” Phys. Rev. B 102, 064202 (2020).
  * Nahum  _et al._ (2021) Adam Nahum, Sthitadhi Roy, Brian Skinner, and Jonathan Ruhman, “Measurement and entanglement phase transitions in all-to-all quantum circuits, on quantum trees, and in landau-ginsburg theory,” PRX Quantum 2, 010352 (2021).
  * Nahum  _et al._ (2017) Adam Nahum, Jonathan Ruhman, Sagar Vijay, and Jeongwan Haah, “Quantum entanglement growth under random unitary dynamics,” Phys. Rev. X 7, 031016 (2017).
  * Nahum  _et al._ (2018) Adam Nahum, Sagar Vijay, and Jeongwan Haah, “Operator spreading in random unitary circuits,” Phys. Rev. X 8, 021014 (2018).
  * Zhou and Nahum (2019) Tianci Zhou and Adam Nahum, “Emergent statistical mechanics of entanglement in random unitary circuits,” Phys. Rev. B 99, 174205 (2019).
  * Skinner  _et al._ (2019) Brian Skinner, Jonathan Ruhman, and Adam Nahum, “Measurement-induced phase transitions in the dynamics of entanglement,” Phys. Rev. X 9, 031009 (2019).
  * Li  _et al._ (2018) Yaodong Li, Xiao Chen, and Matthew P. A. Fisher, “Quantum zeno effect and the many-body entanglement transition,” Phys. Rev. B 98, 205136 (2018).
  * Jian  _et al._ (2020) Chao-Ming Jian, Yi-Zhuang You, Romain Vasseur, and Andreas W. W. Ludwig, “Measurement-induced criticality in random quantum circuits,” Phys. Rev. B 101, 104302 (2020).
  * Bao  _et al._ (2020) Yimu Bao, Soonwon Choi, and Ehud Altman, “Theory of the phase transition in random unitary circuits with measurements,” Phys. Rev. B 101, 104301 (2020).
  * Potter and Vasseur (2022) Andrew C. Potter and Romain Vasseur, “Entanglement dynamics in hybrid quantum circuits,” in  _Entanglement in Spin Chains: From Theory to Quantum Technology Applications_, edited by Abolfazl Bayat, Sougato Bose, and Henrik Johannesson (Springer International Publishing, Cham, 2022) pp. 211–249.
  * Fisher  _et al._ (2023) Matthew P.A. Fisher, Vedika Khemani, Adam Nahum, and Sagar Vijay, “Random quantum circuits,” Annual Review of Condensed Matter Physics 14, 335–379 (2023), https://doi.org/10.1146/annurev-conmatphys-031720-030658 .
  * Napp  _et al._ (2022) John C. Napp, Rolando L. La Placa, Alexander M. Dalzell, Fernando G. S. L. Brandão, and Aram W. Harrow, “Efficient classical simulation of random shallow 2d quantum circuits,” Phys. Rev. X 12, 021021 (2022).
  * Martinis et al’ (2019) John M. Martinis et al’, “Quantum supremacy using a programmable superconducting processor,” Nature 574, 505–510 (2019).
  * Nishino and Okunishi (1997) Tomotoshi Nishino and Kouichi Okunishi, “Corner transfer matrix algorithm for classical renormalization group,” Journal of the Physical Society of Japan 66, 3040–3047 (1997).
  * Levin and Nave (2007) Michael Levin and Cody P. Nave, “Tensor renormalization group approach to two-dimensional classical lattice models,” Physical Review Letters 99 (2007), 10.1103/physrevlett.99.120601.
  * Orús and Vidal (2009) Román Orús and Guifré Vidal, “Simulation of two-dimensional quantum systems on an infinite lattice revisited: Corner transfer matrix for tensor contraction,” Physical Review B 80 (2009), 10.1103/physrevb.80.094403.
  * Vidal (2003) Guifré Vidal, “Efficient classical simulation of slightly entangled quantum computations,” Physical Review Letters 91 (2003), 10.1103/physrevlett.91.147902.
  * Orús and Vidal (2008) R. Orús and G. Vidal, “Infinite time-evolving block decimation algorithm beyond unitary evolution,” Physical Review B 78 (2008), 10.1103/physrevb.78.155117.
  * Li  _et al._ (2021) Yaodong Li, Romain Vasseur, Matthew Fisher, and Andreas WW Ludwig, “Statistical mechanics model for clifford random tensor networks and monitored quantum circuits,” arXiv preprint arXiv:2110.02988 (2021).
  * Aaronson and Gottesman (2004) Scott Aaronson and Daniel Gottesman, “Improved simulation of stabilizer circuits,” Physical Review A 70, 052328 (2004).
  * Fattal  _et al._ (2004) David Fattal, Toby S Cubitt, Yoshihisa Yamamoto, Sergey Bravyi, and Isaac L Chuang, “Entanglement in the stabilizer formalism,” arXiv preprint quant-ph/0406168 (2004).
  * Anand  _et al._ (2022) Sajant Anand, Johannes Hauschild, Yuxuan Zhang, Andrew C Potter, and Michael P Zaletel, “Holographic quantum simulation of entanglement renormalization circuits,” arXiv preprint arXiv:2203.00886 (2022).

##  Appendix A Numerical Methods: iPEPS

Consider a 2d square lattice made of $L_{x}\times L_{y}\equiv N$ sites,
$L_{x}\leq L_{y}$, where each site is represented by a $d$-dimensional vector
space $\mathbb{C}^{d}$, and a state $\ket{\Psi}$ on the lattice that can be
expressed as a PEPS as introduced in the main text. The calculation of the
norm $\mathcal{N}=\braket{\Psi}{\Psi}$ is a common computation in PEPS
algorithms, closely related to the computation of more complicated quantities
of interest, e.g. expectation values $\braket{\Psi}{O}{\Psi}$ or correlation
functions $\braket{\Psi}{O_{1}O_{2}\cdots O_{k}}{\Psi}$ of local operators.
The exact contraction of the norm of $\mathcal{N}$ has a $O(L_{y}\exp(L_{x}))$
computational cost. The boundary MPS method is one of many related approaches
Nishino and Okunishi (1997); Baxter (2007); Levin and Nave (2007); Jordan  _et
al._ (2008); Orús and Vidal (2009); Lubasch  _et al._ (2014) for approximately
computing $\mathcal{N}$, while keeping the computational cost to scale
linearly in $L_{x}$ and $L_{y}$. In this Appendix we consider an infinite PEPS
(iPEPS), with $L_{x},L_{y}\to\infty$, by leveraging the translational
invariance of the state, and where the computational cost is independent of
$L_{x}$ and $L_{y}$.

The tensor network for $\mathcal{N}$ consists of an infinite square lattice
with a pair of tensors $TT^{*}$ on each site, where $T$ and $T^{*}$ are
connected over their physical index (see Figure 5a). One can show that the
norm $\mathcal{N}$ diverges/vanishes with $L_{x},L_{y}$ as
$\mathcal{N}=\nu^{L_{x}L_{y}}$ where $\nu$ is a non-negative real number known
as the norm-per-site. Although it is possible to compute $\nu$, in practice
this is not needed if we are interested in normalized quantities such as
-$\braket{\Psi}{O}{\Psi}/\mathcal{N}$ or the reduced density matrix $\varrho$,
given by the ratio of two 2d tensor networks that only differ by the insertion
of (one or more) local operators (see Figure 4). In this case, we need to
compute the fixed point boundary states $\left|\psi_{\text{\tiny\rm
top}}\right\rrangle,\left\llangle\psi_{\text{\tiny\rm bottom}}\right|$
(assumed to be unique) of the 1d transfer matrix given by an infinite row of
$TT^{*}$ (see Figure 5b). From now on we drop the double braket notation
$\left|\right\rrangle,\left\llangle\right|$ for boundary states, in favor of
single braket notation $\ket{},\bra{}$. The key idea behind the boundary MPS
(bMPS) algorithm is to approximately represent the infinite top and bottom
boundary states $\ket{\psi_{\text{\tiny\rm top}}}$ and
$\bra{\psi_{\text{\tiny\rm bottom}}}$ with two infinite 1d MPS of bond
dimension $\chi$. The algorithm detailed below describes how to obtain the top
fixed point boundary MPS (the generalisation to the bottom boundary is
straightforward).

(a) (b)

Figure 4:  Computation of observable and reduced density matrix: (a) Tensor
network for the computation of observables
$\braket{\Psi}{O}{\Psi}/\braket{\Psi}{\Psi}$. (b) Tensor netwok for the
computation of normalized $\varrho(\chi)$. More detailed calculation of these
objects is described below.

(a) (b) (c) (d)

Figure 5:  iPEPS norm computation: (a) Infinite square lattice tensor network
of the norm $\mathcal{N}$ of an iPEPS. (b) Norm boundary states
$\ket{\psi_{\text{\tiny\rm top}}}$ and $\bra{\psi_{\text{\tiny\rm bottom}}}$
with a 1d transfer matrix in between, which is proportional to the iPEPS norm.
(c) The bMPS algorithm obtains an approximation to $\ket{\psi_{\text{\tiny\rm
top}}}$ and $\bra{\psi_{\text{\tiny\rm bottom}}}$ with an infinite bMPS. (d)
Contracting from the left and from the right of the infinite chain in (c) we
obtain the left and right boundary states
$\bra{l_{1}[\mathcal{T}]},\ket{r_{1}[\mathcal{T}]}$ (see later for details on
their calculation).

(a) (b) (c) (d)

Figure 6:  iPEPS reduced density matrix computation: Computation of the iPEPS
2d reduced density matrix for one site $\mathcal{\varrho}$. It is closely
related to $N$ in Figure 5 but a physical index is left open. Note that to
obtain the correct $\varrho(\chi)$ we need to normalize the object in (d) as
in Figure 4a.

###  A.1 Boundary MPS: Method

We build an initial boundary MPS by multiplying pairs of tensors $TT^{*}$
after removing the vertical upper index (by fixing their value to 1)222Because
we are interested in the fixed point (assumed to be unique) of the 1d transfer
matrix, the choice of initial $\left|\psi(y=0)\right\rrangle$ is not
important.. This translation-invariant bMPS is characterized by a tensor $B$
with 4 indices represented graphically as “legs”: 2 bond legs connecting the
tensors with each other, and other 2 legs, one joining $B$ with $T$ and the
other with $T^{*}$. This initial bMPS has bond dimension $D^{2}$. If
$D^{2}\leq\chi$, we act with the transfer matrix $n$ times until the first
truncation is required (namely, when $D^{2n}>\chi$). After this point, we
iteratively perform two operations: i) truncation of the boundary MPS bond
dimension down to $\chi$ and ii) application of the 1d $TT^{*}$ transfer
matrix which increases the bMPS bond dimension to $\chi D^{2}$. This algorithm
ends once the bMPS has converged (see A.1.2). The converged bMPS
$\ket{\tilde{\psi}}$, characterized by tensor $\tilde{B}$, is our
approximation to the true fixed point boundary state $\ket{\psi_{\text{\tiny
top}}}$.

####  A.1.1 Infinite MPS Bond Dimension Truncation

A key component of the algorithm described above is to perform a bond
dimension truncation of the infinite bMPS from $\chi^{\prime}\to\chi$, where
$\chi^{\prime}>\chi$, while trying to keep a high fidelity between the states
represented by $B(\chi)$ and the original $B(\chi^{\prime})$. In order to do
this, we consider a bond bipartition into left and right halves of the
infinite spin chain represented by the bMPS. We bring this bond into the
canonical form Vidal (2003), such that the bond index, $i$, corresponds to the
labeling of Schmidt vectors in the Schmidt decomposition of the bMPS
representation of the state $|\psi\rangle$ across that index. More explicitly:

$|\psi\rangle=\sum_{i=1}^{\chi^{\prime}}\sigma_{i}\ket{\Phi^{L}_{i}}\otimes\ket{\Phi^{R}_{i}},$
(6)

where $\sigma_{i}$ label the Schmidt coefficients (normalized such that
$\sum_{i}\sigma_{i}^{2}=1$ and arranged in decreasing value
$\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{\chi^{\prime}}\geq 0$) and
$\ket{\Phi^{L}_{i}},\ket{\Phi^{R}_{i}}$ are the Schmidt vectors, which form an
orthonormal set
$\braket{\Phi^{L}_{i}}{\Phi^{L}_{j}}=\braket{\Phi^{R}_{i}}{\Phi^{R}_{j}}=\delta_{ij}$.
The basis $\ket{\Phi^{L}_{i}},\ket{\Phi^{R}_{i}}$ are also the eigenvectors of
the reduced density matrices $\rho^{A}$ and $\rho^{B}$, with eigenvalues
$\sigma_{i}^{2}$. Once we perform a gauge transformation to bring the bond
into this basis, we can truncate the bond dimension of the bMPS by keeping
only the largest $\chi$ Schmidt coefficients. We apply this truncation for
every bond in the infinite chain. This truncation procedure will keep
truncation errors low as long as the entanglement spectrum eigenvalues
$\\{\sigma_{i}^{2}\\}$ decay sufficiently fast with $i$ Orús and Vidal (2008);
Vidal (2003).

We summarize the truncation algorithm in steps i) to vi) as depicted in Figure
7a. We need to obtain the invertible matrices $M$ and $P$ that bring the
translationally invariant bMPS tensor $B$ into the left and right canonical
form $B_{l}$ and $B_{r}$ respectively. This gauge requirements are such that
the left and right dominant eigenvectors of $B_{l}B_{l}^{*}$ and
$B_{r}B_{r}^{*}$ respectively (contracted over the physical bMPS index as
shown in Figure 7b) should be the identity and its corresponding dominant
eigenvalue, $\lambda_{1}$, should be 1 for normalisation. These conditions are
depicted in Figure 7b. Starting from the large bond dimension infinite bMPS
that we wish to truncate (Fig 7a.i)), we insert two resolutions of the
identity $M^{-1}M$ and $PP^{-1}$ at every bond (Fig 7a.ii)). This yields an
infinite chain of $P^{-1}BM^{-1}$ tensors denoted $B_{c}$ (Fig 7a.iii)) with a
matrix $MP$ on the bonds (Fig 7a.iv)). We consider the singular value
decomposition of the matrix $MP=USV^{\dagger}$. The diagonal entries of the
matrix $S$ and the isometries $U,V^{\dagger}$ correspond to the Schmidt
decomposition presented in Equation 6 if we had considered only one bond
partitioning the infinite chain in half. We truncate these objects and keep
only the largest $\chi$ singular values and corresponding vectors therefore
approximating the bond with smaller sized matrices $MP\approx
U^{\prime}S^{\prime}V^{\prime\dagger}$ (Fig 7a.v)). Putting everything
together, as in Fig 7a.vi), we obtain a final expression for
$B(\chi)=S^{\prime
1/2}V^{\prime\dagger}P^{-1}B(\chi^{\prime})M^{-1}U^{\prime}S^{\prime 1/2}$.

(a) (b)

Figure 7:  Truncation of bMPS bond dimension: (a) Steps, i)-iv), in the
algorithm for truncating the bond dimension of the bMPS. (b) Finding left and
right canonical form gauge transformation matrices $M$ (top) and $P$ (bottom).

####  A.1.2 Convergence of the Boundary MPS to a Fixed Point

After a sufficient number of iterations, $y_{*}$, of the 1d $TT^{*}$ transfer
matrix followed by the truncation as described in A.1.1 we observe that our
bounday MPS has converged, meaning that $B(y_{*})\approx B(y_{*}+1)$ within a
set tolerance, in which case we define $\tilde{B}=B(y_{*}+1)$. This
equivalence can be assessed by looking at the fidelity-per-site (see A.2.2 for
definition) between $B(y)$ and $B(y+1)$ (which should be approximately 1 if
convergence has been reached) or by comparing the Schmidt spectrum of the
boundaries in successive iterations (which should be equal). In our work we
use both metrics to ensure convergence.

###  A.2 Infinite MPS and Infinite PEPS: Computation of Metrics

(a) (b) (c)

Figure 8:  Transfer matrices: (a) Norm $\braket{\psi}{\psi}$ of an infinite
MPS (left) and the transfer matrix $\mathcal{T}_{\braket{\psi}{\psi}}$
involved in its calculation. (b) Fidelity $F=\braket{\psi_{1}}{\psi_{2}}$
between two infinite MPS and the transfer matrix
$\mathcal{T}_{\braket{\psi_{1}}{\psi_{2}}}$. (c) Transfer matrix $\mathcal{T}$
used for computing $\mathcal{N},\varrho$ and $\braket{\psi}{O}{\psi}$.

####  A.2.1 Infinite MPS: Norm and Correlation Length

Consider an infinite MPS $\ket{\psi}$ given by tensor $B$ of bond dimension
$\chi$. Its transfer matrix $\mathcal{T}_{\braket{\psi}{\psi}}$ (see Figure
8a) has eigenvalue decomposition:

$\mathcal{T}_{\braket{\psi}{\psi}}=\sum_{i=1}^{\chi^{2}}\lambda_{i}\ket{r_{i}}\bra{l_{i}},$
(7)

where $\\{\lambda_{i}\\}$ is the set of eigenvalues, organized such that
$|\lambda_{1}|>|\lambda_{2}|\geq\cdots\geq|\lambda_{\chi^{2}}|$, and
$\ket{r_{i}},\bra{l_{i}}$ are the right and left eigenvectors respectively.
One can show that: i) each eigenvalue $\lambda$ is either real or part of a
complex conjugate pair $\lambda,\lambda^{*}$ and ii) the dominant eigenvalue
$\lambda_{1}$ is real and non-negative. We have also assumed that the dominant
eigenvalue $\lambda_{1}$ is not degenerate. The norm $\braket{\psi}{\psi}$
(8a) is given by $\lambda_{1}^{L_{x}}$, where $L_{x}$ is the number of tensors
in the MPS, and thus $L_{x}\to\infty$ for an infinite MPS. If $\lambda_{1}\neq
1$, implying that $\braket{\psi}{\psi}\neq 1$, we may obtain a normalized MPS
$\ket{\psi^{\prime}}$ with tensor $B^{\prime}=B/\sqrt{\lambda_{1}}$ so that
$\lambda_{1}^{\prime}=1$ and $\braket{\psi^{\prime}}{\psi^{\prime}}=1$.

The correlation length $\xi$ is calculated by considering the dominant and
second largest eigenvalue $\lambda_{1},\lambda_{2}$ and it is given by:

$\xi=\frac{-1}{\log(|\lambda_{2}[\mathcal{T}_{\braket{\psi}{\psi}}]/\lambda_{1}[\mathcal{T}_{\braket{\psi}{\psi}}]|)},$
(8)

where we introduce the notation
$\lambda_{i}[\mathcal{T}_{\braket{\psi}{\psi}}]$ to refer to the
$i^{\text{th}}$ eigenvalue of the transfer matrix
$\mathcal{T}_{\braket{\psi}{\psi}}$.

####  A.2.2 Infinite MPS: Fidelity

Consider two iMPS $\ket{\psi_{1}},\ket{\psi_{2}}$ given by tensors
$B_{1},B_{2}$ and bond dimensions $\chi_{1},\chi_{2}$ respectively, normalized
such that $\braket{\psi_{1}}{\psi_{1}}=\braket{\psi_{2}}{\psi_{2}}=1$, that is
$\lambda_{1}[\mathcal{T}_{\braket{\psi_{1}}{\psi_{1}}}]=\lambda_{1}[\mathcal{T}_{\braket{\psi_{2}}{\psi_{2}}}]=1$
(see previous section). Their fidelity,
$F[\psi_{1},\psi_{2}]=\braket{\psi_{1}}{\psi_{2}}$ can be expressed as
$F=F_{s}^{L_{x}}$, where the fidelity-per-site $F_{s}$ is such that
$0\leq|F_{s}|\leq 1$. One can show that $|F_{s}|=1$ if and only if $B_{2}$
equals $B_{1}$ up to a so-called gauge transformation, in which case both
tensors give rise to the same state. To calculate $F_{s}$ consider the
transfer matrix given by $B_{1}B_{2}^{*}$, as depicted in Figure 8b. This
transfer matrix can be decomposed in the standard way:

$\mathcal{T}_{\braket{\psi_{1}}{\psi_{2}}}=\sum_{i=1}^{\chi_{1}\chi_{2}}\lambda_{i}\ket{r_{i}}\bra{l_{i}},$
(9)

where $\\{\lambda_{i}\\}$ is the set of eigenvalues, organized such that
$|\lambda_{1}|>|\lambda_{2}|\geq\cdots\geq|\lambda_{\chi_{1}\chi_{2}}|$, and
$\ket{r_{i}},\bra{l_{i}}$ the right and left eigenvectors respectively. Here,
we have assumed that the dominant eigenvalue $\lambda_{1}$ is not degenerate.
Then, the fidelity-per-site $F_{s}$ is simply given by this dominant
eigenvalue, $F_{s}=\lambda_{1}[\mathcal{T}_{\braket{\psi_{1}}{\psi_{2}}}]$.

####  A.2.3 Infinite PEPS: Reduced Density Matrix

The reduced density matrix for a local region $R$ of the 2d lattice is denoted
$\varrho_{R}=\text{tr}_{\bar{R}}\ket{\Psi}\bra{\Psi}/\mathcal{N}$, where the
trace is over all sites except those in $R$. For simplicity, we consider the
case where $R$ consists of a single site and refer to this single-site reduced
density matrix simply as $\varrho$. In a PEPS calculation with boundary MPS
with bond dimension $\chi$, we obtain an approximation $\varrho(\chi)$ to the
exact $\varrho$. The calculation of $\varrho(\chi)$ is shown in Figure 4b,
where the dominant left and right eigenvectors $\ket{r_{1}},\bra{l_{1}}$ are
those of the transfer matrix $\mathcal{T}$, depicted in Figure 8.

The error in $\varrho(\chi)$ due to finite $\chi$ in the bMPS depends directly
on the square of the Schmidt coefficients of the bMPS, as seen in Figure 12a,
and therefore on the $n=1$ (von Neumann) Renyi entropy.

###  A.3 iPEPS and Boundary MPS: Additional Numerical Results

In this appendix we provide additional numerical results for the contraction
of 2d infinite PEPS with the boundary MPS formalism.

####  A.3.1 Averaging Over Several Realisations

In the main text we have mostly displayed results for a single instance of a
random tensor $T$ for each PEPS bond dimension $D$. Here, we examine several
random instances to build confidence that the reported behavior is typical.
Specifically, we explore average properties of $N_{s}=15,15,10,5,1$ instances
of clean, random iPEPS with bond dimension $D=2,3,4,5,6$ respectively. In all
cases, we again observe rapid convergence of the entanglement spectrum (Figure
9) and Renyi entropies (Figure 10) with the truncated bMPS bond dimension,
$\chi$. Moreover, the variance in these quantities decreases rapidly with
increasing PEPS bond dimension $D$. This suppression qualitatively agrees with
the statistical mechanical model, for which fluctuations around the ground
state are suppressed at large $D$.

(a) (b)

(c) (d)

(e) (f)

Figure 9:  bMPS average Schmidt entanglement spectrum, several 2d PEPS
instances $T$: (a) and (b) $D=2$ and $D=3$ (respectively) of the average and
standard deviation (over 15 instances) of the Schmidt spectrum of
$\tilde{B}(\chi)$ for a range of $\chi$ between $\chi_{\rm min}=2,\chi_{\rm
max}=72$. The average and standard deviation is calculated for each eigenvalue
index. (c)-(f) $D=2,3,4,5$ (respectively) Schmidt spectrum for
$\tilde{B}(\chi_{\rm max}=72)$ for varying number of instances.

(a) (b) (c)

Figure 10:  Renyi entropy order dependence: (a) and (b) $D=2$ and $D=3$
(respectively) system Renyi entropy as a function of order $n$, for all
simulated bMPS bond dimension $2\leq\chi\leq 72$. Average over 15 random
initializations of $T$. (b) Renyi entropy for $\tilde{B}(\chi_{\rm max}=72)$
as a function of order $n$ for $D=2,3,4,5,6$ averaged over 15,15,10,5,1
iterations respectively. We note that the $n=1$ (von Neumann) entropy is
expected to control the error of the truncated bMPS approximation, and that
the $n=0$ entropy simply counts the number of Schmidt weights and trivially
saturates to $\chi_{\rm max}$.

####  A.3.2 Fidelities of Converged bMPS

In order to quantify how close $\ket{\tilde{\psi}(\chi)}$ is to the exact
fixed-point boundary state we calculate the fidelity-per-site $F_{s}$ between
$\ket{\tilde{\psi}(\chi)}$ and $\ket{\tilde{\psi}(\chi_{max})}$, where we use
$\ket{\tilde{\psi}(\chi_{max})}$ as a proxy for $\ket{\tilde{\psi}}$. We see
in Figure 11a that this fidelity is very close to 1 already for a very small
value of $\chi$, characteristic of weakly entangled states. This small value
of $\chi$ increases slightly as a function of the 2d PEPS bond dimension $D$.
The sharp transitions in fidelity match the Schmidt spectrum decay in Figure
1b. We also calculate the required $\chi$ to obtain a certain fidelity-per-
site as a function of $D$ (Figure 11b) and we notice a $\chi\sim D$ relation.
This is in contrast with the $\chi\sim D^{2}$ relation often quoted in the
literature for ground states of local Hamiltonians. This discrepancy is likely
due to the weakly entangled character of the random PEPS.

(a) (b)

Figure 11:  bMPS fidelity numerical results: (a) Convergence of the fidelity-
per-site of the largest simulated bMPS network ($\chi_{\rm max}=72$) with all
other smaller systems. (b) Bond dimension required for the bMPS to obtain a
certain fidelity-per-site with the largest simulated boundary.

####  A.3.3 iPEPS Reduced Density Matrix

In Figure 1 of the main text we presented the change $\Delta\varrho(\chi)$ in
the physical single-site density matrix $\varrho(\chi)$ as a function of
$\chi$, where we took $\varrho(\chi_{\rm max}=72)$ as a reference. We observed
that $\Delta\varrho(\chi)$ is exponentially suppressed with $\chi$. We expect
this result to apply also to the density matrix $\varrho(\chi)$ on a larger
local region. As an example, Figure 12b shows $\Delta\varrho(\chi)$ for
$\varrho(\chi)$ for a region of two contiguous sites, where indeed we again
observe exponential suppression with $\chi$.

(a) (b)

Figure 12:  2d PEPS reduced density matrix: (a) Largest eigenvalue of the
difference between the largest and smallest bMPS $\rho(L)$ divided by the
squared Schmidt spectrum of the boundary for $\chi_{\rm max}$. This plot shows
a constant dependence, indicating that the order-1 Renyi entropy is to be
considered when analysing observable accuracy of the 2d PEPS contraction. (b)
Change $\Delta\varrho(\chi)$ in the physical $L=2$ density matrix
$\varrho(\chi)$ as a function of $\chi$

####  A.3.4 Correlation Length

The correlation length that we obtain from the converged bMPS at maximum bond
dimension are low: averaging at order $\sim 1$ for all bond dimensions $D$
(Figure 13b). The convergence of the correlation length as a function of bond
dimension $\chi$ (Figure 13a) also matches the gaps and exponential decay in
the spectrum as per Figure 1. This also implies that for relatively low bond
dimensions we obtain correlation lengths that are comparable to our best
estimate given by $\chi_{\rm max}$.

(a) (b)

Figure 13:  bMPS correlation length numerical results: (a) Convergence of the
correlation length as a function of bMPS bond dimension $\chi$ and for varying
2d PEPS sizes $D$. (b) Average correlation length for $\tilde{B}(\chi_{\rm
max}=72)$ as a function of $D$ for $D=2,3,4,5,6$, averaged over 15,15,10,5,1
initialisation respectively.

####  A.3.5 bMPS Entanglement Entropy for the 2d Tensor Network Corresponding
to the Overlap of two iPEPS

We have numerically and analytically confirmed that for 2d random PEPS the
approximate computation of the norm $\mathcal{N}$ and expectation values of
local observables can be performed efficiently. However, other types of
observables, such as individual wave-function components, $\langle s_{1}\dots
s_{N}|\Psi\rangle$, or overlaps between distinct random PEPS,
$\langle\Psi^{\prime}|\Psi\rangle$, are expected to be hard since they involve
a 1d transfer matrix map $T{T^{\prime}}^{*}$, where $T\neq T^{\prime}$, which
is not positive.
We consider the computation of the overlap of two distinct 2d random PEPS
$\braket{\Psi}{\Psi^{\prime}}$ given by tensors $T$ and $T^{\prime}$, where we
build $T^{\prime}$ from $T$ with noise tuned by a parameter $\eta\in[0,1]$,
such that $T^{\prime}=(1-\eta)T+\eta\mathcal{C_{N}}$, where $\mathcal{C_{N}}$
is the normal complex ditribution with zero mean and standard deviation 1. In
the results in Figure 14 we observe an increase in the Schmidt entanglement
entropy at the boundary for the contraction evolution of the overlap
$\braket{\Psi}{\Psi^{\prime}}$, where, for large $\eta$ the entanglement
saturates to a constant value given by the bond dimension $\chi$ of the bMPS.

(a) (b)

Figure 14: Global vs local observables: (a) and (b) show the evolution of the
entanglement entropy at the boundary MPS for the overlap
$\braket{\Psi}{\Psi^{\prime}}$, where both 2d PEPS have bond dimension $D=2$
and $D=3$ (respectively) and the bMPS has bond dimension $\chi=35$.

##  Appendix B Stabilizer PEPS formalism and additional numerical results

In this appendix, we provide details about the stabilizer PEPS formalism and
its simulation.

A stabilizer PEPS is a PEPS whose composing tensors $\\{T_{\textbf{[r]}}\\}$
are stabilizer tensors, i.e. each of the states defined through

$\ket{T_{\textbf{[r]}}}=\sum_{i,j,k,l,s}\left(T_{\textbf{[r]}}\right)_{s,ijkl}|s,ijkl\rangle,$
(10)

is a quantum stabilizer state over $p-$qudits, for some prime number $p$. In
the expression above, $i,j,k,l$ are virtual bond indices, and $s$ is the
physical bond index. Each bond represents a Hilbert space containing an
integer number of $p$-qudits, thus both physical and virtual bond dimensions
need to be some integer power of $p$: $D=p^{k_{D}},d=p^{k_{d}}$.

The contraction of two stabilizer tensors can be simulated efficiently by
performing several Bell measurements on the contracted bonds, as detailed
below. Let us assume $T^{1}_{ab_{1}}$ and $T^{2}_{b_{2}c}$ are two stabilizer
tensors where the two $b$ indices are of the same dimension $p^{k}$. Let
$T^{12}_{ac}$ be the tensor obtained by contacting the two tensors over the
$b$ indices:

$T^{12}_{ac}:=\sum_{b}T^{1}_{ab}T^{2}_{bc}.$ (11)

The contraction can be realized by performing several (forced) Bell
measurements on qudits within $b_{1}$ and $b_{2}$ Li  _et al._ (2021). If we
label the qudits in the $b_{1}$ with $\\{1_{1},...k_{1}\\}$, and those in the
$b_{2}$ with $\\{1_{2},...k_{2}\\}$, then:

$|T^{12}\rangle\propto\prod_{i=1}^{k}\left(\mathcal{M}[X_{i_{1}}X_{i_{2}}]\mathcal{M}[Z_{i_{1}}Z^{-1}_{i_{2}}]\right)~{}|T^{1}\otimes
T^{2}\rangle,$ (12)

where $\mathcal{M}[P]$ is the projector to the $+1$ subspace of the Pauli
operator $P$. Using the Gottesman-Knill theorem Aaronson and Gottesman (2004),
one can show that the state $|T^{12}\rangle$ is still a stabilizer state.
Further, its stabilizers can be obtained from those of $T^{1}$ and $T^{2}$
with $O(\log^{3}D_{\text{total}})$ time complexity, where
$D_{\text{total}}=\dim a\cdot\dim c\cdot(\dim b)^{2}$ is the product of all
the open bonds’ dimensions.

It is worth noting that the complexity of contracting stabilizer tensors is
independent of the entanglement property of the underlining states. Thus for a
given stabilizer PEPS, we are able to compute the boundary state
$|\psi(t)\rrangle$’s evolution exactly without any truncation and study its
entanglement properties. The latter can be computed from the state’s
stabilizers using the algorithms introduced in Fattal  _et al._ (2004).

We consider both disordered and clean (i.e. translational invariant)
stabilizer PEPS. In both cases the PEPS is finite and takes periodic boundary
conditions along the $x$ direction, thus the boundary state $|\psi(y)\rrangle$
is also finite and has the periodic boundary condition. To reach the large $D$
regime, we take $p=173$. The simulated von-Neumann entropy $S_{A}(y)$ 333For a
stabilizer state (tensor) and any bi-partition of bonds, the non-zero singular
values associated with the bi-partition are all of the same value. Thus the
Rényi entanglement entropy across the bi-partition is independent of the Rényi
index $\alpha$, in particular it equals the von Neumann entropy ($\alpha=1$).
of $|\psi(y)\rrangle$ is presented in the Figure 3 in the main text as well as
the Figure 15.

Figure 15:  Further numerical results for stabilizer PEPS – (left) Profile of
$|\psi(y)\rrangle$’s entanglement entropy in a disordered PEPS. The $L$\- and
$|A|$\- independent plateau suggests area-law entanglement of
$|\psi(y)\rrangle$. (mid) The simulated ‘entanglement barrier’ for various
choices of $D$ and $d$ in the disordered PEPS. (right) The simulated
‘entanglement barrier’ in the clean PEPS.

We start by focusing on the disordered case, where each unit stabilizer tensor
$T_{\textbf{r}}$ is sampled independently. The simulation shows that at any
given layer number $y$ and when $|A|$ is far from $1$ or $L$, the von Neumann
entropy $S_{A}(y)$ takes a constant value that is independent of $|A|$ or the
system size $L$ (Figure 15, left). Further, the constant value first increases
linearly and then drops with the increase of the number of contracted layers
$y$. Both the peak value and the turning point are dependent on the bond
dimensions $D$ and $d$ (Figure 15, mid). The two plots together suggest that
the boundary state $|\psi(t)\rrangle$ is at most area-law entangled at any
time (layer number) $y$. Further, the simulated entanglement barrier’s
dependence on $(d,D)$ matches with the prediction of the stat-mech mapping in
Equation 5, as is shown in Figure 3 in the main text.

Next we come to the clean PEPS case, where all unit tensors $T_{[{\bm{r}}]}$
are identical and taken to be a randomly sampled stabilizer tensor. The
simulation results suggest that the behavior of $S_{A}(y)$ is almost identical
to that in the disordered case. Namely, the $|\psi(y)\rrangle$ is also area-
law entangled in the clean case, with an area-law value following the
expression Equation (5) (Figure 15(right) and Figure 3(right)).

##  Appendix C Details of Statistical-Mechanics Mapping

In this appendix, we review the derivation of the stat-mech description of
entanglement features of random PEPS.

In the gaussian random PEPS ensemble, the tensor, $T_{r}$, for each site, $r$,
is chosen independently and identically distributed from a Gaussian
distribution:

$\displaystyle\mathbb{E}\left[\left(T^{*}_{[r^{\prime}]}\right)^{s^{\prime}}_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}\left(T_{[r]}\right)^{s}_{ijkl}\right]=\delta_{ss^{\prime}}\delta_{ii^{\prime}}\delta_{jj^{\prime}}\delta_{kk^{\prime}}\delta_{ll^{\prime}}\delta_{r,r^{\prime}},$
(13)

where $\mathbb{E}[\dots]$ denotes averaging over the ensemble, $s=1\dots d$ is
the physical index, and $i,j,k,l=1\dots D$ are bond indices, and $r$ label
sites of the lattice. In the following, we drop the indices $s,i,j,k,l$ on the
tensors.

###  C.1 Mapping RTNs to Replica-Magnets

Consider the tensor network contraction to compute the norm of the PEPS:
$\mathcal{N}=\langle\Psi|\Psi\rangle$ using the MPS method outlined in Section
A above. Denote the (unnormalized) density matrix of the evolved boundary
state as: $\rho(y)=|\psi(y)\rrangle\llangle\psi(y)|$. Our aim will be to
compute the evolution of the ensemble-averaged Renyi entanglement entropy of a
region $A$ of the evolved boundary-state, $|\psi(y)\rrangle$.

$S_{A}^{(n)}=\frac{1}{1-n}\log\frac{{\rm tr}\rho_{A}^{n}}{({\rm
tr}\rho)^{n}},$ (14)

where and $\rho_{A}$ is the reduced density matrix in some contiguous interval
$A$ of size $L_{A}$ obtained from tracing out the complement of $A$ in
$\rho=\left|\psi\right\rrangle\left\llangle\psi\right|$.

Since the wavefunction $\left|\psi\right\rrangle$ is not necessarily
normalized, so the denominator in Equation (14) is crucial to obtain a
meaningful entanglement entropy. Directly computing the disorder average of
this non-linear quantity is challenging. To avoid this difficulty, we employ a
standard replica trick based on the identity:

$\displaystyle\log{\rm tr}\rho_{A}^{n}=\lim_{m\to 0}(({\rm
tr}\rho_{A}^{n})^{m}-1)/m.$ (15)

This allows us to express the disorder average of eq. (14) as

$\mathbb{E}\left[S_{A}^{(n)}\right]=\frac{1}{n-1}\lim_{m\to
0}\frac{1}{m}\left({\cal F}_{A}-{\cal F}_{0}\right),$ (16)

with $\mathcal{F}_{A,0}=-\log\mathcal{Z}_{A,0}$ and
$\mathcal{Z}_{0}\equiv\mathbb{E}\left[(\text{tr}\rho^{n})^{m}\right]$,
$\mathcal{Z}_{A}\equiv\mathbb{E}\left[(\text{tr}\rho_{A}^{n})^{m}\right]$.
Using this identity, the calculation of the Renyi entropies reduces to
computing $\mathcal{Z}_{0}$ and $\mathcal{Z}_{A}$, and to evaluate the replica
limit (16).

When $m$ and $n$ are integers, the averages in $\mathcal{Z}_{0}$ and
$\mathcal{Z}_{A}$ can be evaluated analytically using Wick’s theorem. One can
then express the partition functions $\mathcal{Z}_{A}$ and $\mathcal{Z}_{0}$
in terms of a classical statistical mechanics model, whose degrees of freedom
are permutations labelling different Wick contractions at each vertex of the
tensor networks: at each site, each tensor $T_{r}$ must be paired with a
$T_{r}^{*}$ possibly belonging to a different replica. Let $Q=nm$ be the
number of copies of $\rho$. Then, the partition function $\mathcal{Z}$
involves computing quantities like $\mathbb{E}\left[\rho^{\otimes Q}\right]$.
Note that $|\psi\rrangle$ contains both $T_{r}$ with $T_{r}^{*}$ at each site
$R$. Then, in the replicated theory there are $2Q$ copies of $T$, and $2Q$
copies of $T^{*}$ for each site, which we label, $T^{\alpha}_{r}$ with a
replica index $\alpha=1\dots 2Q$. We adopt the following ordering for the $2Q$
copies:

$\\{1,\overline{1},2,\overline{2},\dots Q,\overline{Q}\\}$ (17)

where $k$ labels the state (“ket”) in replica $Q$, and $\overline{k}$ denotes
the dual state (“bra”) in replica $Q$. To label permutations we use cycle
notation, for example $(124)(35)\in S_{6}$ denotes the permutation
$123456\rightarrow 245136$, i.e. with separate cyclic permutations of elements
$(124)$ and of elements $(35)$ [for convenience, we only list the cycles with
more than one element]. We will also need to define the cycle counting
function

$\displaystyle C(g,h)\equiv C(g^{-1}h)=\text{\\# of cycles in $g^{-1}h$},$
(18)

where $C(g)$ also includes single-element cycles that are not listed
explicitly in our notation e.g. for the above example, $C[(124)(35)]=3$.

According to Wick’s theorem, upon averaging over the Gaussian random tensors,
a non-zero contribution is obtained only if each tensor $T^{\alpha}_{r}$ is
paired with a permuted copy $\left(T^{g_{[r]}(\alpha)}_{r}\right)^{*}$, where
$g_{[r]}\in S_{2Q}$ labels a permutation of the replicas, and $S_{2Q}$ is the
symmetric group on $2Q$ elements.

The partition function $\mathcal{Z}$ corresponding to the tensor network
contraction, can then be written as a sum over replica-permutation “spins” for
each site:

$\displaystyle\mathcal{Z}=\sum_{\\{g_{[r]}\\}}W[\\{g_{[r]}\\}],$ (19)

where $W$ is the weight of the Wick contraction for the corresponding spin
configuration. The weight can be computed analytically for each bond in the
tensor network. There are three distinct types of contractions to consider:

  1. 1.

Bulk bonds connecting different nearest-neighbor nodes $i$ and $j$ with
permutation “spins” $g_{[r]}$ and $g_{r^{\prime}}$, and bond dimension $D$.
The same contraction occurs in the $2Q$ layers of the replicated tensor
network. Since Wick contractions force indices to be the same, the resulting
weight is equal to $D$:

$D^{C(g_{[r]},g_{[r^{\prime}]})}={\rm e}^{(\log
D)C(g_{[r]},g_{[r^{\prime}]})},$ (20)

on each link of the square lattice, since the number of independent bond
indices is equal to $C(g)$, the number of cycles in the permutation $g$. This
is most easily seen by a graphical representation: each cycle in
$g_{[r]}^{-1}g_{[r^{\prime}]}$ leads to a “loop” where indices have to be the
same, with weight $\sum_{\alpha=1}^{D}\delta_{\alpha\alpha}=D$. Interpreting
this positive weight as a Boltzmann weight, this terms leads to a
ferromagnetic interaction (favoring $g_{[r]}=g_{r^{\prime}}$ for neighboring
$i,j$ to maximize the number of cycles to $2Q$) with interaction strength
$\log D$. This Boltzmann weight has a left/right symmetry $\left(S_{2Q}\times
S_{2Q}\right)\rtimes\mathbb{Z}_{2}$ (where the extra $\mathbb{Z}_{2}$ symmetry
corresponds to $g\to g^{-1}$).

  2. 2.

Bulk contraction of $T_{r}$ with $T_{r}^{*}$ along the physical leg with
dimension $d$. This contraction can be implemented by adding a site with fixed
permutation equal to identity $e=()$: in each replica $k$, we pair $k$ with
itself (corresponding to gluing $T$ with $T^{*}$ in the ket), and
$\overline{k}$ with itself (corresponding to gluing $T$ with $T^{*}$ in the
bra). This leads to a factor

$d^{C(eg_{[r]})}={\rm e}^{(\log d)C(g_{[r]})},$ (21)

on each site. This can be seen as a $S_{2Q}\times S_{2Q}\to S_{2Q}$ symmetry-
breaking field favoring the identity permutation. This bulk field on every
site prevents any phase transition by creating an energy costs for domains of
spins with $g_{[r]}\neq e$ that scales as the volume of the domain.

  3. 3.

Boundary contractions at the top layer to compute
$\mathcal{Z}_{0}\equiv\overline{(\text{tr}\rho^{n})^{m}}$,
$\mathcal{Z}_{A}\equiv\overline{(\text{tr}\rho_{A}^{n})^{m}}$. At the top
layers we have dangling legs, that should be contracted to implement the trace
and partial trace operations to compute entanglement. In $\mathcal{Z}_{0}$, we
want to compute $\text{tr}\rho$ in each replica. This means that in each
replica (and at each site at the boundary), we want to contract $T_{r}$ (resp.
$T_{r}^{*}$) in the ket with $T^{*}_{i}$ (resp. $T_{r}$) in the bra. In our
language this corresponds to the permutation

$g_{0}=(1\overline{1})(2\overline{2})\dots(Q\overline{Q}),$ (22)

Note that this permutation is not identity, it has $Q$ cycles while $e=()$ has
$2Q$ cycles. At the end of each leg, we fixed the permutation to $g_{0}$:

$D^{C(g_{0}^{-1}g_{[r]})}={\rm e}^{(\log D)C(g_{0}g_{[r]})},$ (23)

for $i=1,\dots,L$ at the top boundary. To implement the partial trace in
$\mathcal{Z}_{A}\equiv\overline{(\text{tr}\rho_{A}^{n})^{m}}$, we fixed the
permutation to $g_{0}$ if $i$ is in $\overline{A}$, and to $g_{A}$ is $r\in
A$, with

$g_{A}=\left(\left(1\overline{2}\right)\left(2\overline{3}\right)\dots\left(n\overline{1}\right)\right)^{\otimes
m}.$ (24)

Assembling these ingredients results in the effective Hamiltonian of the main
text.

###  C.2 Comparison to related stat-mech models

A nearly identical replica-spin model was derived in Vasseur  _et al._ (2019)
for holographic random tensor network states (rTNS), i.e. whose tensors had
physical legs only at the boundary, and only virtual bond legs in the bulk. In
that work, the key difference was that the holographic rTNS did not have
positive tensors. As a result the permutation spins were $S_{Q}$\- rather than
$S_{2Q}$\- valued, and the bulk retained the $S_{Q}$ symmetry since there was
no field along $e$. This crucial difference led to two possible phases of the
stat-mech model: a disordered (paramagnetic) phase at weak coupling (low-$D$)
in which the permutation spins are short-range correlated, and an ordered
(ferromagnetic) phase at strong coupling (large-$D$) in which the $S_{Q}$
symmetry was spontaneously broken and the permutation spins have long range
order. In the ordered phase, domain walls had a non-vanishing surface-tension,
resulting in an extensive free-energy for the boundary twist in the
entanglement region size, resulting in volume-law entanglement scaling. By
contrast, in the disordered phase, there is only a local free-energy cost at
the edge of the boundary-domains, resulting in area-law entanglement.

Coming back to the stat-mech model for the PEPS norm computation: the key
difference is that the tensors are completely positive, i.e. are composite
tensors made up of $T$ and $T^{*}$ with physical legs contracted. This results
in a bulk field along the identity ($e$) permutation that explicitly breaks
the $S_{2Q}$ symmetry. A similar statistical mechanics model emerges in the
context of random quantum channels Li  _et al._ (2023). In that language, bond
dimension corresponds to physical dimension in the channel, and our physical
dimension $d$ maps to the strength of channel (effectively the number of Kraus
operators).

Intuitively, this explicit symmetry breaking destroys the ordering transition,
such that the entire phase diagram is effectively ordered (in the sense that
domain walls have a non-zero surface tension). Naively, one might expect this
to result in a volume-law entanglement throughout for any $D$. However, as we
show next, there is an exact cancellation of the volume-law contribution to
the free-energy with twisted boundary conditions, generically resulting in
area-law entanglement for the evolved boundary state.

Figure 16:  Graphical notation for permutations and cycle counting –  (Left)
In each of the $Q=mn$ replicas of $\rho=\ket{\psi}\bra{\psi}$, there are $2Q$
$T$ tensors, and $2Q$ $T^{*}$ tensors. We label $T,T^{*}$ tensors from the bra
$\bra{\psi}$ with an over-bar. It is convenient to draw the $T$ and $T^{*}$’s
in $m$ groups of $n$ (shown here for $m=1$). Averaging over tensors forces a
wick contraction between $T_{i}$ and $T^{*}_{g(i)}$ where $g\in S_{2Q}$. (Top
right) The three permutations corresponding to the bulk ($e$ = identity) and
boundary ($g_{A},g_{0}$) fields. (Bottom right) Graphical calculation of cycle
counting for various fields, $C(\sigma)$ is given by counting the number of
independent loops in the picture (note, for $m>1$, not shown, there would
simply be $m$ independent copies of this picture).

###  C.3 Minimum cut picture of random PEPS contraction

At large-$D$, the permtuation spins are strongly locked to each other by their
ferromagnetic interactions, and pinned to the bulk $e$ fields. Here, domain
walls have a non-zero line tension, and we can approximately compute the
domain wall free-energy for the entanglement entropy, by minimizing this line-
tension.

Let us focus on the thermodynamic limit, $L,y\to\infty$ and on half-system
entanglement. The statistical mechanics model has a bulk symmetry-breaking
field that prevents an entanglement phase transition which would be associated
with a spontaneous breaking of $S_{2Q}$ symmetry. Specifically, the $e$ field
produces an energy cost for domains with $g_{[r]}=e$ that scales like the
volume of the domain. This field favors the identity permutation $e$ in the
bulk (with fluctuations suppressed if $D\gg 1$), while the boundary fields
favor $g_{0}$ or $g_{A}$. However, the energy cost of the domain walls between
those permutations and $e$ are the same, since $g_{A}$ and $g_{0}$ each have
$Q=nm$ cycles. Therefore, ${\cal F}_{A}$ and ${\cal F}_{0}$ each contain an
extensive term $L\log D$, but importantly this extensive term is the same for
both ${\cal F}_{A,0}$ and cancels in the difference. This cancellation can be
traced back to the $S_{Q}\times S_{Q}$ symmetry of the $e$-field in the bulk,
and that $g_{0}$ and $g_{A}$ differ by a transformation in this symmetry group
so that the two types of boundary conditions are locally equivalent.
Consequently, the only difference between $A$ and $0$ arises from local energy
cost associated with the domain wall between $g_{0}$ and $g_{A}$ BCs, which in
2d has constant size independent of $L$. In general, we thus have:

$S_{A}^{(n)}\sim\text{constant}.$ (25)

That is, the top boundary is always area law. In fact, as $D\to\infty$, both
partition functions are dominated by the ground-state configuration where all
spins are $g_{[r]}=e$, and we have ${\cal F}_{A}-{\cal F}_{0}=0$ corresponding
to a disentangled state. We note, in passing, that an identical calculation
for a 3D PEPS shows that the operator entanglement would scale linearly in the
system size. It is plausible that a boundary 2d PEPS, which can account for
such linear scaling of operator entanglement, would again enable an efficient
contraction.

###  C.4 Fluctuation Corrections

As we now briefly discuss, fluctuation corrections to the $D=\infty$ limit can
be viewed as an expansion in dilute gas of spin flips on top of the
$e$-polarized ground state.

For a finite number of replicas, the minimal-energy excitations are single
spin flips (1SF’s), where we replace $g_{[r]}:e\rightarrow\sigma\neq e$ at
some site $r$. Denoting $J=\log D$ and $h=\log d$, and
$c(g,g^{\prime})=C(g,g^{\prime})-C(e,e)$, the 1SF costs energy:

$\displaystyle E_{1}(g)=\begin{cases}(4J+h)c(g,e)&{\rm bulk}\\\
(3J+h)c(g,e)+h\left[c(g,g_{A/0})-c(e,g_{A,0)}\right]&\text{r $\in A/\bar{A}$
boundary}\end{cases}$ (26)

The lowest-energy spin-flips correspond to transposing a single pair of
replicas: $g=(ab)$, which have bulk energy: $E_{1}\left((ab)\right)_{\rm
bulk}=4J+h$. There are $D_{1}^{\rm bulk})=Q(2Q-1)$ different single-
transposes, leading to a corresponding degeneracy of the single SF excitations
in the bulk. Near an $\bar{A}$ boundary, the cheapest spin flip is is of the
form $(k\bar{k})$ for some $k=1,\dots 2Q$, and costs energy
$E_{1}\left((k\bar{k})\right)_{\bar{A}-{\rm bdry}}=2J+h$, and degeneracy
$D_{1}^{\rm edge}=Q$. The minimal-energy spin flips and corresponding
degeneracy near $A$ boundary are related to those near the $\bar{A}$ boundary
by the symmetry generator: $g_{A}g_{0}^{-1}$, which commutes with the bulk $e$
fields. At large $D$, we can use these excitations to approximate the free-
energy by a dilute gas of spin-flip excitations. This expansion is however
subtle in the replica limit, as permutations with a number of cycles
proportional to $m\to 0$ become dominant in the replica limit. While this
caveat prevents us from systematically computing the free energy in a
controllable way, this simple counting of low energy excitations predicts that
the coefficient of the area-law coefficient scales as

$\displaystyle\mathcal{F}_{A}-\mathcal{F}_{0}=-\log\frac{\mathcal{Z}_{A}}{\mathcal{Z}_{0}}\approx
D^{-2}d^{-2}+\mathcal{O}(D^{-3}).$ (27)

Though we are unable to explore a large enough range of $D$ with sufficient
precision in the iPEPS numerics to test this asympotic prediction in detail,
we note that the large-$D$ suppression in sample-to-sample variance of
entanglement features observed in the iPEPS numerics is in qualitative
agreement with the suppression of fluctuation contributions to the stat-mech
model at large-$D$.

###  C.5 Correlation length of random PEPS

The stat-mech mapping also enables one to estimate the correlation length-
scale for observables in random PEPS. Namely, consider computing the typical
amplitude of a correlation function:

$\displaystyle\log C_{\rm
typ}\equiv{\mathbb{E}}\left[\log\left[\frac{\langle\Psi|O_{1}O_{2}|\Psi\rangle\langle\Psi|\Psi\rangle}{\langle\Psi|O_{1}|\Psi\rangle\langle\Psi|O_{2}|\Psi\rangle}\right]\right],$
(28)

where $O_{1,2}$ are local observables on sites $1,2$, and we have explicitly
normalized the wave-function, and also divided by the (normalized) one-point
correlators:
($\frac{\langle\Psi|O_{1,2}\|\Psi\rangle}{\langle\Psi|\Psi\rangle}$) to remove
dependence on the operator norm of $O_{1,2}$.

For concreteness, and without loss of generality, let us specialize to the
case where $O_{[r]}=|s\rangle_{[r]}\langle s|$ is a projector onto physical
state $|s\rangle$ at site $r$ (and identity elsewhere). Generic observables
can be written as linear combinations of such projectors (up to a basis
transformation that can be absorbed into the randomly-drawn tensor on site
$[r]$).

One can evaluate the average of the log in $C_{\rm typ}$ via a replica trick
as outlined above for the bMPS entanglement. The result is that:

$\displaystyle C_{\rm
typ}=\exp\left[-\left(F_{O_{1},O_{2}}-F_{O_{1}}-F_{O_{2}}+F_{0}\right)\right]=\frac{\mathcal{Z}_{O_{1}O_{2}}\mathcal{Z}_{0}}{\mathcal{Z}_{O_{1}}\mathcal{Z}_{O_{2}}},$
(29)

where $F_{0}$ is the free-energy associated with the stat-mech Hamiltonian
(3), $F_{O_{1},O_{2},\dots O_{k}}$ is the same free-energy except with the
projectors inserted at sites $1,2,\dots k$, and $\mathcal{Z}=e^{-F}$ is the
associated partition function. The projectors restrict the sum over the
physical index values to $s$, which is equivalent to removing the replica
symmetry-breaking: $e$-field on that site. Equivalently, $F_{O_{1},O_{2},\dots
O_{k}}$ is given by the free-energy of the Hamiltonian discussed in the main
text, but perturbed by a term: $\Delta H=+\log d\sum_{[r]=1\dots
k}~{}C(e,g_{r})$.

At large-$D$, we can estimate the leading contribution to $C_{\rm typ}$ as
follows. The leading contribution to the stat-mech partition function is from
a uniformly $e$-polarized replica-spin configuration. Fluctuations about this
come in the form of small domains of non-$e$ polarized spins. By inspection,
only domains that include both sites $1,2$ make a non-cancelling contribution
to the ratio in (29). At large-$D$, this contribution is dominated by the
smallest such spanning domain, which is a line of flipped replica-spins,
$g_{i}\neq e$, along a short path connecting points $1$ and $2$. This domain
has a line tension $\mathcal{F}\approx\left(\log D^{4}d\right)r_{12}$, where
$r_{12}$ is the length of the shortest path through the network connecting
points $1,2$. This contributes exponentially decaying correlations:

$\displaystyle C_{\rm typ}\approx e^{-r_{12}/\xi},$ (30)

with characteristic correlation length: $\xi\approx\left(\log
dD^{4}\right)^{-1}$. Notice that the correlation length decreases with
increasing $D$. However, note that random PEPS states at large $D$ _are not_
close to product states, but, in fact saturate the maximal entanglement
allowed for the given bond-dimension PEPS.

This shows that large-$D$ random square lattice PEPS actually have rather
short range correlations, in accordance with previous studies Lancien and
Pérez-García (2021), and our numerical observations for clean random iPEPS.

##  Appendix D Random PEPS vs physically relevant ground states

An important question is to what extent the results presenting in this paper,
strongly indicating that random PEPS can be efficiently approximately
contracted, can be extended to PEPS representing ground states of physically
relevant Hamiltonians, e.g. in the context of condensed matter, materials
science and quantum chemistry. That is, can our results shed some light into
the performance of PEPS algorithms when simulating such systems? Here we
restrict our considerations to two-dimensional ground states that obey the
entanglement area law (2d ground states that violate the entanglement area
law, such as ground states in the presence of a 1d Fermi surface, are expected
to be harder to contract).

We have seen below that random PEPS have a very short correlation length $\xi$
on the order of one lattice site or less. In contrast, the correlation length
$\xi$ in a physically relevant ground state can be arbitrarily large (for
instance, the correlation length diverges as we approach a quantum critical
point). Relatedly, we have numerically seen that the boundary MPS for a random
PEPS has very limited amount of entanglement whereas, as PEPS practitioners
have learned over the last 15 years, the entanglement in the boundary MPS for
a physically relevant 2d ground state can again be arbitrarily large (even in
those cases where the boundary MPS obeys an area law).

We have therefore identified two structural differences between random PEPS
and physically relevant ground states, namely differences in correlation
lengths $\xi$ and in the amount of boundary MPS entanglement. How fundamental
are these structural differences? Based on experience with renormalization
group, random-circuit dynamics, and random matrix theory it is tempting to
conjecture that the large-$D$ random PEPS might represent a sort of coarse-
grained “fixed-point” representation of physically relevant ground states.
However, below we will see that while coarse-graining the PEPS for a
physically relevant ground state would indeed effectively remove the
difference in correlation length, it would not change the difference in
boundary MPS entanglement. Since boundary MPS entanglement determines the
computation cost in approximate PEPS contractions, we cannot conclude that our
results for random PEPS apply to such PEPS.

That is not to say that the stat-mech approach used in this paper to
successfully characterize the boundary MPS entanglement for random PEPS is
restricted to studying states with a short correlation length $\xi$. On the
contrary, as discussed below, we will see that the same approach can be used
for 2d random tensor network states (which are not 2d PEPS) with arbitrarily
large correlation length $\xi$.

###  D.1 Coarse-graining by blocking tensors

First, note that, any finitely correlated PEPS, i.e. with finite correlation
length, $\xi$, can be transformed into a PEPS with shorter correlation length
$\xi^{\prime}\sim 1$ by “blocking” together sites in $\xi\times\xi$ blocks of
physical sites. This blocking adds constant overhead to the bond-dimension of
each tensor, $D^{\prime}\sim D^{\xi}$. While this cost may be severe in
practice, from an asymptotic complexity standpoint, it is merely a constant
overhead. This argument suggests that one can perhaps think of a random PEPS
as reflecting a block-spin renormalization group (RG) style coarse-graining of
a physical PEPS with longer-range correlations.

However, the following observation suggests that there is no connection
between a coarse-grained PEPS for physical models with large $\xi$, and a
random PEPS with bond-dimension $D^{\prime}$. Years of numerical experience
Verstraete and Cirac (2004); Corboz and Mila (2014); Corboz  _et al._ (2014);
Niesen and Corboz (2017); Zheng  _et al._ (2017); Ponsioen  _et al._ (2019);
Chen  _et al._ (2020) show that ground-states with large correlation length
have corresponding large entanglement both for the physical state, and the
bMPS for its norm and correlation functions. While blocking reduces the
correlation length, it does not change the entanglement spectrum of the bMPS
(see Fig. 17,18). Hence, for physical states, the bMPS entanglement should
grow with $\xi$, whereas for random PEPS with bond dimension
$D^{\prime}=D^{\xi}$, the method of Appendix C.4 above predict bMPS
entanglement decreasing exponentially with $\xi$ as $\sim 1/D^{2\xi}$. This
argument shows that large-$D$ random PEPS do not have the correct entanglement
structure to capture block coarse-grainings of long-range correlated states
encountered in simulation of physical systems.

However, in standard renormalization group approaches, coarse-graining does
not simply involve merely blocking sites together, but of hierarchically
decomposing the state via layers of coarse-graining steps that act on
different distance scales. Inspired by this, in the next section, we construct
a class of hierarchical random tensor network states that have arbitrarily-
long correlation length, $\xi$, and which have bMPS entanglement that grows
with $\xi$ in a manner qualitatively consistent with that found in simulations
of physical systems (though it remains an open question whether variants of
such hierarchical tensors networks reflect all the important structure found
in physical states).

(a) (b)

Figure 17:  (a) Square lattice of the 2D PEPS norm. (b) Boundary MPS obtained
from the contraction of the norm in (a).

(a) (b)

Figure 18:  (a) Coarse grained square lattice of the 2D PEPS norm. (b)
Boundary MPS obtained from the contraction of the coarse grained norm.

###  D.2 Random Tensors Networks with large correlation lengths

In this section, we construct an ensemble of random 2d tensor networks that:

  1. 1.

have arbitrarily long correlation length, $\xi$,

  2. 2.

can be viewed as PEPS with effective bond-dimension $D_{\rm eff}\sim{\rm
poly}(\xi)$, and

  3. 3.

can be reliably analyzed by stat-mech mapping in a large-$D$ limit, which
predicts that their physical properties can be efficiently computed via an
area-law bMPS.

Specifically, inspired by expectation that an RG coarse-graining can reduce a
PEPS with any finite correlation length $\xi$ to one with $\xi\lesssim 1$, we
consider a shallow generalized multi-scale entanglement renormalization ansatz
(gMERA) architecture introduced in Anand  _et al._ (2022). A $1d$ version of
this structure is shown in Fig. 19, with obvious generalizations to higher-d.
It consists of a depth, $R$, layers of tensors, in which at layer $1\leq j\leq
R$, the tensors are connected only at distance $2^{j}$ 444We note that, while
this gMERA structure was originally introduced in the context of quantum
circuit tensor networks, and considered unitary or isometric tensors, the
isometry constraint will have little impact on the stat-mech description, and
we drop it for simplicity.. This shallow gMERA geometry allows one to neatly
interpolate between finitely-correlated states ($R$-finite) and critical
states ($R\rightarrow\infty$). Heuristically, we can view this as a discrete
version of the holographic AdS/CFT correspondence, where physical legs live
only at the boundary of a (short) extra “scale” dimension, which runs from
short-distance (UV) at the physical boundary, to longer-distance (IR). Here,
we consider the case where each tensor in this network has Gaussian random
entries, and all internal legs have bond-dimension $D$, and physical legs have
dimension $d$. Adapting the discussion of typical two-point correlation
functions above to this network, one again concludes that the correlations
decay exponentially with the size of the smallest domain that includes both
sites $1,2$. In this network, the smallest domain will run along the IR edge
of the network, resulting in:

$\displaystyle C_{\rm typ}(r)\sim\begin{cases}1/r^{p}~{}~{}{\rm
w/}~{}~{}p\sim\log_{2}D&;r<2^{R}\\\ e^{-r/\xi}~{}~{}\rm{w/}~{}~{}\xi\sim
2^{R}&;r\geq 2^{R}\end{cases}$ (31)

where $\xi\sim 2^{R}$. We note that, the functional forms listed merely
reflect an overall asymptotic decay of the envelop of correlations. In
addition, there is a complicated fractal/self-similar modulation inherited
from the geometry of the network.

From this expression, we see that the correlation length, $\xi$ can be made as
large as desired by controlling the depth of the shallow gMERA. At the same
time, this shallow gMERA can be viewed as a PEPS with effective bond dimension
$D_{\rm eff}=D^{R}=\xi^{\log_{2}D}$, which scales polynomially with the
correlation length.

Figure 19:  Random Tensor Networks with Arbitrarily Long Correlation Length –
Schematic of $1d$ gMERA (with obvious generalization to higher-dimensional
versions) which has arbitrarily-long correlation length $\xi\sim 2^{R}$ even
in the large-$D$ limit, showing that correlation length of random tensor
networks is not necessarily short. Blue dots are tensors with Gaussian-random
entries. Internal bonds have dimension $D$, physical bonds (sticking up at the
top) have dimension $d$. The holographic “scale” dimension has size $R$. The
correlation length is $\xi\sim 2^{R}$. The $2d$ version of this network can be
contracted by a boundary MPS (dashed box) with effective bond dimension is
$D_{\rm eff}\sim D^{R}$, which scales polynomially with $\xi$.

The stat-mech mapping of the bMPS entanglement for these shallow gMERA
proceeds similarly to that for the 2d square PEPS, except that the physial
legs arise only at the UV layers. Therefore, for contracting networks
representing norms and correlations, the replica-symmetry breaking $e$-field
is only present in the UV. Nevertheless, this is still sufficient to
explicitly break the replica symmetry, and give an area-law for the bMPS for
any $D_{\rm eff}$, i.e. for any correlation length, $\xi$.

The area-law coefficient may be estimated in the large-$D$ expansion as
outlined above for the square-lattice PEPS. The leading contribution again
comes form two-site domains of flipped replica spins ($g\neq e$) that straddle
the entanglement cut. A new feature is that the straddling domain can occur at
any layer in the (shallow) holographic dimension without effecting its free-
energy cost, giving rise to an entropic factor that scales as $\sim R$. The
resulting bMPS entanglement in the large-$D,R$ limit is then:

$\displaystyle S\sim\frac{1}{D^{4}}R\sim\frac{1}{D^{4}}\log\xi.$ (32)

We note that the physical dimension, $d$ does not appear in this expression
because the bulk tensors have only virtual legs. At large $\xi$, the bMPS
entanglement in (32) scales like that of a nearly-critical $1d$ system with
effective “central charge” $\sim 1/D^{4}$, with critical scaling cut off by a
finite correlation length $\xi$. In particular,

This example shows that short correlation length is not an intrinsic
limitation of random large-$D$ tensor networks. However, it remains unclear
whether this example fully capture the structure relevant to physical ground-
states.