mechanics problem with very little required structure. Similarly the complex Ising model discussed above (which is an example of a random tensor network) is close to being the most general Ising model that one could write down for this lattice geometry [123]. On the other hand, the true measurement dynamics in the MPT does have some structure (for example, structure associated with causality), which is not present in a generic tensor network. In writing down the field theory in the next section we are assuming that the only aspect of the structure of the MPT that is important for the critical theory is the shift in the number of replicas from $N = 0$ to $N = 1$ that is induced by the Born probability. This assumption should certainly be examined further.

We note that the unitary limit, $r = 0$, is a case where additional structure due to unitarity certainly is important. There the appropriate effective “spin model” has hard constraints on the allowed spin configurations, which, for example, enforce causality [34,42,44] (these are relaxed when projection operators are included [13,14]). As a result, the unitary models do not possess invariance under $O(d+1)$ rotations in space time, even in the scaling limit, and are not described by the field theories below, which do possess this symmetry. However, the unitary models do share some features with the ordered phases of these theories, such as a positive domain wall tension.

### E. A field theory for the measurement transition

With this motivation, let us write the simplest Lagrangian for $X_{ab}$, which can represent a coarse graining either of a permutation matrix or of the composite field above. We see that this simplest Lagrangian passes a basic consistency check for the MPT. (In the next section we see that we need to extend it for the RTN and the FMPT.)

Let us make subtractions so that the row and column sums of the matrix give zero:

$$ \sum_a \hat{X}_{ab} = 0, \quad \sum_b \hat{X}_{ab} = 0. \qquad (153) $$

In the case where $X$ is microscopically a permutation, this simply requires us to subtract a constant:

$$ \hat{X}_{ab} = X_{ab} - \frac{1}{N}. \qquad (154) $$

As a result of these linear constraints, which are preserved under coarse graining, $\hat{X}$ has $(N-1)^2$ independent components, and forms an irreducible representation of $G_N$. Below we omit the caret on $\hat{X}$.

Including terms in the potential only up to cubic order in $X$, and imposing $G_N$ symmetry gives a relatively simple

Lagrangian. The theory we propose for the MPT is

$$ \mathcal{L} = \sum_{ab} \left[ \frac{1}{2} (\partial X_{ab})^2 + \frac{\mu}{2} X_{ab}^2 + g X_{ab}^3 \right]. \qquad (155) $$

We include both time and space derivatives in the first term with the same coefficient, i.e., we set a nonuniversal speed to 1. This field theory has emergent Euclidean rotational invariance (not Lorentz invariance) in space time if this is not broken by boundary conditions. The components of the matrix $X$ are not independent, because of the constraints in Eq. (153). Note that as a result, in contrast to the theory discussed in the next section, the only linear term $\sum_{ab} X_{ab}$ that would be allowed by symmetry is in fact zero. The replica limit $N \to 1$ is also implied. The renormalized squared mass vanishes at the critical point, $\mu^2 \propto (r - r_c)$.

Alternatively, we may write $X$ in terms of an unconstrained $(N-1) \times (N-1)$ matrix field $\phi_{\alpha\beta}$ [124],

$$ \mathcal{L} = \frac{1}{2} \sum_{\alpha\beta} [(\partial\phi_{\alpha\beta})^2 + \mu\phi_{\alpha\beta}^2] + g \sum_{\alpha\beta\gamma} D_{\alpha\beta\gamma}^{\mu\nu\lambda} \phi_{\alpha\mu} \phi_{\beta\nu} \phi_{\gamma\lambda}. \qquad (156) $$

The tensor $D$ is a tensor product of that appearing in the cubic term of the Potts model [80,81]:

$$ D_{\alpha\beta\gamma}^{\mu\nu\lambda} = d_{\alpha\beta\gamma} d_{\mu\nu\lambda}, \quad d_{\alpha\beta\gamma} = \sum_{a=1}^{N} e_{\alpha}^{a} e_{\beta}^{a} e_{\gamma}^{a}. \qquad (157) $$

The theory with the cubic term can only make sense for the replica limit—for $N > 2$ we have an unstable potential and for $N = 2$ the cubic term vanishes. This is also the case for the Landau-Ginsburg-Wilson-like theory for percolation, which we already discuss in Sec. III B. Like that theory, the upper critical space-time dimension of Eq. (155) is $D = 6$.

A basic consistency check on our picture is that this theory indeed sustains a stable ordered phase, with the simple pattern of symmetry breaking described in previous sections, when $\mu^2 < 0$. That is, the masses of fluctuations about the ordered state should remain positive in the replica limit $N \to 1$: otherwise some more complex pattern of symmetry breaking might be required [125–127]. To check this we put $X_{ab} = f(\delta_{ab} - 1/N) + W_{ab}$, where $f$ is the magnitude of the order parameter, and $W$ represents fluctuations (with $\sum_a W_{ab} = 0$, etc.). The saddle-point equation requires $f = \frac{-\mu^2}{3g} \frac{N}{N-2}$. The mass terms in the Lagrangian for $W$ are then $\mathcal{L} = \frac{-\mu^2 N}{2(2-N)} (\sum_{ab} W_{ab}^2 - 2 \sum_a W_{aa}^2)$. We may check that the eigenvalues of the mass matrix appearing here are indeed positive when $\mu^2 < 0$ and $N \to 1$ (Appendix 1), so this consistency check is satisfied.

Now we consider another important consistency check.
## F. Counting fields

Above we start with an $N \times N$ matrix $X_{ab}$ transforming under $G_N$ symmetry. For integer $N > 1$ we may split a general such matrix into four distinct fields, transforming under distinct representations of $S_N \times S_N$:

$$ 
\begin{aligned}
S &\equiv \sum_{ab} X_{ab}, \\
R_a &\equiv \sum_{b} X_{ab} - \frac{S}{N}, \\
C_b &\equiv \sum_{a} X_{ab} - \frac{S}{N}, \\
\hat{X}_{ab} &\equiv X_{ab} - \frac{C_a + R_b}{N} - \frac{S}{N^2}.
\end{aligned} 
\quad (158) $$

The last of these, $\hat{X}$, is in the fundamental (standard) representation for both $S_N$ factors. It lives in an irreducible representation of $G_N$ of dimension $(N-1)^2$. R and C each transform under only one of the $S_N$ factors. Since they are exchanged by the $\mathbb{Z}_2$ generator, together they form a single representation of $G_N$ of dimension $2(N-1)$. S is a singlet.

For the MPT we construct a Landau theory that contained only the field $\hat{X}$. This is the obvious thing to do for various reasons (for example, if we think of $\hat{X}$ microscopically as a permutation matrix, then R, C, and S are trivial constants). We conjecture that for the MPT $\hat{X}$ is the only field that becomes massless at the critical point.

However, the group theory at $N \to 0$ [128] gives additional constraints, which strongly suggest that all of the representations in Eq. (158) become simultaneously massless at the critical point, so that we cannot throw away the representations R, C, and S. Therefore, we have to work with a general matrix X in which the row and column sums are not fixed to zero. The first indication of this is that the subtractions in Eqs. (154) and (158) diverge when $N \to 0$.

As with many other replica field theories, the partition functions that we are interested in become trivial—exactly equal to 1—in the replica limit, for certain choices of boundary conditions. An unusual feature of the circuit models with measurements or forced measurements is that this occurs at two values of N. When $N \to 0$ (FMPT) it occurs for the usual reason—because the partition function is the average of something raised to the power zero. When $N \to 1$ (MPT) the partition function is the sum of the probabilities of all the measurement outcomes—again giving 1 but for a different reason.

The fact that the microscopic partition function is equal to 1 implies constraints on the spectrum of operators in the continuum theory [118,128–130]. Here a minimal heuristic point will be sufficient: there should not be any massless fields left when N is set equal to $N_*$, the desired number of replicas, otherwise we have a nontrivial free energy, contradicting $\mathcal{Z} = 1$.

The Lagrangian, Eq. (155), for the MPT satisfies this condition, since the field is in a representation of dimension $(N-1)^2$, which tends to zero when $N \to 1$. Therefore, it passes this basic consistency check.

At first we might assume that the same field theory can also be continued to $N=0$ in order to describe the RTN and FMPT. However, this is not the case. Since $(N-1)^2$ is equal to one in this limit, rather than zero, this is not consistent.

However, the total multiplicity of all the representations in Eq. (158) is just $N^2$ (the number of components of the matrix), which does tend to zero in the replica limit $N \to 0$. This suggests that we should write a Lagrangian for a matrix X without imposing any condition on its row or column sums [131]. This is what we do next.

## G. Field theory for random tensor network and FMPT

Let us denote the unconstrained real $N \times N$ matrix by Y, to distinguish it from the matrix X above, which obeyed linear constraints. Assuming only $G_N$ symmetry and no constraints on Y, we argue below that the most relevant terms as $N \to 0$ are contained in

$$ \mathcal{L} = \sum_{ab} \left[ \frac{1}{2} (\partial Y_{ab})^2 + r Y_{ab} + g Y_{ab}^3 \right] + \frac{m_F^2}{2} \sum_{abcd} Y_{ab} F_{ab,cd} Y_{cd}, \quad (159) $$

where F is the tensor

$$ F_{ab,cd} = \delta_{bd} + \delta_{ac}. \quad (160) $$

The parameter that drives this theory off criticality is r, the coefficient of the linear term (not to be confused with the measurement rate in previous sections, also denoted r). Since no constraint is imposed on Y, this linear term does not vanish (contrast Sec. E). The term $\sum_{ab} Y_{ab}^2$ is absent because its coupling can be set to zero by a shift $Y_{ab} \to Y_{ab} + C$ with a constant C, i.e., it is redundant [132]. Surprisingly, we find below that for this theory the upper critical dimensionality of space time is $D=10$.

A peculiar feature of the $N \to 0$ limit of Eq. (159), which is shared with some other replica field theories such as the Landau-Ginsburg formulation of the random-field Ising model [82], is the presence of a quadratic coupling, which is not zero at the critical point and which cannot be removed. This is the term $m_F^2 Y.F.Y$.

If we instead study the above theory for a larger value of N, for example in the $N \to 1$ limit, then the effect of $m_F^2$ is simply to give a mass to certain representations in the decomposition of Y. The corresponding fields can therefore be eliminated at large scales and low momenta. Doing so returns us to the critical theory proposed in Sec. E for the MPT, with $\mu^2 \sim -r$. This is shown explicitly in Appendix 2. However, writing the propagator explicitly shows that
the limit $N \to 0$ that is of interest to us in this section does not commute with the limit of small momentum [82]. Therefore, we have to retain the YFY term explicitly.

Note that this term, which can be written

$$ \sum_{abcd} Y_{ab} F_{ab,cd} Y_{cd} = \sum_{a} (\sum_{b} Y_{ab})^2 + \sum_{b} (\sum_{a} Y_{ab})^2, $$

includes contributions such as $Y_{12}Y_{13}$: this is consistent with the “repulsion” that is discussed heuristically towards the end of Sec. VID, between pairing patterns involving a given layer.

$G_N$ symmetry allows many other terms at order $Y^3$ but we argue that in the $N \to 0$ limit they contribute only less relevant couplings. The dimensional analysis may be simplified using an approach [133–136] introduced by Cardy for the field theories of the random-field Ising model [137,138] and the branched polymer [139–141]. Since decomposition into representations of $S_N$ fails in the $N \to 0$ limit, the next best thing is to exploit a decomposition into representations of an $S_{N-1}$ subgroup acting on indices 2, ..., N. Here we must do this for both the row and column indices of Y.

We make a linear transformation to rewrite the field $Y_{ab}$ as a field $y_{\alpha\beta}$ whose indices $\alpha$ and $\beta$ take values in the set $\{+, -, 2, \dots, N\}$:

$$ y_{\alpha\beta} = \frac{1}{2} \mathbf{v}^{\alpha} \cdot \mathbf{Y} \cdot \mathbf{v}^{\beta}. \qquad (161) $$

The index values + and – denote two distinct linear combinations that are invariant under $S_{N-1}$, while the values 2, ..., N are permuted by $S_{N-1}$. The vectors $\mathbf{v}^\alpha$ are

$$ \begin{aligned} \mathbf{v}^+ &= \frac{1}{2} \left( 1, \frac{1}{N-1}, \dots, \frac{1}{N-1} \right), \\ \mathbf{v}^- &= \frac{1}{2} \left( 1, \frac{-1}{N-1}, \dots, \frac{-1}{N-1} \right), \\ \mathbf{v}^i &= (0, \dots, 0, 1, 0, \dots, 0) - \frac{1}{N-1} (0, 1, \dots, 1), \end{aligned} \qquad (162) $$

where in the last line the extra “1” is in the i'th place. The $N-1$ vectors $\mathbf{v}^2, \dots, \mathbf{v}^N$ add up to zero, so span only an $N-2$ dimensional space. Below, indices $i, j, k$ always run over 2, ..., N. Technical details are in Appendix 3.

After this rewriting, the terms in Eq. (159) up to quadratic order in Y become

$$ \mathcal{L}_{\text{quadratic}} = \mathcal{L}^{(1)} + \sum_{j} \mathcal{L}_{j}^{(2)} + \sum_{k} \mathcal{L}_{k}^{(3)} + \sum_{jk} \mathcal{L}_{jk}^{(4)}, \quad (163) $$

with

$$ \begin{aligned} \mathcal{L}^{(1)} ={}& (\partial y_{++}) (\partial y_{--}) + (\partial y_{+-}) (\partial y_{-+}) \\ & + 2ry_{--} + 2m_F^2 y_{--} (y_{+-} + y_{-+}), \end{aligned} \qquad (164) $$

$$ \mathcal{L}_j^{(2)} = \frac{1}{2} (\partial y_{j+}) (\partial y_{j-}) + \frac{m_F^2}{2} y_{j-}^2, \qquad (165) $$

$$ \mathcal{L}_k^{(3)} = \frac{1}{2} (\partial y_{+k}) (\partial y_{-k}) + \frac{m_F^2}{2} y_{-k}^2, \qquad (166) $$

$$ \mathcal{L}^{(3)} = \frac{1}{4} (\partial y_{jk})^2. \qquad (167) $$

Because of the linear constraints $\sum_j y_{j+} = 0$ etc., sectors (2) and (3) each contain $N-2$ copies of the same theory, and sector (4) contains $(N-2)^2$ copies of the same theory [142]. In the above rewriting, terms with couplings that vanish as $N \to 0$ are dropped [133].

Before writing the interaction terms, we use the quadratic terms to assign engineering dimensions to the various fields (see Appendix 3 for details). We assign dimensions $x_{\alpha\beta}$ such that all the quadratic terms in the Lagrangian Eq. (163) are marginal at the critical point $r=0$. This gives

$$ x_{\alpha\beta} = w_{\alpha} + w_{\beta}, \qquad (168) $$

with (recall that in the case of a circuit $D=d+1$ is the space-time dimension)

$$ w_+ = \frac{D-6}{4}, \quad w_i = \frac{D-2}{4}, \quad w_- = \frac{D+2}{4}. \qquad (169) $$

The RG eigenvalue of a cubic interaction term $y_{\alpha\beta}y_{\alpha'}\beta'Y_{\alpha''\beta''}$ is then determined by the difference in the number of + indices it contains and the number of – indices it contains among $\alpha, \dots, \beta''$ (Appendix 3).

However, the terms that can appear are constrained by the $G_N$ symmetry (whose effects are less obvious in the new representation). We confirm in Appendix 3 that the cubic term $g \sum_{ab} Y_{ab}^3$ shown in Eq. (159) is strictly more relevant than the other symmetry-allowed cubic terms (at least for large enough $D$) and is of the form

$$ g \sum_{ab} Y_{ab}^3 = \frac{g}{2} [6y_{++} (y_{++}y_{--} + 2y_{+-}y_{-+}) + 6y_{++} \\ \times (y_{-k}y_{+k} + y_{j-k}y_{j+k} + y_{jk}y_{jk})/4 \\ +3 (y_{+-}y_{j+}y_{j-} + y_{-+}y_{+k}y_{k-}) + 3y_{j+}y_{+k}y_{jk}] \\ + \text{less relevant terms.} \qquad (170) $$

The RG eigenvalue of $g$ is $(10-D)/2$, so the upper critical dimension for this theory is $D=10$.
## H. Consequences of the MPT field theory

We discuss some simple consequences of the putative field theory for the MPT, deferring a detailed analysis, and a discussion of the more complicated theory in the previous section, to another time. However, first we note an important caveat to the discussion.

Our initial hope is that the large-$D$ limits of these field theories would give exact results both for the all-to-all circuits and for tree tensor networks. For example, this is what we find for the classical minimal cut toy model (because all-to-all percolation could be understood using the field theory for percolation in high dimensions, Sec. III E). But the class of tree tensor networks that we understand best, including those derived from the all-to-all FMPT circuit with Haar-random gates, seems not to be described by the field theory of Sec. G, simply because it is hard to imagine the exponential scaling of the order parameter in Eq. (52) being reproduced by a mean-field treatment of Eq. (159). Therefore, it seems unlikely that the all-to-all circuits studied in this paper are described by the $d \to \infty$ limit of the above field theories. We do not yet understand the reason for this difference.

It is not ruled out that our Lagrangians overlook some crucial structure, and that as a result they do not capture any models of measurement circuits or random tensor networks, even in finite dimensions. For present purposes we assume this pessimistic scenario does not hold, and that the two field theories in Secs. E and G do capture at least some class of models for the MPT and for the FMPT and RTN. We explore these issues further elsewhere.

The simpler of the two field theories is that in Sec. E for the MPT, involving a field $X_{ab}$ with vanishing row and column sums. As a result of the cubic term, this theory has upper critical space-time dimension $d+1=6$. Interestingly, the logic of Sec. III E for the percolation problem above five spatial dimensions applies in this case too, since it relies only on the engineering dimensions of the fields. We can therefore carry over the exact exponent values so that (neglecting physics on timescales shorter than $L=N^{1/d}$, see Sec. III E) the natural scaling variables in high dimensions are again

$$t/N^{1/5}, \quad N^{2/5}\delta r, \quad (171)$$

where $\delta r = r - r_c$ is the parameter driving the transition (and the number $N$ of spins should not be confused with the replica number in the preceding sections).

Let us consider the operator entanglement in the ordered phase, still above the upper critical dimension. The plateau value of the operator entanglement, $S_2 \sim sN$, is proportional to the energy cost of a domain wall in $X$ that spans the system in the spatial directions, as discussed in the following section. In high dimensions the scaling of $s$ follows from dimensional analysis, giving $s \propto \mu^5 g^{-2}$, Eq. (155),

or in terms of the deviation $\delta r$ from criticality,

$$s \sim |\delta r|^{5/2}, \quad (172)$$

which is the same exponent as for the classical problem in high dimensions. A similar scaling form will again apply, $\bar{S}_2 = H[t/N^{1/5}, N^{2/5}\delta r]$, but with a different scaling function $H$. The size of the order parameter $X$ itself, which may be measured using appropriate correlation functions, grows linearly with the distance from the critical point, $X \sim |\delta r|$.

Again we have a characteristic timescale $\tau = N^{1/5}W$ ($N^{2/5}\delta r$), for an appropriate scaling function $W$. In the entangled phase this timescale grows exponentially in $N$ (Sec. I), with

$$\ln \tau \sim sN \sim |\delta r|^{5/2}N \quad (173)$$

close to the transition.

Below $5+1$ dimensions the scaling is different, because the cubic term is no longer dangerously irrelevant. The appropriate scaling variables are as usual

$$t/L, \quad L^{1/\nu}\delta r, \quad (174)$$

where $\nu$ is the correlation length exponent for the field theory, Eq. (155). Exponents could be computed in a $6-\epsilon$ expansion and will differ from percolation exponents (since the structure of the field theory is different, despite sharing the same upper critical dimension). In the ordered phase there is still an exponentially long timescale, with

$$\ln \tau \sim sN \sim |\delta r|^{\nu d}N \quad (175)$$

close to the transition (Sec. I).

## I. Long timescale in the entangled phase

So far in this section we focus on the continuum description close to the transition. Here we discuss something simpler, namely the emergence of a timescale that (in the entangled phase) is exponentially large in the number of spins, and the contrasting short timescale in the disentangled phase. We may consider either a model in $d$ spatial dimensions with $N = L^d$ spins, or the all-to-all model. The results in this section are independent of the conjectural field theories above, as they rely only on more basic features of the effective spin model (pairing field) descriptions.

The appearance of a long timescale may be understood in analogy to standard 1D or quasi-1D classical models. Here the 1D coordinate is time: see Fig. 30.

In the ordered phase the pairing field (either $\sigma$ on the lattice or $X$ in the field theory) has long-range order across a temporal slice and, after coarse graining sufficiently, we may think of it as a function only of time. There is then
a competition between the free-energy cost of imposing a domain wall at a particular time, which scales as $sN$ with $s > 0$, and the entropy $\ln t$ associated with translating the domain wall in the time direction. At a timescale $\tau$ with

$$ \ln \tau \propto sN \qquad (176) $$

the translational entropy wins, and domain walls proliferate. Long-range order then no longer extends from the initial to the final time. By the identification of the entanglement with a free energy, this also means that the entanglement begins to decay exponentially with time.

Recently the exponentially long timescale in the entangled phase has been discussed from several points of view. References [12] and [56] consider a limit where the unitary evolution during a unit time can be treated as a $2^N \times 2^N$ Haar random unitary (see also Appendix E here for related considerations). Reference [27] has also given an analysis in terms of Ising domain walls that is similar to our considerations below.

The proportionality in Eq. (176) allows for an order-1 constant: however we expect that $N^{-1} \ln \tau$ vanishes in the same manner as $s$ when the critical point is approached from the entangled side (for example, with the same power of the tuning parameter when this dependence is a power law).

At times sufficiently shorter than $\tau$ the operator entanglement entropy has a plateau at an extensive value. The plateau value is corrected by a negative subleading term whose magnitude grows logarithmically with time. In terms of the pairing field, the plateau regime is that where the number of domain walls is the minimal number allowed by the boundary conditions.

For $S_2$, in the plateau regime, it is in fact sufficient to think about an Ising domain wall in a system with (Ising symmetric) disorder, for reason discussed towards the end of Sec. IV B. That is, we expect that the replica trick can be avoided in the strongly entangled regime. It is also possible to argue for the Ising picture using the replica treatment, by arguing that in this regime the replica theory is equivalent to the replica representation of a disordered Ising model [34]. In the replica treatment we have $m \to 0$ “elementary” domain walls, each associated with one of the transpositions in $\tau_{2,m}$ [Eq. (137)]. These may either bind together, forming a composite domain wall, or may separate for entropic reasons. See Ref. [34] for a discussion in the unitary case. However, this system of $m \to 0$ domain walls can really just be thought of as the replica description of a single Ising-like domain wall in a disordered environment. See Refs. [34,43] for details in the unitary case.

If we neglect quenched disorder, then we obtain

$$ \bar{S}_2 \sim sN - \ln t, \qquad (177) $$

in the plateau region. The second term is the contribution from translational entropy, arising because the centre

of mass temporal coordinate $t_{DW}$ of the domain wall can be located anywhere in $(0, t)$. The form in Eq. (177) was obtained in Ref. [12] in a limit of very dilute measurements, where the system can be viewed as completely scrambled by a random unitary between each measurement. Reference [27] gave a picture in terms of Ising domain walls equivalent to the one presented above. Here we also suggest how the effective Ising model can be justified (in an appropriate regime and at the level of universal properties) rather than being only a heuristic model. Our consideration also implies that we should take into account quenched disorder, as discussed below. (For another application of domain wall entropy in an effective 1D model to quantum chaos, see Refs. [117,143].)

As a check on the replica picture, we also consider a toy model for the entangled phase that involves multiplying large random matrices. A crude treatment in Appendix E (which neglects spatial structure, random fluctuations, and also the $n$ dependence of the Rényi entropies) reduces to computing the singular values of a sub-block of a large Haar-random unitary. This treatment also yields Eq. (177), and shows that the plateau value $sN$ determines the timescale for exponential decay of $S_n$ in the regime of much later times, as expected from the above. An analysis of related random matrix models has recently been presented in Ref. [56].

Equation (177) is the simplest picture, neglecting quenched disorder. In reality there will be more complex crossovers. For example, in the all-to-all model there may be a regime of timescales where the subleading correction is not $\ln t$ but instead proportional to $\sqrt{N} \ln t$ as in the classical minimal cut problem (Sec. III D).

This is because the conditional free energy $F(t_{DW})$, given by fixing $t_{DW}$, will vary with $t_{DW}$ due to randomness: $F(t_{DW}) = sN + \eta(t_{DW})$. In high enough dimensions, and therefore presumably also in the all-to-all model, the typical fluctuations $\eta(t_{DW})$ will be Gaussian with a scale $\sqrt{N}$. Although these fluctuations are much smaller than $N$, they are in principle much larger than 1. Therefore, at early enough times the free energy will be dominated by the optimal (most negative) value of $\eta(t_{DW})$, rather than by translational entropy (see endnote [144]). But at larger times, there may be a regime where $\ln t$ entropy again dominates, giving the functional form in Eq. (177). At still larger times multiple domain walls will proliferate (and the full replica treatment is required) and eventually $S_2$ decays exponentially in time.

The fact that only a single domain wall plays a role in the plateau regime means that there is an approximate factorization property for $S_2$ in a given realization of the circuit. If we divide $V$ into two parts, $V^{(1)}$ corresponding to evolution from 0 to $t'$ and $V^{(2)}$ from $t'$ to $t$, then $e^{-S_2} \simeq e^{-S_2^{(1)}} + e^{-S_2^{(2)}}$. The first term includes configurations with $t_{DW} \in (0, t')$ and the second those with $t_{DW} \in$
(t', t). (This is approximate not only because it neglects configurations with multiple domain walls but also because it does not correctly treat domain walls with $t_{DW}$ close to t'.)

We now contrast the properties of the disentangled phase with those of the entangled phase. Let us take the limit $N \to \infty$ first, so that as usual we can define the operator entanglement per spin at a given time:

$$s_2(r, t) = \lim_{N \to \infty} \frac{\overline{S_2(r, t, N)}}{N} \qquad (178)$$

(we write this equation for $S_2$, but the choice of Rényi index $n \ge 1$ should not be crucial). In contrast to the quasidiag limit discussed above, this is the free-energy cost, in an infinite slab of finite thickness, of imposing the domain-wall boundary conditions described in Sec. IV B. In the disentangled phase the free-energy cost per unit transverse area decays exponentially with the thickness of the slab, so that $s_2(r, t)$ decays exponentially to zero with time.

## J. Variants and comments

In this subsection we discuss a few extensions of the field-theory approach we present, as well as some open questions.

The measurement problems and random tensor networks that we discuss so far have no internal global symmetries. One could also consider, say, measurement dynamics with an Ising symmetry [38,41]. The definition of the pairing field in Sec. VID allows such symmetries to be incorporated, and suggests that in many cases they will change the universality class of the entanglement transition.

For example, if the tensor network in Sec. VID has a $\mathbb{Z}_2$ Ising symmetry that changes the sign of $S^a$ (and if we assume that the field whose mass vanishes at the transition is still $X_{ab} \sim S^a \overline{S^b}$) then odd powers of $X$ are forbidden by symmetry in the continuum Lagrangian, which completely changes its structure in the limits of both $N \to 0$ and $N \to 1$.

This symmetry consideration highlights a feature of the discussion in Sec. VID, which is that the definition of $X_{ab}$ involves choosing a local basis. In many cases this choice may not seem natural: for example, in many random models, the statistical invariance property emphasized in Sec. I ensures that any choice of local basis is equivalent to any other. (The exact mappings to models of permutations avoid having to choose a basis, but on the other hand it is less obvious how to coarse grain them.) An open question is whether this necessity of choosing a basis is just an aesthetic issue, or a fundamental one. Is it possible, for example, that the statistical invariance property imposes constraints on the continuum theory that we neglect to take into account?

Other restrictions on the unitaries, not related to conventional symmetries, can also change the symmetries of the replica theory. For example, if all the unitaries are real valued [145] then there is no distinction in the bulk between forward and backward layers. The symmetry group $G_N$ is then enlarged to $S_{2N}$. In this case we can introduce a pairing field in a similar manner to Sec. VID, now with a replica symmetry action like that in standard disordered magnets and spin glasses. (The restriction to Clifford unitaries [7,12,15,17,22] is a more drastic change, which may require a different theoretical approach.)

The picture in Sec. VID relates random tensor networks (for which the limit $N \to 0$ is the appropriate one) to the language typically used to discuss spin glasses. This relation raises the question of whether other types of replica symmetry breaking, or other types of glass transition [127], are relevant to natural choices of circuit or tensor network. (Of course we could always engineer, say, a glassy phase if we specifically design a tensor network with this in mind.) For example, one could imagine a second transition taking place inside the entangled phase for some choices of tensor network. At the entanglement transition, the $2N$ layers form a collection of $N$ pairs, breaking $G_N$ symmetry down to $S_N \times \mathbb{Z}_2$. Can the residual $S_N$ symmetry be broken in a subsequent transition? What are the entanglement properties of the resulting (presumably glassy) phase?

A statistical mechanics problem that provides a possible analogy for some of these phenomena is the directed polymer with random complex (or random sign) weights [146–151]. The replica formulation of this problem involves $N$ copies of the polymer’s partition function and $N$ copies of its complex conjugate. Averaging over random phases forces the copies to form pairs in order to avoid phase cancellation [148], in analogy to the pairing phenomenon in the circuits. Further, the paired object—a bound state of polymers from different copies—may itself undergo phase transitions due to disorder. Perhaps this simpler problem can provide lessons for the circuit.

## K. Free-fermion measurement dynamics

Models of free fermions subjected to stochastic dynamics [37,38,52–54,152–156] can also show a transition in $d > 1$ between two phases with differing amounts of entanglement [37]. However, instead of an area law and a volume law phase (for states in finite dimensions), we instead have an area law phase and a phase with a logarithmic violation of the area law [37,53,54].

We may also characterize the two phases by transmission of information between initial and final time, which gives a distinction that makes sense in any dimension or for the all-to-all setup. For concreteness we may consider the latter case. The model of Ref. [37], which used the language of Majorana fermions, has a simple field-theory description that is related to a model of classical
<figure><img src="image_34.png" /><figcaption>FIG. 31. The field-theory description of the Majorana measurement model of Ref. [37] has a continuous replicalike symmetry, allowing smooth domain walls that give a more rapid decay of $S_n$ than in the interacting case where replica symmetry is discrete (Appendix B2). This exhibits a more general feature of free-fermion models.</figcaption></figure>

loops (random walks) representing Majorana worldlines. The quasi-one-dimensional regime, which is relevant here has been studied in depth in Ref. [157], which also characterizes the statistical properties of random samples. Here we consider only some more basic average quantities.

The characteristic timescale for the operator entanglement to decay is of order $N$, where $N$ is the number of lattice sites, rather than being exponentially large in $N$ as we find in the interacting case. This is a generic feature of free-fermion models, as discussed below. Within the “more entangled” of the two phases, the scaling of the operator entanglement is

$$S_n \propto \begin{cases} K(r)N/t & t \ll K(r)N \\ \exp\{-ct/[K(r)N]\} & t \gg K(r)N. \end{cases} \quad (179)$$

Here $K$ (the $\sigma$ model stiffness) is an order-1 constant deep in the phase, and vanishes as $K(r) \sim (\delta r)^2$ upon approaching the transition at $r = r_c$ to the disentangled phase ($c$ is a fixed order-1 constant). Note that the scaling in Eq. (179) is identical to the conductivity of a disordered $N$-channel wire, showing the crossover from Ohm’s law to localization on a timesale of order $K(r)N$ [157,158].

The reason for the reduced timescale in the entangled phase (of order $N$ compared to the exponential timescale in interacting models) is that the appropriate replica field theory has continuous, rather than discrete, replica symmetry. In the ordered regime, a nonlinear $\sigma$ model description may be used. Domain walls are smooth objects whose free-energy cost decreases with their thickness, which in the case of interest is the temporal duration $t$ of the evolution: see Fig. 31.

For this reason, we anticipate that the scaling in Eq. (179) applies to more general free-fermion models with measurement. [The scaling of $K(r)$ close to $r_c$ will depend on symmetries and dimensionality. The constant $c$ may also depend on $n$ in general.] General free-fermion models can be formulated using the replica trick, in close analogy to replica $\sigma$ models for Anderson localization [159], leading to continuous replica symmetries. However,

in addition to the $N \to 0$ limit familiar from localization, the $N \to 1$ limit is now also of interest. We discuss this elsewhere [119].

The timescale of order $N$ for free fermions agrees with the recent results of Ref. [56], which studied a model in which measurements of a single fermionic mode were alternated with Gaussian unitaries acting on the entire system. This model has even less locality structure than the all-to-all circuit. In this limit also, the authors found that $O(N^2)$ measurements were required to forget the initial state: this corresponds to $t = O(N)$ in our conventions.

## VII. OUTLOOK

It remains an open question to what extent the properties of the MPT, in various settings, will turn out to be tractable (either analytically or numerically). In this paper, however, we show that exact results are possible in certain regimes. We close by summarizing the regimes we study, and some of the outstanding questions.

We begin our analysis by considering the “classical limit” of the MPT in the all-to-all setting. We show that a fairly complete picture is possible, including an analytical derivation of the critical point, critical exponents, and scaling forms for the entanglement.

Our results for quantum trees, including those obtained from a spin-1/2 all-to-all circuit, show that exact results are also possible even far from this classical limit. In this setting it is possible to demonstrate that an entanglement transition occurs at a definite nonzero measurement rate that is distinct from the classical value. (It may even be possible to obtain rigorous results on the phase diagram using the recursion relation approach.) The critical scaling on the tree is qualitatively different from a simple percolation picture.

We argue that the critical point on the tree is the same as the critical point of the FMPT in the all-to-all circuit (which is locally treelike). Since the location of the critical point in the circuit is difficult to check numerically, this equivalence has not yet been demonstrated clearly by our numerics. In the future we would like to have a clearer demonstration (or disproof) of this relationship between the tree and the all-to-all circuit. Our results based on the tree are also restricted to the FMPT; it would be interesting to understand to what extent they are relevant to the MPT.

The scaling on the tree raises several questions that we hope to return to elsewhere. First, it will be worthwhile to examine the relationship between the random recursion relation studied here and approaches to tree tensor networks based on replicas [68]. Second, we raise the question of whether there are multiple universality classes on the tree. This question remains to be settled, and could perhaps be addressed by generalizing our approach to a broader class of trees (with more general distributions of tensors or with larger bond dimension). Finally, it remains to be
understood how to reconcile the scaling that we find on the tree with field theory.

In our numerical study of the MPT we propose observables that have benefits over the state entanglement, in that they do not require one to specify a spatial subregion. (Constructing such observables is crucial in the all-to-all setting, for which there is no meaningful distinction between area-law and volume-law scaling, but they are also useful in 1+1D, where significant finite-size effects make it important to avoid introducing length scales that are smaller than the system size.) We demonstrate numerically that there is a long timescale in the entangling phase over which some aspects of unitarity are retained; for example, two initially orthogonal states remain approximately orthogonal.

The optimal numerical protocol for studying critical properties in the all-to-all circuit remains to be settled. One complication is the lack of a priori knowledge of how the characteristic timescale scales with N when r = r_c. In the 1+1D problem, establishing that the dynamical exponent is equal to unity [6,15] allows one to reduce the number of independent variables in scaling collapses by fixing t/L to a constant. Our candidate field theory for the MPT suggests that in high dimensions the appropriate scaling variable is t/N^{1/5}, but it is unclear whether this theory applies to the all-to-all circuit.

The proximity of the classical critical point ($r_c^{cl} = 0.8$) to the quantum one (e.g., $r_c = 0.749$ for the FMPT with Haar-random gates) in the ensembles we study may also complicate the numerical analysis. For this reason it might be useful to study an all-to-all model (for example, involving weak measurements) in which the classical transition is eliminated entirely. It will also be interesting to relax the unitary invariance property of the gate distribution: the strong constraints imposed by this invariance are a surprising feature of our analysis of the quantum tree.

Finally, we address the replica approach to the MPT and to random tensor networks, both in the two phases and near the critical point, and we make concrete proposals for field theories for these problems. Determining the domain of applicability of these theories will be the subject of further work. A basic ingredient in Sec. VI is the construction of an “overlap” order parameter for the MPT and random tensor networks that brings these problems closer to the language we use for disordered magnets and spin glasses (in comparison with the more abstract language of permutation group elements used so far): this point of view may shed light on new possibilities for ordering.

## ACKNOWLEDGMENTS

We thank Patrick Draper, John Chalker, Andrea De Luca, Tianci Zhou, Toby Swann, Michael Gullans, David Huse, Ehud Altman, Andreas Ludwig, Romain Vasseur,

Yi-Zhuang You, and David Kesler for useful discussions. A.N. is supported by a Royal Society University Research Fellowship. S.R. is supported by EPSRC Grant No. EP/S020527/1. J.R. acknowledges the funding of the Israeli Science Foundation under Grant No. 994/19 and the support of the Alon fellowship by the Israeli Council of Higher Education.

# APPENDIX A: MORE ON CLASSICAL PROBLEM

## 1. Density of infinite cluster

Here we briefly derive Eq. (18), which describes the probability $f_\infty$ that a given node in the interior of the classical graph is connected to an infinite number of other nodes in the limit of infinite N and T. In other words, $f_\infty$ describes the density of the infinite cluster.

Consider the process of building a tree starting with an arbitrarily chosen node, as depicted in Fig. 4(c). The starting node has four possible edges, each of which may be severed by a measurement. If we denote by $e_\infty$ the probability that following a given edge will lead to a subtree with an infinite number of nodes, then

$$f_\infty = 1 - (1 - e_\infty)^4. \quad (A1)$$

The quantity $(1 - e_\infty)^4$ denotes the probability that none of the four edges connected to the starting node leads to an infinite number of other nodes.

Following a particular edge, one may next encounter either a measurement [with probability $p = r/(2-r)$] or a node (with probability $1-p$). The probability that this node is connected to an infinite number of other nodes at later generations is given by $1 - (1 - e_\infty)^3$. Thus we can write a self-consistency relation for $e_\infty$, given by

$$e_\infty = (1-p)[1 - (1 - e_\infty)^3]. \quad (A2)$$

Near the critical point, $p = 2/3 + \delta p$, where $\delta p = (25/18)\delta r$ and $\delta r = r - r_c \ll 1$. On the disconnected side of the transition, $e_\infty = 0$, while just on the connected side (small negative $\delta r$) $0 < e_\infty \ll 1$. Expanding Eq. (A2) for small $\delta r$ gives $e_\infty \simeq -(25/6)\delta r$. A similar expansion of Eq. (A1) gives Eq. (18) of the main text.

## 2. Effective 1D field theory

Here we derive the mapping between the “layered Erdős-Rényi” percolation model and the one-dimensional field theory that is described in Sec. III B.

This is a bond percolation model with sites labeled $(i, t)$ with $i = 1, \dots, N$ and $t = 1, \dots, T$. Generalizing slightly from the case in the text let a bond between sites $(i, t)$ and $(j, t)$ on the same timeslice be present with probability $b/N$, and a bond between sites $(i, t)$ and $(j, t+1)$ on the next slice be present with probability $b'/2N$. The average degree of a bulk node is $z = b + b'$, and from
considerations like those in Sec. A the critical case is $z = 1$.

Bond percolation can be simply mapped to the Potts model with $Q \to 1$ states (see Ref. [82] for a review). We introduce a Potts spin $\sigma(i, t) = 1, \dots, Q$ on each site $(i, t)$, and couplings for pairs of sites that are allowed to be connected by a bond. For each pair of spins that is allowed to be connected there is a term

$$[(1-p) + p \delta_{\sigma(i,t),\sigma(j,t)}] \quad (A3)$$

in the Boltzmann weight, where $p$ is the bond probability. The two terms correspond, in a diagrammatic expansion, to the presence and absence of the bond, respectively. Sites in the same percolation cluster have the same Potts spin state because of the Kronecker $\delta$s on the bonds. Summing over spin states gives a factor of $Q^{\text{clusters}}$, which becomes 1 in the replica limit. Spin correlation functions can be used to diagnose connectivity. The probability that two sites $(i, t)$ and $(j, t')$ are in the same cluster is [82]

$$p_{\text{conn}}(i,t; j, t') = \lim_{Q \to 1} \frac{\langle \delta_{\sigma(i,t),\sigma(j,t')} - 1/Q \rangle}{1 - 1/Q}. \quad (A4)$$

Below, the limit $Q \to 1$ will be left implicit.

Using the fact that the bond probabilities are of order $1/N \ll 1$, the partition function may be written

$$Z = \sum_{\{\sigma\}} \exp \left[ \frac{b}{2N} \sum_{t=1}^{T} \sum_{i,j} \delta_{\sigma(i,t),\sigma(j,t)} + \frac{b'}{2N} \sum_{t=1}^{T-1} \sum_{i,j} \delta_{\sigma(i,t),\sigma(j,t+1)} \right]. \quad (A5)$$

As is standard in the field-theory formulation of the Potts model [80,81], it is convenient to use a set of $(Q-1)$-component vectors $\mathbf{e}_\sigma$, for $\sigma = 1, \dots, Q$ to represent the spin states, with the vectors satisfying

$$\mathbf{e}^{\sigma} \cdot \mathbf{e}^{\sigma'} = \delta_{\sigma,\sigma'} - Q^{-1}. \quad (A6)$$

For $Q=2$ we can take $e^{1,2} = \pm 1/\sqrt{2}$. For $Q=3$ the three vectors point to the three corners of an equilateral triangle. For $Q=4$ they point to the vertices of a regular tetrahedron, etc. Note that

$$\sum_{\sigma} \mathbf{e}^{\sigma} = 0, \quad \sum_{\sigma} e_{\mu}^{\sigma} e_{\nu}^{\sigma} = \delta_{\mu,\nu}, \quad (A7)$$

as we see by considering $(\sum_{\sigma} e_{\mu}^{\sigma} e_{\nu}^{\sigma}) e_{\nu}^{\tau}$ and applying Eq. (A6). Writing $\mathbf{e}(i, t) = \mathbf{e}^{\sigma(i,t)}$, and denoting the sum of

the spins in a layer by

$$\mathbf{E}_t = \sum_i \mathbf{e}(i, t), \quad (A8)$$

the partition function is (we drop an unimportant multiplicative constant)

$$\begin{aligned} Z &= \sum_{\{\sigma\}} \exp \left( \frac{b}{2N} \sum_{t=1}^{T} \mathbf{E}_t^2 + \frac{b'}{2N} \sum_{t=1}^{T-1} \mathbf{E}_t \cdot \mathbf{E}_{t+1} \right) \\ &= \sum_{\{\sigma\}} \exp \left[ \frac{b-b'}{2N} \sum_{t=1}^{T} \mathbf{E}_t^2 + \frac{b'}{4N} \sum_{t=1}^{T-1} (\mathbf{E}_t + \mathbf{E}_{t+1})^2 \right. \\ & \qquad \left. + \frac{b'}{4N} (\mathbf{E}_1^2 + \mathbf{E}_T^2) \right]. \end{aligned} \quad (A9)$$

We can use two sets of Hubbard-Stratonovich fields, one set located at half-integer times, denoted $\mathbf{f}_{t+1/2}$, to decouple the $b'$ term, and one set located at integer times, denoted $\mathbf{g}_t$, to decouple the $b-b'$ term. Each has $Q-1$ components. Once the $\mathbf{E}$ appear linearly in the exponent we can sum over the spins in a given timeslice $t$ (the prime indicates that the sum is only over these spins) via

$$\sum_{\{\sigma\}}' \exp [\mathbf{E}_t \cdot \mathbf{y}] = \left( \sum_{\sigma=1}^{Q} e^{\mathbf{e}^{\sigma} \cdot \mathbf{y}} \right)^N = Q^N \exp[NV(\mathbf{y})], \quad (A10)$$

which defines $V(\mathbf{y})$. Expanding in $\mathbf{y}$ for small $\mathbf{y}$ and using the identities mentioned above for the set of vectors $\{\mathbf{e}^\sigma\}$,

$$V(\mathbf{y}) = \frac{\mathbf{y}^2}{2Q} + \frac{d_{\mu\nu\lambda} y_\mu y_\nu y_\lambda}{6Q} + \mathcal{O}(y^4). \quad (A11)$$

The tensor $d$ is [80,81]

$$d_{\mu\nu\lambda} = \sum_{\sigma=1}^{Q} e_{\mu}^{\sigma} e_{\nu}^{\sigma} e_{\lambda}^{\sigma}. \quad (A12)$$

After integrating out the spins,

$$Z = \int \mathcal{D}(f,g) e^{-\frac{1}{2}(\sum f^2 + \sum g^2) + \sum_{t=1}^{T} NV(y_t)}, \quad (A13)$$

where the final sum is over integer $t$. For $2 \le t \le T-1$,

$$\mathbf{y}_t = \sqrt{\frac{b'}{2N}} (f_{t-1/2} + f_{t+1/2}) + \sqrt{\frac{b-b'}{N}} g_t. \quad (A14)$$

At the boundaries we have, e.g.,

$$\mathbf{y}_1 = \sqrt{\frac{b'}{2N}} f_{1+1/2} + \sqrt{\frac{b-b'}{N}} g_1. \quad (A15)$$
For the present we neglect the boundary terms. The boundary condition on the field theory is important but below we fix it on physical grounds.

The negative power of N in y will allow us to truncate the action at cubic order. Let us combine **f** and **g** into a field **h** labeled by both integer and half-integer values, **h**<sub>t</sub> = **g**<sub>t</sub>, **h**<sub>t+1/2</sub> = **f**<sub>t+1/2</sub>. The lattice field theory is then (with τ, τ' ∈ Z/2)

$$Z = \int \mathcal{D}\mathbf{h} \exp \left[ -\frac{1}{2} \sum_{\tau, \tau'} h_{\tau} T_{\tau, \tau'} h_{\tau'} - \frac{1}{\sqrt{N}} A_3(\mathbf{h}) \right], \quad (A16)$$

where $A_3$ contains the cubic terms. To avoid clutter, let us immediately set $Q = 1$ in the dispersion relation. The matrix $T$ is then (the first row and column shown correspond to a half-odd-integer index value):

$$T = \begin{pmatrix} \dots & & & & & \\ 1-b' & -\alpha/2 & -b'/2 & & & \\ -\alpha/2 & 1-b+b' & -\alpha/2 & 0 & & \\ -b'/2 & -\alpha/2 & 1-b' & -\alpha/2 & -b'/2 & \\ & 0 & -\alpha/2 & 1-b+b' & -\alpha/2 & \dots \end{pmatrix} \quad (A17)$$

Here we define

$$\alpha = \sqrt{2b'(b-b')}. \quad (A18)$$

This becomes imaginary when $b' > b$—which includes the line $b = 0$ on which we do simulations of this model—but this does not present a problem in the formal derivation below [160]. Let us write

$$b = \frac{z + \Delta}{2}, \quad b' = \frac{z - \Delta}{2}, \quad (A19)$$

where $z$ is the mean degree of a site, and the location of the critical point is $z = 1$ for any value of $\Delta$. The dispersion relation has one “massive” mode, and one mode that becomes massless at the critical point $z = 1$, at frequency $\omega = 0$, with the eigenvalue of $T$ being $(1-z) + [(z-\Delta)/4]\omega^2 + \mathcal{O}(\omega^4)$. At $z = 1$ and $\omega = 0$ the eigenvector of this mode is $(g,f) \propto (\sqrt{\Delta}, \sqrt{1-\Delta})$. For the low-frequency theory we make the coefficient of this mode a slowly varying field, $\phi(t)$. Let us write

$$\delta z = z - 1 \quad (A20)$$

for the parameter that vanishes at the phase transition. Let us drop the small parameter $\delta z$ except in the mass term, where it is the leading factor:

$$Z = \int \mathcal{D}\boldsymbol{\phi} \exp \left\{ -\int dt \left[ \frac{1-\Delta}{8} (\partial_t \boldsymbol{\phi})^2 - \frac{\delta z}{2} \boldsymbol{\phi}^2 \right] - \frac{A_3}{\sqrt{N}} \right\}. \quad (A21)$$

Thanks to the small prefactor $1/\sqrt{N}$ of the cubic term, we may take the continuum limit in a controlled manner. The cubic term is negligible for frequencies $\omega$ of order 1 due to

the small prefactor, but important at parametrically small frequencies (since it is RG relevant). Since only small frequencies are important we can simply insert the form of the low-lying mode at $k = 0$ into the cubic term without any need to explicitly integrate out high-frequency modes. The final result is

$$Z = \int \mathcal{D}\mathbf{h} \exp \left( -\int dt \mathcal{L} \right), \quad (A22)$$

with the “Lagrangian” (again a factor of $Q$ has been set to 1 in the denominator of the final term)

$$\mathcal{L} = \frac{1-\Delta}{8} (\partial_t \boldsymbol{\phi})^2 - \frac{\delta z}{2} \boldsymbol{\phi}^2 - \frac{d_{\mu\nu\lambda}}{6\sqrt{N}} \phi_\mu \phi_\nu \phi_\lambda. \quad (A23)$$

Above, all the ($Q-1$) components of the field $\boldsymbol{\phi}$ are independent. We can write a more explicit form at the cost of using $Q$ fields that obey a linear constraint (summing to zero). For notational convenience we write them as the components of a $Q \times Q$ diagonal matrix with components

$$\Phi_{\sigma, \sigma} = \boldsymbol{\phi} \cdot \mathbf{e}^{\sigma}. \quad (A24)$$

The constraint is tracelessness

$$\mathrm{tr} \, \Phi = 0. \quad (A25)$$

Equation (A23) becomes

$$\mathcal{L} = \frac{1-\Delta}{8} \mathrm{tr}(\partial_t \Phi)^2 - \frac{\delta z}{2} \mathrm{tr} \Phi^2 - \frac{1}{6\sqrt{N}} \mathrm{tr} \Phi^3. \quad (A26)$$

This is the result given in the main text for the special case $b = 0$ (i.e., $\Delta = -z \simeq -1$ close to the critical point).
Close to the critical point, the connectedness correlation function for sites at distinct times is

$$p_{\text{conn}}(i, t; j, t') = \frac{\langle \vec{\phi}(t) \cdot \vec{\phi}(t') \rangle}{N(1 - 1/Q)} \rightarrow \frac{\langle \phi_1(t) \phi_1(t') \rangle}{N}, \quad (\text{A27})$$

where $\phi_1$ is an arbitrarily chosen component. We can also write this as

$$p_{\text{conn}}(i, t; j, t') = \frac{\langle \text{tr} \, \Phi(t) \Phi(t') \rangle}{N(Q-1)}. \quad (\text{A28})$$

In Appendix A4 we present results for the connectedness correlation functions of boundary points. Since the boundary conditions on the Potts spins are free, this corresponds to the “ordinary” surface transition (discussed for percolation in Refs. [161–163]) where the boundary spin operator is $\propto \partial_t \Phi$ in the continuum theory rather than $\Phi$ as in the previous equation [82]. This gives the scaling forms in Appendix A4.

Finally let us consider the percolation probability $P_{\text{perc}}$. This can be used to define a characteristic timescale $t_*(r, N)$ for the classical problem, and it is much simpler to formulate in field theory than the minimal cut cost. As discussed in Sec. III E, all this carries over to finite spatial dimensions $d > 5$ by setting $N = L^d$. $P_{\text{perc}}$ is equal to $1 - e^{-\Delta F}$, where $\Delta F$ is the free-energy cost of imposing twisted boundary conditions [164] on the Potts spins [84].

The scaling form is

$$P_{\text{perc}} = F(t/N^{1/5}, \delta r t^2), \quad (\text{A29})$$

where $t$ is now the total time, and we use the notation $r$ for the parameter driving the transition to match the circuit. Let us consider a few different regimes.

By a rescaling of the field and the time coordinate we can choose to write the action in the form (suppressing order-1 constants)

$$S_{\text{eff}} = \frac{N}{t^5} \text{tr} \int_0^1 du [(\partial_u \tilde{\Phi})^2 + (\delta r t^2) \tilde{\Phi}^2 + \tilde{\Phi}^3]. \quad (\text{A30})$$

This rewriting suggests that if we take the limit of large $N$ and $t$ (and small $\delta r$) in such a way that the scaling variable $\delta r t^2$ is fixed while $N/t^5$ becomes large, $\Delta F$ is given by a saddle-point action (we do not try to make this precise in this replica theory),

$$\ln(1 - P_{\text{perc}}) \sim -\frac{N}{t^5} A (\delta r t^2), \quad (\text{A31})$$

In particular, at the critical point

$$\ln(1 - P_{\text{perc}}) \propto -\frac{N}{t^5} \quad \text{for } r = r_c \text{ and } t \ll N^{1/5}. \quad (\text{A32})$$

Note that for a system with finite $d > 5$, unlike the case $d < 5$, $P_{\text{perc}}$ is parametrically close to 1 at the critical point

of a system with $t \sim L$, i.e., with $t \sim N^{1/d}$. This is because, for percolation above the upper critical dimension, there are many percolating clusters in a large hypercubic sample at $p_c$ [165–170].

Next, let us take $N$ large with $\delta r < 0$ small but fixed, in order to examine the exponential growth of the timescale with $N$ inside the percolating phase. By an alternative rescaling of the field,

$$S_{\text{eff}} = N \text{tr} \int_0^t dt' [(\partial_{t'} \hat{\Phi})^2 + (\delta r t^2) \hat{\Phi}^2 + \hat{\Phi}^3]. \quad (\text{A33})$$

Assuming again that we can make an analogy with saddle-point solutions in more conventional theories with discrete symmetry, we anticipate a localized domain wall or “instanton” solution interpolating between the two boundary condition values of the spin, with a classical action $N \times c(r)$. The scaling form will require $c(r) \sim (\delta r)^{5/2}$ close to the critical point. At sufficiently early times there is at most one such instanton, which can be placed at any time in between 0 and $t$:

$$1 - P_{\text{perc}} \sim t |\delta r|^{1/2} \exp[-\text{const}(\delta r)^{5/2} N]. \quad (\text{A34})$$

The $\delta r$ dependence of the prefactor has been fixed by requiring consistency with the scaling form [171]. Therefore, the plateau at $P_{\text{perc}} \simeq 1$ lasts for an exponentially long time

$$t_* \sim \frac{1}{\sqrt{|\delta r|}} \exp[\text{const}(\delta r)^{5/2} N]. \quad (\text{A35})$$

The interpretation is just that the probability of having a disconnection event at a given time is $p_{\text{break}} \sim |\delta r|^{1/2} e^{-\text{const}(\delta r)^{5/2}}$. Since the probabilities of such events are independent, at long times we have exponential decay of $P_{\text{perc}}$, with a timescale also given by Eq. (A35).

### 3. Criticality in layered Erdős-Rényi graphs

As introduced in Sec. III B, the layered Erdős-Rényi model is a simplification of the classical random graph depicted in Fig. 4(b), in which a large number $N$ of nodes are arranged in discrete layers with time index $t$. Each node may be connected only to nodes in adjacent time layers; there are no connections between nodes within the same time layer. Edges between time layers $t$ and $t+1$ are randomly chosen, such that a total number of edges $cN$ are created between adjacent layers. The connectivity $c$ is the major parameter of the model ($c$ is equal to $b'/2$, in the notation of Sec. III B), and plays a similar role as the complement of the measurement rate, $1-r$. The critical value of $c$ is $c_{\text{crit}} = 1/2$, since a given node at time $t$ has connections to both $t-1$ and $t+1$ and its expected number of total connections is $2c$.
<figure><img src="image_35.png" /><figcaption>FIG. 32. The percolation probability $P_{\text{perc}}$ for the layered Erdős-Rényi graph at the critical connectivity $b'_{\text{crit}} = 1$ is plotted as a function of time. The inset shows the raw data for different system sizes $N$, and the main figure shows the same data plotted as a function of $(t+c_0)/N^{1/5}$, where $c_0 = 6$ is a constant. Compare Fig. 5 of the main text.</figcaption></figure>

Since the layered Erdős-Rényi description is the basis for the theoretical derivation of the scaling forms in Sec. III B, we numerically simulate the layered Erdős-Rényi graph to verify these scaling forms and ensure that they yield the same behavior as the data presented in Figs. 5 and 6 for the full classical graph.

Figure 32 shows the percolation probability at the critical point, plotted as a function of time for different system sizes $N$. A good scaling is observed as a function of the variable $t/N^{1/5}$, as suggested by Eq. (19). At a fixed value of $c$ and $N$, the percolation probability $P_{\text{perc}}$ is observed to decay exponentially with time. As shown in Fig. 33, near the critical point the scaled exponential decay time $\tau(c, N)/N^{1/5}$ is a function only of the variable $(c - c_{\text{crit}})N^{2/5}$. This is consistent with the scaling forms in Eq. (19).

## 4. Two-point correlation functions

In Sec. III C we show that the probability of percolation in the classical graph can be described in terms of the scaling variables $t/N^{1/5}$ and $(r - r_c)N^{2/5}$ [Eq. (19)]. Various correlation functions can also be understood in terms of these same scaling variables. The simplest correlation function, which we denote by $C$, is the probability that two distinct nodes in the graph belong to the same cluster. The nodes may be on a temporal boundary (either the same boundary or different ones) or in the bulk of the graph.

In the Potts language, the correlator is the spin two-point function. The bulk and boundary operators have different scaling dimensions, with the former scaling like $N^{-2/5}$ (or equivalently like $t^{-2}$) and the latter like $N^{-3/5}$ (or $t^{-3}$); see Sec. III E and Appendix A 2.

As a result of this scaling, the probability $C_{\text{oppo}}$ that two nodes on opposite temporal boundaries are connected by a cluster has the form

$$C_{\text{oppo}}(T, N) = \frac{1}{N^{6/5}} F_{\text{oppo}}(T/N^{1/5}, \delta r N^{2/5}). \quad (A36)$$

In Fig. 34 we test this scaling for case $r = r_c$, in order to confirm the theoretical value for the operator’s scaling dimension. This data is for the simplified multilayer Erdős-Rényi model described in Sec. III B.

The probability $C_{\text{same}}$ for two sites on the same temporal boundary to be connected has, in addition to the scaling term, a noncritical contribution of order $1/N$, which is in fact dominant at $r_c$. This $1/N$ factor is on the order of the probability for the two sites to be connected by a “microscopic” path (for example, by a single bond). For simplicity, consider the limit $T \to \infty$, when only one of the arguments of the scaling function remains:

$$C_{\text{same}}(N) = \frac{1}{N^{6/5}} F_{\text{same}}(\delta r N^{2/5}) + \frac{A(r)}{N}. \quad (A37)$$

The noncritical term may be eliminated by a subtraction: $\tilde{C}_{\text{same}}(N) = C_{\text{same}}(N) - 2^{-1}C_{\text{same}}(N/2)$. Figure 35 demonstrates a reasonable scaling collapse for this quantity. This plot constitutes a second check that $\delta r N^{2/5}$ is the appropriate off-critical scaling variable.

If $r \lesssim r_c$ is fixed and $N \to \infty$, Eq. (A37) shows that $C_{\text{same}}(N)$ scales like $(r_c - r)^3$; this is the square of the surface order parameter (the probability that a boundary site lies in the infinite cluster), which is parametrically

<figure><img src="image_36.png" /><figcaption>FIG. 33. The timescale $\tau$ for exponential decay of the percolation probability is plotted as a function of the connectivity $b'$ and system size $N$ for the layered Erdős-Rényi graph (inset). This timescale exhibits critical scaling when plotted as a function of the scaled variables $\tau/N^{1/5}$ and $(b' - b'_{\text{crit}})N^{2/5}$, as suggested by Eq. (19) with critical connectivity $b'_{\text{crit}} = 1$. Compare Fig. 6 of the main text.</figcaption></figure>
smaller than the bulk order parameter [which scales as in Eq. (18)] when $r_c$ is small. In the language of surface critical phenomena, this is the “ordinary” transition [82,161–163].

As an aside, let us make a distinction between the correlation function $C_{\text{same}}$ above and the mutual information between spins in the final state. $C_{\text{same}}$ indicates whether in the final state two spins lie in the same connected tensor network. However, being in the same connected component does not imply that spins’ mutual information, which can be detected with appropriate physical two-point functions, is large [6]. The zeroth Rényi mutual information $I_0$ is given by a different classical correlation function to the one above, which in the finite-dimensional problem behaves as a power law at $r_c$ [6]. A related observable is the distribution of entanglement entropy $S_0$ for a single spin [17]. These observables again map to boundary correlation functions of the Potts spins [172], but in a Potts system with a magnetized boundary condition, rather than with the free boundary conditions used above for $C_{\text{same}}$. (In the classification of surface criticality this is the “extraordinary” transition [173,174].) The critical contribution to these quantities at $r_c$ is smaller than a trivial analytic contribution similar to that mentioned above, we have not been able to see it numerically.

### 5. Extrapolating min-cut tension to $N \to \infty$

As mentioned in Sec. III D, the behavior of the classical minimum cut value $S_0(t, N)$ within the percolating phase,

<figure><img src="image_37.png" /><figcaption>FIG. 34. The correlation function $C_{\text{oppo}}(t)$, which denotes the probability of two randomly chosen nodes on opposite temporal boundaries being connected to the same cluster after a total time evolution $t$. Data here corresponds to the layered Erdős-Rényi model at the critical point, $c = 1/2$. The inset shows $C_{\text{oppo}}(t)N^{6/5}$ as a function of $t$ for different system sizes $N$, while the main figure shows that this same data scales onto a single curve when plotted as a function of $(t + c_0)N^{-1/5}$. Here $c_0 \approx 3.1$ is a constant.</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>FIG. 35. The critical contribution $\tilde{C}_{\text{same}}$ to the probability for two sites on the same temporal boundary to be connected, in the limit of long simulation time $T$. Data corresponds to the layered Erdős-Rényi model, and is plotted against $b'$, twice the expected number of connections between a node at time $t$ and all other nodes at time $t + 1$. The inset shows the raw data for different system size $N$, while the main figure shows the data plotted in terms of scaled variables.</figcaption></figure>

$r < r_c$, has the functional form

$$ \frac{S_0(t, N, r)}{N} \approx s(r) - \frac{d(t, r)}{\sqrt{N}}, \qquad (A38) $$

at large $N$, and for all times $t$ that are larger than an initial transient but short enough to satisfy $\ln t \ll \ln \tau$. Here, $d(t, r)$ is a constant for fixed $t$ and $r$ and $\tau$ denotes the decay time of the percolation probability; $\tau$ grows exponentially with $N$ (see Fig. 6). An example of this scaling is shown in Fig. 36. One can estimate the value of the minimum cut per spin, $s(r)$, by extrapolating this relation to $1/\sqrt{N} = 0$. As illustrated in Fig. 36, this extrapolated value is relatively insensitive to the time $t$, so long as $t$ is larger than a short-time transient ($t \gtrsim 5$ is sufficient for the $r$ values plotted here and in Fig. 7) and shorter than $\tau$.

The data in Fig. 7 corresponds to $t = 10$, and comes from an extrapolation using system sizes $N = 100, 200, 400, 600$, and $800$. As shown in Fig. 36, performing an extrapolation at $t = 20$ yields essentially identical results. The extrapolation procedure becomes numerically difficult at very small $r_c - r$, since the decay time $\tau$ becomes short for all but very large system sizes. The extrapolated data in Fig. 7 is therefore limited to $r_c - r \ge 0.02$.
<figure><img src="image_39.png" /><figcaption>FIG. 36. An example of the extrapolation of $S(t, N)/N$ to the limit of $N \to \infty$ in order to estimate the classical minimal cut per spin, $s(r)$. Points correspond to $S_0/N$ for different simulated system sizes at two different times $t$, and the dashed lines are linear best fit lines to $S_0/N$ as a function of $1/\sqrt{N}$. The $y$ intercept of these lines gives the estimated value of $s(r)$, plotted in Fig. 7. The extrapolated values of $s(r)$ for two values of $t$ agree to within 3%. The top plot corresponds to $r = 0.4$, and the bottom plot corresponds to $r = 0.7$.</figcaption></figure>

## APPENDIX B: MAJORANAS WITH PAIRWISE MEASUREMENT

Reference [37] describes measurement-only dynamics for an even number of Majorana modes $\gamma_i$ for $i = 1, \dots, N \gg 1$. Each measurement is of a fermion parity $i\gamma_i\gamma_j$ for some $i$ and some $j$. (If desired we can impose a bipartite structure, so that measurements are allowed only for $i \in A$ and $j \in B$, where $A$ and $B$ are say two sublattices of a bipartite lattice.) In finite dimensions, such a model allows a phase transition between an area-law phase and a phase with a logarithmic violation of the area law. We briefly note the fate of this transition in the all-to-all setting.

We first select a single preferred grouping of the Majoranas: for $i$ odd, $i$ is grouped with $i+1$. We can think of the two Majoranas within a group as forming a single complex fermion operator. With probability $r$, a measurement is an intragroup measurement, and with probability $1-r$ it involves two Majoranas chosen uniformly at random, so that it is an intergroup measurement with probability $1$ in the limit $N \to \infty$.

The trivial phase arises at large $r$, when intragroup measurements predominate, and the nontrivial phase at small $r$ when intergroup measurements predominate. The dynamics can be thought of as the dynamics of an evolving pairing between Majorana modes [37]. In space time it maps to a loop model [175,176]. The loops represent Majorana worldlines, and the operator entanglement is proportional

to the number of worldlines connecting the initial to the final time. The relevant case for the all-to-all model is a quasi-1D version of such a loop model [157].

The phase transition is described by the $RP^{n-1}$ $\sigma$ model in the limit $n \to 1$ (or $CP^{n-1}$ if we impose a bipartite structure). The relation between the all-to-all case and the finite-dimensional case is similar to that in Sec. III E. For the critical point it is natural to use a soft-spin formulation of the $RP^{n-1}$ or $CP^{n-1}$ model, with a cubic term and upper critical dimension 6 [176,177], so that by the logic of Sec. III E we expect the scaling variables in Eq. (19) to apply. However, the associated scaling functions, and some of the other exponents, will be very different. The basic difference from the field theory discussed in Sec. III E is that the $RP^{n-1}$ and $CP^{n-1}$ models have continuous replicalike symmetries [respectively, $SU(n)$ and $SO(n)$], unlike the discrete symmetries in both the field theory for percolation (Sec. III E) and the replica descriptions of the generic quantum problem (Sec. VI). This difference leads to much smaller entanglement in the nontrivial phase.

The entanglement in the present model is related to the free-energy cost of twisted boundary conditions in the $\sigma$ model. The possibility of continuous twisting of the order parameter makes this smaller than in a system with discrete symmetry. In finite dimensions, the nontrivial phase has only a logarithmic violation of the area law, rather than volume-law entanglement [37].

In the all-to-all setting we again characterize the phases by the operator entanglement between initial and final times. In the nontrivial phase the $RP^{n-1}$ field is “ordered” (on timescales much shorter than $\tau$ below) and so we use a nonlinear $\sigma$ model formulation for an $n \times n$ matrix $Q_{ab}$ that parameterizes $RP^{n-1}$. This has the form $S = K \int dt \operatorname{tr}(\partial_i Q)^2$, with a stiffness that scales as

$$K \propto (\delta r)^2 N \qquad (B1)$$

close to the transition on the nontrivial side. A rescaling of $t$ shows that this leads to a characteristic timescale of order

$$\tau \sim (\delta r)^2 N. \qquad (B2)$$

For $t \gg \tau$, the simple one-dimensional field theory has exponentially decaying correlations, and the operator entanglement (computed from the cost of imposing twisted boundary conditions between the initial and final time [178]) decays exponentially. For $t \ll \tau$, the operator entanglement is extensive in $N$ and scales as $S \sim \delta r^2 N/t$.

Note the large difference between the timescale in Eq. (B2), which is linear in $N$, and the timescale in the entangled phase (which is accessed in more generic dynamics) that is exponentially large in $N$. The model outlined in this appendix has only pairwise Majorana correlations, and a very restricted entanglement structure. However, this particular feature of the present problem is
likely to carry over to a larger set of models involving unitaries and measurements or projections for free fermions [53,156], since the key feature is having continuous, rather than discrete, replica symmetry [156,159].

# APPENDIX C: CALCULATIONS FOR QUANTUM TREE

## 1. Numerical recursion for quantum tree

In this section we briefly describe the numerical procedure used to simulate the quantum tree. Due to the single-site Haar rotations the Schmidt basis of the tree Eq. (56) becomes uniformly distributed and we can thus characterize the wave function using only the two Schmidt values (that is, using a single real positive number between 0 and 1/2, which is the minimal Schmidt value squared $Z = \lambda_{min}^2$).

The recursive procedure to generate the tree at generation $k+1$ consists of using three singular values at generation $k$, $Z_k$, $Z_k'$, and $Z_k''$ and connecting them to a node using Eq. (57). Thus, the number of eigenvalues required to describe a certain instance of the tree *exactly* grows exponentially as $3^k$, which is clearly not simulable at large $k$. Here we take advantage of the fact that in the case of the FMPT the nodes of the tree are statistically independent. Thus, at a certain level $k$ we can generate a large constant pool of $N$ singular values, where $1 \ll N \ll 3^k$. Assuming the pool spans the distribution function of $Z_k$ faithfully we can then draw randomly three singular values from this pool to generate a member in the pool of the next generation. This is known as the “pool method” [105–107].

To verify that the pool spans the distribution of $Z_k$ faithfully we test the convergence of the evolution of $Z_k$ with the generation number $k$ as a function of $N$. It is known that convergence in $N$ can be very slow for the pool method [105,106]. For example, in Fig. 37 we present $Z_k^{typ}$ at $\ln(r_c - r) = -5.3$ as a function of $k$ for different values of $N$, which is the point closest to the critical point in Fig. 14. The origin of the strong $N$ dependence lies in the exponent $\lambda$, which approaches 1/2 at the transition causing the distribution of $Z$ to become broad. Upon tuning farther away from the critical point the minimal $k$ and $N$ required for convergence is found to decay rapidly (not shown in the figure).

Finally we also note that due to the forced measurements the distribution function of the singular values has a $\delta$ function at $Z = 0$ with a known prefactor. Namely, the probability of a singular value at generation $k$ to be exactly zero is given by the recursive relation

$$f_k(p) = p + (1-p)[f_{k-1}(p)]^3. \quad (C1)$$

In principle we can keep these zeros in our pool. However, this is highly inefficient, especially when $p$ starts

<figure><img src="image_40.png" /><figcaption>FIG. 37. $Z_k^{typ}$ vs $k$ at $\ln(r_c - r) = -5.3$ for Haar ensemble and for various pool sizes $N$. This plot demonstrates the convergence with pool size $N$ and $k$ for the point closest to the phase transition in Fig. 14. The minimal $N$ and $k$ required for convergence diminish for points farther away from the critical point. In the case of half-integer powers we round $N$ to the closest integer.</figcaption></figure>

to get close to the classical transition, where $f_\infty(p) = [\sqrt{1 + (2 - 3p)p/(1 - p) - 1]/2$ becomes unity. Thus, in our simulations we keep only nonzero eigenvalues and, if needed, account for the zeros using Eq. (C1).

## 2. Averages of tree recursion constants

In Sec. E we describe the recursion relation for the singular values of the tree. The linearized recursion relation involves the random multiplicative constants in Eq. (E),

<figure><img src="image_41.png" /><figcaption>FIG. 38. The distribution function of the parameters $V_j \equiv \ln A_j$ for a Haar-random two-site unitary (left) and the ensemble in Eqs. (39) and (38) with $\Delta t = 0.3$. The data here is collected from $N = 10^7$ random unitaries.</figcaption></figure>
which we repeat for convenience:

$$
\begin{align*}
A_1 &= \frac{|U_{11}^{11} U_{21}^{21} - U_{21}^{11} U_{11}^{21}|^2}{(|U_{11}^{11}|^2 + |U_{11}^{21}|^2)^2}, \\
A_2 &= \frac{|U_{11}^{11} U_{12}^{21} - U_{12}^{11} U_{11}^{21}|^2}{(|U_{11}^{11}|^2 + |U_{11}^{21}|^2)^2}, \tag{C2} \\
A_3 &= \frac{|U_{11}^{11} U_{11}^{22} - U_{11}^{12} U_{11}^{21}|^2}{(|U_{11}^{11}|^2 + |U_{11}^{21}|^2)^2}.
\end{align*}
$$

In this appendix we derive the facts about the distribution of these quantities that are given in Sec. E. Some of these facts also hold for general node tensors, not necessarily expressed in terms of unitaries.

As in the text, we assume that the distribution of $U$ is invariant under multiplication by single-site U(2) matrices on any of its four legs. Initially however we do not assume that it is invariant under 2-site unitaries, i.e., we do not assume that $U$ is Haar-distributed in $U(4)$.

First we show that, so long as the unitary $U$ is non-trivially entangling (defined below) with probability 1, the average of any of the above quantities is exactly equal to 1,

$$ \langle A_i \rangle = 1, \qquad (C3) $$

and also

$$ \langle A_i^{1/2} \ln A_i \rangle = 0. \qquad (C4) $$

The argument is the same for any of the three $A_i$, so consider $A_3$ for definiteness. The argument only relies on the property of U(2) invariance on a leg mentioned in Sec. B (together with the assumption that certain singular values are not fine tuned to zero), so they hold for more general choices of $t$ satisfying this requirement. (We could also consider unitaries acting on more sites and trees with a larger branching number.)

The expression for $A_3$ involves only the matrix elements $U_{11}^{ad}$ for $a, d = 1, 2$. Regarding this as a $2 \times 2$ matrix with row index $a$ and column index $d$, we make a singular value decomposition, with positive singular values $\eta_1$ and $\eta_2$:

$$ U_{11}^{ad} = \sum_{\mu=1,2} w_{a\mu} \eta_{\mu} v_{\mu d}. \qquad (C5) $$

Here, $w$ and $v$ are U(2) matrices. Now we note that, for any given $U$, the singular values of the tree we are considering are invariant under unitary basis transformations for the bond at the top of the tree. This implies that $A_3$ must be invariant if $U_{11}^{ad}$ is multiplied by an arbitrary single-site unitary acting on the $a$ index [see Eq. (41)]. We choose

this unitary to be the inverse of $w$, so that $U_{11}^{ad}$ is replaced in Eq. (C2) by

$$ U_{11}^{ad} \rightarrow \eta_a v_{ad}. \qquad (C6) $$

Together with $|\det v|^2 = 1$, this gives the expression

$$ A_3 = \frac{\eta_1^2 \eta_2^2}{(\eta_1^2 |v_{11}|^2 + \eta_2^2 |v_{21}|^2)^2}. \qquad (C7) $$

For some trivial, nonentangling two-site unitaries, such as the identity or swap, one of the singular values $\eta$ is exactly zero, and $A_3$ vanishes. We assume that the distribution of $U$ is such that, with probability 1, both singular values $\eta$ are nonzero. This is our definition of “nontrivially entangling” above.

The above expression involves a single column, $v_{a1}$, of the U(2) matrix $v_{ad}$. Since we assume that the distribution of $U$ is invariant under single-site rotations, $v_{ad}$ is Haar distributed, and $v_{a1}$ is just a unit vector (with two complex or four real components) that must be averaged uniformly over the sphere $S^3$.

This can be done in a standard way by relating the average over the sphere to a Gaussian average. Let us write the four real components of the unit vector ($v_{11}, v_{21}$) as $V = (w, x, y, z)$. If $\langle \cdots \rangle_{\mu}$ is the Gaussian average with weight proportional to $e^{-\mu V^2}$, then

$$ \langle \cdots \rangle_{\mu} = 2\mu^2 \int dR R^3 e^{-\mu R^2} \langle \cdots \rangle_{|V|=R}, \qquad (C8) $$

as we see by splitting the Gaussian integral on the rhs into radial and angular parts. The latter gives the integral over a sphere of fixed radius, which is the last expression on the rhs. We are interested in $\langle \cdots \rangle_{|V|=1}$, the average over the unit sphere.

The Gaussian average of $A_3$, which is of order $|V|^{-4}$, diverges at small $R$, so we instead first consider

$$ f(\mu) = \left\langle e^{-(\eta_1^2(w^2+x^2)+\eta_2^2(y^2+z^2))} \right\rangle_{\mu} = \frac{\mu^2}{(\mu + \eta_1^2)(\mu + \eta_2^2)}. $$

Using Eq. (C8) we may alternately write $f(\mu)$ as an average over the sphere. By scaling out a factor of $R$ from the components of $V$ we obtain an average over the unit sphere, and performing the $R$ integral gives

$$ f(\mu) = \left\langle \frac{\mu^2}{\{\mu + [\eta_1^2(w^2+x^2) + \eta_2^2(y^2+z^2)]\}^2} \right\rangle_{|V|=1}. $$

Equating the expressions for $f(\mu)$, and taking the limit $\mu \to 0$, the U(2) Haar average is

$$ \left\langle \frac{1}{(\eta_1^2 |v_{11}|^2 + \eta_2^2 |v_{21}|^2)^2} \right\rangle_v = \frac{1}{\eta_1^2 \eta_2^2}, \qquad (C9) $$
where we restore the previous notation for the vector. Plugging this into Eq. (C7) gives $\langle A_3 \rangle = 1$, as stated above, regardless of the precise distribution of $\eta$.

The same argument applies for $A_1$ and $A_2$, using the appropriate singular value decomposition. Note that each of the $A_i$ involves only a subset of the components of $U$, so that it effectively reduces to a $2 \times 2$ matrix, as above for the matrix $U_{11}^{ad}$. That is, in each case two of the four legs of $U$ are set to index value “1.”

By multiplying $f(\mu)$ with $\mu^a$ and integrating over $\mu$, we find

$$ \langle A_i^\lambda \rangle = \frac{1}{2\lambda - 1} \left\langle \frac{1/H^{\lambda-1} - H^\lambda}{1-H} \right\rangle, \quad (C10) $$

with $H = \eta_1^2/\eta_2^2$. The remaining average on the rhs is over these singular values, again for the appropriate singular value decomposition of $U$.

Differentiating with respect to $\lambda$ at $\lambda = 1/2$ gives Eq. (C4), irrespective of the distribution of $H$.

While Eqs. (C3) and (C4) simplified, more general moments depend on the detailed distribution of $U$. For the two-site Haar case we may write analytical formulas for $\langle A_i^\lambda \rangle$ [given in Eqs. (64) and (65) of the main text].

First consider $A_3$, as above. This simplifies because, for a Haar-distributed unitary, $U_{11}^{ad}$ can be viewed as a normalized and uniformly random vector in a Hilbert space of dimension $2 \times 2$ (a Page-random state). We are interested in the singular values when this state is split into two equal subsystems. Writing $s_i = \eta_i^2$, this distribution is [179],

$$ P(s_1, s_2) = 3\delta(s_1 + s_2 - 1)(s_1 - s_2)^2 \\ \rightarrow P(H) = \frac{3(1-H)^2}{(1+H)^4}, $$

for $0 < H < \infty$. Applying this to Eq. (C10) gives

$$ \langle A_3^\lambda \rangle = \frac{\pi \lambda (1 - \lambda)}{\sin(\pi \lambda)}. \quad (C11) $$

This is finite for $\lambda \in (-1, 2)$.

By the right-invariance of the Haar measure, the distribution of $A_1$ is the same as that of $A_2$. We must consider the singular values for a decomposition of the matrix $U_{b1}^{ad}$, which is the upper left $2 \times 2$ block of a $4 \times 4$ Haar matrix. The distribution of singular values for such a sub-block of a Haar unitary may be found in Ref. [180]. Writing again $s_i = \eta_i^2$,

$$ P(s_1, s_2) = 6(s_1 - s_2)^2, \quad (C12) $$

with the constraint $0 < s_i < 1$ but, unlike in the previous case, no constraint on $s_1 + s_2$. Since there is a relabeling symmetry under $s_1 \leftrightarrow s_2$, or equivalently under $H \leftrightarrow$

$1/H$, we may insist $0 < H < 1$. Then

$$ P(H) = 3(1 - H)^2. \quad (C13) $$

Applying this to Eq. (C10) gives

$$ \langle A_1^\lambda \rangle = \langle A_2^\lambda \rangle = \frac{12}{12 + 8\lambda - 7\lambda^2 - 2\lambda^3 + \lambda^4}, \quad (C14) $$

as stated in the text. Again this is finite for $-1 < \lambda < 2$.

Finally, let us discuss the asymptotics of the distributions of the $A_i$. This will clarify the following point. In the main text we describe an ensemble of two-site unitaries with an entangling strength parameterized by $\Delta t$. In the limit of small $\Delta t$, these unitaries become closer and closer to the identity. For the identity, $A_2 = A_3 = 0$ exactly. But we show above that, for any nonzero value of $\Delta t$, no matter how small, $\langle A_2 \rangle = \langle A_3 \rangle = 1$. Therefore, the limit $\Delta t \to 0$ does not commute with the average over unitaries. This is because the distribution of $A_i$ develops a long tail when $\Delta t$ becomes small (see for example the intermediate tail in Fig. 38 for $\Delta t = 0.3$).

Define the two-component vectors

$$ \phi_a = U_{11}^{a1} \quad \psi_a = U_{21}^{a1}, \quad \chi_a' = U_{11}^{a2} \quad \psi_a' = U_{12}^{a1}. \quad (C15) $$

Then we may write Eqs. (C2) as

$$ \begin{aligned} A_1 &= \frac{|\psi|^2 |\phi|^2 - |\psi^\dagger \phi|^2}{|\phi|^4}, \\ A_2 &= \frac{|\psi'|^2 |\phi|^2 - |\psi'^\dagger \phi|^2}{|\phi|^4}, \\ A_3 &= \frac{|\chi|^2 |\phi|^2 - |\chi^\dagger \phi|^2}{|\phi|^4}. \end{aligned} \quad (C16) $$

We see that $A_i$ can become arbitrarily large if $|\phi|$ becomes small, with $A_i$ scaling like $|\phi|^{-2}$ in this limit. Since $\phi$ has two complex (or four real) components, we expect that for a generic distribution of unitaries, the cumulative probability distribution of $|\phi|$ scales like $|\phi|^4$ at small $|\phi|$. This gives

$$ P_{\ln}(\ln A_i) d \ln A_i \sim \frac{d \ln A_i}{A_i^2}, \quad A_i \gg 1 \quad (C17) $$

at large $A_i$, as stated in the text. At small $A$, similar considerations for the numerators in Eq. (C2) show that generically

$$ P(\ln A_i) d \ln A_i \sim A_i d \ln A_i, \quad A_i \ll 1. \quad (C18) $$

These power laws are consistent with numerics and also with the fact that the moments $\langle A_i^\lambda \rangle$ for the Haar case diverge at $\lambda = 2$ and at $\lambda = -1$. The $A_i$ are of course
correlated, but we do not consider their joint distribution here.

Now consider the case of a weakly entangling unitary with random single-site scramblers,

$$U = (u_1 \otimes u_2)e^{-i\Delta tH}(u_3 \otimes u_4), \quad (C19)$$

with $\Delta t$ small but fixed. (H may either be fixed, as in an ensemble discussed in the main text, or random.) We focus on $A_2$ and $A_3$, whose distributions become broad at small $\Delta t$ (that of $A_1$ does not). The two cases are similar, so consider $A_3$, which is given by Eq. (C7). We expect that for small $\Delta t$ we typically have $\eta_1^2 \sim \Delta t^2$ (here we keep only the scaling with $\Delta t$). Therefore, so long as $|v_{22}|^2 \gg \Delta t^2$, i.e., in the regime $A_3 \ll \Delta t^{-2}$, we have

$$A_3 = \frac{\eta_1^2 \eta_2^2}{(\eta_1^2 |v_{12}|^2 + \eta_2^2 |v_{22}|^2)^2} \sim \frac{\Delta t^2}{|v_{22}|^4}. \quad (C20)$$

For small $|v_{22}|$, the cumulative distribution of $|v_{22}|$ scales like $|v_{22}|^2$, which gives

$$P(\ln A_3) d \ln A_3 \sim \frac{\Delta t}{\sqrt{A_3}} d \ln A_3, \quad \Delta t^2 \ll A_3 \ll \Delta t^{-2}. \quad (C21)$$

On the other hand, for $A_3 \gg \Delta t^{-2}$, we expect to recover the generic exponent $-2$ [Eq. (C17)] for the distribution, suggesting

$$P(\ln A_3) d \ln A_3 \sim \frac{(\Delta t)^{-2}}{A_3^2} d \ln A_3, \quad A_3 \gg \Delta t^{-2}. \quad (C22)$$

Similarly for $A_3 \ll \Delta t^2$ we expect to recover the exponent 1 for the generic case [Eq. (C18)]. The existence of three regimes with different power law exponents, 1, -1/2, and -2, when $\Delta t$ is relatively small, is in good agreement with the numerical data as shown in Fig. 38. Note that the scaling above, with the -1/2 tail being cut off at $A_3 \sim \Delta t^{-2}$, is also consistent with $\langle A_3 \rangle$ being of order 1 at small $\Delta t$.

### 3. Continuum recursion relation at $\Delta > 1$

We start with the FKPP equation with the spatially varying diffusion constant $D(x) = 1 + e^x$,

$$\partial_{\tau} H = \partial_{x} [D(x) \partial_{x} - a] H + \Delta H (1 - H), \quad (C23)$$

in the regime $\Delta > 1$, such that the linearized problem is in the paramagnetic phase. The entanglement transition is at $a = -(\Delta + 1)$ so we write

$$a = -(\Delta + 1) + \sigma, \quad (C24)$$

with $0 < \sigma \ll 1$. When $\sigma > 0$, we converge at late times to a stationary solution satisfying $\partial_{\tau} H = 0$. We would like

to determine the position $x_f$ of the front in this solution, or equivalently the value of $Z^{\text{typ}}$ (recall that $x_f \sim \ln Z^{\text{typ}}$).

Since H varies by an exponentially large factor over the relevant range of $x$, it is useful to look instead at the local exponential decay rate, which is order 1:

$$R(x) \equiv \partial_x \ln H(x). \quad (C25)$$

From the stationary version of Eq. (C23), this satisfies

$$\partial_x R = -\Delta f(x) - [1 - (\Delta - \sigma)f(x)]R - R^2 + \Delta f(x)H(x), \quad (C26)$$

where

$$f(x) = \frac{1}{e^x + 1}. \quad (C27)$$

$R$ tends to 1 in the limit $x \to \infty$. [To see this note that in this limit Eq. (C23) becomes $\partial_x e^x \partial_x H = 0$. Together with $\lim_{x \to \infty} H = 0$ this gives $H \propto e^{-x}$.] Let us define

$$R(x) = -1 + \sigma S(x). \quad (C28)$$

Then

$$\partial_x S = -f(x) + [1 - (\Delta - \sigma)f(x)]S - \sigma S^2 + \frac{\Delta}{\sigma}f(x)H(x). \quad (C29)$$

We will see in a moment that we must treat the cases $1 < \Delta < 2$ and $\Delta > 2$ separately. For now let us just note that if $x$ is sufficiently large, the final term above will be subleading since both $H$ and $f$ tend to zero at large $x$. Assuming we are in the range of $x$ where this term is negligible, and dropping terms that are subleading in $\sigma$,

$$\partial_x S \approx -f(x) + [1 - \Delta f(x)]S. \quad (C30)$$

Solving this equation, and fixing the integration constant by demanding that $S$ does not blow up as $x \to \infty$,

$$S = \frac{1}{\Delta(\Delta - 1)} \left[ e^{-(\Delta - 1)x} (e^x + 1)^\Delta - (e^x + \Delta) \right]. \quad (C31)$$

Now consider this solution for large negative $x$

$$S = \frac{1}{\Delta(\Delta - 1)} \left[ e^{(\Delta - 1)|x|} + \dots \right]. \quad (C32)$$

We check below that for $1 < \Delta < 2$, there is a range of $x$ where this expansion is valid, i.e., where the final term in Eq. (C29) can indeed be neglected.
According to this expansion, as we increase $-x$, the value of $R$ begins to increase significantly from $-1$ once

$$x \sim -\frac{1}{\Delta - 1} \ln \frac{1}{\sigma}. \qquad (C33)$$

This suggests that the front is at $x_f \sim -\frac{1}{\Delta-1} \ln \frac{1}{\sigma}$, and that here we can match onto the stationary solution of the traveling-wave equation, Eq. (C23), with $\sigma = 0$ and with a spatially constant diffusion coefficient. (This matching makes sense since the forward part of this solution has $R = -1$ when $\sigma = 0$.)

However, we must check the self-consistency of our neglect of the final term in Eq. (C29). Assuming the above scaling for $x_f$, we have

$$\frac{1}{\sigma} f(x) H(x) \sim \sigma^{(2-\Delta)/(\Delta-1)} \frac{e^{-2x}}{1+e^{-x}} \qquad (C34)$$

for $x \gg x_f$. If $1 < \Delta < 2$, this term is indeed much smaller than the rhs of Eq. (C30) for $x \gg x_f$. Therefore, for this range of $\Delta$ the above analysis, giving

$$x_f \sim -\frac{1}{\Delta - 1} \ln \frac{1}{\sigma} \quad (1 < \Delta < 2), \qquad (C35)$$

is self-consistent, though not rigorous.

On the other hand, for $\Delta > 2$, when the power of $\sigma$ in Eq. (C34) is negative, this term cannot be dropped from Eq. (C29) for $x \gg x_f$. Numerically solving Eq. (C29) suggests that instead the final term in Eq. (C29) contributes at leading order if we fix $x$ and take $\sigma \to 0$. Indeed, the alternative would be to have some $x_c$, with $x_f \ll x_c \ll 0$, such that the term is negligible for $x \gg x_c$ but not for $x_f < x < x_c$, and this may be seen to be inconsistent by examining the ratio of this term to the right-hand side of Eq. (C30). Using this fact, that $H/\sigma$ should be of order 1 when $x \sim 0$, and assuming that $H \sim e^{-(x-x_f)}$ for $x \gg x_f$, we find that

$$x_f \sim -\ln \frac{1}{\sigma} \quad (\Delta > 2). \qquad (C36)$$

Equations (C35) and (C36) give the power laws $Z^{\text{typ}} \sim \sigma^{1/(\Delta-1)}$ and $Z^{\text{typ}} \sim \sigma$ stated in Sec. I.

These power laws can be checked directly by making a numerical solution of Eq. (C23) for $H(x, \tau)$ at different values of $\Delta$ and $a$. We use a numeric differential equation solver, solved over a wide domain of discrete values $x \in (x_L, x_R)$, with boundary conditions such that $H(x_L) = 1$ and $H(x_R) = 0$. These are exponentially close in $x_L, x_R$, respectively, to the true values for the solution on the infinite domain. An initial guess is used for $H(x, 0)$ and then evolved until a very long time $\tau = \tau_f$ in order to arrive at the steady-state solution, from which we can read off the position of the front. We define this as the value of

$x = x_{\text{front}}$ such that $H(x_{\text{front}}, \tau_f) = 1/2$. For the data presented in Fig. 17, $x_L = -70$, $x_R = 40$, $\tau_f = 10^8$, and the domain of $x$ is discretized into 8001 points.

Making a linear fit of $x_{\text{front}}$ against $\ln \sigma$ for a given value of $\Delta$, and recalling $x_{\text{front}} \sim \ln Z^{\text{typ}}$, allows one to determine the value of the exponent $\gamma$, defined by $Z^{\text{typ}} \propto \sigma^\gamma$. We make this linear fit over the range of $\sigma$ such that $-10 < \ln \sigma < -6$. Smaller $\sigma$ requires very high numerical accuracy (a dense discretization of the domain of $x$), while at larger $\sigma$ the critical behavior may not be apparent. The results are shown in the main text in Fig. 17. The error bars in this figure are defined by the difference in slope obtained from fits using only the left half of this range, $-10 < \ln \sigma < -8$, as compared to only the right half of the range, $-8 < \ln \sigma < -6$.

### 4. Minimal cut formula on tree

Assume that the minimal cut chops out $m$ subtrees from the full tree, with each subtree being cut only once, at its apex. The singular values of a given subtree $a \in \{1, \dots, m\}$ are $\{\lambda_1^{(a)}, \lambda_2^{(a)}\}$, with $[\lambda_1^{(a)}]^2 + [\lambda_2^{(a)}]^2 = 1$ and $[\lambda_2^{(a)}]^2 = Z^{(a)}$. The full state may be written

$$|\psi\rangle = \sum_{i_1, \dots, i_m} \lambda_{i_1}^{(1)} \cdots \lambda_{i_m}^{(m)} |i_1, \dots, i_m\rangle_{\text{subtrees}} |i_1, \dots, i_m\rangle_{\text{rest}}. \qquad (C37)$$

The states $|i_1, \dots, i_m\rangle_{\text{subtrees}}$ are products of the Schmidt states at the base of the subtrees, as in Sec. E, and are orthonormal. The states $|i_1, \dots, i_m\rangle_{\text{rest}}$, for the remaining spins, are neither normalized nor orthogonal. They are obtained by contracting the tensor network on the other side of the min cut with Schmidt states at the tops of the subtrees. The Rényi entropies are determined by the singular values of the (unnormalized) density matrix

$$\rho = \sum_{i_1, \dots, i_m} [\lambda_{i_1}^{(1)}]^2 \cdots [\lambda_{i_m}^{(m)}]^2 |i_1, \dots, i_m\rangle_{\text{rest}} \langle i_1, \dots, i_m|_{\text{rest}}. \qquad (C38)$$

In the limit of small $Z^{(a)}$, with everything else fixed,

$$\rho \approx |e_0\rangle\langle e_0| + \sum_{a=1}^{m} Z^{(a)} (|e_a\rangle\langle e_a| - |e_0\rangle\langle e_0|), \qquad (C39)$$

with $e_0 = |1, \dots, 1\rangle_{\text{rest}}$ and $e_a = |1, \dots, 2, \dots, 1\rangle_{\text{rest}}$, with the “2” in slot $a$. For $n > 1$ this gives

$$S_n \approx \frac{n}{n-1} \sum_{a=1}^{m} \left( \frac{|e_a|^2 |e_0|^2 - |e_0^\dagger e_a|^2}{|e_0|^4} \right) Z^{(a)} + \mathcal{O}(Z^2). \qquad (C40)$$

Each term is associated with a bond lying on the minimal cut. Let the height of this bond above the base be $k(a)$.
For a given term, the coefficient in brackets will vanish if the states $e_a$ and $e_0$ become parallel. Each of these is a state in the “rest” Hilbert space, given by the tensor network made up of the tensors on one side of the minimal cut. The bonds lying on the minimal cut have fixed states attached to them: one of these is changed in going from $e_0$ to $e_1$. Exploiting the fact that these truncated tensor networks are still trees, we can write $e_0$ and $e_1$ in terms of the singular value decompositions of two subtrees. We see that the term in brackets is of the same order as the singular value for a tree of depth $k(a)$. (This assumes that the singular values are not growing with $k$: this can occur if we fine tune the boundary conditions, but is not relevant to the case we are discussing.) Therefore, in the present limit of small $Z$ we confirm the minimal cut conjecture in the main text, according to which each bond $a$ on the minimal cut contributes an entanglement of order $Z^{(a)}Z'^{(a)}$, where $Z^{(a)}$ and $Z'^{(a)}$ are two random variables each distributed like $Z_k$ for $k = k(a)$. (The $Z^{(a)}$ values are independent of each other and of the $Z'^{(a)}$ values, but the $Z'^{(a)}$ are correlated among themselves.)

## APPENDIX D: MORE ON CIRCUIT SIMULATIONS

We present some additional details regarding the simulations of quantum circuits in Sec. V. The nonunitary time-evolution operator, $V$, for a given realization of the circuit is built by choosing at every step a measurement with probability $r$ on a randomly chosen site or entangling two randomly chosen sites with the unitary with probability $1-r$. For a system of size $N$, time progresses by one unit for every $N$ unitaries applied. The Haar-random unitaries are generated using Mezzadri’s algorithm [181].

Note that the $\tau_{\text{eff}}$ described in the main text is obtained by taking a log derivative of $D(t)$. Since, numerical derivatives are notoriously noisy, we smooth the data for $D(t)$ using a Savitzky-Golay filter [182] and then take the derivative. Representative examples of $\tau_{\text{eff}}$ and the plateaux therein are shown in Figs. 39 and 40.

From the $\tau$ so obtained, we extract $a(r)$ by fitting the data to a form $a(r)N + b(r) + c(r)/N$. The fits are shown in Fig. 41.

## APPENDIX E: TOY MODEL FOR THE SLOW DECAY OF $S$

In the entangled phase the operator entanglement has a plateau at an extensive value $S_n/N \approx s_n > 0$, which persists for a time that scales exponentially in the number of spins, $N$. In this appendix we describe the crudest toy model for this plateau, which neglects both locality of the interactions and, for the most part, distinctions between different Rényi entropies. We consider the

<figure><img src="image_42.png" /><figcaption>FIG. 39. $\tau_{\text{eff}}$ from time derivatives of $\ln D^{\text{typ}}$. The plateaux denoted by the dashed lines show the $\tau(r, N)$ values, which are used to extract $a(r)$ (see Fig. 41). Different colors correspond to the different $N$ following the same convention as in Fig. 22. Results are for the Haar circuit, with the top and bottom panels corresponding to the forced-measurement protocol and measurement protocol respectively.</figcaption></figure>

forced-measurement case, where the circuit is made up of uncorrelated random pieces.

Imagine dividing up the circuit up into blocks of temporal duration $\Delta t$, which corresponds to writing the time-evolution operator $V(t)$ as a product of random matrices $W_i$, with each matrix of size $2^N \times 2^N$ and $i = 1, \dots, t/\Delta t$. $\Delta t$ is chosen to be much larger than 1 but much smaller than the timescale $\tau$ that will emerge below. Each block has a singular value decomposition $W_i = U_i^{(1)} D U_i^{(2)\dagger}$. As a toy model, we treat the unitaries $U^{(1)}$ and $U^{(2)}$ as Haar random (neglecting locality) and we make the simplest choice of $D$ that yields a given value of $S_1 = sN$, which is a flat

<figure><img src="image_43.png" /><figcaption>FIG. 40. Same as in Fig. 39 but for the Haar circuit with measurements.</figcaption></figure>
<figure><img src="image_44.png" /><figcaption>FIG. 41. Fits of $\ln \tau$ to a function of the form $a(r)N + b(r) + c(r)/N$ for different values of $r$ for the Haar circuit with forced measurements (top) and measurements (bottom). We use $N \ge 14$ for the fits shown, which are used to extract $a(r)$ shown in Fig. 24. Fits including smaller $N$ (not shown) are used to estimate the errorbars on $a(r)$.</figcaption></figure>

entanglement spectrum:

$$D = \frac{1}{\sqrt{B}} \begin{pmatrix} \mathbb{I}_{B \times B} & 0 \\ 0 & 0 \end{pmatrix}, \quad B = e^{sN}, \quad 0 < s < \ln 2.$$

Note that the nonzero block is a small fraction of the size of the matrix for large $N$. Finally, we make an uncontrolled simplification by also treating the entanglement spectrum of $V(t)$ as flat:

$$V(t) = U_L(t)D(t)U_R(t)^\dagger,$$

$$D(t) = \frac{1}{\sqrt{B(t)}} \begin{pmatrix} \mathbb{I}_{B(t) \times B(t)} & 0 \\ 0 & 0 \end{pmatrix},$$

with $B(t) = e^{S(t)}$. Note that $B$ and $D$ refer to quantities for a single slice of width $\Delta t$, while $B(t)$ and $D(t)$ refer to the complete evolution operator up to time $t$.

We are interested in the singular values of the new evolution matrix $V(t + \Delta t)$, which may be written

$$V(t + \Delta t) = U' D(t) U D U^{\dagger} \quad (E1)$$

for Haar unitaries $U'$, $U$, $U''$. These values are also the singular values of

$$\tilde{V} = D(t)UD. \quad (E2)$$

Up to a normalization factor, $\tilde{V}$ is just a rectangular block, of size $B(t) \times B$, taken from a Haar unitary of exponentially larger size $(2^N \times 2^N)$. Correlations between unitary matrix elements become weaker as the size of the matrix increases, so we expect that we can treat them as Gaussian, with $E U_{ab} = 0$ and $E U_{ab} (U_{a'b'}^*) = (1/2^N) \delta_{aa'} \delta_{bb'}$. (Higher cumulants are suppressed by powers of $2^N$.) The singular values squared of $\tilde{V}$, denoted $v_i = \eta_i^2$, are eigenvalues of the $B(t) \times B(t)$ matrix

$$M = \tilde{V}\tilde{V}^{\dagger}. \quad (E3)$$

When the matrix elements are of $\tilde{V}$ are complex Gaussian random numbers, this is as a Wishart random matrix (see, e.g., Ref. [179] for an application in a related context). The distribution of its eigenvalues depends on $B$ as well as on $B(t)$. We assume that $1 \ll B(t) \ll B$. Normalizing the matrix so the $s_i$ sum to one, the eigenvalue density for $\tilde{V}\tilde{V}^{\dagger}$ is the Marcenko-Pastur distribution,

$$\rho(v) = \frac{8\sqrt{v-v_{-}}\sqrt{v_{+}-v}}{\pi(v_{+}-v_{-})^{2}v}, \quad v_{\pm} = \frac{[\sqrt{B} \pm \sqrt{B(t)}]^{2}}{BB(t)}.$$

From this distribution the operator Rényi entropies can be calculated as

$$e^{-(n-1)S_n(t+\Delta t)} = \int_{v_{-}}^{v_{+}} dv v^n \rho(v) \\ \approx B(t)^{-(n-1)} \left[ 1 + \frac{n(n-1)}{2} \frac{B(t)}{B} + \dots \right] \\ S_n(t + \Delta t) - S(t) \approx - \exp \left\{ - \left[ sN - S(t) - \ln \frac{n}{2} \right] \right\}.$$

As expected, the entanglement spectrum does not remain flat. In our crude approximation, however, we neglect this, and apply the above transformation iteratively, so that in the continuum limit for times $\gg \Delta t$:

$$\partial_t S(t) \sim -C \exp\{-[sN - S(t)]\}. \quad (E4)$$

Here $C$ is an order-1 constant. At times larger than $\Delta t$, but short enough such that $S(t) \gg 1$ (which we assume above), this equation gives a solution:

$$S(t) = sN - \ln t + \dots \quad (E5)$$

independently of the value of $C$. Note that we neglect random fluctuations (see Sec. I).
This analysis suggests a characteristic timescale $\tau$ with $\ln \tau \simeq sN$. We may verify this dependence directly in the opposite limit of asymptotically late times, where (as usual for a product of random matrices) there is a separation of scale between the largest singular value, the second largest, and so on. This separation allows us to consider only the two largest singular values. At a given time $t$, let them be normalized as

$$ \{\eta_1, \eta_2\} = \{1, \epsilon(t)\}. \qquad (E6) $$

Then in place of $M$ in Eq. (E3) we have a $2 \times 2$ matrix

$$ M_{ik} = \epsilon^{i+k-2} \sum_{j=1}^{B} U_{ij} U_{kj}^* . \qquad (E7) $$

Each of the elements of $M_{ik}$ is a different sum of many random variables, so we assume that $M$ can be approximated as Gaussian:

$$ \langle M_{ik} \rangle \approx \frac{B \epsilon^{2(i-1)} \delta_{ik}}{2^N}, \quad \langle \langle M_{ik} M_{i'k'} \rangle \rangle \approx \frac{B \epsilon^{2(i+i'-2)} \delta_{ik'} \delta_{i'k}}{2^{2N}}. $$

After absorbing a normalization constant into $M$,

$$ M = \begin{pmatrix} 1 & 0 \\ 0 & \epsilon^2 \end{pmatrix} + \frac{1}{\sqrt{B}} \begin{pmatrix} a & \epsilon\beta \\ \epsilon\beta^* & \epsilon^2 b \end{pmatrix}, \qquad (E8) $$

where $a, b, \beta$ have mean zero and $\langle a^2 \rangle = \langle b^2 \rangle = |\langle \beta \rangle|^2 = 1$. The new (small) singular value squared is

$$ \epsilon_{\text{new}}^2 = \epsilon^2 \left( 1 - \frac{a-b}{\sqrt{B}} + \frac{a^2 - ab - |\beta|^2}{B} + \dots \right). \qquad (E9) $$

Let us study the typical value of this exponentially small quantity, defined by $(\epsilon_{\text{new}}^2)_{\text{typ}} = \exp\langle \ln \epsilon_{\text{new}}^2 \rangle$. We have

$$ \ln \epsilon_{\text{new}}^2 = \ln \epsilon^2 + \frac{b-a}{\sqrt{B}} + \frac{a^2 - b^2 - 2|\beta|^2}{2B} + \dots \qquad (E10) $$

Here the average is taken over $a, b, \beta$. Note that the leading fluctuation term, of order $1/\sqrt{B}$, averages to zero. The next term, however, gives a negative drift under the recursion. Recalling that $B = e^{sN}$, and applying this map iteratively with each added block,

$$ \epsilon_{\text{typ}}^2(t) \sim \exp\left(-\frac{t}{\tau}\right), \quad \tau = \Delta t \times \exp(sN). \qquad (E11) $$

Therefore, at the latest times,

$$ S_n(t) \sim \exp\left(-\frac{t}{\tau}\right) \times \begin{cases} [n/(n-1)] & n > 1 \\ (2t/\tau) & n = 1 \end{cases}. \qquad (E12) $$

Thus, the main conclusion is that in this toy model $sN$ sets both the value of the early time plateau and also the timescale for the late-time exponential decay, in agreement with the picture from the replica treatment.

## APPENDIX F: FIELD THEORY: FURTHER DETAILS

### 1. $N \to 1$ ordered phase

In Sec. H we state that for $\mu^2 < 0$ the field theory of Sec. E has an ordered phase with $X_{ab} = f(\delta_{ab} - 1/N) + W_{ab}$, where $f = (-\mu^2/3g)[N/(N-2)]$ is the order parameter, and $W_{ab}$ represents fluctuations around the saddle-point value whose quadratic Lagrangian is

$$ \mathcal{L} = \frac{1}{2} \sum_{ab} (\partial W_{ab})^2 + \frac{-\mu^2 N}{2(2-N)} \left( \sum_{ab} W_{ab}^2 - 2 \sum_a W_{aa}^2 \right). \qquad (F1) $$

We want to check that if we compute the masses of the fluctuation modes using this expression, and then take $N \to 1$, these masses remain positive. Viewing $W_{ab}$ as a vector, the term $\sum_a W_{aa}^2$ is $W.M.W$ where the matrix $M_{ab,cd}$ (with row index $a, b$ and column index $c, d$) is 1 if $a=b=c=d$ and zero otherwise. We want the eigenvalues of $M$ when projected onto the subspace of $W$ satisfying $\sum_a W_{ab} = 0$ and $\sum_b W_{ba} = 0$, i.e., the eigenvalues of $\tilde{M}_{ab,cd} = P_{aa'}P_{bb'}M_{a'b'c'd'}P_{c'c}P_{d'd}$ where $P$ is the projector $P_{aa'} = \delta_{aa'} - 1/N$. Since $M$ is nonzero only when its four indices are equal, drawing a diagram shows that $\text{tr} \tilde{M}^k = \text{tr} S^k$, where $S_{a,b} = (1-2/N)\delta_{ab} + 1/N^2$. This gives the nonzero eigenvalues of $\tilde{M}$ as $(N-1)/N \to 0$ with multiplicity 1 and $(N-2)/N \to -1$ with multiplicity $N-1$. Altogether, the eigenvalues of the matrix $(\mathbb{I} - 2\tilde{M})$, which appears in Eq. (F1) are either 1 or 3 in the limit and are positive.

### 2. Effect of YFY coupling when $N > 0$

Here we show that for $N > 0$ (for example, in the replica limit $N \to 1$, but not in the limit $N \to 0$), and in the vicinity of the critical point $r = 0$, the theory

$$ \mathcal{L} = \sum_{ab} \left[ \frac{1}{2} (\partial Y_{ab})^2 + r Y_{ab} + g Y_{ab}^3 \right] + \frac{m_F^2}{2} \sum_{abcd} Y_{ab} F_{ab,cd} Y_{cd} \qquad (F2) $$

can be reduced to the theory

$$ \mathcal{L} = \sum_{ab} \left[ \frac{1}{2} (\partial X_{ab})^2 + \frac{\mu^2}{2} X_{ab}^2 + g X_{ab}^3 \right], \qquad (F3) $$

for a matrix $X$ with vanishing row and column sums, by discarding massive modes.

First, shifting the field by $Y_{ab} \to Y_{ab} + c$ with $c = -r/(2Nm_F^2) + \mathcal{O}(r^2)$ removes the linear term and generates a
mass term. The quadratic part of the Lagrangian is then

$$ \mathcal{L}_2 = \frac{1}{2} \sum_{ab,cd} Y_{ab} \left[ (k^2 + \mu^2) \delta_{ac} \delta_{bd} + m_F^2 F_{ab,cd} \right] Y_{cd}, \quad (F4) $$

with $\mu^2 = -3gr/(Nm_F^2) + \mathcal{O}(r^2)$. As a matrix, $F = \mathbb{I} \otimes E + E \otimes \mathbb{I}$, where $E$ is the $N \times N$ matrix with unit elements: $E_{ab} = 1$, so the matrix appearing in the brackets in Eq. (F4) is

$$ (k^2 + \mu^2)\mathbb{I} \otimes \mathbb{I} + m_F^2 (\mathbb{I} \otimes E + E \otimes \mathbb{I}). \quad (F5) $$

This can be decomposed in terms of the projection matrices $P_1 \equiv N^{-1}E$ and $P_{N-1} \equiv \mathbb{I} - P_1$ as

$$ (k^2 + \mu^2) (P_{N-1} \otimes P_{N-1}) \quad (F6) $$

$$ + (k^2 + \mu^2 + m_F^2 N) (P_1 \otimes P_{N-1} + P_{N-1} \otimes P_1) \quad (F7) $$

$$ + (k^2 + \mu^2 + 2m_F^2 N) (P_1 \otimes P_1). \quad (F8) $$

This is a decomposition into three representations of $G_N$ of dimensions $(N-1)^2$, $2(N-1)$, and 1. We see that, at the critical point (where $r$ and therefore $\mu^2$ vanish) the second and the third representations remain massive and only the first becomes massless. Retaining only this representation is equivalent to fixing the row and column sums of $Y$ to zero. Doing so and renaming the resulting field $X$ yields precisely Eq. (F3). More precisely the two massive representations ought to be integrated out, renormalizing the values of the couplings in Eq. (F3).

These manipulations manifestly require $N > 0$: for example, they are appropriate for the replica limit $N \to 1$ relevant to the MPT. The replica limit $N \to 0$ must be handled separately.

### 3. Rewriting Lagrangian in $N \to 0$ limit

We give details of field redefinitions in Sec. G. First define the $N$-component vectors $\mathbf{v}^+$, $\mathbf{v}^-$, and $\mathbf{v}^i$ for $i = 2, \dots, N$ [133–136]:

$$ \mathbf{v}^+ = \frac{1}{2} (1, 0, \dots, 0) + \frac{1}{2} \frac{(0, 1, \dots, 1)}{N-1}, \quad (F9) $$

$$ \mathbf{v}^- = \frac{1}{2} (1, 0, \dots, 0) - \frac{1}{2} \frac{(0, 1, \dots, 1)}{N-1}, \quad (F10) $$

$$ \mathbf{v}^i = (0, \dots, 0, 1, 0, \dots, 0) - \frac{(0, 1, \dots, 1)}{N-1}, \quad (F11) $$

where the extra “1” in the third line is in the $i$th place. There are $N+1$ of these vectors but the $\mathbf{v}^i$ are not linearly independent: $\sum_{i>1} \mathbf{v}^i = 0$.

We use these vectors to rewrite $Y_{ab}$ in terms of $y_{\alpha\beta}$:

$$ y_{\alpha\beta} = \mathbf{v}^{\alpha} \cdot \mathbf{Y} \cdot \mathbf{v}^{\beta} = \sum_{ab} v_a^{\alpha} Y_{ab} v_b^{\beta}. \quad (F12) $$

In this appendix we use $a, b, a', b', \dots$ to denote indices that run from 1 to $N$, and $\alpha, \beta, \dots$ to denote indices that take the $N+1$ values $\{+, -, 2, \dots, N\}$. We use $i, j, k$ to denote indices that run only over 2 to $N$. Note that

$$ \sum_{j=2}^{N} y_{j, \beta} = 0, \quad \sum_{k=2}^{N} y_{\alpha, k} = 0. \quad (F13) $$

(Alternately we could define an $N \times N$ matrix without this redundancy.)

To invert the transformation defining $y$, define $N+1$-component vectors

$$ \mathbf{x}^1 = (1, 1, 0, \dots, 0), \quad (F14) $$

$$ \mathbf{x}^{i>1} = (1, -1, 0, \dots, 1, \dots, 0). \quad (F15) $$

The final 1 in the second line is in $(i+1)$st place, which in our labeling convention corresponds to the component $x_{\beta}^i$ with $\beta = i$. We have $\sum_{\alpha} x_{\alpha}^a v_{a'}^{\alpha} = \delta_{a'}^a$, so that

$$ Y_{ab} = \mathbf{x}^a \cdot \mathbf{y} \cdot \mathbf{x}^b = \sum_{\alpha\beta} x_{\alpha}^a y_{\alpha\beta} x_{\beta}^b. \quad (F16) $$

Explicitly (for $j, k > 1$):

$$ Y_{11} = (y_{++} + y_{+-} + y_{-+} + y_{--}), \quad (F17) $$

$$ Y_{1k} = (y_{++} - y_{+-} + y_{-+} - y_{--}) + (y_{+k} + y_{-k}), \quad (F18) $$

$$ Y_{j1} = (y_{++} + y_{+-} - y_{-+} - y_{--}) + (y_{j+} + y_{j-}), \quad (F19) $$

$$ Y_{jk} = (y_{++} - y_{+-} - y_{-+} + y_{--}) \quad (F20) $$

$$ + (y_{j+} - y_{j-}) + (y_{+k} - y_{-k}) + y_{jk}. \quad (F21) $$

Inserting these relations into the derivative term,

$$ \begin{aligned} \sum_{ab} (\partial Y_{ab})^2 = & 8 (\partial y_{+-} \partial y_{-+} + \partial y_{++} \partial y_{--}) \\ & + 4 \sum_j \partial y_{j+} \partial y_{j-} \\ & + 4 \sum_k \partial y_{+k} \partial y_{-k} \\ & + \sum_{jk} (\partial y_{jk})^2 \\ & + \mathcal{O}(N), \end{aligned} \quad (F22) $$
where the final line contains terms whose coefficients contain an explicit factor of N, which we assume can be neglected in the limit $N \to 0$.

Next consider the F term in the Lagrangian, which has the form

$$ \sum_{abcd} Y_{ab}F_{ab,cd}Y_{cd} = \sum_{a} \left( \sum_{b} Y_{ab} \right)^2 + \sum_{b} \left( \sum_{a} Y_{ab} \right)^2 . \qquad (F23) $$

This simplifies to

$$ \begin{aligned} \sum_{abcd} Y_{ab}F_{ab,cd}Y_{cd} = & \ 16y_{--}(y_{+-} + y_{-+}) \\ & + 4\sum_{j}y_{j-}^2 \\ & + 4\sum_{k}y_{-k}^2 \\ & + \mathcal{O}(N). \end{aligned} \qquad (F24) $$

The linear term in the Lagrangian is simply

$$ \sum_{ab} Y_{ab} = 4y_{--} + \mathcal{O}(n). \qquad (F25) $$

Since the coefficient of the linear term is zero at the critical point, we assign “engineering” dimensions to the fields such that the quadratic terms in Eqs. (F22) and (F24) are marginal. Denoting the space-time dimension by D, the inverse length dimensions of the fields are, in order of increasing scaling dimension,

$$ [y_{++}] = (D - 6)/2, \qquad (F26) $$

$$ [y_{j+}] = [y_{+k}] = (D - 4)/2, \qquad (F27) $$

$$ [y_{+-}] = [y_{-+}] = [y_{jk}] = (D - 2)/2, \qquad (F28) $$

$$ [y_{j-}] = [y_{-k}] = (D + 0)/2, \qquad (F29) $$

$$ [y_{--}] = (D + 2)/2, \qquad (F30) $$

so that if the number of “+” indices for a component of y is $n_+$, and the number of “-” indices is $n_-$, the engineering dimension is

$$ x(n_{+}, n_{-}) = \frac{(D - 2) - 2(n_{+} - n_{-})}{2}. \qquad (F31) $$

This formula implies that, at a given order in y, terms with the largest number of + indices and the smallest number of - indices are most relevant.

Now consider the additional terms in the Lagrangian beyond those in Eqs. (F22), (F25), (F24), order by order in Y.

$G_N$ symmetry does not allow any linear terms other than Eq. (F25). At quadratic order the remaining possibilities are $\sum_{ab} Y_{ab}^2$ and $(\sum_{ab} Y_{ab})^2$. The former is redundant—it can be cancelled by a shift in Y because of the presence of $\sum_{ab} Y_{ab}^3$. We also check this below in the new parameterization. The latter, $(\sum_{ab} Y_{ab})^2$, is proportional to $(y_{--})^2$ by Eq. (F25), so is less relevant than the terms on the rhs of Eq. (F24) and can be neglected.

Cubic terms are obtained by contracting indices in

$$ Y_{ab}Y_{a'b'}Y_{a''b''}, \qquad (F32) $$

i.e., by setting some indices equal to others and summing. Left indices may only be set equal to other left indices, and similarly for right indices. This allows three different types of index contraction: an index (e.g., a) can be summed without being set equal to any other index; two indices can be set equal and then summed (e.g., a = a'); or three indices can be set equal and summed (a = a' = a''). When we rewrite the contracted expression in terms of $y_{\alpha\beta}y_{\alpha'\beta'}y_{\alpha''\beta''}$, the single, double, and triple index contractions lead to contractions with the tensors

$$ d_{\alpha}^{(1)} = \sum_{a} x_{\alpha}^{a}, \quad d_{\alpha\alpha'}^{(2)} = \sum_{a} x_{\alpha}^{a}x_{\alpha'}^{a}, \quad d_{\alpha\alpha'\alpha''}^{(3)} = \sum_{a} x_{\alpha}^{a}x_{\alpha'}^{a}x_{\alpha''}^{a}, \qquad (F33) $$

respectively, for the indices $\alpha, \alpha', \alpha''' \in \{+, -, 2, \dots, N\}$. For each of these, we may check the maximal value of

$$ \Delta \equiv n_{+} - n_{-}, \qquad (F34) $$

the difference in the number of + and - indices, that may appear on the right-hand side. ($\Delta$ should not be confused with a scaling dimension.) Since the difference in the total number of + and - indices is what determines the engineering dimension of a field or a product of fields [Eq. (F31)], identifying $\Delta_{max}$ for each type of index contraction allows us to say which types of contraction will give the most relevant cubic terms: they are those for which the sum of $\Delta_{max}$, over all contractions, is largest.

Explicitly,

$$ d_{\alpha}^{(1)} = (N, 2 - N, 1, \dots, 1)_{\alpha}. \qquad (F35) $$

However, y vanishes when contracted with $(0, 1, \dots, 1)$ [Eq. (F13)]. After dropping this part, and taking the limit
$N \to 0$ directly in the coefficients,

$$d_{\alpha}^{(1)} \to 2\delta_{\alpha, -\alpha} \qquad (\text{F36})$$

Therefore, this pattern of index contraction contributes $\Delta = -1$. Using similar simplifications,

$$d_{\alpha, \alpha'}^{(2)} \to 2(\delta_{+\alpha} \delta_{-\alpha'} + \delta_{-\alpha} \delta_{+\alpha'}) + \sum_{j=2}^{N} \delta_{j\alpha} \delta_{j\alpha'}. \qquad (\text{F37})$$

This pattern of index contraction contributes $\Delta = 0$. Finally, $d^{(3)}$ contains various patterns of index contraction with different values for $\Delta$:

$$d_{\alpha, \alpha', \alpha''}^{(3)} \to 2(\delta_{-\alpha} \delta_{+\alpha'} \delta_{+\alpha''} + \dots) \qquad (\text{F38})$$

$$+ \sum_{i=2}^{N} (\delta_{+\alpha} \delta_{i\alpha'} \delta_{i\alpha''} + \dots) \qquad (\text{F39})$$

$$+ \sum_{i} \delta_{i\alpha} \delta_{i\alpha'} \delta_{i\alpha''} \qquad (\text{F40})$$

$$- \sum_{i} (\delta_{-\alpha} \delta_{i\alpha'} \delta_{i\alpha''} + \dots) \qquad (\text{F41})$$

$$+ 2\delta_{-\alpha} \delta_{-\alpha'} \delta_{-\alpha''}. \qquad (\text{F42})$$

(Ellipses indicate terms related to those shown by cyclic permutations of $\alpha, \alpha', \alpha''$.) The first two lines on the rhs have $\Delta = 1$; all the others have smaller values of $\Delta$. If we are interested only in keeping the most relevant terms in a given expression, we can truncate to only the first two lines:

$$\tilde{d}_{\alpha, \alpha', \alpha''}^{(3)} \equiv 2 (\delta_{-\alpha} \delta_{+\alpha'} \delta_{+\alpha''} + \dots) \qquad (\text{F43})$$

$$+ \sum_{i=2}^{N} (\delta_{+\alpha} \delta_{i\alpha'} \delta_{i\alpha''} + \dots). \qquad (\text{F44})$$

Therefore, the most relevant cubic terms allowed by $G_N$ symmetry are those with $\Delta = 2$ (in total, i.e., counting both row and column indices) that arise by discarding the less-relevant parts of $\sum_{ab} Y_{ab}^3$:

$$\sum_{ab} Y_{ab}^3 \to \sum_{ab} \tilde{d}_{\alpha, \alpha', \alpha''}^{(3)} \tilde{d}_{\beta, \beta', \beta''}^{(3)} y_{\alpha\beta} y_{\alpha'\beta'} y_{\alpha''\beta''}, \qquad (\text{F45})$$

which is equal to

$$\left(\sum_{ab} Y_{ab}^3\right)_{\Delta=2} = 12y_{++}(y_{++}y_{--} + 2y_{+-}y_{-+}) \\ + 12y_{++}\left(y_{-k}y_{+k} + y_{j-y_j+} + \frac{y_{jk}y_{jk}}{4}\right) \\ + 6\left(y_{+-}y_{j+y_j+} + y_{-+}y_{+k}y_{+k}\right) \\ + 6y_j+y_{+k}y_{jk}. \qquad (\text{F46})$$

(Repeated $j$ or $k$ indices are summed from 2 to $N$.) All other contractions of $YYY$, such as $\sum_{aa'bb'} Y_{ab}Y_{ab'}Y_{a'b'}$, give terms that are strictly less relevant according to the engineering dimensions.

Let us check that, having dropped less relevant terms from the Langrangian in the process of rewriting it in terms of $y$, the coupling of the quadratic term

$$\frac{1}{4} \sum_{ab} Y_{ab}^2 = 2(y_{+-}y_{-+} + y_{++}y_{--}) + y_{j+y_j-} + y_{+k}y_{-k} \\ + \frac{y_{jk}y_{jk}}{4}, \qquad (\text{F47})$$

is still redundant. The shift $Y_{ab} \to Y_{ab} + c$ for each element corresponds to shifting $y_{++} \to y_{++} + c$ and leaving other elements of $y$ unchanged [by Eq. (F12)]. Under this shift, the cubic term, Eq. (F46), generates precisely Eq. (F47), with a coefficient of order $c$ [and the linear term in Eq. (F25) with a coefficient of order $c^2$]. Therefore, by an appropriate shift of $y_{++}$ we may eliminate the quadratic term, Eq. (F47).

[1] T. Basché, S. Kummer, and C. Bräuchle, Direct spectroscopic observation of quantum jumps of a single molecule, Nature 373, 132 (1995).

[2] S. Gleyzes, S. Kuhr, C. Guerlin, J. Bernu, S. Deléglise, U. Busk Hoff, M. Brune, J.-M. Raimond, and S. Haroche, Quantum jumps of light recording the birth and death of a photon in a cavity, Nature 446, 297 (2007).

[3] R. Vijay, D. H. Slichter, and I. Siddiqi, Observation of Quantum Jumps in a Superconducting Artificial Atom, Phys. Rev. Lett. 106, 110502 (2011).

[4] L. Robledo, L. Childress, H. Bernien, B. Hensen, P. F. A. Alkemade, and R. Hanson, High-fidelity projective read-out of a solid-state spin quantum register, Nature 477, 574 (2011).

[5] Z. K. Minev, S. O. Mundhada, S. Shankar, P. Reinhold, R. Gutiérrez-Jáuregui, R. J. Schoelkopf, M. Mirrahimi, H. J. Carmichael, and M. H. Devoret, To catch and reverse a quantum jump mid-flight, Nature 570, 200 (2019).

[6] B. Skinner, J. Ruhman, and A. Nahum, Measurement-Induced Phase Transitions in the Dynamics of Entanglement, Phys. Rev. X 9, 031009 (2019).

[7] Y. Li, X. Chen, and M. P. A. Fisher, Quantum Zeno effect and the many-body entanglement transition, Phys. Rev. B 98, 205136 (2018).

[8] A. Chan, R. M. Nandkishore, M. Pretko, and G. Smith, Unitary-projective entanglement dynamics, Phys. Rev. B 99, 224307 (2019).

[9] Y. Li, X. Chen, and M. P. A. Fisher, Measurement-driven entanglement transition in hybrid quantum circuits, Phys. Rev. B 100, 134306 (2019).

[10] M. Szyniszewski, A. Romito, and H. Schomerus, Entanglement transition from variable-strength weak measurements, Phys. Rev. B 100, 064204 (2019).
[11] S. Choi, Y. Bao, X.-L. Qi, and E. Altman, Quantum Error Correction in Scrambling Dynamics and Measurement-Induced Phase Transition, Phys. Rev. Lett. 125, 030505 (2020).

[12] M. J. Gullans and D. A. Huse, Dynamical purification phase transitions induced by quantum measurements, arXiv:1905.05195 [quant-ph] (2019).

[13] Y. Bao, S. Choi, and E. Altman, Theory of the phase transition in random unitary circuits with measurements, Phys. Rev. B 101, 104301 (2020).

[14] C.-M. Jian, Y.-Z. You, R. Vasseur, and A. W. W. Ludwig, Measurement-induced criticality in random quantum circuits, Phys. Rev. B 101, 104302 (2020).

[15] Y. Li, X. Chen, A. W. W. Ludwig, and M. P. A. Fisher, Conformal invariance and quantum non-locality in hybrid quantum circuits, arXiv:2003.12721 [quant-ph] (2020).

[16] A. Zabalo, M. J. Gullans, J. H. Wilson, S. Gopalakrishnan, D. A. Huse, and J. H. Pixley, Critical properties of the measurement-induced transition in random quantum circuits, Phys. Rev. B 101, 060301 (2020).

[17] M. J. Gullans and D. A. Huse, Scalable Probes of Measurement-Induced Criticality, Phys. Rev. Lett. 125, 070606 (2020).

[18] Q. Tang and W. Zhu, Measurement-induced phase transition: A case study in the nonintegrable model by density-matrix renormalization group calculations, Phys. Rev. Res. 2, 013022 (2020).

[19] Y. Fuji and Y. Ashida, Measurement-induced quantum criticality under continuous monitoring, Phys. Rev. B 102, 054302 (2020).

[20] O. Lunt and A. Pal, Measurement-induced entanglement transitions in many-body localized systems, arXiv:2005.13603 [quant-ph] (2020).

[21] M. Szyniszewski, A. Romito, and H. Schomerus, Universality of entanglement transitions from stroboscopic to continuous measurements, arXiv:2005.01863 [cond-mat.stat-mech] (2020).

[22] X. Turkeshi, R. Fazio, and M. Dalmonte, Measurement-induced criticality in (2 + 1)-dimensional hybrid quantum circuits, Phys. Rev. B 102, 014315 (2020).

[23] M. V. Regemortel, Z.-P. Cian, A. Seif, H. Dehghani, and M. Hafezi, Entanglement entropy scaling transition under competing monitoring protocols, arXiv:2008.08619 [quant-ph] (2020).

[24] R. Fan, S. Vijay, A. Vishwanath, and Y.-Z. You, Self-organized error correction in random unitary circuits with measurement, arXiv:2002.12385 [cond-mat.stat-mech] (2020).

[25] O. Shtanko, Y. A. Kharkov, L. P. García-Pintos, and A. V. Gorshkov, Classical models of entanglement in monitored random circuits, arXiv:2004.06736 [cond-mat.dis-nn] (2020).

[26] S. Vijay, Measurement-driven phase transition within a volume-law entangled phase, arXiv:2005.03052 [quant-ph] (2020).

[27] Y. Li and M. P. A. Fisher, Statistical Mechanics of Quantum Error-Correcting Codes, arXiv:2007.03822 [quant-ph] (2020).

[28] D. Aharonov, Quantum to classical phase transition in noisy quantum computers, Phys. Rev. A 62, 62311 (2000).

[29] D. Gottesman, An Introduction to Quantum Error Correction and Fault-Tolerant Quantum Computation, arXiv:0904.2557 [quant-ph] (2009).

[30] L. Bonnes and A. M. Läuchli, Superoperators vs. Trajectories for matrix product state simulations of open quantum system: A case study, arXiv:1411.4831 [cond-mat.quant-gas] (2014).

[31] J. Napp, R. L. L. Placa, A. M. Dalzell, F. G. S. L. Brandao, and A. W. Harrow, Efficient classical simulation of random shallow 2D quantum circuits, arXiv:2001.00021 [quant-ph] (2019).

[32] M. B. Plenio and P. L. Knight, The quantum-jump approach to dissipative dynamics in quantum optics, Rev. Mod. Phys. 70, 101 (1998).

[33] A. J. Daley, Quantum trajectories and open many-body quantum systems, Adv. Phys. 63, 77 (2014).

[34] T. Zhou and A. Nahum, Emergent statistical mechanics of entanglement in random unitary circuits, Phys. Rev. B 99, 174205 (2019).

[35] P. Hayden, S. Nezami, X.-L. Qi, N. Thomas, M. Walter, and Z. Yang, Holographic duality from random tensor networks, J. High Energy Phys. 2016, 1 (2016).

[36] R. Vasseur, A. C. Potter, Y.-Z. You, and A. W. W. Ludwig, Entanglement transitions from holographic random tensor networks, Phys. Rev. B 100, 134203 (2019).

[37] A. Nahum and B. Skinner, Entanglement and dynamics of diffusion-annihilation processes with majorana defects, Phys. Rev. Res. 2, 023288 (2020).

[38] S. Sang and T. H. Hsieh, Measurement protected quantum phases, arXiv:2004.09509 [cond-mat.stat-mech] (2020).

[39] A. Lavasani, Y. Alavirad, and M. Barkeshli, Measurement-induced topological entanglement transitions in symmetric random quantum circuits, arXiv:2004.07243 [quant-ph] (2020).

[40] M. Ippoliti, M. J. Gullans, S. Gopalakrishnan, D. A. Huse, and V. Khemani, Entanglement phase transitions in measurement-only dynamics, arXiv:2004.09560 [quant-ph] (2020).

[41] N. Lang and H. P. Büchler, Entanglement Transition in the Projective Transverse Field Ising Model, arXiv:2006.09748 [cond-mat.dis-nn] (2020).

[42] A. Nahum, S. Vijay, and J. Haah, Operator Spreading in Random Unitary Circuits, Phys. Rev. X 8, 021014 (2018).

[43] T. Zhou and A. Nahum, Entanglement Membrane in Chaotic Many-Body Systems, Phys. Rev. X 10, 031066 (2020).

[44] N. Hunter-Jones, Unitary designs from statistical mechanics in random quantum circuits, arXiv:1905.12053 [quant-ph] (2019).

[45] H. Liu and S. Vardhan, Entanglement entropies of equilibrated pure states in quantum many-body systems and gravity, arXiv:2008.01089 [hep-th] (2020).

[46] P. Zanardi, Entanglement of quantum evolutions, Phys. Rev. A 63, 40304 (2001).

[47] T. Prosen and I. Pižorn, Operator space entanglement entropy in a transverse Ising chain, Phys. Rev. A 76, 032316 (2007).

[48] I. Pižorn and T. Prosen, Operator space entanglement entropy in $xy$ spin chains, Phys. Rev. B 79, 184416 (2009).
[49] J. Dubail, Entanglement scaling of operators: A conformal field theory approach, with a glimpse of simulability of long-time dynamics in 1+1D, J. Phys. A 50, 234001 (2017).

[50] T. Zhou and D. J. Luitz, Operator entanglement entropy of the time evolution operator in chaotic systems, Phys. Rev. B 95, 094206 (2017).

[51] C. Jonay, D. A. Huse, and A. Nahum, Coarse-grained dynamics of operator and state entanglement, arXiv:1803.00089 [cond-mat.stat-mech] (2018).

[52] X. Cao, A. Tilloy, and A. D. Luca, Entanglement in a fermion chain under continuous monitoring, SciPost Phys. 7, 24 (2019).

[53] X. Chen, Y. Li, M. P. A. Fisher, and A. Lucas, Emergent conformal symmetry in nonunitary random dynamics of free fermions, Phys. Rev. Res. 2, 033017 (2020).

[54] O. Alberton, M. Buchhold, and S. Diehl, Trajectory dependent entanglement transition in a free fermion chain – from extended criticality to area law, arXiv:2005.09722 [cond-mat.stat-mech] (2020).

[55] In general the membrane tension can depend on the Rényi index $n$ [34]. It can also depend on the orientation of the membrane [51], but here we are interested only in membranes that are “horizontal” on large scales.

[56] L. Fidkowski, J. Haah, and M. B. Hastings, How dynamical quantum memories forget, arXiv:2008.10611 [quant-ph] (2020).

[57] M. J. Gullans and D. A. Huse (Private communication).

[58] B. Derrida and H. Spohn, Polymers on disordered trees, spin glasses, and traveling waves, J. Stat. Phys. 51, 817 (1988).

[59] Y.-Y. Shi, L.-M. Duan, and G. Vidal, Classical simulation of quantum many-body systems with a tree tensor network, Phys. Rev. A 74, 022320 (2006).

[60] L. Tagliacozzo, G. Evenbly, and G. Vidal, Simulation of two-dimensional quantum systems using a tree tensor network that exploits the entropic area law, Phys. Rev. B 80, 235127 (2009).

[61] V. Murg, F. Verstraete, O. Legeza, and R. M. Noack, Simulating strongly correlated quantum systems with tree tensor networks, Phys. Rev. B 82, 205105 (2010).

[62] P. Silvi, V. Giovannetti, S. Montangero, M. Rizzi, J. I. Cirac, and R. Fazio, Homogeneous binary trees as ground states of quantum critical hamiltonians, Phys. Rev. A 81, 062335 (2010).

[63] W. Li, J. von Delft, and T. Xiang, Efficient simulation of infinite tree tensor network states on the bethe lattice, Phys. Rev. B 86, 195137 (2012).

[64] N. Nakatani and G. K.-L. Chan, Efficient tree tensor network states (TTNS) for quantum chemistry: Generalizations of the density matrix renormalization group algorithm, J. Chem. Phys. 138, 134113 (2013).

[65] V. Murg, F. Verstraete, R. Schneider, P. R. Nagy, and O. Legeza, Tree tensor network state with variable tensor order: An efficient multireference method for strongly correlated systems, J. Chem. Theory Comput. 11, 1027 (2015).

[66] G. Vidal, Entanglement Renormalization, Phys. Rev. Lett. 99, 220405 (2007).

[67] B. Swingle, Entanglement renormalization and holography, Phys. Rev. D 86, 065007 (2012).

[68] J. Lopez-Piqueres, B. Ware, and R. Vasseur, Mean-field entanglement transitions in random tree tensor networks, Phys. Rev. B 102, 064202 (2020).

[69] R. N. C. Pfeifer, G. Evenbly, and G. Vidal, Entanglement renormalization, scale invariance, and quantum criticality, Phys. Rev. A 79, 040301 (2009).

[70] The term “random tensor network” allows for almost any structure, so infinite numbers of universality classes can in principle be accessed, most of them extremely fine tuned.

[71] K. H. Fischer and J. A. Hertz, Spin Glasses (Cambridge University Press, Cambridge, UK, 1993), Vol. 1.

[72] For generic versions of the MPT and FMPT, the emergent space-time symmetry is of course partly a conjecture. It is perhaps made more plausible by the existence of such symmetries in some simpler limiting models. The minimal cut problem [6], and also some alternative $q \to \infty$ limits [13,14], map to percolation problems, which have this symmetry. Some measurement induced critical points with a free-fermion structure also map to conformally invariant models [37–39,41,54]. Conformal invariance in $1+1D$ Clifford measurement circuits has been demonstrated numerically [9,15]. There is numerical evidence that the dynamical exponent is unity for the Haar-random MPT [6]. Finally, we may use dual-unitary circuits to set up measurement circuits that have 90° rotational invariance in space time even microscopically (we discuss this elsewhere).

[73] $S_0$ counts the (logarithm of the) number of nonzero Schmidt values in the singular value decomposition of a state or, as we focus on here, an operator.

[74] F. Pastawski, B. Yoshida, D. Harlow, and J. Preskill, Holographic quantum error-correcting codes: Toy models for the bulk/boundary correspondence, J. High Energy Phys. 2015, 149 (2015).

[75] H. Casini, H. Liu, and M. Mezei, Spread of entanglement and causality, J. High Energy Phys. 2016, 77 (2016).

[76] A. Nahum, J. Ruhman, S. Vijay, and J. Haah, Quantum Entanglement Growth under Random Unitary Dynamics, Phys. Rev. X 7, 031016 (2017).

[77] J. T. Chayes, L. Chayes, and R. Durrett, Critical behavior of the two-dimensional first passage time, J. Stat. Phys. 45, 933 (1986).

[78] B. Bollobás, Random Graphs, Cambridge Studies in Advanced Mathematics (Cambridge University Press, 2001), 2nd ed.

[79] C. M. Fortuin and P. W. Kasteleyn, On the random-cluster model: I. Introduction and relation to other models, Physica 57, 536 (1972).

[80] R. K. P. Zia and D. J. Wallace, Critical behaviour of the continuous $n$-component Potts model, J. Phys. A 8, 1495 (1975).

[81] D. J. Amit, Renormalization of the Potts model, J. Phys. A 9, 1441 (1976).

[82] J. Cardy, Scaling and Renormalization in Statistical Physics, Cambridge Lecture Notes in Physics (Cambridge University Press, 1996).
[83] A. Aharony, Y. Gefen, and A. Kapitulnik, Scaling at the percolation threshold above six dimensions, *J. Phys. A* **17**, L197 (1984).

[84] J. L. Cardy, Critical percolation in finite geometries, *J. Phys. A* **25**, L201 (1992).

[85] A conclusive numerical check of the exponent in Eq. (27) would require larger system sizes since it requires the scaling function in Fig. 6 to have large negative argument.

[86] We write this formula for the case where $r_c - r$ is of order 1. Otherwise the lower limit on the range may involve a critical timescale that is larger than 1 but much smaller than $\tau$.

[87] For a naive picture of the scaling of $S_0$ with $N$, consider a minimal cut that has zero temporal width within the entangled phase. There are $O(t)$ choices for the time at which to place the cut. Each choice has a random cost, which for the present illustration we assume to be Gaussian with variance $N$, arising from a sum of $O(N)$ random contributions. Taking the minimum gives the formula above. The second term is subleading so long as $\ln t \ll N$.

[88] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein, *Introduction to Algorithms* (MIT Press, Cambridge, MA, USA, 2009).

[89] As usual, this geometrical disconnection is reflected in the vanishing of the zeroth Rényi entropy $S_0$ between apex and base.

[90] This is a weaker condition than that satisfied by the Haar ensemble of two-site unitaries, which is invariant under left or right multiplication by any two-site unitary.

[91] This is because when we redraw the section of the circuit as a tree, we turn some of the unitaries upside down (which amounts to taking a transpose) or reflect them left-right. The discrete invariances of the distributions mean that we do not need to keep track of this.

[92] The link may contain multiple projections, but this reduces to the case with a single projection. In the case where all the projection operators are identical, this is immediate. More generally, it holds because the distribution of unitaries is invariant under single-site rotations (and because the overall normalization factor for the tree is not important).

[93] In linear algebra terms, the tree represents a $2 \times 2^\ell$ matrix: 2 is the number of states associated with the bond at the apex, and $2^\ell$ is the number of states associated with the $\ell$ spins at the base. We are interested in the singular values in the SVD for this matrix.

[94] The trees we encounter in the circuit have a nonzero probability of terminating before they achieve $k$ generations, in which case $Z_k$ is identically zero. So long as we are below the classical transition, an order-1 fraction of the probability distribution for $Z_k$ is supported on nonzero values even when $k$ is finite but large. For other tree tensor networks $Z_k$ may be nonzero with probability 1 (for any finite $k$), in which case we simply define $\ln Z_k^{\text{typ}} = \langle \ln Z_k \rangle$.

[95] Considering $T \equiv T_{k+1}$ as a state, then $\lambda_1^2 + \lambda_2^2 = 1$ and $\lambda_1^4 + \lambda_2^4 = \text{Tr} \rho^2 / (\text{Tr} \rho)^2$, where $\rho$ is the un-normalized reduced density matrix for the spin at the apex: $\rho_{a,a'} = \sum_{bcd} T_{bcd}^a T_{bcd}^{*a'}$.

[96] This can be seen from the fact that $Z_{k+1}$ is manifestly independent of $Z_k^\dagger$ and $Z_k'$ in this limit.

[97] The latter converges at large $k$ to a constant below 1, so long as we are in the classically percolating phase for the tree. The value of this plateau can be set to zero by a linear transformation of $G$. See Sec. H.

[98] R. A. Fisher, The wave of advance of advantageous genes, *Ann. Eugen.* **7**, 355 (1937).

[99] The value $\lambda = 1$ appears here because it is the exponential decay constant for the initial condition at $k = 0$: $G_0(x) \sim 1 - \langle Z_0 \rangle e^{-x}$ at large $x$. At late times the solution converges to a traveling wave in which the decay constant is $\min\{1, \lambda_*\}$. That is, the decay constant is either that of the initial condition, or that of the minimal-speed wave, whichever decays more slowly as $x \to \infty$.

[100] Recall that we assume various invariances of the distribution of unitaries to simplify the treatment (Sec. B).

[101] Equation (83) is the exact growth rate of $\langle Z_k \rangle$ for any $p$ in the linear problem, but it is only for $p < p_{\text{glass}}$ that $Z_k^{\text{typ}}$ has the same growth rate as $\langle Z_k \rangle$. It is $c_p$ as defined by $Z_k^{\text{typ}}$ that will be relevant when we include nonlinearity. The failure of the annealed approximation when $\lambda_* < 1$ is because, in this regime the distribution of $Z$ becomes broad in the sense that $\lim_{k \to \infty} \langle Z/Z_k^{\text{typ}} \rangle = \infty$. [In this regime the tail in the probability distribution for $\ln(Z/Z_k^{\text{typ}})$ decays as $e^{-\lambda_* \ln(Z/Z_k^{\text{typ}})}$.]

[102] These have $Z=0$, so do not contribute to the recursion.

[103] This is because, close to the classical transition, only a small fraction, of order $(2/3-p)$, of subtrees survive. As a result most nodes in a “pruned” tree have a single descendant, with a fraction of order $(2/3-p)^{m-1}$ having $m$ descendants ($m \le 3$).

[104] We note that due to the forced measurements there is a finite probability for $Z$ to be exactly zero, i.e., the distribution function has a $\delta$ function with a finite weight. These are instances where the tree is classically disconnected. When we compute $Z_k^{\text{typ}}$ and present distribution functions we do not include these trivial instances.

[105] J. D. Miller and B. Derrida, Weak-disorder expansion for the Anderson model on a tree, *J. Stat. Phys.* **75**, 357 (1994).

[106] C. Monthus and T. Garel, Anderson transition on the Cayley tree as a traveling wave critical point for various probability distributions, *J. Phys. A* **42**, 75002 (2008).

[107] I. García-Mata, O. Giraud, B. Georgeot, J. Martin, R. Dubertrand, and G. Lemarié, Scaling Theory of the Anderson Transition in Random Graphs: Ergodicity and Universality, *Phys. Rev. Lett.* **118**, 166801 (2017).

[108] The multiplicative factor appearing here is the probability that a tree never terminates. In the limit of large $k$ this is also the probability that $Z$ is nonzero.

[109] E. Brunet and B. Derrida, Shift in the velocity of a front due to a cutoff, *Phys. Rev. E* **56**, 2597 (1997).

[110] By examining the first-order equation satisfied by $R = \partial_x \ln H$ after dropping the $O(H^2)$ and $O(\sigma)$ terms, we see that the value of this slope at some arbitrarily chosen value $x$ in the vicinity of the origin must be tuned to the correct $O(1)$ value, $R_x$, in order to match onto the correct solution at large positive $x$, where the slope is $-1$. Further this value $R_x$ is necessarily less than $-\sqrt{\Delta}$, which is the approximate slope of Eq. (96) in the region where the sine is of order 1. To achieve this, we must take advantage
of the negative divergence of the cotangent in Eq. (97) at argument $\pi$. This means that to leading order in $\sigma$ we must have $\phi = \pi$.

[111] Defining $\mathbf{t}_{b_1 \dots b_l}$ as the vector with components $t_{b_1 \dots b_l}^a$, then
$$A_1 = \begin{aligned}[t]
& |\mathbf{t}_{11\dots1}|^{-4} (|\mathbf{t}_{21\dots1}|^2 |\mathbf{t}_{11\dots1}|^2 - |\mathbf{t}_{21\dots1}^\dagger \mathbf{t}_{11\dots1}|^2) = \\
& |\mathbf{t}_{11\dots1}|^{-4} |\mathbf{t}_{21\dots1}^1 \mathbf{t}_{11\dots1}^2 - \mathbf{t}_{11\dots1}^1 \mathbf{t}_{21\dots1}^2|^2 \quad \text{(for } A_m \text{ the "2" is the } m\text{th subscript). The argument in Appendix 2 making use of the invariance of } t \text{ under single-leg rotations also applies here.}
\end{aligned}$$

[112] The prefactor depends on the singular vectors at the top of the subtrees, but is of order 1.

[113] Here we use the U(2) invariance property of the unitaries to show that lower terms vanished by phase cancellation.

[114] We use a naive linear extrapolation: the detailed functional form of the subleading corrections to $S/N$ may depend on the precise regime of $N$ and $t$ (Sec. I).

[115] Note that the latter value of $r$ is right at the putative critical point. This value is chosen to avoid being above, or to close to, the classical critical point $r = 0.8$ where the network is trivially disconnected.

[116] These finite size effects, together with the fact that extracting $a(r)$ requires two separate fits [to go from $\mathcal{D}_{\text{typ}}$ to $\tau(r, N)$ and then to $a(r)$] also make it hard to estimate error bars on $a(r)$ accurately.

[117] S. J. Garratt and J. T. Chalker, Many-body quantum chaos and the local pairing of feynman histories, arXiv:2008.01697 [cond-mat.stat-mech] (2020).

[118] J. Cardy, Logarithmic conformal field theories as limits of ordinary CFTs and some physical applications, J. Phys. A 46, 494001 (2013).

[119] A. Nahum, Field theories for measurement and entanglement phase transitions (to be published).

[120] Simplifications arise in the fully unitary case. Even there, for general $N$ it is possible to have negative Boltzmann weights $W_\sigma$. Simplifications also arise at large local Hilbert space dimension [13,14,34,36].

[121] The dependence is only on the cycle structure of $\sigma^{-1}\tau$ thanks to $G_N$ symmetry.

[122] M. Stone and P. Goldbart, Mathematics for Physics: A Guided Tour for Graduate Students (Cambridge University Press, Cambridge, UK, 2009).

[123] A given realization has no symmetry. Because of the distribution we choose for the disorder there are statistical symmetries, i.e., symmetries of the disorder distribution. For example the distribution is invariant under $S \to -S$ on a given site. However, this symmetry does not act nontrivially on the effective field $X$, so we expect it could be broken without changing the universality class of the transition under discussion.

[124] We use the set of $N$ vectors $\mathbf{e}^1, \dots, \mathbf{e}^N$, each of $N-1$ components, that are familiar in the context of the Potts field theory [80,81] (see Appendix A 2) and satisfy $\sum_a e_a^a e_b^a = \delta_{\alpha\beta} : X_{ab} = \sum_{\alpha,\beta=1}^{N-1} (e_a^a e_b^b) \phi_{\alpha\beta}.$

[125] E. Pytte and J. Rudnick, Scaling, equation of state, and the instability of the spin-glass phase, Phys. Rev. B 19, 3603 (1979).

[126] J. R. L. de Almeida and D. J. Thouless, Stability of the sherrington-kirkpatrick solution of a spin glass model, J. Phys. A: Math. General 11, 983 (1978).

[127] G. Parisi, A sequence of approximated solutions to the SK model for spin glasses, J. Phys. A: Math. General 13, L115 (1980).

[128] J. Cardy, Logarithmic correlations in quenched random magnets and polymers, arXiv:cond-mat/9911024 [cond-mat.stat-mech] (1999).

[129] R. Vasseur, J. L. Jacobsen, and H. Saleur, Logarithmic observables in critical percolation, J. Stat. Mech. 2012, L07001 (2012).

[130] R. Vasseur and J. L. Jacobsen, Operator content of the critical Potts model in $d$ dimensions and logarithmic correlations, Nucl. Phys. B 880, 435 (2014).

[131] At first sight the interpretation of these additional fields may seem obscure, given that for a permutation matrix they are trivial constants. This may be more transparent in the approach of Sec. VI D.

[132] M. E. Fisher, Yang-Lee Edge Singularity and $\phi^3$ Field Theory, Phys. Rev. Lett. 40, 1610 (1978).

[133] J. L. Cardy, Nonperturbative aspects of supersymmetry in statistical mechanics, Phys. D 15, 123 (1985).

[134] J. L. Cardy and A. J. McKane, Field theoretic approach to the study of Yang-Lee and Griffiths singularities in the randomly diluted Ising model, Nucl. Phys. B 257, 383 (1985).

[135] J. Cardy, Exact scaling functions for self-avoiding loops and branched polymers, J. Phys. A 34, L665 (2001).

[136] A. Kaviraj, S. Rychkov, and E. Trevisani, Random field Ising model and Parisi-Sourlas supersymmetry. Part I. Supersymmetric CFT, J. High Energy Phys. 2020, 1 (2020).

[137] G. Parisi and N. Sourlas, Random Magnetic Fields, Supersymmetry, and Negative Dimensions, Phys. Rev. Lett. 43, 744 (1979).

[138] A. J. Bray and M. A. Moore, Scaling theory of the random-field Ising model, J. Phys. C: Solid State Phys. 18, L927 (1985).

[139] T. C. Lubensky and J. Isaacson, Field Theory for the Statistics of Branched Polymers, Gelation, and Vulcanization, Phys. Rev. Lett. 41, 829 (1978).

[140] G. Parisi and N. Sourlas, Critical Behavior of Branched Polymers and the Lee-Yang Edge Singularity, Phys. Rev. Lett. 46, 871 (1981).

[141] D. C. Brydges and J. Z. Imbrie, Branched polymers and dimensional reduction, Ann. Math. 1019 (2003).

[142] At the quadratic level, these various sectors can be replaced by a theory with fermions but without a replica limit [133,137] (for example, $N-2 \to -2$ copies of a bosonic theory can be replaced by a fermionic version of the theory).

[143] A. Chan, A. De Luca, and J. T. Chalker, Spectral Statistics in Spatially Extended Chaotic Quantum Many-Body Systems, Phys. Rev. Lett. 121, 060601 (2018).

[144] This is similar to what happens for the classical minimal cut. If $\eta$ may be treated as Gaussian, the correction to the entanglement is of order $\sqrt{N \ln t}$ in this regime.

[145] N. Hunter-Jones, Operator growth in random quantum circuits with symmetry, arXiv:1812.08219 [quant-ph] (2018).
[146] V. L. Nguen, B. Z. Spivak, and B. I. Shkovskii, Tunnel hopping in disordered systems, Sov. Phys. JETP **89**, 1770 (1985).

[147] Y.-C. Zhang, Directed polymers with complex amplitudes, *Europhys. Lett.* **9**, 113 (1989).

[148] E. Medina, M. Kardar, Y. Shapir, and X. R. Wang, Interference of Directed Paths in Disordered Systems, Phys. Rev. Lett. **62**, 941 (1989).

[149] J. Cook and B. Derrida, Lyapunov exponents of large, sparse random matrices and the problem of directed polymers with complex random weights, J. Stat. Phys. **61**, 961 (1990).

[150] Y. Y. Goldschmidt and T. Blum, Directed walks with complex random weights: Phase diagram and replica symmetry breaking, J. Phys. **12**, 1607 (1992).

[151] B. Derrida, M. R. Evans, and E. R. Speer, Mean field theory of directed polymers with random complex weights, Comm. Math. Phys. **156**, 221 (1993).

[152] D. Bernard and T. Jin, Open Quantum Symmetric Simple Exclusion Process, Phys. Rev. Lett. **123**, 080601 (2019).

[153] D. Bernard and T. Jin, Solution to the Quantum Symmetric Simple Exclusion Process: The Continuous Case, arXiv:2006.12222 [math-ph] (2020).

[154] R. Frassek, C. Giardina, and J. Kurchan, Duality in quantum transport models, arXiv:2008.03476 [cond-mat.stat-mech] (2020).

[155] F. Thiel, E. Barkai, and D. A. Kessler, First Detected Arrival of a Quantum Walker on an Infinite Line, Phys. Rev. Lett. **120**, 040502 (2018).

[156] T. Swann, J. T. Chalker, and A. Nahum, (to be published).

[157] A. Chan and J. T. Chalker, Loop soups in quasi-one dimension: The Anderson tower of states and the Poisson-Dirichlet distribution (to be published).

[158] C. W. J. Beenakker, Random-matrix theory of quantum transport, Rev. Mod. Phys. **69**, 731 (1997).

[159] F. Evers and A. D. Mirlin, Anderson transitions, Rev. Mod. Phys. **80**, 1355 (2008).

[160] Formally, when $b' > b$, we use as our modes the coefficients of the right eigenvectors of $T$, which is equivalent to analytically continuing the formula for the $b' < b$ case. The square root nonanalyticity does not appear in the physical quantities below, only in the intermediate formulas.

[161] A. Theumann, Bond percolation problem in a semi-infinite medium. Landau-Ginzburg theory, Phys. Rev. B **19**, 6295 (1979).

[162] J. P. Carton, Surface effects and percolation: The repulsive case, J. Phys. Lett. **41**, 175 (1980).

[163] K. De'Bell and J. W. Essam, Mean field theory of percolation with application to surface effects, J. Phys. A **14**, 1993 (1981).

[164] Translating these BCs into the continuum field theory Eq. (A30) gives boundary magnetic field terms in Eq. (A23) of the form $h\delta(t)e_\sigma\cdot\phi(t)$ (and similarly at the final time boundary) with $h \propto N^{1/2}$. At first sight one worries that the $N$ dependence of $h$ introduces another exponent that could appear in scaling forms. However, we believe that the basic point is just that $h$ diverges with $N$, so that the asymptotic scaling forms are those

of the $h \to \infty$ limit and the detailed $N$ dependence of $h$ determines only subleading corrections.

[165] A. Coniglio, in *Physics of Finely Divided Matter: Proceedings of the Winter School, Les Houches, France, March 25–April 5, 1985*, edited by N. Boccara and M. Daoud (Springer Science & Business Media, 2012), Vol. 5.

[166] M. Aizenman, On the number of incipient spanning clusters, Nucl. Phys. **B 485**, 551 (1997).

[167] L. de Arcangelis, Multiplicity of infinite clusters in percolation above six dimensions, J. Phys. A **20**, 3057 (1987).

[168] S. Fortunato, D. Stauffer, and A. Coniglio, Percolation in high dimensions is not understood, Phys. A **334**, 307 (2004).

[169] S. Fortunato, A. Aharony, A. Coniglio, and D. Stauffer, Number of spanning clusters at the high-dimensional percolation thresholds, Phys. Rev. E **70**, 056116 (2004).

[170] R. Kenna and B. Berche, Universal finite-size scaling for percolation theory in high dimensions, J. Phys. A **50**, 235001 (2017).

[171] In the saddle-point language this factor would come from a fluctuation integral (cf. e.g., [183,184]): it is consistent with $(Q - 2) \to -1$ modes with eigenvalue $\propto |\delta r|$.

[172] $S_0$ for a single spin is either 0 or 1 bits, with the former holding if the spin is connected to no other spins on the final time boundary. This may be written as a one-point function of the spin operator in a system where other boundary spins are fixed. $I_0$ for two spins is either 0, 1, or 2 bits, and is given by a minimal cut formula. If we assume that the probabilities for both $I_0 = 1$ and $I_0 = 2$ have the same scaling form (this can be demonstrated in $1+1D$ [6]) then we can focus on the simpler case $I_0 = 2$, which occurs only if the two spins are connected to each other, but not to any other spins on the final time boundary. This is the two-point function of the same Potts operator.

[173] T. C. Lubensky and M. H. Rubin, Critical phenomena in semi-infinite systems. II. Mean-field theory, Phys. Rev. B **12**, 3885 (1975).

[174] A. J. Bray and M. A. Moore, Critical behaviour of semi-infinite systems, J. Phys. A **10**, 1927 (1977).

[175] C. Candu, J. L. Jacobsen, N. Read, and H. Saleur, Universality classes of polymer melts and conformal sigma models, J. Phys. A **43**, 142001 (2010).

[176] A. Nahum, J. T. Chalker, P. Serna, M. Ortuño, and A. M. Somoza, Phase transitions in three-dimensional loop models and the $cP^{n-1}$ sigma model, Phys. Rev. B **88**, 134411 (2013).

[177] A. Nahum and J. T. Chalker, Universal statistics of vortex lines, Phys. Rev. E **85**, 031141 (2012).

[178] A. Nahum, P. Serna, A. M. Somoza, and M. Ortuño, Loop models with crossings, Phys. Rev. B **87**, 184204 (2013).

[179] C. Nadal, S. N. Majumdar, and M. Vergassola, Statistical distribution of quantum entanglement for a random bipartite state, J. Stat. Phys. **142**, 403 (2011).

[180] M. Kieburg, A. B. J. Kuijlaars, and D. Stivigny, Singular value statistics of matrix products with truncated unitary matrices, Int. Math. Res. Not. **2016**, 3392 (2015).
[181] F. Mezzadri, How to generate random matrices from the classical compact groups, arXiv:math-ph/0609050 [math-ph] (2006).

[182] A. Savitzky and M. J. E. Golay, Smoothing and differentiation of data by simplified least squares procedures, *Anal. Chem.* **36**, 1627 (1964).

[183] J. L. Cardy, Electron localization in disordered systems and classical solutions in Ginzburg-Landau field theory, *J. Phys. C* **11**, L321 (1978).

[184] R. Hayn and W. John, Instanton approach to the conductivity of a disordered solid, *Nucl. Phys. B* **348**, 766 (1991).