Many-body quantum teleportation via operator spreading in the traversable wormhole protocol

By leveraging shared entanglement between a pair of qubits, one can teleport a quantum state from one particle to another. Recent advances have uncovered an intrinsically many-body generalization of quantum teleportation, with an elegant and surprising connection to gravity. In particular, the teleportation of quantum information relies on many-body dynamics, which originate from strongly-interacting systems that are holographically dual to gravity; from the gravitational perspective, such quantum teleportation can be understood as the transmission of information through a traversable wormhole. Here, we propose and analyze a new mechanism for many-body quantum teleportation -- dubbed peaked-size teleportation. Intriguingly, peaked-size teleportation utilizes precisely the same type of quantum circuit as traversable wormhole teleportation, yet has a completely distinct microscopic origin: it relies upon the spreading of local operators under generic thermalizing dynamics and not gravitational physics. We demonstrate the ubiquity of peaked-size teleportation, both analytically and numerically, across a diverse landscape of physical systems, including random unitary circuits, the Sachdev-Ye-Kitaev model (at high temperatures), one-dimensional spin chains and a bulk theory of gravity with stringy corrections. Our results pave the way towards using many-body quantum teleportation as a powerful experimental tool for: (i) characterizing the size distributions of operators in strongly-correlated systems and (ii) distinguishing between generic and intrinsically gravitational scrambling dynamics. To this end, we provide a detailed experimental blueprint for realizing many-body quantum teleportation in both trapped ions and Rydberg atom arrays; effects of decoherence and experimental imperfections are analyzed.

By leveraging shared entanglement between a pair of qubits, one can teleport a quantum state from one particle to another. Recent advances have uncovered an intrinsically many-body generalization of quantum teleportation, with an elegant and surprising connection to gravity. In particular, the teleportation of quantum information relies on many-body dynamics, which originate from stronglyinteracting systems that are holographically dual to gravity; from the gravitational perspective, such quantum teleportation can be understood as the transmission of information through a traversable wormhole. Here, we propose and analyze a new mechanism for many-body quantum teleportationdubbed peaked-size teleportation. Intriguingly, peaked-size teleportation utilizes precisely the same type of quantum circuit as traversable wormhole teleportation, yet has a completely distinct microscopic origin: it relies upon the spreading of local operators under generic thermalizing dynamics and not gravitational physics. We demonstrate the ubiquity of peaked-size teleportation, both analytically and numerically, across a diverse landscape of physical systems, including random unitary circuits, the Sachdev-Ye-Kitaev model (at high temperatures), one-dimensional spin chains and a bulk theory of gravity with stringy corrections. Our results pave the way towards using many-body quantum teleportation as a powerful experimental tool for: (i) characterizing the size distributions of operators in strongly-correlated systems and (ii) distinguishing between generic and intrinsically gravitational scrambling dynamics. To this end, we provide a detailed experimental blueprint for realizing many-body quantum teleportation in both trapped ions and Rydberg atom arrays; effects of decoherence and experimental imperfections are analyzed.
Quantum teleportation leverages entanglement to transmit quantum information between distant locations [2][3][4][5][6]. Typically, one thinks about teleportation in the context of a few, well-controlled degrees of freedom. For example, two distant observers might share a pair of maximally entangled qubits (i.e. an Einstein-Podolski-Rosen (EPR) pair [7]), enabling a measurement by one observer to teleport an unknown quantum state to the other.
Recently, a confluence of seminal results has unveiled several novel instances of teleportation in stronglyinteracting, many-body systems [1,[8][9][10][11][12][13][14][15][16][17]. Similar to conventional quantum teleportation, these protocols utilize shared entanglement as well as measurement and classical communication. However, they differ from conventional quantum teleportation in a few key aspects. Most notably, prior to teleportation, the initial quantum state is scrambled by the application of a many-body unitary. At first glance, this coexistence of scramblingbroadly speaking, the increasing complexity of initially simple quantum information under many-body time dynamics [18][19][20][21][22]-and teleportation might seem counter- * These authors contributed equally to this work.
intuitive. Indeed, one often thinks of teleportation as a directed quantum channel moving information between two specific locations; in contrast, scrambling disperses quantum information across all of the degrees of freedom in a system. The most natural way to reconcile these two perspectives is through the language of quantum error correction [23]: by encoding, via scrambling, one observer's local information into non-local correlations across a many-body system, one can in fact teleport this information with access only to any few of the system's qubits.
The most notable example of many-body teleportation is the so-called traversable wormhole (TW) protocol, discovered in the context of quantum gravity [1,8,9,[15][16][17]. From the bulk gravitational perspective, this protocol consists of a particle traveling from one side of a wormhole geometry to the other; the wormhole is rendered traversable by the application of a coupling between the two sides. In the boundary theory, the wormhole geometry corresponds to a highly entangled thermofield double (TFD) state shared between two copies of a many-body system, and the coupling is implemented via measurement and feed-forward operations [ Fig. 1(a)]. Crucially, for this bulk-boundary correspondence to hold, the Hamiltonian describing the boundary system must exhibit "coherent", gravitational scrambling dynamics- The protocol hosts two mechanisms of teleportation: peaked-size (red) and gravitational (blue). The channel capacity of peaked-size teleportation decreases with increasing time (dark to light red), while its fidelity decreases with decreasing temperature (dark to light red, again). At high temperature and late times, it is equivalent to teleportation in the HPR protocol (red diamond). Gravitational teleportation occurs at low temperatures in systems dual to semiclassical gravity (e.g. the SYK model), and exhibits the same channel capacity but higher fidelity compared to peaked-size teleportation. Increasing the strength of stringy corrections to the gravity theory interpolates between gravitational and peaked-size teleportation. (c) The two mechanisms display distinct time profiles for the teleportation fidelity at fixed coupling strength, g. In systems dual to gravity (top), the fidelity features a single O(1) peak near the scrambling time (gravitational, blue), and a late time revival (peaked-size, red) to a fidelity suppressed by the two-point function G β [1]. In generic thermalizing systems (bottom), the fidelity oscillates between 0 and G β with phase proportional to the operator size, may subsequently decay if sizes become not peaked, and revives at late times.
this is realized, most notably, in the Sachdev-Ye-Kitaev (SYK) model at low temperatures [24,25]. Interestingly, recent work has uncovered a number of instances of many-body teleportation without gravitational dynamics. For example, teleportation in the TW protocol was recently demonstrated analytically in the SYK model at high temperatures [17], and numerically in chaotic spin chains at late times [15,16]; in both cases, the microscopic mechanism for teleportation remains an outstanding puzzle. In addition to the TW protocol, an alternate many-body teleportation protocol was introduced in the context of the Hayden-Preskill variant of the black hole information paradox [11,23]. This socalled Hayden-Preskill recovery (HPR) protocol allows for many-body teleportation via generic scrambling dynamics. Although the two protocols bear some structural similarity, the HPR protocol is exponentially less efficient for teleporting multiple qubits. To this end, understanding the precise relationship between these protocols remains an essential open question.
In this work, we present a unified framework for manybody teleportation from the perspective of the growth of operators under scrambling time-evolution. Most significantly, this framework leads to the identification of a new teleportation mechanism-dubbed peaked-size teleportation-which succeeds for a wide variety of physical systems and encapsulates all known examples of manybody teleportation outside of the gravitational regime.
We emphasize that peaked-size teleportation represents a distinct teleportation mechanism compared to "gravitational" teleportation. Although the same TW protocol can host either mechanism, the features of peaked-size teleportation differ markedly from those of gravitational teleportation [ Fig. 1(c), Table I]. Crucially, this distinction implies that the TW protocol can act as a litmus test for identifying intrinsically gravitational dynamics. More broadly, our results pave the way towards utilizing the TW protocol as a powerful experimental tool for characterizing the growth of operators in strongly interacting systems.

I. SUMMARY OF RESULTS
We now provide a technical overview of our main results and the organization of our manuscript. This summary is intended to introduce the overarching concepts of our work, such that the remaining sections are selfcontained and can be read according to individual preference. A more detailed, section-by-section guide to the reader is included at the end of this summary.
In Section II, we begin with a general description of the TW protocol [ Fig. 1(a)]. In this protocol, locally encoded quantum information is inserted into one side of an entangled thermofield double (TFD) state and teleported to the other side through a combination of (i) unitary evolution of each side individually, and (ii) a simple twosided coupling that acts on a large subsystem of each side. The coupling is quite flexible in form, and corresponds to unitary evolution, e igV , under a two-sided interaction where {O i } are any set of K local, non-identity operators applied to the left (l) and right (r) side of the system. This coupling can be performed as either a quantum gate, or through local measurements of O i on the left side, followed by classical communication and feed-forward operations on the right side [ Fig. 1(a)].
In Section III, we discuss the general requirements for successful teleportation in the TW circuit. In particular, we relate the teleportation fidelity to the following correlation functions of the two-sided system [1]: where Q(±t) is a time-evolved operator initially acting on the qubit(s) to be teleported. Our analysis leads to two conditions on these correlators that, when combined, are necessary and sufficient for teleportation to succeed with unit fidelity: 1. The magnitudes of the correlators must be maximal for every Q.
2. The phases of the correlators must be the same for every Q.
Here, Q runs over a complete basis of operators on the qubits to be teleported. We show that Condition 1 is naturally satisfied, even without the coupling V , if the TFD state is at infinite temperature, in which case it reduces to an extensive set of maximally entangled EPR pairs. On the other hand, Condition 2 requires that the coupling acts non-trivially on the operators Q.
In Section IV, we describe the relation between the coupling, V , and the growth of time-evolved operators, Q(t). For the purposes of teleportation, this growth is characterized by the size distribution of the operators [26][27][28], which provides a finer-grained measure of quantum information scrambling compared to more conventional quantities such as out-of-time-ordered correlators (OTOCs) [19,21,29]. Specifically, writing Q(t) as a sum over Pauli strings, Q(t) = R c R (t)R, we define the size distribution as: where the sum is over Pauli strings, R, of size, S (equal to the string's number of non-identity components). By probing correlations between the two sides of the doubled Hilbert space, the coupling V directly measures the operator size [27].
In Section V, we introduce the peaked-size mechanism for many-body teleportation. This mechanism is possible whenever the size distributions of time-evolved operators, Q(t), are tightly peaked about their average size. In this scenario, the exponentiated coupling, e igV , applies approximately the same phase, proportional to the size, to each coefficient, c R , and therefore to the entire operator, Q(t). We show that these applied phases are sufficient to align the correlators' phases for all operators Q, thereby achieving Condition 2. We also demonstrate that the magnitudes of the correlators are unchanged by the coupling when size distributions are tightly peaked. This implies that peaked-size teleportation achieves perfect fidelity at infinite temperature, where Condition 1 is automatically satisfied; at finite temperature, peaked-size teleportation can still occur, but with a reduced fidelity (Table I).
In Sections VI-VII, we analyze examples of peakedsize teleportation across a wide variety of interacting, many-body dynamics. We demonstrate that the capabilities of peaked-size teleportation-most notably, the fidelity and the number of qubits that can be sent (i.e. the channel capacity)-depend on the temperature, coupling strength, evolution time, and the specific scrambling dynamics of the model under study (Table I).
More specifically, in Section VI, we provide general arguments that all scrambling systems exhibit peakedsize teleportation at late times, after the system's scrambling time (t t s ). In this regime, operators have become fully delocalized across the system, so that their size distributions are peaked about a typical, extensive value. We also show that late time peaked-size teleportation is limited to transmitting only a single qubit.
In Section VII, we show that many scrambling quantum systems also feature peaked-size teleportation at intermediate times, i.e. after the local thermalization time but before the scrambling time (t th t t s ). We begin with ergodic short-range interacting systems in ≥1D, which we show naturally possess peaked-size distributions due to thermalization within the bulk of a timeevolved operator's light cone. In contrast, the size distributions of operators in all-to-all coupled (0D) systems are not intrinsically peaked; nevertheless, peaked sizes can be engineered by non-locally encoding the quantum information before insertion into the teleportation circuit. Interestingly, in both of these classes of dynamics, we find that multiple (∼ O(K)) qubits can be teleported simultaneously via the peaked-size mechanism, in contrast with the unit channel capacity of late time teleportation. We substantiate these claims through extensive numerical and analytic studies on a variety of physical models: random unitary circuits (RUCs) in dimensions d = 0, 1, and 2 [30], the SYK model, and (in Section IX D) experimentally relevant spin chain Hamiltonians [31].
In Section VIII, we discuss the interplay between peaked-size and gravitational teleportation. Notably, we expect gravitational teleportation to occur only at low temperatures, where certain quantum mechanical mod-

Model
Teleportation mechanism Protocol parameters Maximum per qubit fidelity

Channel capacity
All scrambling systems at late times (Refs. [1,11], Section VI) peaked-size g = π mod 2π ∼ G β 1 qubit ≥ 1D RUCs & chaotic spin systems (Sections V B,VII B,IX D) peaked-size η d gS(t)/N = π mod 2π ∼ G β ∼ K qubits 0D RUCs, with encoding (Section VII C) peaked-size η d gS(t)/N = π mod 2π ∼ 1 ∼ K qubits High-T SYK, with encoding (Ref. [17], Sections VII D) peaked-size η d gS(t)/N = π mod 2π ∼ 1 ∼ K qubits Low-T SYK / AdS2 gravity (Refs. [1,8,15,17], Fig. 1) AdS2 gravity with strong stringy corrections, with encoding (Section VIII D) peaked-size gS(t)/N ∼ π mod 2π ∼ G β - Table I. Summary of our expectations for teleportation in a variety of physical models. For each model, we specify the associated teleportation mechanism, the optimal value of the coupling strength g, the optimal teleportation fidelity, and the channel capacity. Here G β is the imaginary time two-point function (Section V B), S(t) is the size of a time-evolved operator, K is the number of measured qubits [ Fig. 1 els (e.g. the SYK model) are known to possess a dual description in terms of conformal matter coupled to gravitational dynamics. From the perspective of operator growth, the unique feature of gravitational teleportation is the presence of non-trivial phase winding in a variant of the size distribution [15]. Crucially, this effect enables gravitational teleportation to satisfy Condition 1, and thereby achieve high teleportation fidelity at low temperatures, in sharp contrast with peaked-size teleportation (Table I). Intriguingly, while it may seem that there is a sharp distinction between peaked-size and gravitational teleportation, we find that this is not always this case. In particular, we show that varying the temperature of the SYK model provides a continuous interpolation between gravitational teleportation at low temperature and peaked-size teleportation at high temperature. In the dual picture, perturbing away from the low temperature limit corresponds to adding stringy corrections to the gravity theory [25,32,33]. Following this intuition, we show that teleportation in a gravity theory with strong stringy corrections [1] bears a remarkable qualitative similarity to peaked-size teleportation, thus providing a first step towards a bulk understanding of this phenomenon.
Finally, in Section IX, we discuss experimental applications of the TW protocol for probing many-body dynamics. In particular, we demonstrate that the protocol can function as a diagnostic tool for scrambling dynamics in near-term quantum simulators, enabling one to starkly distinguish between generic thermalizing systems and gravitational dynamics. To this end, we pro-vide detailed blueprints for realizing the protocol in two complementary experimental platforms-Rydberg atom arrays [31,[34][35][36][37][38] and trapped ions [39][40][41][42][43]. Specifically, the observation of a high teleportation fidelity at low temperatures would be a tantalizing experimental indicator of gravitational scrambling dynamics. In addition, gravitational dynamics exhibit unique qualitative features as a function of both evolution time and protocol parameters [ Fig. 1(c), Table I]. More broadly, our analysis suggests that the TW protocol can provide insights into manybody dynamics outside the gravitational regime. In particular, we demonstrate that the fidelity of peaked-size teleportation probes higher moments of operator size distributions [28].
Guide to the reader -Considering the wide scope of results presented in this work, we encourage readers to skip to sections that align with their specific interests and refer to the above summary for context. To this end, we highlight below the nature of each section and provide recommendations for readers of different backgrounds. Sections II-IV introduce the formal tools and derivations necessary for rigorously understanding our results. These sections will be of interest to readers with a background in quantum information who wish to understand the precise connection between teleportation and operator sizes. Sections V-VII introduce peaked-size teleportation and analyze its realization in several example systems. Since many these systems are experimentally accessible, these sections will be most relevant to members of the quantum simulation and many-body physics communities. Section VIII focuses on the interplay of peaked-size teleportation and gravitational physics, both in the SYK model and from a bulk gravitational perspective. For brevity, background material on gravitational physics is relegated to references, making this section best suited for experts at the intersection of quantum information and quantum gravity. Finally, Section IX contains a summary of the experimental signatures of the TW protocol, detailed blueprints for Rydberg atom and trapped ion implementations, and a discussion of the protocol's behavior under experimental error. This section will be of interest to AMO experimentalists and all readers interested in near-term realizations of many-body quantum teleportation [44].

A. Relation to previous works
To further elaborate on the broad context of our results, a brief summary of the relevant prior studies and their relation to our work is provided as follows.
Gravitational teleportation in the TW protocol-Traversable wormhole teleportation was originally introduced in Refs. [1,8] in the context of gravitational physics, where it was realized that a coupling of the form V enables a traversable channel between the boundaries of a two-sided black hole. The explicit quantum mechanical circuit implementing this teleportation [ Fig. 1(a)] was later introduced in Refs. [15,17], alongside exact calculations for the teleportation fidelity in the large-q SYK model [17]. While the emphasis of our work is not on the bulk interpretation of gravitational teleportationindeed, the peaked-size teleportation mechanism is intended to contrast with the gravitational mechanism-it will be helpful to recall the main results from the gravitational perspective.
We focus on the specific case of two-dimensional antide Sitter space, which is the bulk dual of the SYK model at low temperatures [1,45]. In the simplest case (ignoring gravitational backreaction), the two-sided correlator, Eq. 2, can be explicitly calculated and is given by [1]: Here, G N is Newton's constant, β = 1/T is the inverse temperature of the black hole, ∆ O is the conformal dimension of the coupling operators O i [Eq. (1)], and ∆ Q is the conformal dimension of the operator Q. In the context of the SYK model, G N is inversely proportional to the number of Majorana fermions, N , and the black hole temperature is equal to the temperature of the TFD state [1,17]. For our purposes, the most notable feature of the correlator is that it exhibits a sharp peak at time t ≈ G N log(g) [ Fig. 1(c)], corresponding to the moment a particle inserted on one side of the black hole emerges on the other side. While in the above formula [Eq. (4)], the correlator diverges at this time, in the large-q SYK model, this divergence is regularized and the correlator peaks at its maximal value of unity [17]. Thus, at time t ≈ G N log(g), the correlator satisfies Condition 1 for successful teleportation; in Ref. [17], it was shown that Condition 2 is also satisfied for certain conformal dimensions of the operators Q. In combination, this leads to unit teleportation fidelity.
Another notable feature of gravitational teleportation is the ability to teleport multiple qubits simultaneously, as discussed in Ref. [1]. In the gravitational picture, multi-qubit teleportation has an intuitive explanation: particles corresponding to different qubits pass through the black hole in parallel, without interacting with one another. However, for sufficiently many qubits, the effects of gravitational backreaction become important, leading to a predicted channel capacity of O(K).
HPR teleportation-An independent, but closely related, set of protocols for many-body teleportation was introduced in Ref. [11] for the recovery of information in the Hayden-Preskill thought experiment [23]. Unlike previous works on traversable wormholes, in Ref. [11] teleportation succeeds for any fully scrambling unitary dynamics (i.e. at late times, t t s ), with no reliance on gravitational physics. However, the channel capacity of HPR teleportation is fundamentally limited: multi-qubit teleportation requires a protocol whose circuit depth grows exponentially in the number of qubits to be teleported [11].
In Appendix B, we show that a deterministic variant of the HPR protocol (for single-qubit teleportation) is in fact equal to the TW protocol in Fig. 1(a), restricted to infinite temperature and with a particular choice of the coupling operators, O i . Furthermore, in Section VI B we show that teleportation at late times via the peaked-size mechanism is equivalent to this variant of HPR teleportation. However, peaked-size teleportation is more powerful than HPR teleportation in the sense that: (i) it succeeds for a much larger class of couplings, V , (ii) it can succeed at intermediate times, and (iii) at such times, it is capable of sending multiple qubits with no change in the protocol's complexity, an exponential improvement over the HPR protocol.
Previous many-body teleportation experiments-Manybody quantum teleportation has recently been demonstrated in both trapped ion [13] and superconducting qutrit [14] experiments. Both Refs. [13,14] implement a probabilistic variant of the HPR protocol, which differs slightly from the TW protocol, while Ref. [13] also implements the deterministic variant discussed above. In all cases, the scrambling dynamics, U , are generated by digital quantum gates acting on a small number of qubits. Teleportation is performed for a single qubit and a fully scrambling unitary, placing the experiments in the same physical regime as late-time, peaked-size teleportation.
Our work demonstrates that experiments in the TW protocol at intermediate times can access new regimes of many-body quantum teleportation, with the potential to provide more information about the scrambling dy-namics under study. Most notably, such experiments can distinguish between teleportation in generic many-body systems (via the peaked-size mechanism) versus systems with a gravity dual (via the gravitational mechanism), which is not possible in the HPR protocol.
SYK teleportation in the TW protocol-In Ref. [17], the two-sided correlator of the TW protocol [Eq. (2)] was calculated exactly for the large-q SYK model (defined in Section VII D). As anticipated in Ref. [1], the correlator at low temperatures-where the model is dual to gravity-agrees with the gravitational result [Eq. (4)] up to the previously mentioned regularization. More surprisingly, it was shown that teleportation with unit fidelity is also possible at high temperatures-where the model is not dual to gravity. As we will see in Section VII D, all features of high temperature teleportation in the SYK model are in precise agreement with the peaked-size mechanism; our work thus provides a microscopic understanding for this previously unexplained result.
Gravity in the lab-Ref. [15] discusses various instances of teleportation in the TW protocol. The authors distinguish two teleportation mechanisms: (i) an "operator transfer" mechanism, which occurs at intermediate times in gravitational systems and is capable of teleporting multiple qubits, and (ii) a "state transfer" mechanism, which occurs at late times in all scrambling systems, and is capable of sending only a single qubit. Moreover, they introduce a microscopic interpretation for the teleportation mechanism in gravitational systems, termed "size winding", which we connect to in Section VIII B.
In our terminology, the first teleportation mechanism corresponds to gravitational teleportation, while the second mechanism corresponds to peaked-size teleportation at late times 1 . In our work, we provide a microscopic interpretation for late time teleportation (i.e. the peakedsize mechanism) and demonstrate that it is equivalent to teleportation in the HPR protocol. In addition, we demonstrate that peaked-size teleportation is a more general phenomenon that also occurs at intermediate times in many systems, where we show that it is capable of teleporting multiple qubits.
In a follow-up work, Ref. [16], whose pre-print was posted concurrently with that of this work, the same authors elaborate on their previous results and provide more detailed examples and calculations. These agree with our own results in areas of overlap.

II. INTRODUCTION TO DIAGRAMMATIC NOTATION
We begin by introducing a diagrammatic "tensor network" notation for depicting the teleportation circuit. Adapted from Ref. [11], this notation provides a precise visual framework for analyzing teleportation in Section III and will be convenient for deriving rigorous results on the teleportation fidelity in Section V C.
To begin, we represent a quantum ket |ψ and bra ψ| as: (5) Note that time proceeds upwards-an initial state |ψ terminates the bottom of a leg, while a final projection ψ| terminates the top. Similarly, much as in Fig. 1(a), we represent an operator, for instance the many-body unitary U , as a box with input (bottom) and output (top) legs: (6) Here we have decomposed the input and output into two subsystems, A and its complement for the input, C and its complement for the output, in reference to the teleportation protocol. Specifically, comparing to Fig. 1(a), subsystem A consists of the qubits supporting the input state |ψ , while subsystem C consists of the coupled qubits.
The diagrammatic notation is particularly useful when working with EPR states. The EPR state on two qubits is defined as |EPR = (|00 + |11 )/ √ 2; for a system of N d-dimensional qudits, this is generalized to 1 Here {i} is an arbitrary d Ndimensional basis, * denotes time-reversal (i.e. complex conjugation), and l and r denote the left and right system, respectively. In the diagrammatic notation, we represent this as: We have again decomposed each system into two subsystems, A and its complement,Ā, for convenience (subsystem A is chosen to be identical between the left and right sides). Each dot represents a normalization factor given by the inverse square root of the subsystem's dimension.
To see the utility of the diagrammatic notation, recall that a fundamental property of the EPR state is that an operator acting on the left side is equivalent to its transpose acting on the right: where the middle equality swaps the i, j indices of the sum. In diagrammatic notation, this becomes simply (9) i.e. the operator O "slides" from the left to right side of the EPR pairs, with its input and output indices correspondingly transposed. Similarly, expectation values in the EPR state can be easily computed in terms of the trace of (one-sided) operators, e.g.
where the final equality follows from The EPR state is closely related to the thermofield double Here H is a time-reversal symmetric Hamiltonian, H = H * , with eigenstates |E i , and eigenvalues E i . The TFD state is parameterized by an effective "temperature" 1/β. At infinite effective temperature (β = 0), the TFD and EPR states are equal. At finite temperature, the TFD state is obtained by applying the square root of the density matrix, ρ 1/2 ≡ e −βH/2 / tr e −βH 1/2 , to either side of the EPR state, which we represent as: (11) For the finite temperature TFD state, the analog of Eq. (9) holds only for operators that commute with the Hamiltonian. Most notably, such operators include the time-evolution operator, U = e −iHt , which thus obeys: (12) Eq. (12) also holds for backwards time-evolution, replacing U → U † , U T → U * . We note that for time-reversal symmetric Hamiltonians, U = U T . In this case, combining Eqs. (12,9), we have the useful identity: Applying Eq. (10), we can again express 'two-sided' expectation values in the TFD state in terms of 'one-sided' correlation functions, e.g.
Let us now re-draw the full teleportation protocol in Fig. 1(a) using the diagrammatic notation: (15) This circuit proceeds as follows: (i) prepare the TFD state, (ii) insert the state |ψ on subsystem A of the left side, (iii) time-evolve the two sides by U , U * , (iv) couple the two sides via the unitary operator e igV , with V as in Eq. (1), (v) evolve the right side by U T , (vi) apply a 'decoding' operator D, and (vii) measure the output state of subsystem A on the right side. Compared to Fig. 1(a), we have made two modifications. First, we have replaced the measurement and classical communication with a quantum coupling e igV , as described in Section I. Second, we now include a simple decoding operator, D, applied at the end of the circuit before state recovery. We will find that D = Y ⊗ . . . ⊗ Y for peaked-size teleportation of a multi-qubit subsystem, where Y is the single-qubit Pauli Y operator (Section V).
Finally, we note that a straightforward application of Eq. (12) allows us to re-express the circuit as (16) This equivalent version of the protocol was introduced in Refs. [15,17] and will be more convenient for analysis from here on.

III. GENERAL REQUIREMENTS FOR SUCCESSFUL TELEPORTATION
We now introduce heuristic arguments for when teleportation succeeds in this protocol. This will culminate in the two requirements for teleportation listed in Section I. In Section V, we derive these conditions more formally by providing exact relations between the two-sided correlators in Eq. (2) and the teleportation fidelity.
We begin with the protocol in Eq. (16). To proceed, we insert a resolution of the identity 1 = φ |φ φ| on the "swapped out" subsystem A (the output of U † l ) 2 : This reformulation makes it clear that teleportation depends on the action of the coupling on states of the form Teleportation succeeds when the coupling "transfers" |ψ(t) φ(t)| from the left to right side of the TFD state. More precisely, the following identity, if true for all operators Q A on A, would guarantee successful teleportation for all states: (18) Here θ Q is an overall phase and we represent conjugation by the decoding operator asQ A ≡ D † Q A D. One can verify this explicitly by plugging the RHS of the above equality into Eq. (17): the topmost applications of DU T and U * D † cancel, leaving Q A → |ψ φ| as the topmost operator on the right side, i.e. subsystem A is in the state |ψ .
To quantify whether this equality holds, we measure 2 At infinite temperature, using Eq. (9), |φ can be understood as the counterpart of |ψ , to be teleported from right to left instead of left to right. To see this, use Eqs. (9,12) to re-expess φ| l U † l → U † l |φ r , and apply D l U † l after coupling to recover |φ . 3 Traditionally, this would be considered reverse time-evolution, and denoted Q A (−t). For brevity, we have flipped the sign of t throughout the text. the inner product between the two states 4 : (19) This is precisely the two-sided correlation function introduced in Eq. (2), now modified to include the decoding operator. In particular, if the correlation function is maximal for all operators Q A , then Eq. (18) holds and teleportation succeeds with perfect fidelity for all initial states.
In practice, it is sufficient to evaluate the correlators for a complete basis of operators on subsystem A (e.g. the Pauli operators). In this case, we now have two requirements on the operator correlators, as listed in Section I: (i) all correlators must have maximal magnitude, i.e. equal to 1, and (ii) all correlators must have the same phase-if two operators both individually obey Eq. (18) but with different phases, their sum will not.
At infinite temperature, owing to Eq. (9), we will see that the first requirement is satisfied even in the absence of the coupling, for any symmetric or antisymmetric operator. To satisfy the second requirement, the role of the coupling e igV must be to apply a Q A -dependent overall phase. In the following section, we analyze the action of the coupling and show precisely when such an overall phase occurs.

IV. CONNECTION TO OPERATOR SIZE
In this section, we outline the connection between the coupling V and the operator size when V is acted on states of the form: This connection was discovered in a number of previous works, focusing primarily on a specific bilinear coupling in fermionic systems [15,16,26,27,[46][47][48]. In the following, we introduce this connection in the context of bosonic systems and argue that it applies to a good approximation for any generic, local couplings. From this, we then show that the action of the exponentiated coupling, e igV , is particularly simple-it applies an overall phase-whenever operator size distributions are tightly peaked.

A. Coupling measures size
In bosonic qudit systems, we define the size of a Pauli string as its number of non-identity elements [26]. For instance, the Pauli string has size 3. A more general operator can be written as a sum of Pauli strings, R: and possesses a corresponding size distribution [26,27] 5 : The distribution is normalized to 1 if Q A is unitary, One can naturally characterize the size distribution via its moments-for instance, the average size, S[Q A (t)ρ 1/2 ] ≡ S P (S)S (when context is clear, we denote this simply as S), and the size width, δS.
We will now show that the coupling V approximately measures the operator size, in the sense that it acts on states of the form Eq. (20) as: is an order one constant determined by the local qudit dimension, d. Expectation values of V thus measure the average size, while higher powers of V measure higher moments of the size distribution [26,27]. In particular, the exponentiated coupling in the teleportation protocol applies a size-dependent phase to each Pauli string of Q A (t)ρ 1/2 : We derive this connection by first introducing an exact measure of operator size in bosonic qudit systems, generalizing previous measures for Majorana fermionic systems [26,27]. We then argue that successively more generic couplings display approximately the same behavior, when acted on time-evolved operators in generic many-body scrambling dynamics.
In bosonic qudit systems, we find that the operator size is precisely measured by a sum of individual EPR projectors on each qudit i: where d is the local qudit dimension, N is the number of qudits, and P i form a complete basis of single-qudit operators (e.g. for qubits P i ∈ {1, X, Y, Z}). To see this, let us first analyze the action of a single EPR projector, P EPR,i . Writing a given Pauli string as a tensor product of single-qudit Paulis, R = N j=1 R j , we find using Eq. (10) and tr i (R i )/d i = δ Ri,1 . A single EPR projector thus acts as a binary variable, giving eigenvalue 1 or 0 if a given Pauli string is, or is not, the identity on the designated qudit. The full coupling is a sum of these binary variables over all qudits and therefore counts the total number of non-identity elements in the Pauli string, i.e. the operator size. Its eigenvectors are the states R l |EPR with eigenvalues 1 − S[R]/N , as in Eq. (25).
We now turn to more general local couplings. First, as a trivial but useful modification, we can remove the identity operators from V s , since these are not included our original definition of the coupling, V [Eq. (1)]. These constitute a fraction 1/d 2 of the complete basis, P i , summed in Eq. (27). Removing these terms renormalizes the eigenvalues of the coupling: which now match those quoted in Eq. (25). Note that the left side sum is now over N (d 2 −1) non-identity operators and normalized accordingly. Second, we consider omitting some of the non-identity P i at each site. Intuitively, under thermalizing dynamics, if an operator has spread to some qudit i it should not matter which Pauli operator we use to probe the operator's presence. For example, for qubits, we could omit the O j = X i , Y i couplings and keep only O j = Z i . A random Pauli string has equal probability to commute with Z i as it would with X i and Y i ; thus, coupling using only Z i operators is sufficient for measuring a thermalized operator's size.
Third, we expect even more general couplings, composed of O i that are local but not necessarily Pauli operators, to behave similarly. Specifically, each individual coupling, O i,l O i,r , will asymptote to two different expectation values before and after the time-evolved operator has spread to the support of O i . Before, the coupling will maintain its expectation value in the unperturbed TFD state, tr After, the spread of Q A (t) will disrupt the two-sided correlations in the TFD state that give rise to this initial expectation value, and the coupling will instead asymptote to its value in two thermal states, tr(O i ρ) · tr(O i ρ). As before, the sum of many terms, each behaving as above, leads to an approximate measure of operator size.
Lastly, we consider the case where the coupling is restricted to act only on some subsystem C, consisting of K qudits 6 . The coupling now measures the number of nonidentity elements of a Pauli string within C-we denote this as the K-size, S K , of the Pauli string. The eigenvalues of the coupling are the same as those in Eq. (29), with the replacement S/N → S K /K. For a typical Pauli operator, we expect the K-size distribution of an operator to be similar to its full size distribution when K is large and the coupled qubits are distributed randomly. In particular, in this scenario we expect the average K-size, S K , to be proportional to the average size, S, For simplicity, we will make this substitution in the remainder of the work. However, if C is a spatially local subsystem (instead of a random subsystem), then this replacement will be modified depending on the spatial extent of the operator.
As a final remark, we note that the operator size distribution is directly related to out-of-time-order correlators (OTOCs), a more familiar quantity for probing operator growth [19,21,29]. In particular, the average size is equal to a sum of OTOCs between Q A and O i [26,27], (31) using Eqs. (9)(10)(11)(12)(13)(14). Higher moments of the size distribution can also be probed by OTOCs, now between Q A and various products of the O i , e.g. O i O j for the size width. We discuss these relations further, paying particular attention to subtleties that arise at finite temperature, in Section VIII.

B. Peaked-size distributions
The exponentiated coupling [Eq. (26)] has a particularly action when the size distribution of Q A (t)ρ 1/2 is tightly peaked about its average size. In this regime, each Pauli string gains approximately the same phase, and so the action of the coupling reduces to applying a Q Adependent overall phase, where the applied phase is proportional to the average K-size [see Eq. (29,30)], Corrections to this behavior are controlled by higher moments of the size distribution. Focusing on the overlap of the coupled and uncoupled states, the leading order correction is equal to the K-size variance, δS 2 Q , multiplied by g 2 : The K-size variance receives contributions from two sources: the variance of the full size distribution, δS 2 , and a statistical error from sampling only K of N qubits for the K-size. If the K qubits are distributed randomly, these errors scale as δS K ∼ δS · (K/N ) and δS K ∼ √ S K ≈ SK/N , respectively (see Appendix F for a detailed derivation of the latter). These are small compared to the average K-size whenever δS S and 1 S K . In Appendix A, we go beyond these leading order corrections and provide quantitative bounds on when the peaked-size approximation in Eq. (32) is valid. In general, we can strictly prove that this approximation holds whenever there is a parametric separation between an asymptotic size width, defined in the appendix, and the average size.

V. PEAKED-SIZE TELEPORTATION
Having established general conditions for successful teleportation (Section III) as well as the connection between the coupling in the TW protocol and operator size distributions (Section IV), we are now ready to introduce the peaked-size mechanism for teleportation. In this section, we first demonstrate peaked-size teleportation in its simplest context: teleportation of a single qubit at infinite temperature. We then show that the fidelity of peaked-size teleportation is necessarily suppressed at finite temperature. For ease of reading, we relegate rigorous results supporting each of the above arguments to the end of the section. We turn to specific physical systems realizing peaked-size teleportation in the following sections: in Section VI we show that peaked-size teleportation of a single qubit occurs in all scrambling systems at late times, while in Section VII we show that peaked-size teleportation of multiple qubits occurs in certain systems at intermediate times.

A. Single-qubit teleportation
To analyze teleportation of a single qubit, we turn to the two-sided correlators in Eq. (19), with Q A ∈ {1, X, Y, Z} running over the single-qubit Pauli operators. We recall that the requirements for teleportation are for all C Q to have (i) maximal magnitude and (ii) the same phase.
The first requirement is naturally satisfied at infinite temperature even before coupling and decoding but the second requirement is not. In particular, the four correlators with D = 1, g = 0 are: where the left entries are qubit operators, Q A , and the right entries are the correlators, C Q . The correlators have maximal magnitude because each operator can be transferred perfectly from left to right using Eq. (9). However, the Y operator picks up an overall minus sign during this process, since Y T = −Y , and so the correlator phases are not aligned. One can verify the resulting teleportation fidelity is indeed trivial. Our goal will be to show that the action of the coupling in Eq. (32), as well as a simple decoding operation, are sufficient to align the four phases.
To begin, we assume that all time-evolved Pauli operators have a tightly peaked size distribution and that the average size S is the same for all non-identity operators. From Eqs. (32)(33), we have that the coupling applies a total phase difference η d gS/N between the thermofield double state (the identity operator; size zero) and all perturbed states (time-evolved Pauli operators; size S). Our table of correlator phases is thus modified to: We again do not achieve perfect phase alignment. However, we can now correct the misaligned phases using the decoding operator, D = Y . This applies an additional minus sign to the X and Z correlators: The correlator phases are now aligned whenever leading to perfect teleportation at these values.
B. Peaked-size teleportation at finite temperature There are two important modifications to peaked-size teleportation at finite temperature. First, the relevant notion of operator size is modified [27]. In particular, in the peaked-size regime, the difference in phase applied between the identity and non-identity Pauli operators is modified to Second, the maximal fidelity of peaked-size teleportation is reduced at finite temperature. In particular, when sizes are tightly peaked, the two-sided correlators factorize into a constant magnitude multipled by an overall phase: where θ Q combines the effects of transposition, coupling, and decoding, and the correlator magnitude corresponds to an imaginary-time Green's function, This Green's function is unity at infinite temperature and generically decreases at finite temperatures, due to the reduced entanglement of the TFD state. This violates the maximal magnitude requirement for teleportation, and therefore leads to a corresponding decrease in the teleportation fidelity.
The astute reader will recall that finite temperature teleportation is known to succeed with O(1) fidelities (i.e. higher than G β ) in theories with a gravity dual [1,8,17]; this is a signature of physics outside the peaked-size regime, which we connect to in Section VIII.

C. Rigorous expressions for teleportation fidelity
We now derive formal expressions of the teleportation fidelity for n teleported qubits as a function of the correlator phases. To do so, we consider a variant of the protocol where instead of teleporting a quantum state we attempt to distill an EPR pair: (39) Here state insertion is replaced by swapping in one "half" of an EPR pair with a reference subsystem R (far right) into subsystem A of the left side. When subsystem A is teleported from left to right, the circuit results in an EPR pair between the reference subsystem R and subsystem A of the right (top arrows). The fidelity of EPR distillation is precisely related to the average fidelity of state teleportation [12], where d A = 2 n is the dimension of subsystem A when teleporting n qubits.
We calculate the teleportation fidelity by Pauli decomposing the SWAP operator as SWAP This gives: .
where the third equality utilizes the diagrammatic identities Eqs. (9,10), and the fourth equality inserts the identity, 1 = D r U r U † r D † r , in the center of the right side (recall our notationQ 1/2 = D † Q 1/2 D). Writing the rightmost diagram as an equation, we have: Similar expressions for teleportation of quantum states are contained in Appendix C.
In general, the teleportation fidelity and two-sided correlators are related only by a lower bound, 7 This is obtained diagrammatically by inserting the pro- 7 Under special circumstances, namely large-N models, one may be able to factorize the above expression in terms of correlators of the form Eq. (19) [17]. jector, |TFD TFD|, into the center of Eq. (40): (43) A similar bound was obtained in Ref. [15,16], conditional on certain assumptions about operators' size distributions.
At infinite temperature in the peaked-size regime, we have C Q = e iθ Q and the fidelity is equal to the lower bound: The sum is over d 2 A terms, and is unity only when all the operators' phases are the same. In the case of a singlequbit teleportation at infinite temperature in the peakedsize regime, plugging the final table of Section V A into the above equation gives a fidelity: which oscillates between trivial fidelity (F EPR = 1/4) and unity as a function of the operators' size. At finite temperature in the peaked-size regime, we instead find where the maximum fidelity is again achieved when the correlator phases align. However, its value is now less than unity, and instead is equal to a sum of various imaginary time Green's functions, i.e. the correlator magnitudes [Section V B, Eq. (38)].

VI. PEAKED-SIZE TELEPORTATION AT LATE TIMES
We now introduce the simplest physical example of peaked-size teleportation: teleportation in any scrambling system at late times (after the scrambling time). There are two distinguishing features of this regime: (i) the circuit can only teleport a single qubit, i.e. the channel capacity is one, and (ii) as for all peaked-size teleportation, the teleportation fidelity is suppressed at low temperatures. We also demonstrate that this regime of peaked-size teleportation, as well as the full quantum circuit implementing the TW protocol, are equivalent to HPR teleportation of a single qubit. In Section VII, we will demonstrate that the single-qubit late time channel capacity can be overcome at intermediate times in many scrambling systems.

A. Teleportation at late times
At late times, the dynamics of a scrambling system can be approximated by a Haar random unitary 8 [23,52]. In this case, each time-evolved operator, Q A (t), becomes a sum of random Pauli strings, each with probability 1/d 2 to be the identity at any individual site. As a result, time-evolved operators have an average size, and a size width, where the scaling is based on the central limit theorem. The K-size distribution takes the same form, replacing N with K, and is tightly peaked as long as K is large (specifically, gδS K /K ≈ g/ √ K 1). For simplicity, we will focus on late time teleportation at infinite temperature; finite temperature modifications follow according to Section V B. Using Eqs. (32)(33), we find that the coupling applies a relative phase e ig between the identity operator (size zero) and all non-identity Pauli operators (size above) [1]: The lack of an applied phase for non-identity Pauli operators corresponds to the vanishing of V Q at late times, when OTOCs have decayed to zero [see Eq. (33)]. From Section V A, we see that whenever single-qubit teleportation succeeds. A brief argument shows that late time teleportion of higher dimensional quantum states is not possible. Consider teleportation of a d-dimensional qudit, with a basis of states |i , i = 0, . . . , d − 1. The qudit Pauli operators are generated by the 'clock' and 'shift' operators: Z |i = e iω |i , with ω = 2π/d, and X |i = |i + 1 . The two generators obey the commutation relation, XZ = e −iω ZX. After transposition, each Pauli operator, X p Z q , becomes Meanwhile, late time dynamics ensure that the coupling applies an overall phase only to the identity operator. For teleportation to be successful, we would therefore require a decoding operation, D, that acts as DX −p Z q D † ∼ X p Z q . Suppose there was such a unitary operator 9 , and consider its action on the generators: DXD † = X −1 and DZD † = Z. The above action implies that commuting the two generators gives a different phase before and after decoding: DXZD † = e −iω DZXD † = e −iω ZX −1 and DXZD † = X −1 Z = e +iω ZX −1 . This is a contradiction whenever e +iω = e −iω , i.e. whenever d > 2.

B. Equivalence to HPR protocol
We now turn to the equivalence between peaked-size teleportation and teleportation in the HPR protocol. The latter was originally introduced to recover information in the Hayden-Preskill thought experiment [11,23], and is reviewed in detail in Appendix B.
Here, we restrict our attention to teleportation in the deterministic variant of the protocol, of a single qubit 9 The astute reader may note that this operation is in fact implemented by the anti-unitary operator, D |i = |−i mod d * . However, if one decomposes state insertion in terms of Pauli operators as |ψ φ| = Q A c Q Q A (see Section III), one desires that the entire operator |ψ φ| be transferred from left to right for all possible φ|. The preceding anti-unitary operator will complex conjugate the coefficients c Q , thus spoiling teleportation for any φ| where these are complex.
at infinite temperature [11,13]. The protocol takes the form: (52) where P EPR projects onto an EPR pair between subsystems C on the left and right sides. The equivalence between this protocol and the TW protocol [Eq. (15)] is manifest, with the only difference being the locality of the coupling. Specifically, the HPR coupling is of the same general form as the TW coupling [Eq. (1)]: where the sum is over of a complete basis of d 2 C Pauli operators on C. However, the operators P C are typically non-local across C, whereas the coupling considered in the TW protocol was restricted to local operators. As a consequence, the HPR coupling functions as a binary variable measuring whether or not an operator has support on subsystem C (see Section IV). In contrast, the TW coupling measures the operator size within C, which takes an approximately continuous range of values when C is large. Crucially, at late times under scrambling dynamics, the effect of both couplings will be the same: to apply an overall phase to non-identity operators.
A few additional remarks are in order. First, while the leading order effect of the HPR and TW couplings is the same, they lead to different finite-size corrections. In particular, in a fully scrambled system, the variance in the phases applied by the HPR coupling is equal to the probability of a random Pauli string not having support on C, which is suppressed exponentially in the size of C, i.e. 1/d 2 C . On the other hand, the variance in phases applied by the TW coupling is suppressed only polynomially, by ∼ g 2 δS 2 (48) and the discussion below Eq. (34)]. These enhanced phase fluctuations are relevant for finite-size implementations of the TW protocol, as discussed further in Section IX.
Second, it has previously been shown that an extended version of the HPR protocol allows for teleportation of multiple qubits at late times [11]. Because of the equivalence between the protocols, this extension would also allow for multi-qubit teleportation via the peaked-size mechanism. However, the enhanced channel capacity comes with a trade-off: the circuit complexity (measured by the number of applications of the unitary U ) grows exponentially in the number of qubits to be teleported. As we will see in the following section, this limitation can be overcome by peaked-size teleportation in the TW protocol at intermediate times, owing to the locality of the TW coupling.

VII. PEAKED-SIZE TELEPORTATION AT INTERMEDIATE TIMES
We now turn to analyzing the behavior of peaked-size teleportation at intermediate times, i.e. before the scrambling time. In this regime, multiple qubits can be teleported given a certain condition on the growth of timeevolved operators, namely when the overlap of the operators' support is sufficiently small.
We explicitly demonstrate that this condition is satisfied, and multi-qubit teleporation is possible, in a wide variety of physical systems at infinite temperature. These include random unitary circuits (RUCs) in ≥1D, for which peaked sizes naturally occur due to local thermalization within each operator's light cone, and timeevolved operators are non-overlapping due to spatial locality. More surprisingly, we show that multi-qubit peaked-size teleportation can also be achieved in 'fast scrambling', all-to-all coupled systems, including 0D random unitary circuits and the SYK model (at infinite temperature) [18,23]. In this case, operators are not spatially separated at any nonzero time; nonetheless, the overlap of their size distributions remains probabilistically small at sufficiently early times. Furthermore, we demonstrate that while size distributions of local operators are generically not tightly peaked in all-to-all systems, peaked size distributions can be engineered in the TW protocol by encoding one's initial state into large p-body operators.
Finally, we consider the channel capacity-i.e. the maximum number of qubits that can be teleported (allowing both g and t to vary)-of peaked-size teleportation in all-to-all coupled systems. This is an essential question for comparing the capabilities of peaked-size teleportation with those of gravitational teleportation in traversable wormholes [1]. Remarkably, we provide analytic and numerical evidence that the channel capacity of peaked-size teleportation in 0D RUCs, a quite simple microscopic system, is asymptotically equivalent to that of the gravitational mechanism! Namely, the number of qubits n that can be teleported scales with the number of couplings in the protocol, n ∼ K.

A. Multi-qubit teleportation: additive operator sizes
We begin with a few simple examples of multi-qubit teleportation to build intuition. First, consider a unitary U that factorizes as U = U 1 ⊗ · · · ⊗ U n , where each U i acts on a disjoint subsystem. If we insert n qubits individually into the n different subsystems, then the entire protocol decouples into n independent channels and there is no restriction on sending multiple qubits. This trivial example relies on the fact that U does not scramble information across the entire system but only within each disjoint subsystem. We see that full scrambling of information by U in fact inhibits the teleportation protocol's channel capacity (considered for a fixed set of qubits and dynamics).
A similar situation occurs even when the dynamics are not factorizable, as long as the teleported qubits are in causally separated regions. For example, consider a (D ≥ 1)-dimensional system with short-range interactions, where the inserted qubits are spatially separated. At intermediate times, the time-evolved qubit operators will have support within a local 'light cone' about their initial location, but will continue to act on disjoint subsystems. This scenario is therefore no different from the previous example and multi-qubit teleportation remains possible, as long as (i) the size distributions of each operator is tightly peaked, (ii) the coupling V has support within each qubit's light cone, and (iii) the light cones of each qubit are non-overlapping. This final requirement constrains the number of qubits that can be sent at a given time t. In particular, the light cone of each operator will have a radius v B t where v B is the butterfly velocity. The maximum number of non-overlapping light cones-equal to the total number of qubits n that can be teleported-is therefore n N/(v B t) D , where N is the total system volume.
More formally, we can analyze the success of n-qubit teleportation using the two-sided correlators, C Q . We are concerned with n-qubit operators Q(t) = Q 1 (t) . . . Q n (t), where each Q i ∈ {I, X, Y, Z} is a single-qubit Pauli on the i th teleported qubit. We work at infinite temperature and assume that sizes are tightly peaked. Teleportation therefore succeeds whenever all correlators have the same phase.
Inspired by the example of n decoupled protocols, we will take the decoding operator to be the tensor product, D = Y ⊗ . . . ⊗ Y . The combination of transposition and conjugation by D thus applies a minus sign to every single-qubit non-identity Pauli operator. An additional phase is applied by coupling proportional to the size of each operator. For example, for n = 2 qubits, we have: In order for all correlators to have the same phase, we require that η d gS 1 /N = η d gS 2 /N = π mod 2π, and that the operator sizes add, such that e −iη d gS12/N ≈ e −iη d g(S1+S2)/N = e i(π+π) = 1.  In 1D and 2D, sizes grow ballistically in time, while the size width grows with a slower power of t and matches predictions from the KPZ universality class (Section VII B). Because of the separation between the size and size width, the teleportation fidelity for a single qubit exhibits an oscillatory behavior at intermediate times, with nearly perfect maximum fidelity. At late times, the teleportation fidelity saturates close to 1 for odd values of g/π, as expected for any scrambling system (Section VI). (c) In 0D all-to-all coupled RUCs, both the size and size width grow exponentially in time and obtaining a large separation between them requires encoding the initial state into p-body operators. With this encoding, the teleportation fidelity displays a distinct three-regime profile for g 1. In particular, as in 1D and 2D, peaked-size teleportation succeeds (i) at early times, with an oscillating fidelity, and (ii) at late times, where the fidelity saturates close to 1 (for odd g/π). Between these regimes, no teleportation occurs because the size width has grown too large, gδS/N 1.
This requirements generalize straightforwardly to n qubits. Specifically, teleportation succeeds whenever the single-qubit operator sizes obey η d gS i /N = π mod 2π and the multi-qubit operator sizes add under operator multiplication: This latter requirement implies that the phases applied by the coupling, e igV , factorize, and allows the n qubits to be teleported 'in parallel' as in the previous simple examples.
The size addition requirement naturally bounds the channel capacity in terms of the number of couplings, K. Specifically, the K-size takes integer values between 1 and K. However, the requirement that all three singlequbit Pauli operators have the same K-size increases the minimum K-size to 2. From Eq. (54), this implies that an n-qubit operator has a K-size of at least 2n, which is only possible if Indeed, this strict upper bound can also be understood from an information theoretic perspective: teleporting n qubits requires an increase of 2n in the mutual information between the left and right sides of the system. Each of the K classical bits sent from left to right in Fig. 1(a) increases the mutual information by at most 1, so at least 2n bits are required.

B. ≥1D random unitary circuits
As a first concrete example of intermediate time peaked-size teleportation, we consider a random unitary circuit (RUC) applied to a lattice of N qubits in one or higher dimensions. At each time step, pairs of neighboring qubits are evolved via independent Haar random unitaries arranged in a 'brick-layer' fashion, with periodic boundary conditions [ Fig. 2(a,b)]. Operator growth in such systems has been studied at great length, and is believed to be a good model for many aspects of information scrambling under Hamiltonian dynamics [30,49,50,[53][54][55]. We extend these previous studies by demonstrating new results on the behavior of the operator size widthi.e. power-law scaling at intermediate times and suppression at late times-which we show can be detected by the teleportation fidelity (Fig. 3).
A key property of Haar random unitary circuits is that the expectation values of many circuit quantities can be computed by replacing the Haar random unitaries with randomly chosen Clifford unitaries, thereby enabling efficient classical simulation [30,56]. Generally, this equiv-alence holds for any quantity that contains no more than two copies each of U and U † (e.g. the Renyi-2 entropy, or the OTOC); however, for systems of qubits, this property holds for up to three copies [57][58][59]. From Eq. (41), we see that the teleportation fidelity contains three copies of U and U † , so the average fidelity is efficiently simulable 10 . Moreover, by definition, the size distributions of operators under Clifford dynamics are perfectly tightlypeaked, since a Pauli operator Q A evolved under a Clifford unitary remains a single Pauli string. Hence, the teleportation fidelity can be computed using the simplified expression given in Eq. (44).
In more detail, we calculate the average EPR fidelity for teleporting n qubits through the following procedure. First, we choose a particular realization of U by sampling each 2-qubit unitary from a uniform distribution of 2-qubit Clifford unitaries. Second, we determine the K-size of U Q A U † for each n-qubit Pauli operator, Q A , or, if n is large, for a random subset of these operators; such simulations can be performed efficiently with a time cost that scales linearly with the circuit depth. Third, we compute the fidelity for a given coupling g using Eq. (44), where the latter term captures the fact that decoding and transposition apply a minus sign for each non-identity element of the initial Q A . Finally, we average the EPR fidelity over multiple realizations of U .
The results of these simulations for n = 1 qubit in 1D and 2D are shown in Fig. 2(a,b). As expected, the average operator size grows ballistically, S ∝ t D , until the operator's light cone reaches the edge of the system, at which point the size saturates to 3/4N . While the behavior of the size width is more complex, in both dimensionalities it grows more slowly than the average size. This implies that the size distribution is tightly-peaked and the teleportation fidelity can be approximated by (45)]. We verify that the time profile of the fidelity follows this prediction, and nearly perfect fidelity is achieved when η d gS/N = π mod 2π. In Appendix E, we also demonstrate that teleportation of n > 1 qubits is also possible at intermediate times, as long as their light cones do not overlap.
Probing the size width-Let us now turn to the time profile of the size width, which exhibits a peak near the scrambling time in both 1D and 2D. Qualitatively, this behavior arises from fact that the size width receives contributions from two sources: the interior of the light cone, and the boundary of the light cone. Within the light cone, we expect a ≥1D system with a small local Hilbert space to 'locally thermalize' as the operator spreads. This implies that the bulk's contribution to the size width scales as δS bulk ∝ √ S ∝ t D/2 and saturates at the scrambling time. Second, the size width also receives η d gδS/N = 1 The size width initially grows as t 1/2 and reaches a peak at the scrambling time t * ∼ N = 10000. (bottom) We probe this behavior by measuring the teleportation fidelity of a single qubit with a large coupling g = 57π ∼ √ N . The fidelity exhibits a distinct decay-revival profile, controlled by whether the size width has exceeded the threshold gδS/N ≈ 1: nearly perfect fidelity initially, power law decay towards a trivial fidelity at intermediate times, and partial revival at late times.
contributions from the light cone's boundary, which has not yet thermalized. At late times, the boundary of the light cone reaches the edge of the system and these additional contributions subside, leading to the peak in the size width at the scrambling time.
To quantify these effects, we note that the growth of operators in ≥1D RUCs is predicted to fall in the Kardar-Parisi-Zhang (KPZ) universality class [30,60]. In 1D, fluctuations in the light cone boundary have been verified numerically to have a growing width ∼ t α with the KPZ growth exponent α = 1/2 [30]. This implies that the contribution of the boundary to the size width is δS boundary ∝ t 1/2 , and the full width is We note that the maximum size width relative to the late-time size width is a constant set by (α bulk + α boundary )/α bulk . Comparing the size width of multiple system sizes, we observe excellent agreement with predicted scalings over a wide range of system sizes (Appendix E). The time profile of the size width is directly observable in the peaked-size teleportation fidelity if we scale g ∼ t 1/2 scr ∼ N 1/2 . In particular, by setting N/g to lie between the maximum size width and the late time size width, we observe a distinct decay-revival profile for the teleportation fidelity (Fig. 3). At early times, we observe successful teleportation with an oscillating fidelity. The fidelity decays slowly, as a power law in time, as it receives corrections proportional to the growing size variance ∼ g 2 δS 2 /N 2 . After the scrambling time, we see a revival in the teleportation fidelity as the size width narrows. The lack of a parametric separation between the maximum and late time size widths means that late time teleportation will also have some finite error for this value of g.
In 2D, we find that the scaling of the size width also matches predictions from the KPZ universality class. In this case, the width of the boundary scales as ∼ t α , with α = 1/3 [30]. However, to calculate the boundary's contribution to the size width, one must take into account two additional considerations. First, the boundary is 1-dimensional, so its length trivially grows in time as ∼ t. Second, fluctuations of the boundary are expected to have a finite correlation length, ξ ∼ t 1/z , where z = 3/2 is the KPZ dynamic exponent [61]. Thus, the boundary can be modeled as n ξ ∼ t/ξ = t 1/3 uncorrelated regions, each of length ξ. Each region contributes ∼ ξt α to the size width; adding the uncorrelated contributions from all regions yields a total size width δS ∼ √ n ξ ξ t α = t 1/6+2/3+1/3 = t 7/6 .
The time profile of the size width in 2D is thus given by We confirm these scalings in our numerics ( Fig. 2(b) and Appendix E). Notably, the size width is now dominated by the boundary contribution at intermediate times, such that the ratio of the maximum size width to the late time size width scales as t 1/6 scr ∼ N 1/12 . As in 1D, one can probe this behavior using the peaked-size teleportation fidelity, now with g ∼ N/t 7/6 scr ∼ N 5/12 . We emphasize that in 2D, the scaling of the size width is determined by correlations between different points on the light-cone boundary. This goes beyond the behavior studied in previous works on RUCs, which focus on quantities probed by local OTOCs.

C. 0D random unitary circuits
We now turn to random unitary circuits in zero dimensions, a prototypical model for 'fast scramblers' [18,23]. These circuits are constructued as follows: at each timestep, we partition the N qubits into randomly chosen pairs, and apply independent Haar random 2-qubit unitaries to each pair.
Below we analyze such circuits using theoretical arguments, in combination with numerical simulations via Clifford circuits. As the later parts of our analysis are rather technical, we briefly summarize the main results: (i) peaked size teleportation remains possible but only if the input state is initially encoded in non-local, p-body operators; (ii) even though there is no complete separation of operator light cones, size addition still occurs at intermediate times in a probabilistic sense and enables mutli-qubit teleportation; and (iii) the maximum channel capacity is linear in the number of coupled qubits, K. These results are depicted numerically in Fig. 2(c) and 4.
Peaked sizes-In all-to-all coupled systems, operators are generally expected to grow exponentially in time, S ∼ e λt , where λ is the Lyapunov exponent [18]. The reason is simple: at each time step, every term in an operatorrather than just those on a 'light-cone' boundary-has a fixed probability of spreading under random pairwise unitaries. A somewhat less intuitive expectation is that the size width also generally grows exponentially [27]. One way of understanding this is by imagining two realizations of the dynamics: in one realization the initial operator doubles at the first time and in the other it does not. In effect, the latter system now lags behind the former by one time step, ∆t, and the difference in their sizes at later times will be exponentially magnified, to e λt (1 − e −λ∆t ).
The lack of separation between the size and size width seems to preclude the possibility of peaked-size teleportation at intermediate times. Nevertheless, we can engineer such a separation by encoding the information of each input qubit into p-body operators, with p 1 [17]. As an example, consider encoding a single qubit into p = 5 qubit operators via Here, E is a Clifford unitary encoding operation that conjugates state insertion and decoding [explicitly, replacing U → U E, U * → U * E * , and U T → E T U T in Fig. 1(a)]. The success of teleportation is now dependent on the size distributions of time-evolved p-body operators, Q A (t) = U EP E † U † , where P runs over the initial unencoded single-qubit Pauli operators. As we will soon verify explicitly, before the scrambling time the support of each of the p operators composing Q A will be approximately non-overlapping, so that their size distributions will convolve. Thus, the total operator size is multiplied by a factor of p but, through the central limit theorem, the size width is multiplied only by √ p.
In more detail, consider the size growth of an operator, Q A , with initial size S 0 = p. During a single time step, each qubit i in the support of Q A (t) is paired with another random qubit; for simplicity, we assume the second qubit is outside the support of Q A (t), which should be valid at times well before the scrambling time. Under random two-qubit Clifford time-evolution, Q A (t) grows to have support on both qubits with probability ν = 1−2(d 2 −1)/(d 4 −1) (9/15 for qubits). The operator size, S t , therefore grows stochastically in time, according to where each s i is a binary random variable that increases the size by 1 with probability ν and 0 with probability 1− ν, and Bi t (S t , ν) denotes the binomial distribution with S t trials and probability ν, which we can approximate as a normal distribution, N t (νS t , S t ν(1 − ν)). The size at time t can thus be written as a sum of random variables drawn at each time step: from which we see that the average size grows exponentially in time with Lyapunov exponent e λ = 1 + ν. Deviations arise at each time step t , with typical magnitude Since this decays exponentially in t , we can approximate the total variation, δS t , as the largest term in the sum (t = 0), which has magnitude As anticipated, the size width is dominated by early time errors that have exponentially grown in time, so that the ratio of the size width to the size remains constant at ∼ 1/ √ p (after some period of growth from its initial value, 0). To support these claims, we numerically simulate the time-evolved size distribution of operators with an initial size p ≈ 1000 [ Fig. 2(c)]. As expected, we observe that the average size grows exponentially as ∼ pe λt and saturates at a timescale t * ∼ log(N/p). Moreover, the size width grows at the same exponential rate but its magnitude is suppressed by a factor of √ p compared to the average size.
To verify that this allows for teleportation, we next compute the fidelity for teleporting a single qubit, in the regime g 1. As shown in Fig. 2(c), teleportation occurs with near perfect fidelity beginning at t ≈ t * − log(gp), corresponding gS/N ≈ 1. Thereafter, the teleportation fidelity decreases exponentially in time, consistent with the increase of the size width.
teleportation stops succeeding entirely, since the size width has reached the limit δS/N ∼ 1. Finally, at late times t ≈ t * − log(p), the fidelity revives as the system becomes fully scrambled and the operator size width narrows to δS ∼ √ S. Size addition-We now turn to the possibility of teleporting multiple qubits in 0D RUCs. Within the peakedsize regime, this reduces to the question of whether operator sizes add according to Eq. (54). Satisfying this requirement in all-to-all coupled systems is not as trivial as in ≥ 1D, since time-evolved operators typically act on overlapping subsystems at any finite time. Nevertheless, we now provide a simple argument for why size addition holds despite this.
To do so, we model each time-evolved Pauli operator Q i (t) as an independent random Pauli string of size S[Q i ]. Consider two such strings, P 1 and P 2 , with support on regions A 1 and A 2 and sizes S[P 1 ] = |A 1 | and S[P 2 ] = |A 2 |. The size of the product, P 1 P 2 , is the size of the union A 1 ∪ A 2 , minus the number of sites where the two strings overlap and have the same single-qubit Pauli operator. This occurs with probability 1/(d 2 − 1) = 1/3 at each site in the region A 1 ∩ A 2 , giving The deviation from the simple additive rule S[P 1 If the Pauli strings P 1 , P 2 have independently random areas of support, the size of this intersection scales as: which is subleading to S[P i ] at intermediate times (when S/N 1). To derive this, note that the probability for both strings to have support on a given qubit is ∼ (S[P 1 ]/N )(S[P 2 ]/N ); summing over N qubits gives the above result.
For n-qubit teleportation, one must consider the combined size, S[P 1 . . . P m ], of m independent Pauli strings, where m takes a typical value m ≈ 3n/4 (a typical nqubit operator has non-identity support on 3n/4 qubits). In general, this quantity will receive corrections from m k different k-way intersections of the strings, for all 2 ≤ k ≤ m. For random Pauli strings, the expected size of these intersections scales as where S ∼ |A i | is the typical size of a single Pauli string [see Eq. (63) above]. For a given k, the correction to size addition will be the sum of m k ∼ m k different intersections and therefore scales as mS(mS/N ) k−1 . These corrections can be neglected if they are small compared to the total size; this occurs when mS N , which corresponds to a timescale much less than the scrambling time.
To demonstrate this claim, we numerically simulate the teleportation protocol with n > 1 qubits in the regime 1 p, np K [ Fig. 4]. Analogous to single-qubit teleportation, the teleportation fidelity exhibits oscillations beginning at t ≈ t * − log(gp), and vanishes at t ≈ t * −log g √ pn due to the growth of the combined size width. However, in contrast to the single-qubit case, teleportation of multiple qubits is not possible at late times, t t * − log(gpn), as predicted in Section VI. Interestingly, between these two regimes, we observe a partial revival of the fidelity: this indicates that the operator size widths begin to narrow before the additive condition is completely invalidated.
Error analysis-While we have confirmed that multiqubit teleportation can be achieved in certain ideal limits, a key question remains: how does the maximum number of qubits that can be teleported scale as a function of K, i.e. what is the protocol's channel capacity? To answer this question, we now estimate how deviations from these ideal limits lead to errors in peaked-size teleportation and ultimately constrain the channel capacity. Throughout this discussion, we assume that the size, S, is extensive, but K is not; this is the natural regime for probing the channel capacity of the protocol at intermediate times, and is the physical scenario in the context of traversable wormholes [1]. The details of this and the following subsection are quite technical in nature, and may be skipped by most readers.
In summary, we identify four distinct sources of error in the multi-qubit teleportation fidelity, F = 1 − : 1. Errors due to finite p: ∼ ng 2 S 2 K /K 2 p 2. Errors due to finite K: ∼ ng 2 S K /K 2 3. Errors due to imperfect size addition: where ellipses indicate higher orders in (nS K /K) 2 4. Errors due to fluctuations in size addition: where ellipses indicate higher orders in nS K /K We discuss each of these errors in detail below.
The first and second sources of error are due to imperfectly peaked K-size distributions. The K-size width receives contributions from finite-p corrections, ∼ S K / √ p, To translate these into errors in the teleportation fidelity, we multiply the size width by g/K and take the square. This gives fidelity errors ∼ g 2 S 2 K /pK 2 and ∼ g 2 S K /K 2 per teleported qubit.
The third and fourth sources of error arise from imperfect size addition. This leads both to 'systematic' errors, due to the average overlap of operators, as well as 'sampling' errors, due to random fluctuations in this overlap. We begin with the systematic errors: as we recall, the size addition of m time-evolved operators receives corrections from k-way overlaps of the operators, each scaling as ∼ mS K (mS K /K) k−1 , for 2 ≤ k ≤ m (rescaling our previous results to the K-size instead of the size). The nonlinear dependence on m indicates that sizes do not add perfectly. Nevertheless, when teleporting an n-qubit initial state for large n, we can correct for the above effect at leading order by using a linear approximation for m k about its typical value, (3n/4) k . This leads to an effectively smaller operator size, which can be observed in the reduced frequency of the fidelity oscillations for 10-qubit teleportation compared to 1-,3qubit teleportation in Fig. 4(a). The leading errors after this shift are quadratic in δm ≡ m − 3n/4, which has a typical magnitude δm ∼ √ n. Multiplying by g/K and taking the square, we therefore find multi-qubit fidelity errors ∼ (gS K /K) 2 (nS K /K) 2k−2 ; at leading order k = 2, this gives ∼ n 2 g 2 S 4 K /K 4 . Finally, each intersection above is subject to additional random fluctuations about its average value. When operator sizes are much smaller than the system size, we can treat each intersection as arising from a binomial process, in which case fluctuations are proportional to the square root of the intersection's average size (see Appendix F for a detailed accounting). These add in quadrature for ∼ n k overlaps, producing a total fidelity error ∼ (g 2 /K)(nS K /K) k .
Channel capacity-To define the channel capacity of the teleportation protocol, we fix a per qubit error threshold th , and determine the maximum number of qubits that can be sent while maintaining a multi-qubit fidelity above this threshold 11 , i.e. F ≥ 1 − n th . We are interested in how the channel capacity scales with the number of couplings, K, while allowing both g and S K (determined by the evolution time) to vary.
In 0D RUCs, all errors increase with g, so it is optimal to set g to its minimal value, η d gS/N = π. This gives a per qubit error The first term is negligible in the large p limit and so we will neglect it from here on. We minimize the remaining terms with respect to S K . There are two relevant regimes. For n √ K, the minimum is determined entirely by the leading order contributions in nS K /K to the error (i.e. neglecting the ellipses). Taking the derivative and setting to zero, we have the minimum at S (1) K ∼ K 2/3 /n 1/3 . As we increase n, the optimal size approaches the value S (2) K ∼ K/n. At this point, size addition errors of all orders (i.e. the ellipses) become large, and so the true minimum becomes fixed just below S (2) K . This crossover between these two minima occurs at n ∼ √ K, at which S K . The above minima give two distinct scalings for the per qubit error and thus the channel capacity. The first minimum has a per qubit error (1) /n ∼ (n/K 2 ) 1/3 , which gives rise to a superlinear channel capacity, n 3 th K 2 . However, as we increase K, this capacity eventually surpasses the value √ K. Above this, the optimal size is given by the second minimum, which has an error (2) /n ∼ n/K, and thus the channel features an asymptotically linear capacity, This is a stronger instance of the strict general bound Eq. (55). Intuitively, this channel capacity arises because the individual K-sizes must be large, S K 1, for the K-size to be tightly peaked, while at same time the combined K-size must be much smaller than K, nS K K, for the K-sizes to add; hence n K.
We test this scaling numerically by simulating the teleportation protocol and measuring the per qubit fidelity, F (1) EPR , as a function of n and K. Specifically, for each value of K, we sweep the number of qubits n and determine the maximum qubits that can be sent before the infidelity exceeds a threshold, 1 − F (1) EPR = th . These results are shown in Fig. 4(b) and exhibit a clear linear trend across two orders of magnitude, confirming our prediction of a linear channel capacity.
A few final remarks are in order. First, while in principle the per qubit fidelity can be calculated by taking the n th root of the full n-body fidelity, this approach is numerically unstable for large n. Thus, we instead compute the fidelity of a single qubit, while trying to send multiple qubits, using an approach derived in Appendix E. This amounts to performing a sum analogous to Eq. (44), but only including pairs of Q 1 and Q 2 that are equal on all sites except for one.
Second, the range of system parameters that lie within the linear scaling regime is ultimately constrained by the finite total system size, N = 10 8 . In particular, to maximize the linear scaling regime, we choose p = 101 and th = 0.07. The former ensures that finite-p errors are negligible, while the latter allows the number of qubits at the threshold to be large enough to access the n √ K regime but small enough that the operators are initially dilute, i.e. n N/p.

D. Large-q SYK model: infinite temperature
We now demonstrate peaked-size teleportation in a 0D Hamiltonian system, the large-q SYK model, at infinite temperature. While teleportation at low temperatures in the SYK model is known to succeed via the gravitational mechanism, teleportation at infinite temperature was discovered only recently [17]. In addition to showing that this mechanism is in fact peaked-size teleportation, we also find that, remarkably, all qualitative aspects of this teleportation match those of 0D RUCs.
The large-q SYK model is defined by the Hamiltonian [25,27]: where ψ i are Majorana fermions, {ψ i , ψ j } = 2δ ij , and the couplings are drawn independently from a Gaussian distribution with zero mean and a variance J 2 j1,...,jq = J 2 /2q N −1 q−1 . This model is exactly solvable at all temperatures in the large-q, large-N limit [25,27].
To construct the teleportation protocol for the SYK model, we first define the N -fermion EPR state, ψ j,l |FEPR ≡ −iψ j,r |FEPR , ∀ j = 1, . . . , N (67) From this, the TFD state is obtained as before, ∼ K F (1) EPR Figure 4. Teleportation of multiple qubits in 0D RUCs. (a) Many-body teleportation fidelity, FEPR, as a function of time for teleporting n = 1, 3, 10 qubits with fixed coupling strength (g = 177π). Compared to a single qubit, the decayrevival profile for multiple qubits is shifted to earlier times, since multi-qubit operators both have a larger size width and saturate the system size earlier. Moreover, multi-qubit teleportation is not possible at late times, resulting in a trivial late-time fidelity (Sec. VI A). (b) Numerical results for the channel capacity nmax as function of the number of coupled qubits K, which exhibit a clear linear scaling. To determine the channel capacity, we compute the maximum per qubit fidelity F EPR for a fixed number of qubits, n, and couplings, K, while allowing the coupling strength, g, and evolution time to vary. For fixed K, F EPR decreases as the number of qubits n is increased, as depicted in the inset for K = 9000. The channel capacity nmax is defined as the maximum number of qubits for which the fidelity is above a fixed threshold (dashed line).
For the two-sided coupling, we consider the simple bilinear interaction, which measures the size of operators in the Majorana string basis, divided by qN [26,27]. As in 0D RUCs, the size and size width of timeevolved operators in the SYK model increase exponentially in time, and exhibit a large separation only when initially encoded in p-body operators. To see this, we can generalize previous computations of size distributions in the large-q SYK model [27] to initial p-body operators, ψ = ψ 1 ψ 2 . . . ψ p ; this relies on the factorization of SYK correlation functions in the large-N limit [17]. After the relaxation time (t 1/J), but before the scrambling time (t log(N/p)/J), the size and size width are: The scaling δS ∼ S/ √ p matches that found for 0D RUCs; in particular, ensuring a large separation between the size and size width requires p q. Note that our condition for peaked size distributions depends on the (large) parameter q, through the size width.
This large separation suggests that peaked-size teleportation is possible at early times in the large-p limit.
To verify this, we analyze the two-sided correlator, which is given by [8] C ψ (t) = e −igV ψ r (−t)e igV ψ l (t) at infinite temperature before the scrambling time 12 . For large p and early times, we can approximate the correlator as using We refer to this regime as the "early time regime", and analyze its analog in large-N systems at finite temperature in Section VIII A.
Crucially, as expected for peaked-size teleportation, the early time correlator consists of an overall phase equal the average operator size, Eq. (70), multiplied by g/qN . This indicates that teleportation succeeds with nearly maximal fidelity beginning when gS/qN ≈ 1. Based on its similarity with 0D RUCs, we expect that teleportation in this regime is capable of teleporting O(K) qubits (Table I); however, we do not calculate this explicitly. Teleportation continues to succeed until the above approximation breaks down, which occurs when the size width, δS, becomes of order (g/qN ) −1 . As for all scrambling systems, the two-sided correlator is expected to revive at late times, t log(N/p)/J, at which point the sizes saturate the entire system [1,8] (see Section VI); this is not reflected in Eq. (71), which is valid only before the scrambling time.

VIII. INTERPLAY BETWEEN PEAKED-SIZE AND GRAVITATIONAL TELEPORTATION
In this section, we seek to understand the interplay between peaked-size and gravitational teleportation. A central theme in this understanding is a comparison between the size distribution introduced in Section IV, and the winding size distribution introduced in Ref. [15,16].
To illustrate the distinction between these distributions, consider a time-evolved Majorana fermion operator, decomposed in a basis of Majorana strings, χ [26,27]: From this decomposition, one defines the size distribution [26,27], and the winding size distribution [15,16], where S[χ] is the size of the string χ. Note that the size distribution is real-valued, while the winding size distribution may be complex. The teleportation correlators [under coupling Eq. (69)] are, in fact, directly related to the winding size distribution [15,16]: which can be derived by explicitly plugging Eq. (73) into the teleportation correlator. The size distribution, by contrast, is related to "one-sided" correlation functions, e.g. Eq. (31), where both instances of the time-evolved operator appear on the same side of the TFD state [27]. Despite this distinction, we have so far been able to analyze teleportation using the size distribution, as opposed to the winding size distribution, because the two are equal in two circumstances. The first is at infinite temperature, where the coefficients c χ are real because ψ(t) is Hermitian. The second has been precisely our focus: when size distributions are perfectly tightly peaked, in which case both distributions approach a delta function.
In what follows, we describe several scenarios in which the distinction between the two distributions becomes relevant. First we begin in large-N systems, where large-N factorization provides a precise relation between the teleportation correlator and the OTOC at early times. We find that, even in the presence of the large-p encoding, the correlator deviates from the peaked-size prediction whenever the OTOC contains an imaginary part. Large-N systems encompass both peaked-size and gravitational teleportation-our results suggest that the former occurs in systems where the OTOC is real (e.g. at infinite temperature with large-p encoding, see Section VII), while the latter occurs where the OTOC is imaginary (e.g. at low temperature in SYK) [33,45]. Second, we review recent results showing that this deviation eventually leads an O(1) correlator magnitude when the winding size distribution takes a particular form, thereby enabling teleportation with unit fidelity (see Section III). This is conjectured to be the microscopic origin of gravitational teleportation [15,16], and so we expect it to occur only in low temperature models with a gravity dual. Third, we return to teleportation in the large-q SYK model and show that this model interpolates between gravitational teleportation at low temperatures and peaked-size teleportation at high temperatures. Surprisingly, this interpolation occurs despite the fact that the large-p encoding ensures a large separation between the size and size width, i.e. the size distribution naively appears tightly peaked, even at low temperatures. Finally, motivated by this smooth interpolation, we conclude this section by searching for a 'dual' description of peaked-size teleportation in a bulk gravitational theory. In particular, we argue that strong stringy effects lead to the same qualitative features as peaked-size teleportation.

A. Early time teleportation in large-N systems
In Section IV, we saw that for peaked-size operators the teleportation correlator depends only on the first moment of the size distribution, i.e. the average size [Eq. (31)]. We will now show that a more general relationship holds for large-N systems at early times, where we substitute the average size with the first moment of the winding size distribution. Specifically, using Eqs. (9)(10)(11)(12)(13)(14), the first moment of the winding size is given by a two-sided OTOC: (77) using Eqs. (9)(10)(11)(12)(13)(14). This differs from the one-sided OTOC, for probing the average size [Eq. (31)], in terms of the placement of the thermal density matrix.
To relate the OTOC and the teleportation fidelity, we consider two simplifying assumptions. First, we focus on 0D large-N systems, e.g. the SYK model, with a pbody initial encoding. In such systems, the teleportation correlator in fact factorizes into a product of single-body correlators (up to 1/N corrections) [17]: where ψ 1 is a single-body operator. Second, generalizing Eqs. (34,72), we consider sufficiently early times to work at leading order in g 13 : where G β = i ψ 1,r ψ 1,l = tr ρ 1/2 ψ 1 ρ 1/2 ψ 1 is the imaginary time Green's function, and F 2 (t) is the first-order, connected component of the two-sided OTOC [Eq. (77)], Similar to Eq. (34), the leading correction to Eq. (79) is , and the approximation holds when this is small.
Let us now consider the behavior of the teleportation correlator, Eq. (79), under different physical scenarios. We focus on chaotic systems during the so-called Lyapunov regime, which occurs between the thermalization time, t ∼ O(1), and the scrambling time, t ∼ O(log N ). In this regime, the connected OTOC is characterized by a simple exponential F 2 (t) ∼ e λt with a prefactor that is generally complex. As a result, we expect the teleportation correlator to exhibit two distinct effects: (i) the real part of F 2 (t) causes rapid phase oscillations in the teleportation correlator, while (ii) the imaginary part increases/decreases the teleportation correlator magnitude, depending on the sign of the coupling g.
At infinite temperature, F 2 (t) is strictly real and thus only effect (i) can occur. Indeed, in this case, the twosided OTOC directly measures the operator size and Eq. (79) is equivalent to Eq. (34). It follows that peakedsize teleportation can be achieved with perfect fidelity: the teleportation correlator magnitudes are equal to one due to the infinite temperature, and their phases can be aligned by tuning g or t. More generally, at finite temperature, F 2 (t) contains both a real and imaginary part, and the real part-which leads to effect (i)-is formally distinct from the first moment of the size distribution. Rather, recent work has shown that Re{F 2 (t)} is computable via a ladder diagram identity and is physically interpreted as a 'branching time' [33,62]. Here teleportation is similarly possible by tuning g or t to align the correlator phases, however the teleportation fidelity is bounded from above if the correlators do not have magnitude one (Section III).
At the opposite extreme, effect (ii) is dominant in systems with a gravity dual [33,45] (as well as other maximally chaotic systems, e.g. maximally chaotic 2D CFTs with a large central charge [63]). In such cases, F 2 (t) is mostly imaginary and leads to the growth (or decay) in the magnitude of the correlator. This opens the door to magnitudes greater than the two-point function, |C ψ (t)| > G β , which is not possible in peaked-size teleportation (Section V B). Interpolating between the two above limits, it has been conjectured that the prefactor of F 2 (t) is proportional to e iλβ/4π [33,45]. This would imply that the imaginary part is dominant if and only if λ ≈ 2πβ, i.e. the system approaches the bound on chaos [21].

B. Gravitational teleportation and the size-winding mechanism
We now move beyond early times and provide a brief review of how the correlator can achieve its maximal magnitude, 1, at finite temperatures. This occurs via the 'size winding' phenomenon introduced in Ref. [15,16] as the microscopic mechanism for gravitational teleportation [1,8]. We refer the reader to Ref. [16] for a complete discussion of this mechanism, including its connection to physical quantities in the bulk gravity theory. As we emphasize in Section III, maximizing the magnitude of the correlators is necessary for high fidelity teleportation, but it is not sufficient: we must also align the correlator phases, for every operator on the subspace to be teleported.
To begin, note that the winding size distribution is normalized to the two-point function, G β ≤ 1, in contrast to the size distribution, which is normalized to 1. From Eq. (75), we see that this norm being less than one implies that the phases of the coefficients c χ are not perfectly aligned for different strings χ. It is convenient to separate this misalignment into two classes: first, when coefficients of strings of the same size S are misaligned, which manifests in the magnitude of f (S) being less than maximal for a given S, and second, when the phases of f (S) for different sizes S do not align with each other.
We focus on the latter case and, more specifically, consider an ansatz in which the coefficients' phases wind with the size [15,16]: In this case, the coupling of the teleportation protocol, by applying a phase that is also proportional to the size, can serve to unwind the phases of f (S) at the value g/N = −2α [see Eq. (76)]. This increases the teleportation correlator magnitude from its initial value, G β , to unity. Although seemingly artificial, in the following subsection we show that this ansatz holds exactly for the SYK model at low temperatures.
C. Large-q SYK model: finite temperature We now turn to explore the interplay between peakedsize and gravitational teleportation in an explicit example: the large-q SYK model at finite temperature and large-p encoding [27]. Despite the fact that this model features a large separation between the size and size width, we show that teleportation is not governed by the peaked-size mechanism at low temperatures, due to the presence of strong size winding.
To begin, let us consider the finite-temperature teleportation correlator, given by [17]: where (G β ) p = i p ψ r ψ l = (λ/2J) 2p/q is the p-body two-point function, and the Lyapunov exponent λ corresponds to the solution of βλ = 2βJ cos(λβ/4) and interpolates between 2π/β at low temperatures and 2J at high temperatures. At infinite temperature, the correlator reduces to Eq. (71), and follows our expectations for peaked-size teleportation (see Section VII D). At low temperatures, where the model is known to possess a gravitational dual [24,25,45], the correlator behaves substantially differently; most notably, its magnitude increases from G p β at time zero to unity when gJe λt /2λN = 1 [illustrated in Fig. 1(c)].
From this correlator, we can verify the two predictions made in Sections VIII A and VIII B: (i) the early time behavior is governed by the two-sided OTOC, and (ii) the size winding mechanism is responsible for the O(1) peak in the correlator magnitude at low temperatures. To see the former, we expand the correlator in the early time regime: Indeed, the term in the exponent is directly proportional to the connected piece of the two-sided OTOC [33], matching Eq. (79) 14 . At high temperatures this OTOC is equal to two times the operator size [Eq. (70)], resulting in phase oscillations, whereas at low temperatures the OTOC rotates to become predominantly imaginary, leading to an exponential growth in the correlator magnitude. Next, to understand the role of size winding, we must analyze the full winding size distribution. We can derive this distribution by expanding the teleportation correlator in powers of e −ig/qN to match Eq. (76) [15,16,27]. To do so, it is convenient to consider the exact correlator (before a g/N 1 approximation) [17,27]: Rewriting this correlator using Eq. (83) and the Taylor expansion, and identifying the n th coefficient with the winding size distribution, we have: At intermediate times and large p, the distribution takes a particularly simple form, where we define the size decay rate, γ, as 14 More precisely, the correlator in Eq. (84) is missing a factor of G p β compared to Eq. (79). This same mismatch is noted in Ref. [27], and is attributed to the large-q limit utilized for the calculation, since in this limit G β approaches 1.
The fact that the correlator magnitude increases in time, and moreover reaches an O(1) value at low temperatures, is a hallmark of gravitational teleportation and signals physics outside the peaked-size regime. Naively, this result is surprising, as we expect the p-body encoding to ensure a peaked size distribution. Indeed, the average size and size width remain separated by √ p at all temperatures [27]: This demonstrates that our simple intuition, of judging a size distribution to be tightly peaked if the ratio between the size width and average size is small, is not always correct. Rather, in Appendix A, we provide a more precise condition for when peaked-size teleportation holds, and explicitly show that this condition breaks down for the SYK model at finite temperature (but remains satisfied at infinite temperature). Let us now provide intuition for how peaked-size teleportation is modified by size winding at low temperatures. To this end, we express the SYK correlator in terms of the winding size distribution parameters: At early times, this integral can be solved using a saddle-point approximation. At infinite temperature, the saddle point, n s , occurs precisely at the average size, n s = (2p/q)/γ = S/q, giving the peaked-size correlator, C ψ = (−iG β ) p · exp(−igS/qN ). In contrast, at finite temperature, the size winding α shifts the saddle point in the imaginary direction of the complex plane, giving n s = (2p/q)/(γ + 2iα) and a correlator C ψ = (−iG β ) p · exp(−ign s /qN ). From this, we recognize the saddle point as precisely the two-sided OTOC, n s = p 2q F 2 (t).
The inclusion of the size winding in the low temperature saddle point thus has two effects. First, it contributes an imaginary part to the OTOC and thereby increases the magnitude of the teleportation correlator. More subtly, it also alters the real part of the OTOC. At low temperatures, α/γ ≈ βJ 1, and we can approximate the saddle as n s ≈ (2p/q)/(2iα) + (2p/qγ)(γ/2α) 2 . Recognizing S = 2p/γ, we see that the real part of the OTOC now corresponds to the average size suppressed by two factors of the ratio (α/γ) 2 .

D. Gravity with stringy effects
While the bulk of this paper approaches teleportation firmly through the lens of quantum mechanics, we would be remiss not to explore the analog of peaked-size teleportation in gravitational physics. Specifically, we would like to ask: is there a teleportation mechanism in gravitational systems that shares the same features as peakedsize teleportation? Such a connection might seem surprising, given the prevalence of peaked-size teleportation in quantum mechanical models with no apparent connection to gravity. Nonetheless, the smooth blending between gravitational teleportation and peaked-size teleportation in the SYK model suggests a positive answer.
Here, we demonstrate-in a particular gravitational geometry, AdS 2 -that an analog of peaked-size teleportation indeed occurs when strong stringy corrections [1,32] are included in the gravitational theory 15 . Intuitively, our results are consistent with our previous analysis of the SYK model, where, in the dual gravitational theory, increasing the temperature is known to add stringy effects [45].
Our derivation closely follows that of Ref. [1] and assumes a background familiarity with the gravitational description of teleportation in AdS 2 (a thorough summary of which can be found in the seminal works of Refs. [1,8]). In this setting, the teleportation correlator can be calculated explicitly by considering gravitational scattering in a wormhole geometry [ Fig. 5]. We will maintain our SYK notation, so that V consists of K single-body fermion operators, ψ i , and our input operator is a p-body fermion, ψ. The correlator can be solved for by decomposing the fermion operators in a momentum basis and applying the scattering matrix: where Ψ l/r (k, t) is the wavefunction for the p-body operator inserted on the left/right boundary with in-falling 15 We are grateful to Zhenbin Yang and Douglas Stanford for discussions leading to this connection. momentum k (and similarly Ψ 1,l/r (s, 0) for any singlebody operator in V ), and e iδ(k,s) is the scattering matrix element between ψ(t) and ψ 1 (0). In pure gravity, i.e. in the absence of stringy effects, these quantities take the form [1]: where we have set β = 2π for convenience, Θ(x) is the Heavyside function, and ∆ = p/q is the conformal weight of ψ. The single-body wavefunction, Ψ 1 (s, 0), is obtained by setting t = 0 and replacing ∆ → ∆ 1 = 1/q (i.e. the conformal weight of a single fermion).
In the semiclassical limit, we can evaluate the correlator by expanding e iδ to linear order in G N [1]. We find: whereg ≡ g4 −∆1 ∆ 1 /2. This expression is almost identical to the large-q SYK correlator of Eq. (94), setting the size decay rate to zero, γ = 0, and identifying the momentum k in the gravitational calculation with the size n in the SYK model [64]. Notably, the correlator diverges at the teleportation time, 4 =gG N e t . In bulk gravity, this divergence is exactly the light-cone pole between the left and right sides of the traversable wormhole, and is regulated by including higher order terms in G N or stringy corrections [1]. While the full effects of stringy scattering in an AdS background are not known, we will take a phenomenological treatment as in Ref. [1,32]. Here, the total effect of stringy corrections is to change the scattering amplitude to where ε controls the strength of stringy effects, and varies from 1 in pure gravity to 0 in the highly stringy limit.
Again expanding e iδ to leading order in G N , and Wick rotating k → −ik, we can write the correlator as where A ε is a constant of order 1. Note that the kdependence in front of exponential is a Poisson distribution with a saddle point at k s ≈ ∆/2 in the heavy particle limit, ∆ = p/q 1. At early times, e εt G N 1, and for strong stringy effects, ε → 0, the change in this saddle Ψ r (t) Figure 5. Schematic of the teleportation protocol from the bulk gravitational perspective in AdS2, under both (a) semiclassical gravity, and (b) strong stringy corrections. The TFD state corresponds to a two-sided black hole. Local quantum mechanical operators, ψ l/r , create or annihilate particles near the two boundaries, with wavefunctions Ψ l/r (red). The protocol begins by inserting a particle on the left side, with wavefunction Ψ l (red, bottom left), at time −t, which then falls towards the interior of the geometry during time-evolution (red line). The two-sided coupling, g N i ψ i,l ψi,r, is then applied, producing a shock wave (blue) that interacts with the in-falling particle [1,8]. (a) In the semiclassical limit, the shock wave shifts the position the in-falling particle outside of the right horizon (dashed), which enables the particle to reemerge near the right boundary (red, top right) [1,8]. (b) When stringy effects are present, the scattering amplitude between the in-falling particle and the shock wave is modified according to Eq. (99) [1,32]. In the highly stringy limit and at early times, the interaction results in an overall phase shift, θ = gGN Aε(∆/2) ε e εt [Eq. (101)]. The overlap between the in-falling particle and a particle at the right boundary is nevertheless non-zero (red, top right), and is given by the unperturbed two-point function, G β = i ψ l ψr . [Note that stringy effects may also modify the initial wavefunctions of Ψ l/r , as we discuss in the context of Eq. (103).] point from the scattering, g, is negligible. In these limits, the saddle point approximation thus gives the correlator: which has exactly the same form as in peaked-size teleportation [Eq. (37)] 16 ! Specifically, the correlator is equal to the two-point function, G β = i ψ l ψ r , multiplied by a pure phase. Tentatively, this suggests interpreting the phase as the operator size in a dual boundary theory. This size, grows exponentially in time with a non-maximal Lyapunov exponent, 2πε/β. A few remarks are in order. First, while in the above treatment the strength of stringy effects depends on a 'free' parameter ε, we expect that in a UV complete theory ε would in turn depend on the temperature (and other physical parameters). In particular, we expect ε → 1 at low temperature in theories that are dual to pure gravity, and ε → 0 at high temperature, where stringy, UV effects should play an important role. This statement also follows from the point of view of the boundary field theory, since the scattering matrix is proportional to an OTOC of the boundary theory, which is real at infinite temperature. 16 Note that the phase in Eq. (101) becomes order-one within the Lyapunov regime, i.e. t 1/ log(1/G N ), but at sufficiently early times to satisfy G N e εt 1. These conditions are consistent as long as ∆ = p/q is sufficiently large to ensure Aε(∆/2) ε 1.
Second, if we would like to recover the infinite temperature SYK correlator, Eq. (71), from the scattering computation, choosing a proper ε as a function of β is not enough. One also needs to modify the wavefunction of ψ, to: Such a wavefunction modification due to UV data should be model dependent, and it would be interesting to understand how to derive this 'stringy-corrected' wavefunction from the bulk point of view. Nevertheless, one particular feature of the modified wavefunction has a clear motivation from the boundary perspective. Specifically, Wick rotating Eq. (103), k → −ik, leads to a distribution whose width, δk ∼ ∆ 1/ε , broadens as ε → 0. This broadening increases the phase variations in the exponential of Eq. (100) and results in the decay of the correlator at the timescale e εt G N / √ ∆ ≈ 1 for small ε. From the boundary point of view, this decay corresponds to the requirement that the size width must be small, gδS/N 1, for peaked-size teleportation, as we saw for 0D RUCs and infinite temperature SYK (Section VII). We expect this decay to be common to many 0D quantum systems at high temperatures, which suggests that the broadening of the bulk stringy wavefunction as ε → 0 might also be a general feature.
Finally, the most obvious effect of a non-unity ε is to change the scattering phase, δ(k, s), from being realvalued to complex. Indeed, in the strong stringy limit, δ(k, s) becomes purely imaginary. In general scattering theory, a complex δ means that the scattering matrix, e iδ , is no longer normalized, and implies the existence of inelastic scattering [32]. Since peaked-size teleportation is replicated in the limit ε → 0, this suggests a more general relationship between peaked sizes and inelastic scattering. In Appendix H, we demonstrate that these two phenomena also coincide at infinite temperature, for arbitrary wavefunctions and scattering amplitudes.

IX. EXPERIMENTAL PROPOSALS
Having illustrated the wide breadth of physics that enters into the TW protocol, in this section we outline explicitly how one can probe this physics in the laboratory. We begin with a summary of the key signatures of teleportation, and how they can be applied towards (i) characterizing operator size distributions in generic scrambling dynamics, and (ii) distinguishing generic vs. gravitational scrambling dynamics. For (i), we show that the TW protocol can be simplified dramatically at infinite temperature, where an equivalent 'one-sided' protocol eliminates the need to experimentally prepare the thermofield double state. We next present two near-term experimental realizations of the protocol: first with neutral atoms and second with trapped ions. The fundamental requirement is the ability to time-evolve forwards and backwards under many-body scrambling dynamics; recent experimental progress has demonstrated this in a number of quantum simulation platforms [65][66][67][68][69]. We conclude with a discussion of the effect of experimental error, and a comparison of the TW protocol with other diagnostics of scrambling physics.

A. Signatures of the TW protocol
We begin by reviewing the key signatures of the TW protocol, as discussed in the previous sections and summarized in Table I. We first recall that the simplest experimental signal-that is, any non-trivial teleportation fidelity of a single qubit-has already been demonstrated experimentally in the closely-related HPR protocol [13,14]. As discussed in Section VI, this signifies that the implemented unitary is scrambling but does not distinguish between peaked-size or gravitational teleportation. In what follows, we discuss two more refined applications of the TW protocol.
Characterizing size distributions in generic scrambling dynamics-The dynamics of the teleportation fidelity within the TW protocol can be used to probe the size distributions of time-evolved operators. This approach relies on the peaked-size teleportation mechanism and thus applies to generic scrambling systems, including the examples analyzed in Section VII (e.g. RUCs, spin chains, high T SYK).
Specifically, the teleportation fidelity as a function of time exhibits three relevant features. First, since peakedsize teleportation relies on the width of the size distribu-tion being small, gδS/N 1, its success or failure indicates whether the width has surpassed the tunable value, N/g. Depending on the model and the value of g, this leads to a temporal profile that exhibits three regimes: initial teleportation when the size width is small, no teleportation when δS N/g, and late time teleportation once the size width converges to its small final value in a finite-size system [as depicted schematically in Fig. 1(c) and observed numerically in 0D RUCs in Fig. 2(c)].
Second, within the peaked-size regime, oscillations in the teleportation fidelity as a function of time, F =  (45)], provide a direct measurement of the growth in operator size. In particular, setting g = 2πn + π, one expects to see n oscillations in the teleportation fidelity before it reaches its late time plateau. The peaks in these oscillations give the operator size as a function of time: S = (m/n)(1 − 1/d 2 )N at the m th peak.
Third, the teleportation of multiple qubits demonstrates the equivalent channel capacities of peaked-size and gravitational teleportation (Section VII). Formally, multi-qubit teleportation probes whether the sizes of time-evolved operators add under operator composition. While this is trivial when the operators are causally separated, determining the requirements for size addition under more general dynamics-e.g. all-to-all or power-law interactions-remains an open question 17 .
Distinguishing gravitational scrambling dynamics-The TW protocol can also be used as an experimental litmus test for gravitational dynamics. To this end, we propose to use two experimental signatures that distinguish between gravitational and peaked-size teleportation: (i) the teleportation fidelity at low temperature, and (ii) the behavior of the teleportation fidelity as a function of time, t, and the coupling strength, g. For (i), the observation of a high teleportation fidelity, ∼ O(1), at low temperatures strongly suggests the occurrence of gravitational teleportation, since the fidelity of peakedsize teleportation is limited at such temperatures by the (small) two-point function, G β . For (ii), one observes that the qualitative profile of the teleportation fidelity as a function of time differs between the two mechanisms (see Fig. 1(c) for a comparison between the two, and Figs. 2, 3 for additional examples of peaked-size teleportation). Namely, keeping g fixed, the fidelity of gravitational teleportation is expected to display a single peak as a function of time, whereas the fidelity of peaked-size teleportation is highly oscillatory in time. Furthermore, gravitational teleportation works only for a specific sign of the coupling, g > 0, while the peaked-size teleportation fidelity is an even function of g [1,8,15,16].
Contrasting with finite-size effects-Finally, we would like to distinguish many-body teleportation from spuri-ous effects that may be seen in the TW protocol at smallsize systems. The most effective way to avoid such signals is by utilizing a coupling gV [Eq. (1)] whose individual terms have a small magnitude, i.e. g/K 1; this is most naturally achieved by including many couplingswhich requires a sufficiently large system-and setting g ∼ O(1). In this limit, the action of the coupling is negligible unless local operators have grown significantly under many-body dynamics, i.e. S ∼ K/g 1 (see Section IV); any teleportation signal is thus necessarily a result of scrambling dynamics. Furthermore, we expect large-size operators to generically exhibit smooth size distributions, justifying our approximation (Section IV B) that the teleportation fidelity is governed by the distributions' first few moments.
Away from this limit, our general framework relating the teleportation fidelity to operator size distributions remains valid [e.g. Eq. (26)]. However, for g/K 1, we expect the fidelity to be sensitive to the discrete nature of the size distributions, and our predictions based on the first few moments may no longer apply. Fortunately, as we show in the following subsections, none of these complications are evident for experimentally relevant system sizes (e.g. K ∼ N ∼ 20) and g ∼ O(1) coupling strengths; indeed, our finite-size numerical results agree very well with predictions from the peaked-size teleportation framework [ Fig. 7(b) and 8].
Lastly, in the case where g/K ∼ 1, operator growth is no longer necessary for the coupling to have a strong effect, leading to the possibility of a teleportation signal unrelated to scrambling. Indeed, for g/K = π, the coupling effectively 'swaps' the left and right qubits. This is made precise for the coupling V s [Eq. (27)], where exp(iπN V s ) = (SWAP)Y l Y r . In this case, one would observe perfect teleportation fidelity even without manybody time evolution, i.e. U = 1; in fact, if U is perturbed away from the identity via scrambling dynamics, the teleportation fidelity would actually become suppressed. The simplest way to see this is via Fig. 1(a)-in particular, any subsequent time-evolution on the right side of the system is in the wrong direction to refocus the timeevolved state (one would want to apply U † after the coupling, not U T ). To achieve a large teleportation fidelity, the combined time-evolution, U T U , would therefore need to preserve the "teleported" state, ψ| U T U |ψ ∼ 1, a situation that is only likely to occur if the dynamics are non-scrambling (U = 1 is a special case of this) or undergo a late-time, fine-tuned, Poincare-type recurrence.

B. One-sided implementation of teleportation circuit
Before proceeding to the experimental blueprints, we first introduce a simpler implementation of the teleportation protocol that works at infinite temperature (Fig. 6). The outcome of this protocol is equivalent to that of the two-sided protocol (up to experimental errors), yet  9) [replacing U → U T for convenience, compared to Fig. 1(a)]. Blue arrows denote the sequence of operations in the one-sided protocol, the green band marks the teleported qubit and its corresponding component in the one-sided protocol, and the red band marks the initial EPR state and its corresponding component.
it eliminates the need to prepare EPR pairs and requires half as many degrees of freedom. The cost of this simplification is two-fold: (i) it is restricted to simulating an infinite temperature TFD state, and (ii) it requires a higher depth quantum circuit.
We derive the one-sided implementation from the 'twosided' implementation [copied in Fig. 6 from Fig. 1(a)] by sliding all operations from the left side of the manybody EPR pairs to the right side, using Eq. (9). The initial state of the one-sided circuit thus corresponds to the top left of the two-sided implementation. Namely, we initialize the K 'measured' qubits of subsystem C in a definite outcome state, |o 1 · · · o K (purple). These states should be drawn from the distribution of measurement outcomes, but when teleporting an EPR pair at infinite temperature they will be uniformly distributed. For the N − K 'unmeasured' qubits, we use the resolution of the identity 1 ∝ s |s s| to replace the unterminated legs with an initial product state in the computational basis, |o K+1 · · · o N (gray). This state should be sampled from shot-to-shot over all 2 N −K basis states, in effect preparing a maximally mixed state on these qubits. Finally, we include one ancillary qubit for each qubit to be teleported, whose initial state is sampled over a complete basis |φ for the teleported subsystem (i.e. subsystem A in Section II). Similar to the unmeasured qubits, this corresponds to the unterminated leg of the thermofield double state when we insert the teleported qubit |ψ in the two-sided implementation.
Having defined an initial pure state, we now implement the circuit starting from the top left of the twosided implementation and proceeding counter-clockwise (Fig. 6). The circuit consists of three successive applications of U or U † , interspersed with a swap gate exchanging subsystem A with the ancillary qubit(s), and operationsV i = e igoiÔi/K determined by the initial state of the 'measured' qubits. The outcome of the circuit is an EPR measurement between the ancilla qubit and subsystem A (black arrows).
As one can see in Fig. 6, the one-sided implementation no longer performs teleportation, but rather prepares an EPR pair from an otherwise scrambled, many-body system. Specifically, we know that upon swapping out, subsystem A is maximally entangled with the remaining qubits whenever the unitary, U , is scrambling; the one-sided circuit distills this entanglement into an output EPR pair. This connection has been noted in gravity, where similar one-sided protocols can be interpreted as distilling the partner operators of emitted Hawking radiation [71,72] or observing behind the horizon in the SYK model [73].

C. Preparing the thermofield double state
In the previous subsection, we introduced a one-sided protocol that obviates the need to prepare the highly entangled TFD state. However, this approach was restricted to infinite temperature; at finite temperature, one must implement the original two-sided protocol, which necessitates preparing a finite temperature TFD state. A number of recent works have explored the preparation of TFD states variationally using quantum approximate optimization algorithms (QAOA) [74][75][76]; we note that these preparation strategies require no additional experimental capabilities beyond those already necessary for the TW protocol. The optimization step within a QAOAbased TFD preparation relies on a cost function that requires one to measure the entanglement entropy between the two sides [74,75]. While challenging, this can in principle be experimentally realized by either using several copies of the system [77][78][79] or via randomized measurements [80], both of which have been demonstrated in small-scale trapped ion experiments [81,82].

D. Implementation with neutral Rydberg atoms
One particularly promising platform for implementing the traversable wormhole protocol is a quantum simulator based on neutral alkali or alkaline-earth atoms held in a reconfigurable and controllable array of optical dipole traps. Recent experiments have already achieved near-deterministic trapping and loading of atoms into arbitrary geometries in one, two, and three dimensions [34,84,85]. By leveraging the strong dipole coupling between atomic Rydberg states, high-fidelity analog quantum simulations and digital gates have also recently been demonstrated [31,[34][35][36][37][38]. These demonstra-tions have primarily used two natural schemes of encoding qubits into neutral atoms:

A qubit can be encoded by choosing an atomic
ground state |g to be the |0 state, and a highly excited Rydberg state |r with principal quantum number n 1 as the |1 state [see Fig. 7(a)].
2. Alternatively, the qubit states can also be chosen as two long-lived hyperfine ground states (for alkali atoms or fermionic alkaline earth atoms) or a ground state and a metastable clock state (for bosonic alkaline earth atoms), such that the |1 state can be coupled to a Rydberg state to perform entangling gates [see Fig. 7(c)].
We will show how both encodings can be used to realize the teleportation protocol in feasible near-term experiments. We find that the first encoding is naturally suited to 'analog' time-evolution under the native (Isingtype) Hamiltonian for a Rydberg setup, but is limited to system sizes of 30−35 qubits (in one spatial dimension) due to the inability to perfectly time-reverse long-range interactions. On the other hand, the second encoding is more flexible and allows for digital time-evolution including RUCs and Floquet dynamics. This time-evolution can be reversed exactly and is limited only by qubit and gate fidelities. While we will primarily consider realizations of our protocol in experimental setups where the neutral atoms are individually trapped in optical tweezers and undergo (near-)resonant excitation to Rydberg states, we also conclude by discussing how similar physics can be seen in an optical lattice setup where the atoms are primarily in ground states |0 and |1 , but one of these states is 'dressed' by an off-resonant laser field which couples it to a Rydberg state [86][87][88].
Analog implementation-We first consider the encoding where the qubit states |0 and |1 correspond to a ground state |g and a highly excited Rydberg state |r . While neutral atoms are effectively non-interacting in their ground states, nearby atoms interact strongly via van der Waals interactions ∝ n 11 /R 6 if they are both in the Rydberg state, where R is the distance between the atoms. If one drives the transition |g i ↔ |r i at each site i with tunable Rabi frequency Ω i and detuning ∆ i [see Fig. 7(b)], the system will undergo analog time evolution under the Hamiltonian Here, qubits are encoded in two hyperfine ground states. Insets show possible pulse sequences to implement the controlled-phase gate and the single-qubit rotations [83]. The full TW protocol is obtained by inserting this gate sequence (and its Hermitian conjugate) in place of U , U † in Fig. 5.
we reverse the nearest-neighbor interactions by conjugating time-evolution via Pauli operators X i (i.e. applying π-pulses) on every other site. The tunable single-site parameters Ω i and ∆ i are then adjusted to ensure that each single-site term is also reversed. We note that this simple scheme does not reverse the (much weaker) next-nearestneighbor interactions. In a one-dimensional array, the errors in our implementation will arise from two main sources: (i) the finite lifetime of the Rydberg state, which gives rise to a nonzero decoherence rate at each of the N sites, and (ii) the weak next-nearest neighbor interactions ∼ J 0 /2 6 = J 0 /64, which cannot be time-reversed simultaneously with nearest neighbor interactions. To estimate the effect of the former, let us consider the specific case of 87 Rb atoms excited to the 70S Rydberg state [31,35], which has a lifetime τ ≈ 150 µs. Realistically achievable Rabi frequencies and interaction strengths are of order ∼ 2π × 10 − 100 MHz. The total time to implement the three scrambling unitaries of the teleportation protocol is thus ∼ 3N/|Ω i |; when summed over N qubits and compared to the Rydberg lifetime, this gives an estimated many-body error ∼ 3N 2 /|Ω i |τ .
In order to precisely characterize the effects of imperfect backwards time-evolution, we perform large-scale numerical simulations of the teleportation protocol with the Rydberg Hamiltonian, Eq. (104) [89]. Our results are depicted in Fig. 7(b) for a one-dimensional chain of N = 20 atoms and three values of the coupling g. Analogous to our 1D RUC numerics [ Fig. 2(a)], the fidelity increases monotonically in time for g = π; while, for g = 2π and g = 3π, the fidelity oscillates in time, reaching a local maximum whenever the average size satisfies the phase-matching condition [Eq. (35)]. Notably, even with perfect time reversal, the overall fidelity is reduced from unity due to the finite width of the size distribution. This is a general feature of peaked-size teleportation in finite-size systems, since the relative size width scales as δS/S ∼ 1/ √ N (Section VI). Indeed, in Fig. 8, we confirm that the fidelity improves with increasing system size and is consistent with our peaked-size error analysis [e.g. see Eq. (34)].
With imperfect time reversal, we observe an additional ∼ 10% reduction in the fidelity compared to the ideal case at the scrambling time [ Fig. 7(b)]. We can estimate the magnitude of this effect by assuming errors due to the next-nearest-neighbor interactions add coherently over time-intervals δt ∼ 1/J 0 (the local thermalization time), and incoherently at larger time-scales. Within each δt, each atom accumulates an error ∼ (δt J 0 /64) 2 ; summed over N atoms and total time 3t * ≈ 3N δt, this gives a total many-body error ∼ 3N 2 /64 2 . Thus, the error due to imperfect time reversal is magnified at larger system sizes and will eventually outweigh the improvement in fidelity from the narrowing of the size distribution.
Combined with the Rydberg lifetime error, this suggests that near-term experiments should be able to implement peaked-size teleportation in systems of N ∼ 35 qubits. We note that in higher dimensions, the smaller relative distance of next-nearest neighbor atoms gives rise to a larger error contribution from imperfect timereversal.  104), with the same system parameters as in Fig. 7. At late times, the fidelity increases for larger systems but decreases for larger values of g. This is consistent with our error analysis in Section VI; in particular, we expect the error to scale as g 2 δS 2 /N 2 and the size distribution to approach a binomial distribution for which δS ∼ S/ √ N . In contrast, at early times, smaller systems exhibit a larger fidelity not because of the size width but because the acquired phase is η d gS(t)/N , where η d g is fixed and S(t) is initially independent of size. The curves in (a) intersect near the scrambling time due to the transition between the early and late time regimes.
Digital implementation-To implement the protocol in larger systems, higher dimensions, and at finite temperature, we propose a digital scheme, using the second type of qubit encoding (i.e. hyperfine ground states) [ Fig. 7(c)]. In this approach, we envision time-evolution to be formed from alternating layers of nearest-neighbor controlled-phase gates and single-qubit rotations. Here, the controlled-phase gates can be implemented by applying a simple pulse sequence to excite and de-excite qubits from the |1 state to the |r state, so that the wavefunction acquires a phase of −1 if either of the two qubits are in the |1 state, but not if both qubits are in the |0 state [see Fig. 7(c) insets] [83]. As demonstrated in recent experiments [90], these Rydberg-mediated controlled-phase gates can be performed in parallel for sufficiently wellseparated pairs of qubits, and non-nearest neighbor interactions can be avoided by slightly reducing the parallelism within each layer of controlled-phase gates. Singlequbit rotations can be performed with sufficiently high fidelity such that the overall circuit fidelity is primarily limited by the entangling gates [84,91].
For a generic choice of gates, the circuit will be fully scrambling when U is composed of ∼ N layers of controlled-phase gates. The fidelity of the overall implementation is limited by the finite lifetime of the Rydberg state, which is populated for time ∼ 1/J 0 during each controlled-phase gate. Assuming the same experimental parameters as in the analog case, one expects to be able to perform approximately Ωτ ∼ 10 3 − 10 4 controlledphase gates within the decoherence time-scale. Thus, in the digital approach, one expects that the teleportation protocol can naturally be implemented for N ∼ 200 qubits.
The digital approach can also be adapted to experiments using Rydberg-dressed neutral atoms in an optical lattice [86][87][88]. In such a setup, qubits are again encoded in hyperfine ground states and strong Ising-like interactions are generated by coupling the qubit state |1 to a Rydberg state with a far-detuned laser field. In this way, the Rydberg interaction gives rise to an energy shift for two neighboring atoms both in the |1 state. Analogous to our previous discussion, a simple scrambling unitary could consist of alternating layers of Rydbergdressed interactions and single-qubit rotations. While the total accumulated error in the Rydberg-dressing approach is comparable to the gate-based protocol, one potential advantage is an increased tunability of the interactions [92,93].
In addition to scrambling time evolution, there are three ingredients to implement the one-sided teleportation circuit (Fig. 6): (i) the ability to 'swap' in the qubit |φ , (ii) single-qubit rotations, V i = e ±igZi/K , and (iii) the final measurement in the EPR basis. In both digital setups, these are easily accomplished by combining controlled-phase gates, arbitrary single-qubit rotations, and local measurements. In the analog setup, we propose to temporarily 'turn off' the Hamiltonian by transferring each Rydberg state |r to a hyperfine ground state (e.g. the state used as |1 in the digital protocol) using a resonant laser pulse. Once this is done, all of the above operations can be performed identically as in the digital setup. Afterwards, an additional resonant laser pulse returns the system to the analog encoding. The ancillary qubit can be decoupled from the system qubits during Hamiltonian time-evolution in two ways: (i) by physically positioning the ancillary qubit far from the system, or (ii) by encoding the ancillary qubit in the hyperfine subspace throughout time-evolution.
The two-sided, finite temperature TW protocol can be achieved by combining the above techniques with TFD preparation as in Section IX C. A particularly natural geometry for such a realization would be two parallel chains of Rydberg atoms, with each chain forming one side of the TFD state. The coupling between the two sides is naturally realized by the atoms' Ising interactions. This coupling can be applied independently from the one-sided Hamiltonian using either full digital control or by manipulating the inter-vs. intra-chain atomic distance.  Figure 9. (a-b) Chain of atomic ions, with qubit states |0 , |1 represented by hyperfine ground states. The states are coupled by a pair of laser beams, one with individual addressing (with strength g1, purple) and one applied globally (with strength g2). Each beam is strongly detuned from an excited state |e by an amount ∆. The coherent beatnote between the beams, at frequency ω0, drives stimulated Raman transitions between the qubit levels with an effective Rabi frequency g1g2/2∆, and also modulates the Coulomb interaction between qubits to give rise to an effective Ising interaction. (a) A two-qubit entangling gate, XXij(θ), (red) is performed by addressing only ions i and j with the first beam. (b) Half of the qubits are addressed, which leads to analog time-evolution under the Hamiltonian Eq. (105) (blue) for all addressed spins. (c) Quantum circuit implementation of the teleportation protocol at finite temperature. EPR pairs are formed using two-qubit gates. The TFD state is then prepared via a QAOA approach by iterating multiple times between two-qubit gates coupling the sides and analog time-evolution on both sides individually [74,75]. The state |ψ is inserted either by projectively measuring the designated qubit and preparing the state, or by digitally swapping in an additional qubit (not shown). Finally, teleportation is implemented using similar ingredients as well as feed-forward measurements (purple dotted lines).

E. Implementation with trapped ions
A second experimental platform that naturally enables the implementation of the TW protocol is arrays of individual trapped atomic ions [94][95][96]. Trapped ion qubits feature near-perfect replicability, negligible idle errors, and the ability to implement both a universal set of reconfigurable quantum gates [43] as well as analog long-range spin Hamiltonians [39,40]. Entangling quantum gates have been demonstrated between isolated pairs of trapped ions with fidelities exceeding 99.9% [41,42]. Teleportation protocols-including the HPR protocol [13]-involving gate operations, partial measurement and feedforward operations, have been experimentally realized in a number of contexts [4,5,13,97].
Compared to Rydberg atom arrays, trapped ions offer two new regimes for exploring many-body teleportation. First, trapped ions naturally interact via a long-range analog Hamiltonian, whose time-evolution can be fully reversed within certain experimental regimes [98,99]. Implementing the TW protocol in this setting would provide a window into operator spreading and size distributions under such long-range dynamics [100,101]. Second, when operated digitally, the same long-range interaction has already been demonstrated to enable the preparation of thermofield double states [74][75][76]102], a crucial step towards realizing the two-sided TW protocol at finite temperature (see Section IX C).
We begin by outlining the analog and digital forms of time-evolution that are possible in trapped ion systems. Interactions between qubits typically stem from state-dependent optical dipole forces that off-resonantly drive motional sidebands of the qubit [103,104]. These sideband operations mediate entanglement and give rise to an effective Ising coupling. When the optical forces are symmetrically detuned far from the upper and lower sidebands, the motion is only virtually excited, resulting in a long-range Ising Hamiltonian [ Fig. 9(b)]: where J ij ≈ J 0 /|i − j| α , with 0 < α < 3 and J 0 1 kHz, and the effective magnetic field B z can be realized by slightly asymmetrically detuning the driving field [105]. The sign of the couplings can be reversed by changing the detuning of the optical forces from the motional sidebands [98,99]. On the other hand, when the optical dipole forces are closer to resonances of the motional modes, one can mediate interactions significantly faster, allowing for the execution of rapid, entangling quantum gates between pairs of illuminated ion qubits [ Fig. 9(a)] [106,107]. The native entangling gates are based upon Ising interactions between any selected pair of ions with a tunable interaction angle; in particular, both XX ij (θ) = e −iθXiXj /2 and Y Y ij (θ) = e −iθYiYj /2 gates are available and θ = π/2 naturally creates an EPR pair [108,109]. Typical entangling operations have duration 1/J ent ∼ 100 µs, while decoherence time-scales are on the order of τ ∼ 400 ms [110]. Following the estimates of Section IX D and requiring 3N 2 /J ent τ 1, we estimate that near-term state-ofthe-art experiments can support high-fidelity many-body teleportation for up to N ∼ 35 qubits.
Let us now describe an implementation of the onesided TW protocol (Fig. 6). We first focus on the ability to implement both U and its inverse U † . For analog time-evolution [Eq. (105)], U † can be implemented by changing the sign of the detuning [65], while for digital time-evolution, one can directly invert and reverse the ordering of the quantum gates.
The one-sided protocol also requires the ability to locally address a sub-extensive number of individual qubits. In particular, a subset K of the qubits must be initially prepared in a product state, |o 1 , . . . , o K and later rotated byV i = e igoiÔi/K . These rotations can be achieved by takingÔ i =Ẑ i and individually addressing the target ions using an auxiliary "poke" laser beam [96,111].
Following the first application of U , one must swap out the qubit(s) corresponding to the teleported subsystem. This swap can be implemented either digitally by applying a SWAP-gate, or physically, by exchanging the two ions via a modulation of the ion trap's axial fields [40,112,113].
Extending this implementation to the two-sided protocol [ Fig. 1(a)] is straightforward. Initialization into EPR pairs (for infinite temperature) can be accomplished via simple Ising gates at the input of the circuit [ Fig. 9(a,c)], while the TFD state (for finite temperature) can be prepared via variational methods (Section IX C). Timeevolution can again take the form of either digital quantum gates [ Fig. 9(a)] or analog Hamiltonian dynamics. To separately implement analog dynamics on the two sides of the system, one would illuminate only half of the ion chain at any given time [ Fig. 9(b)]; this has the added benefit of avoiding unwanted coupling between the left and right sides, but implies that the time-evolution must be performed serially [ Fig. 9(c)].
Finally, in the two-sided protocol, one must perform projective measurements on K qubits that feed-forward to the conditional rotations,V i . These partial measurements can be accomplished by using multiple ion species (i.e. different elements or isotopes) [97], or alternatively, this entire procedure can be replaced with a specific interaction, e igV , between the two sides; this interaction is naturally realized via an XX ij (θ) gate with θ = 2g/K.

F. Effects of experimental error and relation to quantum error correction
We now turn to the effect of experimental error on the TW protocol. We find that teleportation is robust to nearly all errors that occur on the left side of the TFD state after time-evolution by U , but is strongly sensitive to errors at nearly all other locations in the protocol. These two extremes are emblematic of two different relations between scrambling and error: the former corresponds to interpretations of scrambling as an error-correcting code [23], while the latter reflects recent results showing that the effect of errors on scrambling measurements is enhanced proportional to the size, S, of time-evolved operators [114]. In the following discussion, we demonstrate each of these points through simple but representative examples of experimental error.
We begin with the first case: consider errors occurring on the left side of the TFD state after application of U but before measurement/coupling. Recall that, in the absence of error, one can perform teleportation by using any K ∼ O(1) qubits of the left side. This implies that teleportation is robust to any errors that affect only N −K qubits: as long as one has knowledge of at least K qubits that are unaffected, measuring these qubits performs teleportation identically to the error-free case.
This robustness reflects previously noted connections between scrambling and quantum error correction [23]. In particular, we note that many-body teleportation can be understood as an especially generic example of entanglement-assisted quantum error correction (EAQEC) [115]. Indeed, the setup for EAQEC is identical to that of the teleportation protocol: two parties, Alice and Bob, share entanglement (the TFD state), Alice applies an encoding circuit to her share of qubits (the left unitary, U ), and decoding is achieved by teleporting Alice's quantum state to Bob's share of qubits (via the coupling, V , and unitaries on the right). Previous schemes for EAQEC have focused primarily on encodings via Clifford unitaries. In contrast, many-body teleportation, and more specifically peaked-size teleportation, succeeds for a vastly broader class of encoding procedures-i.e. scrambling many-body time dynamics-indicating that naturally occurring, strongly interacting systems offer novel methods of EAQEC.
On the other hand, errors that occur during encoding or decoding-i.e. during the application of U on the left side or at any point on the right side-strongly inhibit teleportation. As a first example, consider a single local error, W 1 , occurring with probability ε on the right side after coupling but before U T (i.e. just before decoding). If the error, W 1 , grows to have a size, S, after U T is applied, one estimates that it will decrease the teleportation fidelity by an amount, 1 − F ∼ εS/N , proportional to the probability that W 1 has support on the teleported qubit after time-evolution. If we sum over such errors on all N qubits, we have 1 − F ∼ εS.
As a second example, consider a local error, W 2 , occurring with probability ε on the left side simultaneously with state insertion (e.g. a damaged TFD state in Fig. 1). In effect, this error shifts the correlator operators [Eq. (2)], Q → Q ⊗ W 2 ; following the arguments of Section VII, one then requires that the sizes add for teleportation to succeed, S[QW 2 ] = S[Q] + S[W 2 ]. In a 1D short-range system (Section VII B), this condition holds if and only if the light cones of W 2 and Q do not overlap. For O(εN ) randomly distributed errors, we expect this to hold as long as the spacing between errors, 1/ε, is much larger than the size of the light cone, εS 1. A similar scaling holds in 0D (Section VII C). Here, we expect size addition to hold as long as the size of the total error (corresponding to a time-evolved product of ∼ εN initially local operators), is much smaller than the system size, N . Once again, this requires εN S N , or εS 1.
The two previous examples are straightforwardly generalized to errors that accumulate continuously throughout time-evolution. To do so, we replace the error probability with an error rate, ε (now with units of inverse time). The total effect of the error is then given by the integral of the error rate multiplied by the size over time, ε´t 0 dt S(t ) [114]. In one-dimensional systems evolved up to the scrambling time, i.e. S ∼ Jt and t s ∼ N/J for a local interaction strength J, we thus estimate a total error, ε´t s 0 dt Jt ∼ εSt s ∼ εN 2 /J, in agreement with our rough estimates in Sections IX D and IX E.
Finally, we consider a particular form of error that may be relevant for analog time-evolution: mismatches between the evolution times of U , U * , and U T . We denote these three evolution times as t 1 , t 2 , t 3 , respectively, and their mismatches as ∆t 12 = t 2 − t 1 and ∆t 13 = t 3 − t 1 . We can characterize the mismatches' effect on the teleportation fidelity using the correlators, III). From this, we anticipate that the protocol is relatively insensitive to mismatches between t 3 and t 1 , t 2 : teleportation succeeds as long as the mismatch is small compared to the local interaction strength, J, i.e. J∆t 13 1. To estimate this, we set g = 0 and t 1 = t 2 in which case the correlator magnitude is given by an autocorrelation function, C Q = Q(t 1 )Q(t 3 ) = G(∆t 13 ). The teleportation fidelity is bounded above by this expression, which decays on a time-scale ∼ 1/J. On the other hand, teleportation is more strongly sensitive to the mismatch between t 1 and t 2 . To estimate this, we treat the difference in time-evolution between U and U * as a product of ∼ (J∆t 12 ) 2 N local errors occurring simultaneous with state insertion (to motivate this scaling, note that one can approximate U (∆t 12 ) as a product of ∼ N local unitaries for small ∆t 12 , and we expect the error to be an even function of ∆t 12 ). Following our previous analysis, teleportation is successful as long as S(J∆t 12 ) 2 N N , or S(J∆t 12 ) 2 1.

G. Directly measuring the size distribution
In Section IX A, we discussed that the time profile of the teleportation fidelity reveals important features of the operators' size distributions, including the average operator size and the size width. We now demonstrate that a more precise characterization of the operator size distribution can be obtained by sweeping the coupling strength, g, at a fixed time, t.
For simplicity, we restrict to infinite temperature 18 and the coupling V s in Eq. (27), which precisely measures the operator size. In this case, the two-sided correlator [Eq. (2)] is equal to the characteristic function, Φ S (g), of the size: from which the size distribution can be obtained by a Fourier transform in g.
More precisely, to measure the real part of the characteristic function (i.e. the teleportation correlator), we perform the teleportation protocol with two small modifications: (i) we replace state insertion with the specific projection operator, (1 + Q)/2, and (ii) we measure the expectation value of Q applied to the right side, instead of the teleportation fidelity. This yields the quantity: where in the second line we use that the "diagonal" terms between the two copies of (1 + Q)/2 vanish at infinite temperature. The imaginary part of the characteristic function can be obtained similarly, by replacing state insertion, (1 + Q)/2, with an application of the unitary operator, (1 + iQ)/ √ 2. Analogous to Fig. 6, both of these measurement schemes can be adapted into one-sided protocols using Eq. (9) whenever the coupling V is classical (i.e. composed of terms O i,l O * i,r , where {O i } mutually commute). While such couplings do not measure the exact size distribution, we expect their behavior to be similar in most cases (Section IV A).
For completeness, we also note an alternate method to measure the size distribution: one prepares the state Q l (t) |EPR and directly measures the two-sided coupling V s . The probability distribution of the measurement results gives the size distribution [see the discussion below Eq. (27)].
Let us now compare these two protocols to other schemes for characterizing the size distribution of operators. First, we recall that a sum of local OTOCs yields the average operator size [Eq. (31)]. Hence, many existing protocols for measuring local OTOCs [116,117] can be straightforwardly adapted to measuring the average size. Higher order moments of the size distribution can similarly be obtained from local OTOCs, using Eq. (27): where the sum is over every possible combination of n single-qubit Pauli operators P i1 , . . . , P in . Based on this approach, however, the number of measurements required to determine the n th moment scales as O(N n ). In certain situations, this scaling may be reduced through sampling, though this depends on the nature of the size distribution and the desired degree of precision. Furthermore, reconstructing the full profile of the size distribution from a finite number of moments is generally a difficult numerical task [118]. In contrast to these limitations, our proposal directly yields the full size distribution, and can recover its moments with a number of measurements independent of the system size 19 . We can also compare our proposal to an independent protocol for measuring the size distribution introduced in Ref. [28]. The protocol of Ref. [28] is experimentally simpler than our own, and in particular involves only a single application of time-evolution by U (and no backwards time-evolution). However, this simplicity comes at a cost: resolving high-size components of the distribution requires a number of measurements that scales exponentially with size.

X. OUTLOOK
In this work, we developed a unified framework for understanding many-body teleportation from the perspective of operator growth under scrambling dynamics. The unifying concept within this framework is the size distribution of time-evolved operators [15,16,[26][27][28]: these form the backbone of peaked-size teleportation, and provide a more fine-grained measure of operator growth compared to the average operator size (as given by the expectation value of OTOCs).
Our work suggests several future directions for applying and building upon this framework. First, while we have studied the size distributions in 0D and ≥ 1D RUCs, it would be interesting to extend this analysis to a multitude of other physical systems, where one expects to find qualitatively distinct behavior. These include long-range interacting systems [119,120], interacting and non-interacting integrable systems [28], ≥ 1D systems with a large on-site Hilbert space [121], 0D systems with sparse couplings [122], and systems with conserved quantities [50].
Another set of open questions concerns the notion of operator size at finite temperature. In systems with peaked size distributions, we found that the phase of the 19 This is simplest to see in the protocol which measures the two-sided coupling Vs. Here, the error in one's measurement of the n th moment is equal to the expectation value of the moment's variance divided by the number of measurements, (δ S n ) 2 = ( S 2n − S n 2 )/M . If one wishes to resolve the moment to within a relative error ε, i.e. δ S n < ε S n , one requires M ∼ S 2n − S n 2 ε 2 S n 2 measurements. This number does not scale with N since it contains the same powers of S in the numerator and denominator.
two-sided teleportation correlator was directly proportional to the conventional definition of operator size [27]. Surprisingly, we observed that this relationship did not hold in the finite temperature SYK model; rather, the phase was given by the real part of the two-sided OTOC. Unlike the conventional size, this OTOC is not UV divergent, and is thus expected to be inherently independent of the microscopic Hilbert space. Recent work has shown that its real part isolates an incoherent component of operator spreading in large-N models [33]; further work is needed to establish and expand this framework. Related to these considerations, one may hope to better understand the bulk analogue of operator size in theories dual to gravity with strong stringy effects. While we have seen that stringy effects can mimic peaked-size teleportation, developing a physical interpretation of this correspondence would be extremely exciting. Third, we have shown that a promising application of the teleportation protocol is to distinguish between different classes of scrambling dynamics. In particular, we have focused on two classes of scramblers-generic thermalizing systems and those with gravitational dualsand demonstrated that the key distinction between them is their teleportation fidelity at low temperatures. It is intriguing to ask whether the fidelity increase associated with gravitational teleportation may also occur in other systems, without a gravitational dual. For instance, recently the teleportation correlator magnitude was observed to increase slightly above G β in non-local random Hamiltonian systems [15,16]; generalizing this to other physical models would be of tremendous interest.
One may also wonder what role an extensive low temperature entropy-a key feature of the SYK model [25]plays in the teleportation process. In particular, how well can systems with extensive low temperature entropy but no known gravitational dual teleport [123,124]? We conjecture that an extensive entropy would allow one to locally encode each qubit into low-energy degrees of freedom (i.e. operators with an O(1) two-point function), since one would only require O(1) qubits on the left side of the TFD in order to have one qubit of mutual information with the right side. Such an encoding would allow low temperature teleportation with perfect fidelity if operator sizes were peaked, naturally motivating the study of operator size distributions in such models.
Finally, we would like to discuss the relation between our results on the TW protocol and the eternal traversable wormhole (ETW) introduced in Ref. [10]. In the latter, the coupling, V , has an O(1) coefficient and, moreover, is applied simultaneously with single-sided Hamiltonian evolution (i.e. the full system evolves under a Hamiltonian, H l + H r + g j O j,l O * j,r ). Under these conditions, Refs. [10,125] find that the ETW teleportation fidelity oscillates in time under gravitational dynamics, indicating that information is transmitted back and forth between the two boundaries. Intriguingly, unlike the TW protocol, the ETW oscillations occur at a time-scale given by the single-sided thermalization time (∼ β, the inverse effective temperature), and not the scrambling time. Developing a microscopic understanding of the ETW in terms of operator spreading, as well as exploring analogous physics in more generic many-body systems, remains an exciting open direction.
Note added -After this work had been completed, we learned of an independent investigation of gravitational many-body teleportation by Nezami, Lin, Brown, Gharibyan, Leichenauer, Salton, Susskind, Swingle, and Walker, which will appear in the same arXiv posting.
In this Appendix, we provide a precise mathematical bound guaranteeing that the teleportation correlator obeys the peaked-size prediction [Eq. (32), Section IV B] when the size distribution is sufficiently tightly peaked. We apply this bound to two examples where the size distribution is known exactly: late times in all scrambling systems (Section VI), and the large-q SYK model (Sections VII D and VIII C). Notably, in the latter we find that our bound applies only at infinite temperature, despite the profile of the size distribution (e.g. its ratio of size width to average size) behaving similarly at all temperatures. The discrepancy arises instead because the correlator magnitude, (G β ) p , decreases exponentially in the encoding size p at all finite temperatures.

Precise bound
As in the main text, we decompose a time-evolved finite temperature operator into a sum of Pauli strings: In this basis, for qubit systems the correlator takes the form Here we define the winding size distribution [15,16] At finite temperature, this size wavefunction is distinct from the size distribution: which is a real, normalized probability distribution probed by the one-sided correlator [27] TFD| Q Nevertheless, the size distribution bounds the size wavefunction magnitude via the triangle inequality: with equality achieved when all Pauli operators of size n contribute the same phase to f (n).
The average size and size variance are easily found from the size distribution as where we work in the continuum limit replacing sums over the size by integrals for simplicity. We now define the asymptotic size width with error ε as the minimal width W ε about the average size such that i.e. a fraction 1−ε of the size distribution's support is contained in the interval I = [S −W ε , S +W ε ] (the lower limit of the integral should be bounded by zero; for simpler notation we'll deal with this by instead defining P (n) = f (n) = 0 for n < 0). We can now separate the correlator into two pieces, one arising from sizes in the interval I and the other from the interval's complementĪ where the remainder R =´Ī dn f (n)e iη d gn/N is strictly smaller than ε: Peaked size teleportation occurs in the regime where gW ε /N 1. In this limit, we can expand where the deviation for n ∈ I is bounded by which holds as long as gW ε /N ≤ π/2. We then have is the imaginary time two-point function, and the error R = e igS/N´I dn f (n)E(n) is bounded by and the second error R = G β (Q A ) −´I dn f (n) is bounded by We therefore conclude that whenever η d gW ε /N ≤ π/2, the deviation of C Q from the peaked size value is controlled by the upper bound Practically speaking, the lowest value of g for successful peaked-size teleportation is η d gS/N = π. Therefore, for a given size distribution, we can guarantee that peaked-size teleportation is possible if we find ε such that B G β (Q A ), i.e. the error in the correlator is small compared to the correlator magnitude.

Application to late times
We illustrate this with some examples, in the few cases where we can exactly solve for operators' full size distribution. First, consider a thermalized system at late times, which we will approximate by setting the size distribution of Q A (t) to be that of a random Pauli string. For large n, N is a Gaussian distribution with mean S = 3N/4 and variance δS 2 = 3N/16: We therefore have The error function decays exponentially in its argument, so even for exponentially small ε we require only W ε = AδS for some constant A ∼ O(1). Setting g equal to its minimal value, η d gS/N = π, we have both ε 1 and sin(η d gW ε /N ) ≈ AδS/S ∼ 1/ √ N 1, and so peaked size teleportation is guaranteed.

Application to the large-q SYK model
We can also use this method to guarantee peaked-size teleportation in the large-q SYK model at infinite temperature. We begin by writing down the size distribution for the large-q SYK model in detail, quoting the results of Ref. [27]. The generating function for the size distribution is: where we define From this, we can identify the size distribution: The size and size width are Therefore, the ratio of size width to average size is which approaches zero when p → ∞ (∆ → ∞).
To apply the upper bound Eq. (A16), we need to integrate (i.e. sum) the tail of the size distribution in order to compute its asymptotic width [Eq. (A8)]. In this example, the discrete tail can be summed explicitly and we define where B x (a, b) and B(a, b) are incomplete and ordinary beta function respectively. Let us take k =n(1 ± ζ) for some small ζ representing the asymptotic width This width corresponds to an error ε = 1 − I(n(1 − ζ)) + I(n(1 + ζ)).
Taking gS/N = π, the upper bound is At infinite temperature G β (Q A ) = 1, we need to show that the minimum of B tends to zero when ∆ → ∞.
For early time sinh Jt ∼ O(1), 1 − x is an order 1 number, and we take ∆ → ∞ limit to get The bound becomes This basically means that the integrated probability betweenn(1 − ζ) andn(1 + ζ) for any finite ζ is 1. One can thus take ζ → 0 with speed slower than 1/∆ → 0 in order to have the bound vanish. This computation applies for x ∈ (0, 1), which means that the peaked size always holds for early time. This is physically reasonable as the operator has not yet been scrambled extensively. However, since the size is small at such early times, in order for teleportation to work we must choose g ∼ N .
For intermediate times, such that sinh 2 Jt ∼ N and ∆ N ∼ 1/(1 − x), we must take the x → 1 limit first. Using the fact that where F is Gauss hypergeometric function, in x → 1 limit the right portion of Eq. (A30) tends to where Γ(x, a) is incomplete gamma function. Meanwhile, the left portion of the second term of Eq. (A30) gives Combining the two, we have It follows that the upper bound is This function has a unique minimum for ζ ∈ [0, 1/2] and this minimum decreases as ∆ increases. Taking derivative with respect to ζ, we get where in the second step we have taken large ∆ limit. Solving ∂ ζ B = 0 in this limit, we find the minimum at which in turn gives the limit value of B to be zero. This proves that at infinite temperature, teleportation exactly matches the peaked-size prediction for both early and intermediate times. For late times t 1 2J log N the size distribution above breaks down, as can be seen since P (n) is dominated by some n > N , which is unphysical since N is the total number of fermions.
In contrast, we can also show that the above bound does not apply at low temperatures for large-q SYK, as expected from the main text. At low temperature, the upper bound B needs to be much smaller than the two-sided correlation function G β (Q A ) ∼ (βJ) −2∆ in order to guarantee peaked-size teleportation. The low temperature size distribution is essentially the same as at infinite temperature, requiring only the replacement [27]: and adding e −µN δ β to the distribution, which shifts the initial size by a constant amount N δ β (accounting for the size of the thermal density matrix [27]). Following a similar computation to above, one can show that B still asymptotes to zero, but now with a slower speed than G β (Q A ). For example, in the early time and large ∆ limits, ) is exponentially smaller for large βJ. Therefore, the upper bound B fails to guarantee peaked-size teleportation. This is consistent with the fact that the correlation function C Q (t) in Eq. (86) in low temperature is far from being a pure phase.

Appendix B: The Hayden-Preskill recovery protocol
In this appendix we review the HPR protocol following Refs. [11,12] and derive its equivalence to the TW protocol in the case of infinite temperature teleportation of a single qubit (introduced in Section VI B). This single-qubit variant of the HPR protocol was experimentally implemented in Ref. [13], although an explicit derivation of its quantum circuit was not provided.
There are two variants of the HPR protocol: a probabilistic variant, which teleports successfully only with some finite probability, and a deterministic variant, which uses an analog of Grover's search algorithm and succeeds with unit probability, but involves a more complex decoding operation. Both protocols take the general form, shown for teleportation of a quantum state |ψ (the generalization to EPR teleportation is straightforward). We now outline the interpretation of each aspect of the above protocol in the context of the Hayden-Preskill thought experiment. For consistency with past literature, we have used different subsystem labels than introduced in the main text-most notably, subsystem D now denotes the coupled qubits, and subsystem C denotes its complement. Subsystem B represents an eternal black hole that is maximally entangled with its past Hawking radiation subsystem B', as represented by a dimension d B = d B EPR pair between the two subsystems. Subsystem A contains the initial state |ψ of an observer Alice's diary. Upon falling into the black hole, the diary's information is scrambled by the unitary time-evolution U acting on the left subsystem l ≡ AB = CD. Far from destroying the information of Alice's diary, scrambling by U in fact allows an outside observer Bob to decode the diary if he has access to any few qubits of new Hawking radiation D, along with the past Hawking radiation B' and an ancillary EPR pair between A' and R', where d A = d A . This decoding relies on OTOCs between subsystem A and D being minimal, a general feature of thermalizing time-evolution after the scrambling time. We describe each of the decoding protocols of Ref. [11] in detail below.

a. Probabilistic decoding: intuition
Although our main focus will be on the deterministic teleportation protocol, we review the probabilistic protocol here for completeness, and as a convenient platform to introduce the intuition connecting operator spreading to the success of teleportation. The decoding operation of the probabilistic HPR protocol consists of projection onto EPR pairs on a subsystems D, D': Perfect teleportation requires d D ≥ d A , and succeeds with probability 1/d 2 A when U is maximal scrambling. The nonunity success probability signifies that the decoding protocol becomes exponentially more complex with the number of qubits to be teleported.
To provide intuition for the protocol's success, we analyze the action of EPR projection on the initial states Q A,l (t) |EPR . We restrict to infinite temperature, i.e. EPR pairs in place of the TFD state, in keeping with the original introduction of the HPR protocol in Ref. [11]. We write Q A (t) as a sum of Pauli strings S on the entire system: Denoting the EPR projector on subsystems D, D' as P EPR,D and writing each Pauli string as a tensor product R = R C ⊗ R D of Paulis on subsystems D and C, we have since EPR D,D | S D,l |EPR D,D = tr D (R D )/d D = δ R D ,1 . Perfect teleportation is achieved when all input Pauli operators on subsystem A have spread to subsystem D, such that every Pauli string S composing Q A (t) has nonidentity support on subsystem D, for all non-identity Q A . In this situation, the EPR projector has eigenvalue 1 on the thermofield double state and eigenvalue 0 in all perturbed states: However, this is no different than projecting onto EPR pairs between subsystems A and A' before time-evolution by U l U * r ! This projection would, of course, have an action Expressed diagrammatically, this equivalence is: for all initial states ψ. However, performing EPR projection between subsystems A, A' before time-evolution is precisely the standard quantum teleportation protocol, applied to subsystems A, A', and R'. The scrambling dynamics of U allow one to perform this teleportation via coupling any subsystem D of the system's qubits.

b. Deterministic decoding
After scrambling, the probability of successful EPR projection on subsystem D, , is exponentially small in the size of subsystem A, the state to be teleported. In contrast to standard teleportation, non-successful EPR projection (i.e. projection onto a different maximally entangled state, not |EPR D,D ) cannot be corrected via an additional decoding operation. This exponential decrease in success probability is overcome in the deterministic HPR protocol, which uses an analog of Grover's search algorithm to search for an EPR pair between subsystems D, D'. The protocol requires O(d A ) steps for completion, again exponential in the number of qubits to be teleported (albeit with half the exponent of the probabilistic decoding).
Grover's search algorithm involves two operations: the first applies a minus sign to the state one is searching for, and the second applies a minus sign to the system's initial state. We will search for an EPR pair on subsystem D, so for the first step we apply W D ≡ 1 − 2P EPR,D = e iπP EPR,D : In the second step, we flip the sign of the initial state (the time-evolved EPR pair between A' and the reference qubit R') by applying W A ≡ U * W A U T : where W A = 1 − 2P EPR,A acts on A', R' to apply a minus sign if the two are in an EPR pair.
The entire Grover protocol is identical to the probabilistic protocol, but with EPR measurement replaced by repeated applications of the two above steps until the EPR pair is found. Displaying, for instance, only the first two iterations: After O(d A ) iterations, the state |ψ is found on subsystem R'.

c. Single qubit deterministic decoding
Two important simplifications occur to the deterministic HPR protocol in the case of single qubit teleportation, d A = 2. The first is that the Grover operator W A is equal to a SWAP operator composed with single-qubit Y operations. To see this, we expand W A in terms of Pauli operators: where in the final equality we used Y r SWAP = SWAPY l , and in the second equality we used the Pauli decomposition for the swap operator between two d A -dimensional boson systems: Expressed graphically, we have The second simplification is that Grover's search for an EPR pair D, D' succeeds after only one step; this is a general result for Grover's search in a d 2 D = 4-dimensional database [7]. It implies that the Grover protocol can teleport one qubit through the circuit: If we only care about the fidelity of the teleported state, we can neglect the final application of U * . Performing the SWAP gate explicitly, and neglecting the action of the final Y operator on R', we have: This exact circuit has been performed in trapped ion experiment [13]. We now make a small cosmetic adjustment, and move the reference qubit R' from the far right to the far left, Sliding U * to the left side using Eq. (9), we have: This is the same circuit appearing the teleportation protocol of Ref. [15,16], modulo the precise form of the coupling.
In the case of EPR teleportation, we would instead have where subsystems R' and A' are in an EPR pair when teleportation is successful. This is the circuit appearing in Ref. [17], modulo the form of the coupling as well as the Y decoding operation. The lack of a Y decoding operation for fermionic teleportation is discussed in Appendix G.

Appendix C: State teleportation fidelity
In Section V C, we provided a detailed derivation of the teleportation fidelity's relation to the teleportation correlators in the case where one teleports one half of an EPR pair. This allowed us to lower bound the fidelity in Eq. (42) and calculate the peaked-size fidelity in Eq. (44). In this appendix we do the same for teleportation of a quantum state, as shown in Fig. 1(a) and outlined in Section II. Our results provide a rigorous foundation for the arguments of Section III, in particular the insertion of the state φ| and the subsequent replacement of |ψ φ| with a Pauli operator Q A .
We begin with the insertion of φ| into the protocol Eq. (16). We do so by inserting the resolution of the identity 1 d A |φ |φ φ| = 1 into the ancillary qubit leg of the diagram for the state teleportation fidelity. We find: Plugging Eq. (18) into this diagram provides unit teleportation fidelity, as described in the main text. When teleportation is successful each of the d A terms of the sum must succeed individually, so the right input state |φ will not affect the success of the teleportation.
where Q i (0, 0) = I i , Q i (1, 0) = X i , Q i (0, 1) = Z i , and Q i (1, 1) = Y i denote individual Pauli operators within the Pauli string. For example, the operator 1 ⊗ 1 ⊗ Z ⊗ 1 ⊗ 1 is represented as x = 00000 and z = 00100. Normally, one would also keep track of the overall phase of Q, but for our purposes the phase will be irrelevant and is dropped in the above notation. Second, we evolve Q under a random Clifford unitary U to obtain Q(t) = U QU † . We consider the circuits shown in Fig. 2, which are composed of random 2-qubit Clifford unitaries laid out in a "brick-layer" fashion. Each of the 2-qubit unitaries is sampled uniformly from the set of 2-qubit Clifford unitaries. While an algorithm exists to perform this sampling directly [130], in practice we find it more convenient to pre-compute and enumerate the entire 2-qubit Clifford set (which consists of 11520 distinct unitaries) 20 . In the binary notation, each 2-qubit Clifford unitary corresponds to a map which acts on the relevant components v, i.e. a unitary with support on the jth and kth qubits updates the values of (x j , z j , x k , z k ). The time complexity of applying the full circuit thus scales linearly with the number of 2-qubit gates and does not otherwise depend on the number of qubits n. As a reference point, simulating a 0D circuit until the scrambling time with 10 8 qubits for a single realization takes approximately one day on a standard single-core processor.
Third, we compute the average operator size distribution and EPR fidelity of the time-evolved operators. For the former, we simply count the size, i.e. number of non-identity terms, of a time-evolved operator Q(t) for a single circuit realization and determine the distribution of sizes with respect to ∼ 10 3 circuit realizations. Depending on the simulation, we either initialize Q with support on a single site (i.e. p=1) or as a p-body operator. In either case, the specific terms in Q (e.g. whether each site is initialized as X, Y , or Z) is arbitrary since the averaged quantities are basis independent.
Computing the averaged EPR fidelity requires an additional average over the initial operators. In particular, for a single circuit realization U , we compute the EPR fidelity using [Eq. (44)]: and η d ≡ 1/(1 − 1/d 2 ), as defined in Section IV B. Note that the first term in θ Q A corresponds to the phase applied by the coupling, while the second term accounts for minus signs associated with transposition and decoding (see Section VII A). The sum in Eq. (E2) is over the complete basis of Pauli operators in subsystem A. For single-qubit teleportation, this consists of three non-trivial Pauli operators and the identity (for which θ = 0), and the sum can be performed explicitly. However, for teleporting many qubits, the number of terms quickly becomes intractable, and we instead approximate the sum by sampling Q A (e.g. ∼ 100 randomly selected operators). To compute the average EPR fidelity, we repeat this procedure for ∼ 100 realizations of U . Finally, we note that the coupling strength g enters the fidelity calculation in a computationally efficient manner; in particular, upon determining the distribution of sizes for a particular circuit realization, we can compute the teleportation fidelity for arbitrary values of g "offline" with no additional simulation cost.

Extended data for 1D and 2D RUCs
Size distribution-The average size and size width for time-evolved operators in 1D and 2D for various system sizes are shown in Fig. 10. In each case, we apply periodic boundary conditions and begin with a single-qubit operator. These results match the functional forms predicted by the KPZ universality class [Eq. (56) and (57)] and allow us to extract the growth parameters {α bulk , α boundary , β bulk , β boundary } = {0.66, 0.70, 1.2, 4.5}.
Multiple qubits-In Fig. 11, we present numerical evidence to support our claim that multiple qubits can be teleported in ≥ 1D short-range models if their operator light cones are non-overlapping (Section VII B). In particular, we simulate the teleportation of n = 5 qubits that are initially evenly spaced in a 1D RUC with periodic boundary conditions. At early times (t < 1300, Region I), we confirm that high-fidelity teleportation is possible for a wide range of coupling strengths, and by measuring the average operator size we infer that during this time the operator light  56) and (57). In particular, we determine α bulk and β bulk from the saturation values (light gray), and α boundary and β boundary from the initial growth rate (dark gray).

Region I Region II
Region III Figure 11. Teleporting multiple qubits (n = 5) in 1D, where the input qubits are evenly spaced in the system (N = 10000).
(a) Teleportation is achieved with high fidelity for t ≤ 1300 (Region I). This corresponds to the regime in which the light cones of the operators are non-overlapping. Interestingly, order-one fidelity can also occur for 1300 < t < 2600 (Region II), when adjacent light cones have overlapped, but only for certain values of g. No multi qubit teleportation is possible for t ≥ 2600 (Region III), as expected from the lack of size addition. (b) The three Regions can be detected by changes in the slope of the operator size as a function of time. In particular, the growth rate decreases when nearest neighbor light cones, then next nearest neighbor light cones, etc. begin to overlap.
cones have not overlapped. In contrast, after the light cones have overlapped, we generally observe a large suppression in the teleportation fidelity. Interestingly, there is one noticeable exception to this reasoning: When only adjacent light cones have overlapped (i.e. 1300 < t < 2600, Region II), high-fidelity teleportation can still occur for specific values of g. This situation corresponds to when the multi-qubit size is a multiple of 2πK/η d g off from the size addition value, e.g.
where m is an integer value. Therefore, strictly speaking, it is possible to satisfy the conditions for many-body teleportation without size addition; nevertheless, it is a non-generic effect that requires finely tuned values of g and evenly spaced input qubits.

Channel capacity for 0D RUCs
An important result of our numerical simulations is substantiating the claim that 0D RUCs exhibit a channel capacity that scales linearly with the number of coupled qubits K. To this end, we first recall that our working definition for the channel capacity is based on setting a threshold for the per qubit fidelity. The most direct way to compute this fidelity would be to take the n-th root of the many-body EPR fidelity; in practice, however, this approach is numerically unstable for large n. Thus, we instead consider a modified protocol for estimating the per Figure 12. Procedure for determining the channel capacity in 0D RUCs. (a-b) For fixed n and K, we compute the per qubit fidelity while sweeping both the evolution time and coupling strength g. (a) The fidelity as a function of evolution time with coupling strength fixed is optimized at the first local maximum, which corresponds to η d gS/N = π. (b) After optimizing the evolution time, the fidelity as a function of the coupling strength g is maximal when g (and correspondingly the average operator size S) is tuned to balance errors due to size addition and the finite number of couplings (see Section VII C for details). The data shown correspond to n = 38 and K = 9000. (c) The channel capacity is defined as the maximum number of qubits that can be teleported while maintaining the fidelity per qubit above a fixed threshold, i.e. 1 − F (1) EPR ≤ 0.07 (dashed line). To determine this number, we fit the optimal fidelity as a function of n (for each K) with a linear fit in log space and compute the intercept of the fit with the threshold fidelity. The fits approximately collapse with respect to n/K, indicating that the channel capacity is linear in K. qubit fidelity where one attempts to send n qubits but only measures the fidelity of one of the n qubits. At infinite temperature and generalizing from one to m qubits, this fidelity is given by: where Q = Q m ⊗ Q u and d A = d m d u , such that Q m acts on the measured qubit(s), and Q u acts on the unmeasured qubits. This can be derived diagrammatically via (E5) Hence, computing the per qubit fidelity, F EPR , is nearly identical to computing the full many-body fidelity, except we sample only over pairs of Pauli operators (Q 1 , Q 2 ) which are identical on every qubit except for one.
We next discuss how to determine the channel capacity from the teleportation fidelity. Specifically, we compute the maximum number of qubits n max that can be teleported above a certain teleportation fidelity, where we fix the number of coupled qubits K and optimize over the evolution time t and the coupling strength g. We consider each of these parameters in turn. First, when sweeping the evolution time and holding all other parameters fixed, the maximum fidelity occurs during the first peak in the time profile; this corresponds to a size η d gS = π/N . After optimizing the evolution time (but holding n and K fixed), the fidelity is non-monotonic with respect to g, owing to the competition among errors due the size addition and finite K. Finally, after optimizing evolution time and g, we determine the maximum number of qubits that can be teleported while maintaining a per qubit fidelity above a fixed threshold value, i.e. 1 − F 1 EPR ≥ 0.07. Our results from this procedure are shown in Fig. 12 and demonstrate that the channel capacity follows a linear trend in K across two orders of magnitude, in agreement with our analytical predictions.

Appendix F: Random circuit calculations
Here we provide more detailed technical calculations of the size overlap and K-size distribution of random Pauli operators of a fixed size. The former is relevant to 0D RUCs (Section VII C), as the vanishingly small overlap of random Pauli strings with size much less than the system size underlies the circuit's ability to teleport multiple qubits at intermediate times. The latter is applicable to all systems when the K coupled qubits are chosen randomly, and quantifies the width introduced to the K-size by this random sampling (Section IV B). In the appropriate limits, these calculations reproduce the intuitive binomial scalings we argued for in the main text.

Distribution of the overlap of two random Pauli strings
We are interested in the probability distribution of the size of the overlap, p (not to be confused with the large-p encoding, which we do not reference in this Appendix) of two randomly samples Pauli strings of fixed size R 1 , R 2 , in a system of N qubits. We expect this to quantify errors to the size addition formula, Eq. (54) in Section VII, for 0D RUCs with large-p encoding (Section VII C), where the assumption of random Pauli strings of a fixed size is appropriate. Our precise derivation is necessarily quite technical, however our final result matches that obtained by intuitive arguments in Section VII [see beneath Eq. (62)].
This probability distribution is computed exactly as a product of various factorials: The numerator computes the number of distinct configurations with Pauli strings of size R 1 , R 2 and overlap p, while the denominator computes the number of distinct Pauli strings of size R 1 , R 2 regardless of the overlap. We are interested in the case where all variables are extensive (scale with N ), but N R 1 , R 2 p. We will proceed by applying Stirling's approximation to each term above, which holds as long as all quantities are large compared to 1. For instance, for dummy variables n, k, we have: We will apply this to a few pairs of factorials in our original expression for P [p]. For convenience, we only keep track of the p-dependence of the probability, and neglect overall constants which serve to normalize the distribution.
Anticipating that the average p will be R 1 R 2 /N , we expand p = R 1 R 2 /N + δ and work to second order in δ. At the end we will show that this is justified. We have: Expanding the logarithm using This gives The last piece is Exponentiating, The first two terms are precisely a Poisson distribution, which has mean R 1 R 2 /N and width R 1 R 2 /N . The exponential is a Gaussian with the same mean R 1 R 2 /N , and a larger width R 1 R 2 /(R 1 + R 2 ). The smaller width determines the width of the product of the two functions, so we conclude: This is what we would expect for drawing p random sites out of N , where each site has independent probability R i /N of being in either Pauli string (Section VII C). The width is subextensive, δ ∼ ε √ N , justifying the higher order terms we neglected along the way.

Distribution of the K-size
Here we are interested in the probability distribution of the K-size of a Pauli string of fixed total size S, with K randomly chosen couplings. Our results substantiate those obtained by intuitive arguments beneath Eq. (34) in Section IV B in the main text. This objective is in fact an identical problem to calculating the overlap: the K-size is the overlap of the K coupled qubits with the S qubits acted on by the operator of interest. We should just replace R 1 → K, R 2 → S, p → n above, where n is the K-size. This is confirmed by comparing the factorial expressions: where the numerator computes the number of distinct configurations with n qubits overlapping the Pauli operator support of size S and K − n qubits not overlapping, and the denominator computes the number of distinct configurations of the K coupled qubits. There are two regimes of interest: when K and S are both extensive, and when S is extensive but K is not. The former provides a more accurate measure of the full operator size (K → N ), while the latter is relevant for probing the channel capacity. Both regimes share the same mean K-size S K and K-size width δS K : This matches our prediction in Section IV B, which was based on a simple scenario of picking K sites, each with a S/N chance of being in the support of the Pauli operator.
At infinite temperature before coupling, each of the above operators has a correlator equal to −1, which is exactly the result for bosonic systems, but without a need for the decoding operation Y . At late times, the coupling e igV applies a relative phase between the identity and non-identity Paulis, giving correlator phases: When g V = π all correlators have the same phase, and peaked-size teleportation succeeds with perfect fidelity at infinite temperature. At intermediate times, peaked-size teleportation of multiple bosonic qubits will succeed just as in bosonic systems. The second option is to send a fermionic qubit, for instance by inserting half of an ancillary FEPR pair. Here we begin with intermediate times, and discuss a modification necessary for late time teleportation afterwards. We represent a single complex fermion with two Majorana operators ψ 1 , ψ 2 , and suppose that the operators' size distributions are tightly peaked, and the size of iψ 1 ψ 2 is twice that of the individual Majorana sizes, denoted S (this assumption of size addition is appropriate in all-to-all coupled systems, e.g. SYK, but would not necessarily hold for e.g. a 1D Majorana chain). The relevant operator correlators after coupling with a bilinear fermionic interaction, as in Eq. (69), are: At gS/qN = π/2 we have perfect teleportation. This generalizes straightforwardly to multiple fermionic qubits: a p-fermion operator will gain a phase i p from sliding across the FEPR pair, and a phase e igpS/qN from coupling.
At late times, the sizes of initial single-body and two-body Majorana operators are equal, since they have saturated the size of the system, and the above operator correlators do not have the same phase. We now show that an alteration of the encoding procedure can rectify this and lead to perfect late time teleportation. This alteration is explained by the HPR protocol, and we derive it by reexamining the equivalence between the HPR and TW protocols in the case of fermionic qubits. Here, the relevant difference between bosons and fermions is that the fermionic SWAP gate is not related to the Grover search operation 1 − 2P FEPR by single-qubit rotations. Since fermions gain a minus sign upon exchange, the fermionic SWAP gate takes the form The single-fermion exp −i π 2 [iψ 1,l ψ 2,l ] gate occurs on the swapped-out fermion and may be neglected. Inserting this be useful shortly. We find a fidelity: (G10) In the peaked-size regime with correlator phases θ R,j , we have At infinite temperature, late times, and g V = π, we have correlator phases θ R,j = 0 for the identity and two-bosonic operator and θ R,j = π/2 for single-body fermionic operators, and find perfect teleportation fidelity: where we define F j = 1 if S L/R,j is fermionic, and 0 if bosonic. We note that for state, as opposed to EPR, teleportation, the above CZ gate turns out not to be necessary. Since coherent superpositions of different fermion parity cannot be created by physical Hamiltonians, which contain only even combinations of fermionic operators, we should only consider teleporting states of definite fermion parity. The CZ gate applies only an overall phase on these states, and so does not affect the success of teleportation.
We can also briefly analyze the low temperature results of Ref. [17] through the lens of operator correlator phases. Here, state teleportation is found to succeed perfectly at low temperatures only when the initial operators are encoded in p-body Majoranas, with p = q/2 + 1, despite the operator correlators having maximal magnitude for any value of p. At the semiclassical gravity teleportation time, the correlators have phases: i p (i) 2p/q iψ 1 ψ 2 (−1) p (i) 4p/q For single-body Majoranas, p = 1, the correlators clearly do not have the same phase-in fact, their phases are nearly identical to their phases at infinite temperature with no coupling-so state teleportation is not possible. When p = q/2 + 1, in the large-q limit, these phases are 1, ±1, ±1, 1, respectively, where the sign is determined by whether p = 1, 3 mod 4. When the sign is odd, it can be corrected via the decoding operation iψ 1 ψ 2 = (−1) N , which applies a minus sign when conjugating fermionic operators. Either case can therefore achieve perfect teleportation.
Appendix H: Teleportation and inelastic scattering at infinite temperature In Section VIII D, we found that strong stringy corrections to a bulk theory of gravity led to peaked-size teleportation as well as a deeply inelastic scattering amplitude. We will now demonstrate that these two phenomena-peaked-size teleportation and inelastic scattering-also coincide at infinite temperature, for arbitrary functional forms of the wavefunctions and scattering amplitudes. As we argued before, for a UV complete theory of quantum gravity, strong stringy (and in general deep inelastic) effects are expected to dominate only at high temperatures, β → 0.
At infinite temperature, the form of the correlator is constrained by the equality C ψ (t; g) * = −C ψ (t; −g).
This implies that C ψ (t) can be written as a real function of ig multiplied by the two-point function: C ψ (t) = ψ l ψ r e −F (ig,t) .
When g = 0, C ψ (t) is equal to ψ l ψ r , implying where f 1 (t) is a real function. Therefore, at this order in g, the infinite temperature correlator is simply the two-point function multiplied by a pure phase, matching peaked-size teleportation [Eq. (37)].
To justify that higher order terms in g are subleading, we need an additional assumption: that the wavefunction of ψ(t) has a saddle point at some momentum k. This is analogous to the boundary assumption that operator sizes are tightly peaked. At early times, this saddle will not be significantly changed by the coupling, since the derivative of the scattering matrix with respect to k will be suppressed by G N , and at early times the time-dependence of the wavefunction will not be strong enough to compensate for this suppression (for example, in semiclassical AdS 2 , we observed competition between e 2πt/β and 1/G N ). In such cases, it is easy to see that Eq. (95) becomes ψ l ψ r times a pure phase linear in g, with higher powers of g suppressed by G N .