Quantum Error Correction from Complexity in Brownian SYK

We study the robustness of quantum error correction in a one-parameter ensemble of codes generated by the Brownian SYK model, where the parameter quantifies the encoding complexity. The robustness of error correction by a quantum code is upper bounded by the"mutual purity"of a certain entangled state between the code subspace and environment in the isometric extension of the error channel, where the mutual purity of a density matrix $\rho_{AB}$ is the difference $\mathcal{F}_\rho (A:B) \equiv \mathrm{Tr}\;\rho_{AB}^2 - \mathrm{Tr}\;\rho_A^2\;\mathrm{Tr}\;\rho_B^2$. We show that when the encoding complexity is small, the mutual purity is $O(1)$ for the erasure of a small number of qubits (i.e., the encoding is fragile). However, this quantity decays exponentially, becoming $O(1/N)$ for $O(\log N)$ encoding complexity. Further, at polynomial encoding complexity, the mutual purity saturates to a plateau of $O(e^{-N})$. We also find a hierarchy of complexity scales associated to a tower of subleading contributions to the mutual purity that quantitatively, but not qualitatively, adjust our error correction bound as encoding complexity increases. In the AdS/CFT context, our results suggest that any portion of the entanglement wedge of a general boundary subregion $A$ with sufficiently high encoding complexity is robustly protected against low-rank errors acting on $A$ with no prior access to the encoding map. From the bulk point of view, we expect such bulk degrees of freedom to be causally inaccessible from the region $A$ despite being encoded in it.


Introduction
The bulk-to-boundary map in AdS/CFT has a rich structure. For any boundary subregion A, the associated Ryu-Takayanagi surface [1] singles out a certain subregion a of the bulk spacetime called the entanglement wedge of A [2]. The AdS/CFT map then satisfies subregion A A Figure 1: The entanglement wedge of a boundary subregion A can have a rich substructure. The outer white region on the left is the causal wedge of A, while the slightly darker grey region bounded by a non-minimal QES (shown in red) is the simple wedge. Beyond this lies the python's lunch (dark grey). duality: bulk semi-classical degrees of freedom in a are encoded within A and are protected against erasures in A. Furthermore, bulk operators within the entanglement wedge a can be reconstructed as boundary operators localized within the boundary region A, a property sometimes known as entanglement wedge reconstruction [3]. Using the language of [3][4][5][6][7][8][9][10], these properties hold because the Ryu-Takayanagi formula and its quantum generalizations imply that the bulk-to-boundary map in AdS/CFT is a quantum error correcting code with complementary recovery, where the entanglement wedge a of A is protected against the erasure of A, while a is protected against the erasure of A.
Recent progress points towards a sharper characterization of the structure of entanglement wedges that appear in holography and its generalizations [11][12][13][14][15][16]. Given a general boundary subregion, the corresponding entanglement wedge has a layered structure, i.e., it can be broken up into three regions: the causal wedge, the simple wedge [14,15] and the python's lunch [17]. These three regions of the entanglement wedge are defined as follows (see Figure 1): the causal wedge is the region in the bulk which is causally accessible from the boundary, i.e., a boundary observer in the domain of dependence of A can send signals to and receive signals from all points in the causal wedge of A. The simple wedge is defined to be the bulk domain of dependence of the homology region between A and the outermost quantum extremal surface (QES) which need not be minimal among all the QESs associated with A. The simple wedge is generically larger than the causal wedge, and so there are points in the simple wedge which are out of causal contact with the domain of dependence of A. However, it has been argued that the simple wedge can always be brought in causal contact with the boundary by performing backwards and forwards Lorentzian time evolution with sources turned on to de-focus the causal horizons [15]. Finally, the python's lunch region is defined as the portion of the entanglement wedge which lies between the outermost QES and the minimal QES. This region is causally inaccessible from the boundary subregion A, and furthermore since it lies behind an extremal surface it cannot be brought into causal contact with the domain of dependence of A (in contrast with the simple wedge); this follows from the fact that extremal surfaces must always lie behind causal horizons. This seems to lead to a puzzle -on the one hand, bulk operators in the python's lunch are encoded in A and in particular one should be able to create a semi-classical bulk excitation in the python's lunch via an operator acting on the domain of dependence of A. On the other hand, semi-classical gravity seems to forbid this! The evaporating black hole provides a context where this apparent contradiction is particularly sharp. Beyond the Page time, a portion of the black hole interior -the island -lies in the entanglement wedge of the radiation (see [11][12][13][18][19][20][21][22][23][24][25][26] for a partial list of articles discussing this phenomenon for black hole and cosmological horizons). But this portion lies behind a non-minimal QES, namely the empty surface, 1 and therefore constitutes a python's lunch. While one should be able to manipulate operators in the island by quantum operations on the radiation, such operations seem to blatantly violate semi-classical bulk causality. A potential way out is suggested by bounds coming from computational complexity [14,15,17,[27][28][29][30][31][32] -we expect that the encoding map for excitations in the python's lunch region is extremely complex, perhaps exponentially so in the number of qubits, and so any computationally bounded observer (with access only to sub-exponential operations on the radiation) will be unable to manipulate the degrees of freedom in the island. This is how we expect that semi-classical bulk causality will be approximately respected. On the other hand, certain finely tuned, exponentially complex operations on the radiation should be able to manipulate degrees of freedom in the island, but the gravitational mechanism for this involves Euclidean wormholes.
We can get an intuition for why complexity can protect information in this way from an analogy to older results concerning the complexity of black hole microstates and the difficulty of using simple probes to extract information about them [33][34][35]. Consider, for example, a Schwarzschild black hole of mass M in AdS 5 with a length scale . A microstate of this black hole is described in the dual SU (N ) Yang-Mills theory with 16 supersymmetries, by an 1 In more detail, in the toy models where these calculations are possible, the radiation is extracted into an auxiliary reservoir that is not geometrically connected to the island. Even in the absence of a geometric connection, there is still an obvious candidate extremal surface which one can consider as bounding the region dual to the radiation, namely the "empty surface". By this, we mean that the entire black hole spacetime is taken to lie "outside" the would-be entanglement wedge. After the Page time, this surface is no longer the quantum minimal surface which computes the radiation entropy, and the true minimum QES lies in the spacetime near the black hole horizon. Nevertheless, the island is "behind the empty extremal surface" from the point of view of the radiation.
operator O of dimension ∆ = M ∼ N 2 . O is roughly a polynomial of length N 2 built from the elementary fields of the Yang-Mills theory (a gauge field A µ , fermions ψ a and three complex adjoint scalars X, Y, Z) and their derivatives, with indices contracted to make the polynomial gauge and Lorentz invariant. Almost all such long polynomials are random sequences of fields and derivatives up to constraints of gauge and Lorentz invariance. A light probe of the state like the graviton corresponds to an operator of dimension O(1), like P = Tr(XX). The question is whether a measurement, modeled as a correlation function in the state created by O, 0|O † P † P O|0 , can reveal information about the identity of O. The authors of [33,34] argue that the answer is "no" because of the universal statistics of random polynomials, which mean that almost all O will lead to a similar sum of terms from contractions between the fields in the probe and the fields in O in evaluating the correlator. As such, simple (i.e., low-dimension) probes cannot reveal the microstate, but an observer with prior knowledge of the state could construct a fine-tuned, complex probe to check that knowledge, by choosing these probes to match long sequences of the fields composing O. One expects the situation in the python's lunch inside an evaporating black hole to be somewhat analogous: a highly complex encoding map prevents simple operations in the radiation from affecting the black hole interior, but if the encoding map is accessible then finely tuned, complex operations affecting the interior may be performed more easily.
Kim, Preskill and Tang (KPT) have sharpened these expectations [30]. They suggested that the encoding of the black hole interior degrees of freedom in the radiation, thought of as a quantum error correcting code, has robust error correction properties against low-rank, computationally bounded errors on the radiation, or more precisely, errors which effectively see the radiation density matrix as thermal. Note that this is not the standard error correction one encounters in the context of subregion duality; in the KPT formulation, bulk degrees of freedom in the island -while being encoded in the radiation -are nevertheless approximately (up to corrections exponentially small in the black hole entropy) protected against certain errors acting on the radiation itself. KPT then argued that this approximate error correction implies the existence of "ghost logical" operators which act on the radiation to mimic bulk operators in the island and at the same time commute with computationally bounded operators on the radiation -thus realizing the approximate causality of the black hole spacetime. The language of quantum error correction thus enables one to formulate and address the question of bulk causality in a universal manner.
Recently, the novel error correction in evaporating black holes proposed by KPT was tested in a toy model for an evaporating black hole in Jackiw-Teitelboim gravity [31], and it was argued that the bulk degrees of freedom in the island are protected against a large class of low-rank error operations on the radiation which do not have access to the details of the microscopic black hole state. The low-rank criterion can be formalized as a bound on the coherent information of the error in terms of the black hole entropy. In [31], it was also conjectured that this same robust error correction should also work in the python's lunch portions of more general entanglement wedges. The underlying reason is the high complexity of encoding in the python's lunch. The rough picture is the same as KPT -as the encoding map becomes sufficiently complex, any generic, low-rank error operation involving "simple" operations sees only a coarse-grained 2 density matrix on the boundary subregion, with no sign of the encoded subspace. In other words, the encoded subspace gets lost within the exponentially large Hilbert space of the boundary subregion. This "complexity-protected error correction" makes it possible for the semi-classical degrees of freedom in the python's lunch to be encoded in a boundary subregion and yet be causally inaccessible from it using simple probes.
The purpose of this paper is to demonstrate the above phenomenon in a toy model where the behavior of the encoding complexity is known more or less by construction. Such control is difficult to achieve directly in real holography because proving results about the complexity of the bulk-to-boundary map (without resorting to toy models like tensor networks) in different regions of the entanglement wedge is generically a very difficult task. Rather than a single code, here we consider an ensemble of quantum error correcting codes of the type relevant for entanglement wedge reconstruction in AdS/CFT. Since we want control over the complexity of encoding, our ensemble of codes is generated by picking the encoding map from an ensemble of unitaries with fixed circuit complexity.
We accomplish this by taking these unitaries to be time evolution operators U (T ) = T exp[−i T 0 dt H(t)] in the Brownian Sachdev-Ye-Kitaev (SYK) model [37][38][39][40][41], a quantum mechanical theory of N Majorana fermions. The SYK model here is merely a trick to generate an ensemble of unitaries with fixed complexity, parametrized by the number T . When T is small, the corresponding set of unitary operators is clustered around the identity operator, but as T → ∞ this set grows [42][43][44][45][46] to cover (modulo global symmetries) the entire unitary group. When it covers the entire unitary group, the typical complexity of an operator in the set is exponentially large [47]. So, computing the average error correction properties of such sets gives us some insight into the behavior of a family of codes with increasing complexity.
In this paper we will consider typical, low-rank errors with no prior access to the encoding map, and acting on a small fixed fraction of the physical Hilbert space. As a particular instance of such errors, we will consider the erasure of a small fraction of the physical Hilbert space. In quantum information theory, it is standard to model an error in terms of coupling to an external environment and tracing out the environment. The error correction properties of the code can then be studied in terms of the amount of correlation generated by the error between the code subspace and the environment. Error correction works with high accuracy when these correlations are suppressed by a large parameter e.g. the dimension of the physical Hilbert space.
We will study a particular measure of correlation, namely the "mutual purity" between the code subspace and the environment. We define the mutual purity F ρ (A : B) of a density matrix ρ between Hilbert subsystems H A and H B as Tr ρ 2 AB − Tr ρ 2 A Tr ρ 2 B . The fact that this quantity is a good measure of error correction is rigorously justified in Appendix B. Our main result is that for the Brownian SYK ensemble of quantum error correcting codes, there are three complexity regimes of interest.
(i) For T smaller than a scrambling time T ∼ log N (i.e., low encoding complexity) the erasure of a small fraction of the physical qubits generate an O(1) amount of correlation that decays exponentially with T between the code subspace and the environment, and thus there is no robust quantum error correction.
(ii) For T > log N , the mutual purity becomes O(1/N ) but keeps decaying further as the complexity T increases.
(iii) When T ∼ N , the mutual purity becomes exponentially small in N ; at this point, there is an O(e −N ) residual correlation generated by the error which is unavoidable.
The third and final regime corresponds to an exchange of dominance between a leading saddle point and a subleading saddle point 3 in the Brownian SYK calculation, analogously to the exchange of dominance between a disconnected geometry and the Euclidean wormhole in gravity. This quantitative hierarchy of complexity-protected error correction, ranging from a fragile encoding at O(1) complexity, through a logarithmic complexity regime of reasonable protection, and finally an emergent robust error correction at large encoding complexity, is the central result of this paper. We regard this as a step towards understanding the structure of general entanglement wedges ( Figure 1) from the boundary perspective in terms of quantum error correction.
Three sections follow. In Section 2 we review the necessary ideas from quantum error correction. We discuss the class of errors of interest, and show that the mutual purity which is relevant for recovery from these errors can be expressed in terms of the standard purity 4 of a certain density matrix constructed using the encoding map. We also briefly review the Brownian SYK model. In Section 3, we compute this purity in the large N limit using the Brownian SYK time evolution operator to model the encoding map. We conclude with a discussion in Section 4. In Appendix A, we give a Hamiltonian treatment of Brownian SYK to complement the path integral discussion in the main text and in Appendix B we prove that the mutual purity provides a bound on the error correction properties of an encoding map.

Brief review of quantum error correction
The mathematical framework for quantum error correction involves an isometric embedding of a small "code subspace" H code into a larger Hilbert space H phys : It is standard to model the error and recovery operations as completely positive trace-preserving linear maps, or "quantum channels". Any such map E has a representation in terms of its Kraus operators {E m } [48,49]: The minimum number of Kraus operators needed to implement a particular channel is called the rank of the channel. These quantum channels act on physical density matrices, and the goal of error correction is to determine for a given error channel E whether or not there exists a recovery channel R which restores the state ρ code : On the right hand side, we have in mind that the recovery channel has eliminated redundant portions of H phys , leaving behind precisely the matrix ρ code on the remaining subspace of H phys .
A second, convenient description of a quantum channel is given by its isometric extension, also known as its Stinespring dilation: we describe it as coupling the physical system via a unitary operator U E to an auxiliary environment with Hilbert space H env spanned by basis elements {|e m env }, initially in some fiducial state |e 0 . The action of the channel E on ρ is then recovered by tracing out the environment: E(ρ) = Tr env U E (ρ ⊗ |e 0 e 0 | env ) U † E . This implies that E m = e m |U E |e 0 , or, equivalently, where |ψ is any state in the physical Hilbert space.
A standard fact in quantum error correction, sometimes called the decoupling principle, is that there always exists an approximate recovery channel where the error in recovery is bounded in terms of the amount of correlation the error channel generates between the code subspace and the environment. The convenient way to evaluate this correlation is to use the following procedure: (a) introduce a reference system H ref which is isomorphic to and maximally entangled with with code Hilbert space, (b) act with the error quantum channel, (c) trace out the physical Hilbert space space, and (d) evaluate the correlation between the two remaining auxiliary spaces (the environment used to represent the channel and the reference space). Thus, taking |i ref and |i code to be orthonormal bases for the reference space and the code subspace respectively, we construct the state where we have defined d X to be the dimension of a Hilbert space H X . Here the code states are embedded by V into the physical Hilbert space and maximally entangled with the reference, while the error channel acts via E m on the physical states and thus entangles them with the environment. Then we can say that for any error channel E = {E m } there exists a recovery channel R for which the Schatten 1-norm distance between the resulting state and the original encoded state is bounded as [50,51]: Here I Ψ (ref : env) is the mutual information between the environment and the reference space after tracing out the physical Hilbert space. This means that the error E is exactly correctable in the code V if the reference and environment do not share any correlation, hence the term "decoupling".
In this paper, we will be interested in quantum codes with complementary recovery [52], which are the types of codes relevant for entanglement wedge reconstruction in AdS/CFT. For simplicity, consider a code subspace where we have some semi-classical bulk degrees of freedom in the entanglement wedge of a boundary subregion A, but no excitations in the entanglement wedge of the complement region A. Let |i code denote basis states for these bulk degrees of freedom. It was shown by Harlow that the Ryu-Takayanagi formula together with quantum corrections implies the following structure for the encoding map in this situation: where the physical Hilbert space (i.e., the Hilbert space of the dual CFT) is factorized as and |χ is some fixed pure state in the Hilbert space H A 2 ,A . The argument for this involves the decoupling principle applied to the erasure of A. Let us briefly recall how this works (see [52] for details): we introduce an auxiliary system H aux isomorphic to the code subspace, and construct the state Since the bulk degrees of freedom in the code subspace are contained in the entanglement wedge of A, one can show using the RT [1] plus FLM [53] formula that the mutual information I(aux : A) vanishes, which implies that ρ aux,A = ρ aux ⊗ ρ A . Therefore, viewed as a bipartite state on A and aux ∪ A the Schmidt vectors of Ψ should take a factorized form on aux ∪ A.
The canonical purification [36,54] of ρ aux,A will therefore also have factorized states on A. This is why the state inside the parentheses in equation (2.7) takes the factorized form between A 1 and A 2 ; here A 1 is the canonical purifier of aux and A 2 is the canonical purifier of A.
Finally, any two purifications of the same density matrix ρ aux,A should be related by a unitary on A; this is precisely the unitary U A ⊗ 1 A appearing in equation (2.7). It is easy to check that this code subspace is protected against the erasure of A. We will refer to the operator U A as the encoding unitary. 5 There is an important caveat: the bulk-to-boundary map V need not be an exact isometry, and is often an approximate one with corrections of O(e −1/G N ). 6 Relatedly, the quantum generalization of the Ryu-Takayanagi formula, namely the QES formula, is correct to all orders in the G N perturbation theory for appropriate states, 7 but in general there are corrections of 5 The additional Hilbert space component HA 3 in (2.8) is, for our purposes, a bookkeeping device for situations where the physical Hilbert space dimension is not a product of integers, also implying that a part of it does not participate in the code; so we will simply drop HA 3 as it is not pertinent to our considerations. Henceforth, we will focus our attention on codes which have the above structure, but without HA 3 . 6 There are also more extreme situations in which the map is far from an isometry [29,31,32]. O(e −1/G N ). Therefore, the bulk-to-boundary map in AdS/CFT is only approximately of the form (2.7), and has additional exponentially small corrections. As a first pass, we will focus on codes of the type (2.7) in this work. It would be interesting to incorporate the corrections mentioned above in our analysis, but we will not attempt this here.

An error correction bound
Putting together the considerations from above, we first introduce a reference system isomorphic to the code subspace and consider the maximally entangled state: where we included the code subspace structure in (2.7) and an auxiliary environment in some fiducial initial state |e 0 . The error now acts in the form of a joint unitary operator on A ∪ env (we assume the error does not act on the A system): where we applied the error channel as in (2.4) to the state (2.10) in terms of a unitary operator (2.3) entangling the physical system with the environment. Next, we obtain the reduced density on the reference and environment subsystems: where we have performed the trace over A and replaced χ with its reduced density matrix ρ χ A 2 . Finally, following (2.5) we can bound the error in recovery of the original state after action of the error channel in terms of the fourth root of the mutual information between the reference and the environment: The mutual information in (2.13) is difficult to compute directly. A standard approach is to use the replica trick to obtain the von Neumann entropies on the right hand side as analytic continuations of the Rényi entropies which are easier to compute via the relation where S (n) (ρ) = log Tr(ρ n ) is the n th Rényi entropy. We will take a different approach. In Appendix B we study a particular, well-motivated recovery channel and show that the trace distance between the recovered state under this recovery channel and the actual state satisfies where D(ρ, σ) = 1 2 Tr(|ρ − σ|) with |X| = √ X † X is the trace distance between density matrices. As above, d code and d env are dimensions of the code/reference subspace and the environment in the isometric extension of the error channel, respectively. In this work, (2.15) will replace the standard decoupling principle (2.5) due to the ease of evaluating the right hand side. In particular, the expression (2.15) bounds the error in recovery directly in terms of the quantity which we call the mutual purity. If F vanishes, so does the right hand side of the bound (2.15), so that perfect recovery is possible and the error can be corrected. In view of this bound, below we will compute F to quantify the robustness against errors for encoding maps of increasing complexity.

Error correction and maximum complexity encoding
To get a more quantitative understanding of what happens when the encoding unitary becomes complex, as a first pass we can compute the Haar ensemble average with respect to U A of F Ψ . This is because we expect that the typical unitary in the Haar ensemble will be exponentially complex, and that as long as the dimension of H A is large, deviations away from the ensemble average will be exponentially suppressed in the number of qubits. For a Haar random unitary U acting on a Hilbert space H X , a standard formula for Haar integration says: ).
(2.17) This expression has a gravitational analogue: in the Euclidean path integral computation of the radiation purity in the PSSY toy model for an evaporating black hole in JT gravity [10], the two terms displayed above are respectively analogous to the "disconnected" and "wormhole" gravitational saddles.  19) and the associated reduced density matrices (2.20) In the first step, we represented the action of U E in terms of the Kraus operators E m = e m | U E |e 0 and traced out the reference and the environment degrees of freedom. The second step follows from equation (2.17). In the fourth step, we used equations (2.19) and (2.20).
Following similar steps, we can compute the Haar average of Tr ρ 2 env : .

(2.21)
Since ρ ref is maximally mixed, we have where the · · · indicate exponentially small contributions that we have dropped along the way, is the second Rényi entropy of the A subsystem in the state |χ A 2 A . The salient feature of (2.23) is the leading exponential suppression, as we will now describe. The quantity 1 − d −2 code is simply an O(1) prefactor for a nontrivial code subspace.
Two features of (2.23) are worth highlighting. Firstly, note from the final formula that in the typical code drawn from the Haar ensemble, the error channel perceives the state on A as maximally mixed, and gains no access to the microscopic structure of the state. Consequently, as long as the error channel is low-rank, we see that F Ψ (ref : env) Haar is exponentially suppressed by e −S (2) (σ A ) . This is a direct consequence of complexity -a general, complex encoding unitary scrambles the code subspace to a point where generic error channels do not gain any access to it. (A similar coarse-graining picture for apparent horizons and quantum extremal surfaces was advocated in [15,36,56].) Furthermore, there is an additional suppression factor of e −S (2) (χ A ) in equation (2.23) coming from the shared entanglement with A. The combination of these two effects coming from complexity and entanglement thus makes the code robust against generic, low-rank errors.
In the above analysis, we have assumed that the error channel does not have prior access to the encoding unitary U A . This is crucial, because with prior access to the details of the encoding unitary, it is possible to construct low-rank error channels which corrupt the code subspace. For example, consider the error channel: where the unitaries U † A (· · · )U A first undo the encoding, and the partial swap then swaps out the state on the first t qubits with the environment. Since the reference system in Ψ is maximally entangled with the qubits in A 1 , even if the partial SWAP acts on one of the qubits in A 1 , then it will generate an O(1) amount of mutual information between H ref and H env , and thus error correction fails. There is an analogue of this in the JT gravity model [31] -there, one assumes that the error channel does not "generate" additional asymptotic boundaries which can connect up with the bulk geometry and modify the mutual information. Of course, note that this channel is fine-tuned, in that it uses the specific unitary U A which goes into the encoding. Nevertheless, if U A is computationally simple, then the above error channel is also simple. On the other hand, when the encoding unitary U A is exponentially complex, the error channel described above must be equally complex in order to first undo the encoding. Thus, if the python's lunch has an exponentially complex encoding map, then although it will not not be robust against the error channels which are constructed with prior access to the encoding unitary, the channel in question will be exponentially complex. So it will be extremely difficult to implement such errors. This is again a manifestation of the idea that semi-classical causality in the bulk is robust due to complexity.

Random circuit codes: Brownian SYK
Our goal in the rest of the paper is to study in more detail the dependence of the error correction against generic, low-rank errors acting on A relative to the complexity of the encoding unitary. It is convenient, for this purpose, to study the ensemble average of the mutual purity introduced above, but we would like to consider a one-parameter family of ensembles, labelled by the complexity of the typical unitary in the ensemble.
A simple way to generate such an ensemble is to consider the time evolution operators U A = e −iT H , for some ensemble of chaotic Hamiltonians. It is important that the Hamiltonians be chaotic, because for integrable Hamiltonians, the complexity of the time evolution operator is expected to saturate at a sub-exponential time-scale [46]. On the other hand, for chaotic Hamiltonians, it is widely expected that the complexity C(e −iT H ) grows linearly with time T for an exponential amount of time: [44,46,57,58]). Thus, the parameter T is expected to be a good measure of the complexity for chaotic Hamiltonians for exponentially long times. Considering an ensemble of chaotic Hamiltonians then allows us to rely on this property holding only for the typical chaotic Hamiltonian, which is a much weaker assumption than expecting an arbitrary chaotic Hamiltonian to have linearly growing complexity.
More generally, we could consider unitaries U A which are constructed from random circuits. Any unitary can be constructed as a circuit with local quantum gates -in a random circuit, we randomly choose the local gates at each instant of time from some ensemble. The resulting one-parameter family of random circuit ensembles may be thought of as a oneparameter family of measures dµ(T ) on the unitary group U (d A ). To guarantee increasing complexity, we can choose dµ(T ) to be highly concentrated at the identity when T = 0, and as T increases we require that the support of dµ(T ) should expand outward on U (d A ) like a gas, eventually covering the entire group. If we further require that dµ(T ) approaches the Haar measure when T → ∞, we can guarantee that the typical operator selected by averaging with dµ(T ) will have roughly increasing complexity as T increases [59].
To define such a one-parameter ensemble of random circuit codes it is convenient to pick U A (T ) to be the time evolution operator of the Brownian SYK model [37][38][39][40][41]. This model is constructed from N Majorana fermions ψ a , and is defined by a set of random couplings: (2.25) The coupling constants are time-dependent and are chosen to be independently Gaussian at each time point with mean zero and a fixed variance: The associated encoding unitary operator is where T is the time-ordering operator. Note that we are using the Brownian theory not as a model of a holographic boundary theory Hamiltonian (as has been done previously [41]), but rather as the generator of a family of holographic encoding (bulk-to-boundary) maps.
Because H(t) depends on random couplings, U A (T ) is a random variable which has support on certain portions of the unitary group depending on the magnitude of T . The relevant portions are analogous to regions of space covered by a random walk of a certain fixed length.
An subtlety which we will return to later is that the SYK theory obeys certain global symmetries. The presence of these symmetries prevents the effective measure dµ(T ) from covering the entire unitary group as T → ∞. To get around this, we will follow the strategy of [41], where a semi-classical analysis of the SYK theory gave a natural way of extracting results for SYK-like theories which do end up covering the whole unitary group.

Erasure errors
In order to further simplify the problem, we will consider a particular class of errors. It is important that the error channel has no prior access to the encoding unitary, i.e., we want the error to be generic and low-rank. The error channel we consider in this work will be the erasure of some subsystem R.
and let σ R be the reduced density matrix on R, and σ L be the corresponding reduced density matrix on the rest of the system L. We have defined this new σ, which we will use for the rest of this paper, in place of the previous one in (2.19). Then, for such an erasure error, where recall that |Ψ was the state which resulted from applying the error to the maximally entangled state between the reference and the code subspace. For simplicity, we will specialize to the case where R = A 1 and L = A 2 ∪A. In addition, using the fact that |χ A 2 A is maximally entangled, we then arrive at (2.30) Since the dimension of the environment in this case is the same as d code , we find that the robustness of error correction for the above erasure error is bounded by the quantity In the next section, we will turn to the main objective of this work: computing the purity Tr σ 2 L for Brownian SYK codes. In particular, we are interested in the dependence of this quantity on the encoding complexity of the code, which as explained above, is linearly related to the time parameter T . From equation (2.29), we need to compute Tr σ 2 L as a function of T . Actually, since A has no dynamics associated with it (i.e., there is no non-trivial time evolution operator acting on A), this computation boils down to a Lorentzian path integral entirely in the A 1 A 2 subsystem -the relevant time contours are shown in Figure 2. To arrive at Figure 2, we notice that Tr σ 2 L involves two copies of U A and two copies of U † A , and so can be thought of as a matrix element of the operator The matrix element in question is determined by the trace structure: since the R = A 1 system is traced out first to obtain σ L , the adjacent blue A 1 contours are joined in Figure 2, while the secondary trace over L joins the inner and outer red A 2 contours.
When T is small, we expect Tr σ 2 L to be close to 1, and so the mutual purity is non-zero. On the other hand, at very late times, we expect Tr σ 2 L to approach 1/d 2 code and the mutual purity to approach zero. The intuitive argument for this is as follows: let us first purify the density matrix σ by including an auxiliary system aux which is isomorphic to the code subspace: L , the quantity relevant for the mutual purity for the erasure of A 1 . The red contour corresponds to the A 2 fermions while the blue corresponds to the A 1 fermions. The hatched regions denote an application of the time evolution operator U A (T ) or U † A (T ) of Brownian SYK, which couples the A 1 and A 2 systems. We have omitted the contour orientations which determine the forward and backward time evolution, but from left to right the hatched regions alternate between U † A (T ) and U A (T ), beginning with U † A (T ). The arcs at the top and the bottom specify the final and initial conditions respectively; in our calculation, all these arcs are actually infinitesimally small (corresponding to maximal entanglement), but they have been enlarged for visual clarity.
When T = 0, the subsystem A 1 is maximally entangled with aux while L = A 2 ∪ A is in a pure state. When T becomes large (on the order of the scrambling time), we expect U A (T ) to generate nearly maximal entanglement between A 1 and L. By the monogamy of entanglement, therefore, A 1 cannot share much entanglement with aux. However, the unitary operator never acted on aux; thus the reduced density matrix on aux must still be maximally mixed. We therefore conclude that both A 1 and aux are close to being in a maximally entangled state with L, and so the purity of L must approach 1 In what follows, we wish to understand the detailed time-dependence of the mutual purity at late times. In particular, we will demonstrate that the mutual purity becomes O( 1 N ) by the scrambling time T ∼ 1 J log N , but then continues to decay thereafter, approaching its saturation value which is O(e −N ) at a time of order T ∼ 1 J N . The important point is that the mutual purity keeps decaying even beyond the scrambling time, until it reaches an exponentially small plateau, which in the present model happens at polynomial time. 8 We will interpret this phenomenon as "complexity-protected quantum error correction".

Boundary conditions and large N equations
Our task now is to evaluate the path integral of Brownian SYK on the contour in Figure 2. Following [41], we will use the collective-variable description of Brownian SYK. We view the path integral in Figure 2 as an amplitude where we start with an "in" state, then time evolve for a time T and then take the overlap with an "out" state. The boundary conditions relevant for us are as follows. For the in boundary conditions, we have while for the out state we have the adjoint boundary conditions: Here a 1 denotes the index of the N 1 fermions corresponding to the subsystem A 1 , while a 2 denotes the index of the N 2 fermions corresponding to the subsystem A 2 . The superscript index (i) on ψ (i) a (where i = 1, · · · , 4) labels the four contour segments corresponding to real time evolution. From left to right in Figure 2, we label the contours 1, 2, 3, and 4.
In order to evaluate the path integral, it is convenient to define the two matrices: which we can think of as the singlet part of the fermion two-point functions in the A 1 and A 2 sectors respectively. Here we inserted the operators on the right hand side at the specified time into the path integral in Figure 2. We will soon see that these two sets of variables control the classical limit of the Brownian theory on this contour. To solve the classical equations of motion we will obtain in this limit, we require the boundary conditions that are implied by the in and out state relations above. It is also convenient to define the total two-point function (i.e., the summed two-point function of all the fermions): where we have introduced the parameter λ = N 1 N . When evaluating the path integral at large N , it will be convenient to take the double scaling limit: In fact, A 1 has the same dimension as the code subspace, so we would like to take N 1 much smaller than N 2 . This corresponds to taking λ 1. Thus, we can take λ to be a small (but O(N 0 )) parameter and work in perturbation theory in λ. This makes some of the path integral calculations analytically tractable.
Note that both the in and out boundary conditions have the property that (recall that the fermions are normalized such that ψ 2 = 1/2): Since A 1 fermions lie in the same parity sector as the A 2 fermions, and the (effective) Hamiltonian commutes with the fermion parity operator after averaging (see Appendix A, equation (A.2) and discussion below it), the above relations should hold at any time. Equations (3.8) and (3.9) imply the following symmetry properties: 34 , g 23 , g 13 , (3.10) 34 , g 23 , g 13 . (3.11) These should hold at all times because time evolution preserves fermion parity flavor-wise. So the evolution reduces to the six variables x α = 2ig x 2 (0) = −y 2 (0), z 2 (0) = 1. (3.13) x 1 (T ) = 1, y 1 (T ) = −z 1 (T ), (3.14) x 2 (T ) = y 2 (T ), z 2 (T ) = 1. (3.15) It is convenient to also introduce the total variables x = λx 1 + (1 − λ)x 2 , and similarly for y and z. The above boundary conditions imply the following constraints in terms of the (x, y, z) variables: 2λ). (3.16) In order to proceed with the evaluation of the Lorentzian path integral in Figure 2, recall [41] that the action for the Brownian SYK model on the contour in Figure 2 is given by where the flavor indices run over a = 1, · · · , N (i.e., over both A 1 as well as A 2 fermions), and we have introduced the notation ψ (j) aq . The quantity s j is given by 18) and is related to the difference between forward and backward time evolution (see [41] for details). We now wish to perform the average over the couplings. Using the action obtained after ensemble averaging over the couplings is given by 9 (3.20) At this stage, it is convenient to introduce the collective (G, Σ) variables. Since we have two sets of fermions corresponding to A 1 and A 2 , we introduce two collective fields 21) and the corresponding Lagrange multipliers Σ ij (t, t ) and Σ (2) ij (t, t ) to impose the constraints. We can now integrate out the fermions. The action in terms of the collective variables is

22)
9 In the second step, we have made the same imprecise replacement of the Hamiltonian as in [41], discussed in more detail in Appendix A.3 of [40].
where Pf is the Pfaffian, and we have defined ij (t, t ). (3.23) In the large N limit, the path integral over the collective variables can be performed in the saddle point approximation. The equations of motion corresponding to the above action are: where α = 1 corresponds to the fermions in A 1 while α = 2 corresponds to the fermions in A 2 , the repeated k index is summed, and the star product between two bi-local fields is defined as Using the fact that G (α) and Σ (α) are both anti-symmetric, we can rewrite these equations in a more convenient form: Now, a simplification happens in the Brownian SYK model -recall that for Brownian SYK, J 2 (t, t ) = Jδ(t − t ). As a result, Σ is "diagonal" (in time), and only the diagonal components of all the collective variables are relevant; the off-diagonal components drop out of the equations of motion. In fact, it is easy to see from the action that for Brownian SYK, the off-diagonal modes do not have any interesting dynamics and can be integrated out of the full path integral trivially [41].
Let us denote the diagonal components of the collective variables as The resulting equations of motion for the (g, σ) variables are dg dt where the equation on the left is written for the g and σ matrices and As an aside, the above equations of motion have a Hamiltonian structure. To see this, let us denote x α = (x α , y α , z α ) and . Then, equations (3.32) take the succinct formẋ where p α = ( 1 λ , 1 1−λ ). Similarly, the equations for the total variables take the forṁ Thus, these equations take the form of Hamilton's equations of motion -the underlying phase space is that of two copies, labelled by α, of a co-adjoint orbit of sl(2, R) specified by a constant value of the conserved quantity h 2 (x). 10 The Hamiltonian h q (x) couples the two copies, with the effective coupling constant being λ. It may seem unusual that we have an odd number of variables (e.g. (x, y, z)) in Hamiltonian mechanics, but this is simply because we have parametrized the two-dimensional dynamics on the hypersurface h 2 (x) = const. in terms of coordinates in the ambient R 3 .

Normalization
The fermionic path integral depicted in Figure 2 is the result of writing Tr σ 2 L in a Hilbert space form and replacing the maximally entangled state projectors with connections in the contour. However, translating a Hilbert space expression into a fermionic path integral comes with a standard normalization issue since Majorana fermions have a somewhat unusual Hilbert space interpretation (when they admit one at all). So, we must relate the path integral Z(T ) to Tr σ 2 L (T ) with an overall normalization that ensures Tr σ 2 L (T = 0) = 1.
When T = 0, the contour in Figure 2 consists of four disconnected circles which are not coupled by any time evolution. The result of the path integral for a single free Majorana fermion on a circular contour of length T with antiperiodic boundary conditions is actually equal to √ 2, independent of T , so in the limit T → 0 we still have √ 2. As such, when T = 0, the contour in

Solutions: qualitative discussion
We will first qualitatively discuss what the solutions to the equations of motion (3.32) should look like, leaving a quantitative treatment for Sections 3.3, 3.4, and 3.5. When N 1 = 0 (i.e., λ = 0), then the x 2 equations are easy to solve. In this case, the boundary conditions imply that the solution stays at the fixed point x * 2 = (0, 0, 1). When λ is small but non-zero, we expect that this saddle point remains, but with small corrections. In particular, the x 2 variables will stay close to their original fixed point values. The corrections to the x 2 solutions can be obtained in perturbation theory in λ, and we describe them in detail in Section 3.3. (Recall from the discussion under (3.7) that λ 1.) This resulting solution is the dominant saddle point at small times, and is the analogue of the "disconnected" contribution in equation (2.17), or the disconnected geometry in JT gravity [31].
When λ is small, to zeroth order, the solution for x 2 variables will be unaffected by the x 1 variables, and in particular will correspond to a fixed point of the Hamiltonian picked out by the boundary conditions as we have just described. However, the initial backreaction on the x 1 variables will be large. The source of this strong backreaction is the mismatch between the boundary conditions of the x 1 and x 2 variables. As the A 1 system is small compared to A 2 , the Brownian dynamics quickly thermalizes the A 1 system so that the correlation between contours (in Figure 2) can be measured with any subset of all N fermions; at the level of the solutions, this means that we will have x 1 ≈ x 2 at all times except in small neighborhoods around t = 0 and t = T where we expect large transient behaviors for x 1 to arrive at their "thermalized" values. These transient behaviors can be computed analytically in perturbation theory for the disconnected solution and must be treated numerically otherwise. These qualitative properties of the x 1 solutions hold both for the disconnected solution in Section 3.3 as well as the other solutions we now describe.
In addition to generating nontrivial time-dependence for the disconnected solution, turning on a small λ has another important effect -it gives rise to new "tunneling" solutions (i.e., instantons) which are absent at λ = 0. The tunneling allows the x 2 variables to jump between different fixed points; the leading tunneling solution jumps from x 2 (0) ≈ (1, −1, 1) to x 2 (T ) ≈ (1, 1, 1). The x 1 variables again have large transient signals near t = 0 and t = T , but this time both transients are different than the ones which occur for the disconnected solution due to the different fixed points approached by the x 2 variables, and in fact these are the only other two types of transient behavior which can occur. This saddle point, which we describe in Section 3.4, is the analogue of the "connected" saddle point in equation (2.17), or the "wormhole" in JT gravity [31]. While it is suppressed by a factor of e −(1−2λ)N relative to the leading, disconnected solution, the contribution of the disconnected saddle point decays exponentially in time. So, at a time t * ∼ O(N ), there is an exchange of dominance between these two saddle points.
There are also other tunneling saddle points, described in Section 3.5, where the x 2 variables tunnel back and forth multiple times between the two possible initial fixed points x 2 ≈ (0, 0, 1) and x 2 ≈ (1, −1, 1) and the two possible final fixed points x 2 ≈ (0, 0, 1) and x 2 ≈ (1, 1, 1); these are even more subleading in powers of e −N , and occur with all possible combinations of the previously described types of transient behaviors for the x 1 variables. Explicitly, there are two possible behaviors at t = 0 and two at t = T corresponding to the possible initial and final fixed points for x 2 , and all four combinations of initial and final transient behaviors occur in the multiply tunneling solutions. These multiply tunneling solutions show interesting behavior as a function of T . We will see that they become genuine solutions of the equations (3.32) only after certain critical values of T , related to the scrambling time. Before these critical times, these configurations are actually off-shell. Configurations which tunnel more times take longer to become solutions. As the presence of these contributions is important for unitarity of the overall evolution [41], it is intriguing that they can be invisible on-shell for a parametrically (though not polynomially) long time in N .
In summary, we began with the goal of studying the error correction dynamics of a family of unitary operators with increasing average circuit complexity. The specific family which we chose for convenience was the set of time evolution operators in the Brownian SYK model, a family of time-dependent Hamiltonians which are essentially a continuous random circuit. We found that the error correction dynamics are governed by the quantity F Ψ (ref : env), and this quantity in turn depends on a Brownian SYK path integral (Figure 2). What we have just discussed are the saddle point solutions to that path integral. Evaluating the effective action of these solutions will allow us to draw conclusions about the error correction behavior of the family of unitary operators with increasing complexity.

Disconnected solution
We will first solve for the disconnected solution at small, non-zero λ, and evaluate its on-shell action together with the one-loop determinant. We begin by expanding our variables in a power series expansion in λ: 2 . Therefore, the boundary conditions, equations (3.14) and (3.15), imply that at leading order these variables sit at a fixed point of the Hamiltonian: After substituting these solutions in (3.32), we get the following equations for x 1 : We need to solve these equations with the boundary conditions (3.12) and (3.14). The solution is For instance, x 1 starts off at one at t = 0, but after a brief transient behavior, it thermalizes to a small value of about e − JT 2 q−1 owing to its coupling to the x 2 variables, which act like a bath and dynamically force x 1 ≈ x 2 . There is another transient near t = T , where x 1 again deviates from x 2 significantly to reach the final boundary condition.
At order λ, we note from x = λx 1 + (1 − λ)x 2 , that (3.41) Now we use the O(λ 0 ) solutions to find the following boundary conditions up to O(λ): x (1) With these boundary conditions in hand, we can in principle solve for all the variables at O(λ). However, in what follows, we will only need x (1) in order to evaluate the on-shell action up to O(λ). The corresponding differential equations are given bẏ The solution is: (3.46)

Classical on-shell action
Having obtained the classical solutions, we can evaluate the action in (3.22) at leading order in λ. We first compute the Pfaffian for the A 1 fermions: On the right hand side of the above expression, we have used the fact that on-shell, σ ij (t) = σ ij (t); see equation (3.30). Following [41], we can write the Pfaffian in a Hilbert space representation in terms of a single qubit: 11 where h(t) is the qubit Hamiltonian: 12 and X, Y, Z are the Pauli matrices. Further, |+ is an eigenstate of the Pauli X operator with the eigenvalue +1: X|+ = |+ . The initial and final states are fixed by the boundary conditions of the path integral. Since the leading contribution to x q−1 (t), y q−1 (t) is at O(λ q−1 ) we can ignore them in the evaluation of the Pfaffian (assuming q ≥ 4). The Pfaffian is therefore given by (3.50) In the second step, we used the fact that z(t) is a constant at O(λ). Similarly, we can evaluate the Pfaffian for the A 2 fermions: (3.51) 11 Naively, we would need two qubits given that there are four Majorana fermions. However, since time evolution preserves fermion number, we need only use one qubit. 12 The factor of 2 appearing in (3.48) ensures that, when T = 0, the Pfaffian gives 2, as this is the result for a single Majorana fermion path integral on two disjoint circles of any length with antiperiodic boundary conditions.
Having evaluated the Pfaffians to O(λ), we now evaluate the rest of the terms in the action following [41]: is a constant of motion, namely the total energy. Note that this form of the bulk contribution to the action is independent of the particular solution we are considering. Now, we can combine the above terms with equations (3.50) and (3.51) to obtain the full on-shell action for the disconnected saddle point: As the normalization relation (3.36) cancels the leading log 2 in the effective action, the contribution of the disconnected saddle point in the large N limit is given by This formula is consistent with physical expectations; see the discussion around equation (2.31). At small times JT 2 q−1 1, we find that Tr σ 2 L → 1. This is expected, since in this limit, U A (T ) does not introduce much entanglement between L and A 1 . On the other hand, at late times JT , as we anticipated based on monogamy of entanglement. As a further check, we also reproduce the above formula from a "Hamiltonian" point of view in Appendix A. For now, we proceed to evaluate the one-loop determinant around the disconnected solution. But, before doing so, we display the numerical solutions for x 1 and x 2 in Figure 3, where we can clearly see the leading order nontrivial time-dependence of x 1 takes the rough hyperbolic forms we found for x

One-loop determinant
To compute the one-loop determinant in the path integral formalism, we need to expand the action around the saddle point and integrate over small fluctuations. We will follow the notations and conventions of Appendix B in [41]. Recall that the action is given by ij (t)g (2) ij (t) (3.55) We can write the Pfaffian in the Hilbert space representation, as in equations (3.50) and (3.51):

(3.56)
We now expand around the saddle point solution x (α) * found in Section 3.3. We will use hatted variables to denote fluctuations: . The quadratic action for the fluctuations has the following form: (3.60) The matrices K 1 and K 2 can be derived by variation of the Pfaffian terms at quadratic order. S andS are defined as (3.61) To compute the determinant at leading order in λ, it turns out to be sufficient to note that K 2 has non-zero matrix elements only inσ (2) x andσ (2) y . Therefore, K 2 satisfies the relatioñ Moreover, we will only need the (σ We can now compute the determinant of M to leading order in λ: We first compute the trace term in the determinant In the first step, we used the fact that B 22 is the only non-zero entry in B. In the second step, we inserted the A −1 22 and C −1 22 components upto O(λ) corrections. The third and fourth step follow from the relation (3.62). Using relation (3.62) once again, we conclude that det C = 1. We are left with the evaluation of det A.
So, the one-loop determinant does not significantly modify the T dependence of Tr σ 2 L at leading order. The coefficient of the O(λ) term above is bounded by an O(q) number, and qλ is always small in our regime of interest.

Connected solution
When λ is slightly non-zero, x 2 ≈ (0, 0, 1) is not the only fixed point for the x 2 variables which enters the analysis. The leading solution involving more than one fixed point is the tunneling solution between the x 2 ≈ (1, −1, 1) and x 2 ≈ (1, 1, 1) fixed points. This solution has nontrivial time-dependence for the x 2 variables which involves an initial region where x 2 (t) ≈ (1, −1, 1), then a transition to the x 2 (t) ≈ (1, 0, 0) fixed point where the solution remains for a long period, and then a final transition to the x 2 (t) ≈ (1, 1, 1) fixed point. The x 1 solution has large transient behaviors in the initial and final fixed point regions, but matches very closely with x 2 in the long region where x 2 ≈ x 1 ≈ (1, 0, 0).
Because this solution is non-perturbative in λ, we cannot hope to use perturbation theory to evaluate the effective action. We will instead follow the approximate analysis of [41]. Unlike the disconnected solution we described in Section 3.3, the connected solution is suppressed exponentially in N , and the main aim of our approximate analysis will be to demonstrate this suppression quantitatively. The numerical connected solution is shown in Figure 4. We will give an approximate analytical computation of the action for this solution. As argued in Section 3.3, for any solution of the equations of motion the bulk terms in the effective action contribute to the total path integral exp(−JT /2 q−1 ) for r ≈ 1. So what remains is to evaluate the two Pfaffian contributions, again using the qubit Hamiltonian approach.
We see in Figure 4 that there is a long region with x 1 ≈ x 2 ≈ (1, 0, 0). In this region, we may approximate the time-ordered exponential expressions as projectors |+ +|, the lowest energy state of the Hamiltonian −JT X/2 q−1 generating the time evolution in that region. 13 The energy contribution from this ground state exactly cancels the bulk term, so the result of the long region for both x 1 and x 2 is a projector |+ +|. At this point, the analysis splits between x 1 and x 2 . The x 2 variables include an initial region around the (1, −1, 1) fixed point and a final region around (1, 1, 1). Both of these regions share an important property with their adjacent transition regions, namely that the first has y 2 ≈ −z 2 and the second has y 2 ≈ z 2 . Because these regions are adjacent to the long middle region that yields a projector |+ +|, we may use the null state relations +|(iY + Z) = 0 and (iY − Z)|+ = 0 to conclude that the Pfaffian 13 When T is on the order of (1/J) log N and not much larger, there are exponentially suppressed T dependent corrections to this projector which lead to O(1) factors in the wormhole contribution to the Rényi mutual information. While these corrections could be addressed in the path integral formalism we employ here, it is easier to study them in the Hamiltonian picture (Appendix A). We will continue to approximate the long region as a projector because these corrections are highly subleading by the time the wormhole dominates at T ∼ N/J and are therefore unimportant for the qualitative error correction properties of the Brownian circuit.
does not depend on the precise details of these transition regions nor on the initial and final fixed point regions, and the remaining X term in the Hamiltonian simply cancels against the bulk contribution as was the case in the long middle fixed point region. So, the overall Pfaffian for the x 2 variables is determined by the overlap of the initial state |0 and final state 0| with the projector from the middle region: 2 0|+ +|0 = 1.
We may analyze the x 1 variables similarly. Because the long region where x 1 ≈ x 2 again gives a projector |+ +|, and because the transient regions satisfy the same relations between the variables as the transition regions from the x 2 analysis, the same arguments we made about the transition regions for x 2 goes through for the transient behaviors of x 1 , and the Pfaffian does not depend on the precise form of the transient behaviors. The remaining Hamiltonian contribution from X again cancels the bulk term in the transient regions. There are no initial or final fixed point regions for x 1 , so the total contribution is from another projector overlap with the relevant initial and final states: 2 +|+ +|+ = 2. It may seem redundant to analyze the x 1 variables separately as we have done here, since the A 1 Pfaffian term involves the total x variables like the A 2 Pfaffian. However, it was important here to conclude that the transient behaviors do not contribute any time-dependence at O(λ), and we actually obtained a constant result that is independent of T .
Thus, again using the normalization (3.36), from the connected solution we have a contribution Tr (3.69)

Other tunneling solutions
There are also solutions which tunnel from x 2 ≈ (0, 0, 1) to x 2 ≈ (1, 1, 1) and from x 2 ≈ (1, −1, 1) to (0, 0, 1). We name these the "DW" and "WD" solutions, respectively, after the order of transient behavior which occurs for the x 1 variable: the first has "Disk" initial transient behavior and "Wormhole" final transient behavior, while the second has the opposite ordering. The DW solution is shown in Figure 5 while the WD solution is shown in Figure 6. The contribution of these solutions to Tr σ 2 L can be evaluated in the same approximate manner as Section 3.4.
We begin with the DW solution in Figure 5. The x 2 variables has the same long region with x 2 ≈ (1, 0, 0) which appears in the connected solution (Figure 4), and by the same reasoning as in Section 3.4 we conclude that this region yields for the path integral the projector |+ +|. Similarly, the long region with x 2 ≈ (0, 0, 1) gives a projector |0 0|. These two projectors cancel the transition regions and the other constant regions associated with other fixed  points, and we get the overlap 2 0|0 0|+ +|0 = 1. The x 1 variables have a leading contribution determined by simply changing the initial and final states: 2 +|0 0|+ +|+ = 1. Thus, including the normalization (3.36), we have the additional suppression 2 −N for the DW solution: 14 Tr σ 2 We will not bother to compute the O(λ) contribution from the transient behaviors of x 1 in the A 1 Pfaffian (which could lead to nontrivial T dependence), since this solution is already highly suppressed compared to the connected one in Section 3.4.
The WD solution can be analyzed similarly and also has a 2 −N leading suppression. Thus, both the DW and WD solutions are subleading compared to the connected solution from Section 3.4. There are also even more highly suppressed solutions which can be formed by inserting additional periods into any of the four solutions we have discussed up to this point. A full period of the x 2 variables is shown in Figure 7. By the same approximate reasoning, inserting a full period in the solution will suppress the contribution to the path integral by an additional 2 −N .
Interestingly, the long regions of the solution have a minimum length which scales like the scrambling time T s ∼ (1/J) log N . What this means is that they are actually not solutions for all values of T . For instance, the connected solution in Section 3.4 is only a solution for T > (1/J) log N . A configuration with k long regions will not appear as a solution until T > (k/J) log N . This lattice of critical times is interesting from a unitarity perspective.
These subleading saddles are necessary to ensure the total Brownian evolution is unitary, so an experimentalist with access to only on-shell configurations will discover that it is impossible to verify unitarity with accuracy better than 2 −kN until at least T > (k/J) log N .

Summary
We have shown that the leading T dependence of the purity which controls the mutual purity in Brownian SYK is 15 where the first term comes from the disconnected saddle point, the second term from the leading connected saddle point, and the dots represent further subleading solutions that are suppressed in powers of 2 −N . We present a comparison of this saddle point analysis with an exact numerical computation of Tr σ 2 L in Figure 8. The form of (3.71) means F Ψ (ref : env) is initially O(1) and subsequently decays for a polynomial T ∼ N/J amount of time. When T > N/J, the connected solution begins to dominate and leads to an exponentially small 15 In this analysis, we have purposefully ignored the presence of discrete symmetries. The presence of such symmetries generically prevents the time evolution from covering the entire unitary group. Following [41], we can adapt the analysis of Brownian SYK so that the time evolution covers the entire unitary group by only including the saddle points we have discussed. Incorporating the discrete symmetries of the SYK model requires additional saddle points [41]. This means our results are effectively valid for an SYK-like model with no discrete symmetries which does end up covering the whole unitary group. mutual purity. 16 Thus, when the encoding complexity is sufficiently large, the code is robustly protected from the erasure of A 1 . At multiples of the scrambling time T s ∼ (1/J) log N , subleading contributions become genuine on-shell solutions of the equations of motion, though these contributions never dominate F Ψ (ref : env).

Discussion
We have studied the error correction properties of Brownian SYK quantum codes against the erasure of a small number of qubits, but we expect our results to be valid more generally for generic, low-rank errors with no prior access to the encoding map. As a measure of quantum error correction, we computed the mutual purity F Ψ (ref : env), which is related to the purity Tr σ 2 L , where σ = V V † /d code is the density matrix built from the encoding map V , σ L = Tr R σ, and R is a small fraction of the physical Hilbert space which is being erased. In codes defined using Brownian SYK time evolution, which have a linearly growing encoding complexitymimicking the expected behavior of the bulk-to-boundary map for an infalling observer in AdS/CFT -this purity is related to a four-contour Lorentzian (Schwinger-Keldysh) path integral. We found two special saddle point solutions to the large N equations of motion in Brownian SYK -analogous to the disconnected disks and connected wormhole geometries in JT gravity -which dominate this path integral. At early times T N/J, the disconnected solution gives an exponentially decaying value for the mutual purity, while at late times the connected solution dominates and gives a constant, exponentially small mutual purity. Thus, when the encoding complexity is sufficiently large, we find emergent, "complexity-protected" quantum error correction against generic, low-rank errors with no prior access to the encoding map. We should emphasize that it is important that the error does not have access to the encoding map -with prior access, it is possible to violate the above conclusions.

Relation to previous work
Understanding how the complexity of an encoding operator affects certain error correction properties of the code is a problem that has been explored previously from a variety of viewpoints. The most common method of studying codes with increasing complexity is to employ the randomization trick as we have done, where one instead considers a one-parameter family of ensembles of codes with increasing complexity and studies ensemble-averaged properties. 16 Recall that that Tr σ 2 L enters in F Ψ (ref : env) along with a subtraction of a baseline value, and so the contribution which dominates F Ψ (ref : env) is not necessarily the one which makes the largest contribution to Tr σ 2 L .
For instance, [60,61] argued that n-qubit random quantum circuits with O(n log 2 n) twoqubit gates and O(log 3 n) depth can encode k qubits into n while correcting erasure errors on d qubits where with h(x) being the binary entropy function In our analysis, we studied a random error with k/n = d/n = λ, and with these replacements the inequality (4.1) is true for roughly λ ≤ 1/5. It would be interesting to understand whether our analytics can be extended to this rather large value of λ without the need for novel techniques, although the numerical results in Figure 8 suggest we may have an accurate picture of the Brownian theory even when λ = 1/3. At any rate, it appears that the Brownian codes we have studied in this work are able to approximately (with error of order 1/n) correct errors on a fraction λ of the physical qubits with a depth T ∼ (1/J) log n. This polynomial improvement in depth, if true, is likely due to differences in how the random two-qubit quantum circuit theory of [61] and the Brownian SYK theory scramble quantum information.
More recently, [62] studied low depth random circuits with spatial connectivity restrictions in various spatial dimensions D as stabilizer codes. They discovered that such circuits can correct fairly large erasure errors (converging to both the optimal threshold and zero failure probability at large n) with a depth of just O(log n) for D ≥ 2. These results are similar to ours, although we have no restriction on spatial connectivity, but rather a restriction on the number of fermions which can couple in the Hamiltonian. It would be interesting to understand if there is a relation between the universality for D ≥ 2 found in [62] and the expected universality of our results for q ≥ 4. A significant difference of our analysis compared with [62] is that we do not restrict ourselves to stabilizer codes, though we also have not studied the decoding problem in any detail.
Beyond questions of depth, we may also consider the total gate complexity of efficient quantum codes. Several bounds on this complexity exist for stabilizer codes [63][64][65] and their generalizations [65,66]. In particular, for a generic stabilizer code encoding k qubits into n, [65] showed that O(n(n − k)/ log n) gates are sufficient. Entanglement-assisted stabilizer codes were also studied in [65] and were shown to have gate complexity linear in the number of additional entangled qubits c, with O(n(n − k + c)/ log n) gates. As we have not restricted ourselves to stabilizer circuits, our gate complexity is not expected to have such small polynomial asymptotic behavior. 17 However, if we used a sparse SYK model instead [69], we may achieve equal or better gate complexity compared to stabilizer circuits. This issue deserves further study as it would represent an interesting development in efficient random code design.
Our work is also closely related to measurement-induced phase transitions which have recently been studied extensively in the condensed matter community (see for instance [70,71]). In these studies, the quantum circuit usually consists of local unitary gates with some quenched disorder and forms a brickwall pattern. These local unitaries are interspersed with local measurements, which are viewed as "errors". The long range entanglement generated by the random unitary gates is identified as the volume-law phase, suitable for quantum error correction. However, a transition to a short-range entanglement phase can occur when measurement rate is high, i.e. when the error rate is high enough to disentangle different subsystems. The volume-law to area-law transition is identified as a transition in quantum error correction, when the error rate exceeds a critical value [70,71]. In essence, the size of the Hilbert space of the principal quantum system needs to be large enough (spatial) and the time for the unitary gates need be long enough (temporal) to scramble the information so that the entanglement is robust against local disturbances. It would be interesting to compare these results with those presented here.
In another direction, ensembles of encoding maps that satisfy some global symmetry have also been explored [72,73]. The general idea is that there is a tension between the existence of a continuous symmetry leaving the encoding map invariant and strong protection against erasure errors. However, approximate error correction can be achieved in certain circumstances [74]. In Brownian SYK, there are discrete global symmetries (which we did not include in the analysis since we were interested in covering the entire unitary group) but no continuous symmetries, allowing us to avoid these no-go arguments. However, it is easy to implement continuous symmetries in analogues of the SYK model; for instance, SYK with complex fermions satisfies a U (1) global symmetry [75]. It would be interesting to understand the error correction behavior of a complex analogue of Brownian SYK to further elucidate the tension between codes with continuous symmetries and erasure error correction.

Pseudorandom codes
Our results seem to suggest that after a polynomial time, a random quantum circuit, which likely has polynomial circuit complexity, has powerful error correction properties that are essentially as good as a Haar random unitary code, which likely has exponential complexity. relation δ(t − t ) in the variance, along with a sparse query model like the one studied in [67,68]. Because the sparsity of the full SYK Hamiltonian scales with N q , we do not expect simulation to be efficient compared to stabilizer circuits.
One explanation for why this is possible may be that the majority of unitary operators with polynomial complexity are in fact pseudorandom unitary operators, and a simple test of error correction properties cannot distinguish polynomially complex pseudorandom unitary operators from unitary operators of exponential complexity.
A pseudorandom unitary operator is, roughly speaking, an operator which has polynomial complexity but which cannot be distinguished from one with exponential complexity by any sort of simple test which can be implemented efficiently. The transition between disconnected and connected solutions that we found, hints at a sharp transition point where most random circuits with complexity less than some polynomial (N , for the purity transition) are not pseudorandom, while the typical circuit and perhaps the majority of circuits above that critical complexity are in fact pseudorandom, at least for the purposes of error correction. It would be very interesting to understand in more detail what properties of Haar random circuits can be reproduced by such low complexity Brownian circuits.

Complexity and the geometry of the entanglement wedge
From the AdS/CFT point of view, it would be very interesting to understand the bulk significance of our results; indeed, one of our main motivations in this paper was to understand the geometry of Figure 1 in terms of quantum error correction. Following [30,31], we expect that this error correction is a sign of "causal inaccessibility" from the boundary subregion. By this, we mean that including backreaction from turning on simple sources in the asymptotic boundary does not render the relevant degrees of freedom causally accessible from the boundary; the mechanism behind this is that the relevant bulk degrees of freedom lie behind a non-minimal quantum extremal surface. In our calculation, we encountered two significant complexity scales, i.e., the mutual purity becomes O(1/N ) at T ∼ log N , and the mutual purity saturates to an exponentially small plateau at a much larger time-scale. It is tempting to speculate that these thresholds have natural bulk interpretations: the log N time-scale could correspond to the bulk degrees of freedom crossing the causal horizon, while the plateau could correspond to the bulk degrees of freedom crossing over to the python's lunch. In a similar vein, the lattice of subleading solutions we found may also have a geometric meaning, although it is less clear because they do not dominate the calculation of the crucial quantity F Ψ (ref : env).
where the initial and final states are dictated by the boundary conditions in the path integral, and are given by

A.1 Disconnected and connected solutions
In the above expression, we can separate contributions from the ground states and the excited states as follows: Here, g n are all the ground states of H eff (which all have zero energy) and e k are the excited states with energies E k . We first look at the contribution from the ground states. The set of ground states depends on whether we choose q = 4k or q = 4k + 2 but the following two ground states contribute to the leading order independently of q: Therefore, the contribution from the ground states is which reproduces the two leading order terms (i.e., the disconnected and the connected contributions) in the Haar ensemble. The contribution from excitations near the ground states can be approximated in the following manner. The Hamiltonian can be written in terms of the ladder operator as where S x = a X a /2, S y = a Y a /2, and S z = a Z a /2. The matrix elements of the first two terms in the Hamiltonian with excited states near |g 1 are suppressed by a factor of 1/N q/2 and can be ignored for q ≥ 4 at leading order in 1/N . Thus, the Hamiltonian up to O(1/N ) corrections is The contribution from the states near |g 1 is: where |e g 1 k denotes the k th excited state near |g 1 , explicitly given by a choice of k qubits which are flipped to |1 from |0 . These k must come from the first N 1 qubits to give a nonvanishing overlap ψ|e g 1 k .
Similarly, we can compute the correction due to the excited states |e g 2 k near |g 2 . The perturbative Hamiltonian is now and the excited states are formed by flipping k qubits to |− from |+ , where these must come from the last N 2 qubits to give a nonvanishing overlap with the boundary state. We get Note that the contribution from the second and higher excited states (denoted here by ellipsis) is not negligible. Moreover, unlike the case of the disconnected saddle where N 1 /N = λ was a small parameter, we cannot resum all the contributions from higher excited states near |g 2 . Since N 2 /N ∼ 1, one must also take the quantum corrections into account. Nevertheless, the above expression is sufficient to infer that the ground state contribution dominates when T > 2 q−2 J log(N 2 ).
Summing up these contributions we have the following result for Tr σ 2 L : The T dependent term proportional to N 2 /2 N 2 is a contribution from corrections to the projector approximation to the long region we made in Section 3.4. We could have incorporated such terms in the path integral saddle point approximation of Section 3.4 by writing the long region as a projector |+ +| plus an exponentially suppressed correction e −JT |− −|. However, the evaluation of the saddle point including this correction is difficult because the transient regions no longer cancel against the |− −| operator, so this term induces large corrections which depend sensitively on the transient shape. This is the path integral analogue of the Hamiltonian picture difficulty we described under (A.14). Of course, these corrections are only important before T < (1/J) log N , when the connected configuration is not actually a solution of the equations. By the time the connected configuration becomes a genuine saddle point, this T dependence is subleading and the constant 2 −N 2 term dominates up to possible O(1) factors just as T crosses (1/J) log N . Furthermore, by the time the connected solution actually dominates the mutual purity, these terms are suppressed by an even stronger factor of e −N compared to the constant 2 −N 2 term.

A.2 One-loop determinant around disconnected solution
In the Hamiltonian picture, the one-loop determinant is related to corrections in the energy eigenstates and the corresponding eigenvalues near the ground state |g 1 . From equation (A.10), we see that the energy eigenstates |e g 1 k gain corrections from the first two terms related to the ladder operators. However, since they are suppressed by a factor of 1/ √ N q we can ignore these corrections. The correction to energy eigenvalues can be computed by expanding the S q z term to O(1/N 2 ): Thus, the contribution to Tr σ 2 L from the first saddle including corrections at O(1/N ) is To extract the one-loop determinant from the above expression we divide it by the classical saddle point result in equation (A.12) and take the large N 1 limit keeping λ = N 1 N fixed. Define F (T ) as In the large N 1 limit, we use Stirling's approximation for the factorial terms and replace the sum over k by an integral over x ≡ k/N 1 to write F (T ) as follows: where we have defined the functions (A.20) The integrals can be evaluated in the saddle point approximation and we get the following result: Here, x f and x g are saddle points of f (x) and g(x) respectively. Since f (x) and g(x) differ by a term proportional to λ, we can evaluate F (T ) perturbatively in λ. We have the following equations: =⇒ (x f − x g )g (x g ) = λ JT (q − 1) 2 q−3 x g + O(λ 2 ). Using the above relations, we first evaluate the term in the square root. f (A. 25) In a similar manner we can evaluate the expression in the exponential. Finally, we get: The term in the exponential turns out to be equal to the O(λ 2 ) contribution from the classical action while the factor multiplying the exponential piece is the contribution from the one-loop determinant that we computed in (3.68).

B Proof of the error correction bound
We use the two different measures of distance in the proof [49]: the trace distance and fidelity. The trace distance between two states ρ and σ is: In the second expression, we maximize over all possible projectors Q.
The fidelity between two states ρ and σ is defined as: where |ψ ρ and |ψ σ are purification of ρ and σ respectively.
Consider a maximally entangled state |Ψ between the encoded code subspace and a reference system isomorphic to the code subspace: 3) The physical system interacts with the environment initially in some pure state |0 env . This interaction is described by a joint evolution of the physical system and the environment by a unitary U E leading to the following final state: where the set {|φ ij } form an orthonormal basis of the physical Hilbert space. The Schmidt coefficients α j depend only on the environment index because ρ ref is maximally mixed which restricts the form of the Schmidt coefficients in this manner. Indeed, the state ρ env determines the real non-negative coefficients √ α j completely. The condition that |Ψ should be a purification with minimal D(|Ψ , |Ψ ) is hidden in the basis vectors |φ ij phys . Define projection operators Π j as: These projectors satisfy the following relation: Moreover, every subspace corresponding to the projector Π j is isomorphic to the code subspace i.e. for each Π j , there is a unitary operator U j such that U j Π j U † j = Π code , where Π code is a projector onto the code subspace.
Following [76], we construct a recovery channelR which consists of the following two operations: 1. Measurement with some projection operator Π j , and 2. Rotation of the resulting state by the unitary operator U j .
Consider acting withR on |Ψ . Measurement of |Ψ with Π j projects the state |Ψ to the following state with probability α j : The unitary transformation U j acts on |Ψ j as Thus,R acts on |Ψ to give back the original state |Ψ because the unitary U j acts to precisely rotate the basis |φ ij via U j |φ ij = |ψ i . However, we are interested in recovery from the state |Ψ , after the action of the error channel. We will now rephrase the condition for approximate recovery, derived in [76] in terms of the trace distance between the recovered state and the initial state, using the mutual purity between the reference and the environment. We have the following bound on the trace distance between the state obtained by action of the recovery channelR on the actual state |Ψ and the initial state |Ψ . In the second step, we used the monotonicity property of trace distance with respect to the action of a channel (see chapter 9 of [49]). The fourth step follows from the definition of fidelity and the fact that |Ψ is a purification ofρ ref,env that minimizes its trace distance with |Ψ . The sixth step is a standard inequality between fidelity and trace distance [49]. To summarize, we have shown that there exists a set of projection operators Π j , the measurement of which followed by a unitary transformation with U j approximately recovers the maximally entangled state between the reference and the physical system. The accuracy of this recovery in terms of trace distance is bounded by the combination we have found, which is the mutual purity F Ψ (ref : env) from the main text.
The inequality (B.12) was derived for a specific recovery channelR. However, there may exist a better recovery channel R which must also satisfy the inequality: We can use the above inequality to compute a bound on recovery of arbitrary states in the code subspace after the action of the error channel E. We will use the channel-state isomorphism of [77] as follows: Consider a state σ = m,n σ mn |ψ m ψ n | in the code subspace and let σ = R • E(σ). We can write σ in terms of σ ref = m,n σ mn (|m n|) ref  In the fifth step, Π code is the projector on the code subspace. In the final step, we used the inequality in (B.13). The result above is precisely the one quoted in (2.15).