Entanglement of Local Operators and the Butterfly Effect

We study the robustness of quantum and classical information to perturbations implemented by local operator insertions. We do this by computing multipartite entanglement measures in the Hilbert space of local operators in the Heisenberg picture. The sensitivity to initial conditions that we explore is an illuminating manifestation of the butterfly effect in quantum many-body systems. We derive a"membrane theory"in Haar random unitary circuits to compute the mutual information, logarithmic negativity, and reflected entropy in the local operator state by mapping to a classical statistical mechanics problem and find that any local operator insertion delocalizes information as fast as is allowed by causality. Identical behavior is found for conformal field theories admitting holographic duals where the bulk geometry is described by the eternal black hole with a local object situated at the horizon. In contrast to these maximal scramblers, only an $O(1)$ amount of information is found to be delocalized by local operators in integrable systems such as free fermions and Clifford circuits.


References
28 Chaos in classical systems is described by the sensitivity of phase space trajectories to initial conditions. Systems displaying chaos will generically have nearby trajectories diverge exponentially at early times, characterized by a Lyapunov exponent. One can think of this sensitivity to initial conditions as a manifestation of the butterfly effect; a small change, such as a butterfly flapping its wings, can have extraordinary consequences on the state of the system at later times.
Quantum chaos is an old topic with many developments (see e.g. [1,2]) that addresses the connection between classically chaotic systems and their underlying quantum mechanics; how do highly nonlinear classical dynamics emerge from the linear unitary evolution of the Schrödinger equation? Recently, there has been considerable excitement across multiple fields of physics due to a quantum manifestation of the butterfly effect [3][4][5] characterized by out-of-time-ordered four-point correlation functions (OTOCs) where V, W are local operators in the Heisenberg picture and · · · β specifies that the correlator is evaluated in a finite temperature state of inverse temperature β. This directly probes the spreading of a local operator's spatial support. In analogy to classical chaos, (1) can exhibit a quantum Lyapunov exponent, λ L , at early times and tends to zero at late-times, describing the "scrambling" of the quantum information of the initial state. Importantly, this exponential behavior of the OTOC is not a generic feature of quantum chaotic systems. In fact, it has only been found for large-N theories [3][4][5][6][7][8][9][10] and does not generically occur in realistic (finite-N ) quantum chaotic systems 1 (see e.g. [11,12]).
The OTOC also leaves certain information-theoretic questions open about chaos. How close are the quantum states of subsystems with and without the perturbation? How much and how fast is information delocalized (scrambled) by the perturbation? To what extent does the choice of local operator influence the scrambling process? In this paper, we address these questions by studying the local operator entanglement, quantum correlations of local operators 2 . We can study the local operators directly by computing correlation measures not in the original Hilbert space, H, but the doubled Hilbert space of endomorphisms End(H) In practice, this is done by "flipping the bra vector to a ket" |O(x, t) ≡ N m,n e i(En−Em)t O nm (x) |n 1 |m * 2 , (4) where we have expanded in an energy eigenbasis and imposed an appropriate normalization. Entanglement in the Hilbert space of local operators has been considered previously in Refs. [14][15][16][17][18][19][20].
Let us now try to understand how correlations in this state, (4), characterize the butterfly effect by showing how information flows under this time evolution. Consider the trivial case where the operator, O, is the identity, 1. The identity should have no effect on the state. This is shown in Fig. 1 where the time evolution operator moves around quantum information. For integrable systems, information that is initially localized will remained localized at time t in the sense that it will not mix in an ergodic manner in the local Hilbert space. This FIG. 2. The Heisenberg evolution of a (nontrivial) local operator is shown. For integrable channels (left), initially localized information remains largely localized even after the perturbation by operator O. In contrast, for chaotic channels (right), local information delocalizes in time. Then, the perturbation by operator O effects the state such that this information cannot recohere under backwards time evolution, but rather continues to grow (decoheres).
FIG. 3. We show the partitioning used throughout the paper. Left: the BOMI configuration for symmetric intervals A and B. Right: the TOMI configuration for finite interval A and semi-infinite intervals B1 and B2.
is due to the constraints imposed by integrability. This lack of "local ergodicity" does not exclude the possibility of spreading in space as is seen in integrable system exhibiting diffusion. For chaotic systems, the initially localized information becomes spread out at time t and mixed in the local Hilbert space without the severe constraints imposed by integrability. This is when the operator (identity) is inserted. Because the identity acts trivially, the backwards time evolution brings the information back into a localized packet whether or not the system is integrable 3 . This means that the mutual information between subregions is simply proportional to their overlap in the spatial direction.
We progress to nontrivial operators. When these operators are inserted, they may scatter the information (see Fig. 2). For integrable systems, the localized quantum information will remain localized but it may be transmitted to a different location than it started at when it is evolved back to t = 0. For a chaotic system, the butterfly effect implies that the local perturbation created by the operator may ruin its coherence, hence it remains delocalized after the backwards time evolution. Given long enough times, the information will be spread out over the entire system. We propose that an illuminating diagnostic of the amount of quantum information initially in region A that is scrambled by operator O is the tripartite mutual information (TOMI) 4 , defined as where B = B 1 ∪B 2 is the entire output Hilbert space, H 2 . This characterizes how much total (classical + quantum) information from A is lost unless the entire output system is measured. The local operator entanglement allows us to understand how different operators scramble information. Analogously, we also study tripartite operator logarithmic negativity (TOLN) to characterize the purely quantum information that is scrambled. This is defined by replacing the bipartite operator mutual informations (BOMI) on the right hand side of (5) with logarithmic negativities (BOLN). We show the generic setup in Fig. 3. While the operator choice for OTOC may be seen as a disadvantage because it can be misleading (e.g. spin-spin OTOC in the Ising model [4]), it should be seen as an advantage for local operator entanglement because the mutual information probes correlations of all operators; not all butterflies have the same effect.

A. Summary of results
In the rest of the paper, we have many technical results that the casual reader may not wish to sift through. Here, we summarize our central findings. We also present a cartoon summarizing results for I 3 in Fig. 4.
a. Random Unitary Circuits Random unitary circuits are tractable toy models of local Hamiltonians displaying chaotic phenomena. In Section I, we put forward an effective description of the entanglement dynamics of local operators in terms of the free energy of a membrane in spacetime. This effective description is discussed in further detail for Rényi entropy, logarithmic negativity, and reflected entropy by mapping the random unitary circuit to a classical statistical mechanics problem in Appendix A.
We find for Haar random unitary channels that local information is entirely delocalized by the local operator, regardless of the operator chosen. This manifests by the tripartite information increasing in magnitude as fast as is allowed by causality, ultimately saturating to the lower bound on all quantum systems which is proportional to the number of degrees of freedom in subsystem A.
In contrast, when the quantum channel is composed of random unitary elements from the Clifford group instead of the full unitary group, we find that a very small amount of information is scrambled; this value is independent of system sizes but dependent on operator choice. This is notably different than the observed maximal scrambling behavior of Clifford circuits for the unitary time-evolution operator [22]. We explain this discrepancy by emphasizing the importance of Clifford gates being unitary 3-designs. Moreover, we find that depending on the operator, the mutual information and logarithmic negativity behave differently. This demonstrates how quantum and classical information may be scrambled in different ways in quantum channels.
b. Integrable System We study free fermions as an example of an integrable system. In particular, we consider the tight-binding Hamiltonian for simplicity. The local operator is taken to be a fermion parity operator acting on a single site. The resulting local operator state is Gaussian, allowing us to employ the correlator method to compute the local operator entanglement entropy.
When the input and output subsystems are spatially identical, the mutual information starts of at a maximal value before beginning to dip at some time determined by causality. After the wave front of the operator leaves the subsystems, the BOMI begins to relax back to its initial value, though we do not have a proof that the BOMI fully relaxes back to its original value due to finite size effects. This O(1) change (not extensive with system size) in the BOMI is a signature of the (trivially) integrable nature of free fermions. Very little (if any) information is scattered or delocalized. The TOMI is quite similar. It is initially zero but decreases once the operator is within the subregion. Eventually, the TOMI attains its most negative value before it slowly relaxing back to zero. The latetime behaviour of TOMI indicates the lack of scrambling from operators in the free fermion system.
c. Chaotic CFTs The local operator entanglement of holographic 2D CFTs are studied in section III. These are conformal field theories with large central charge and sparse low lying spectra, and are considered maximally chaotic due to their early time exponential behavior in the OTOC [3][4][5]. Another way in which they saturate the fundamental bounds of quantum information scrambling 4 is the decay of the tripartite entanglement of the timeevolution operator [22][23][24].
When the input and output subsystems are symmetric, the BOMI for the local operator begins at its maximum value. After the operator has had time to reach the intervals, it begins to decrease linearly at the maximum rate allowed by causality. Unlike the free fermion BOMI, the BOMI for holographic CFTs decreases all the way to zero. This tells us that the local operator eventually delocalizes the information completely and is consistent with the expectation that these conformal field theories are maximally chaotic. The TOMI for holographic CFTs is also found to decrease from zero to a maximally negative value at the maximal rate where S reg.
A is the UV finite thermodynamic entropy of subregion A at a temperature determined by a regulator i.e. it does not contain the standard UV divergence of von Neumann entropy in continuum theories due to short distance modes near the entangling surface. We stress that these results are significantly stronger than analyses of operator entanglement in the past because this maximal scrambling of information occurs regardless of any details about the operator. The smallest perturbation entirely destroys the quantum information of the state.
An additional notable phenomenon is that the BOMI and TOMI have step function discontinuities associated to when the local operator enters and leaves the associated subregions. The magnitude of these step functions is determined by the conformal weight of the operator. Only for heavy operators (∆ ∼ c) are they discontinuities macroscopic.
Finally, we note that these findings precisely match with the results for the Haar random unitary circuits with bond dimension q in Section I once identifying the bond dimension with the Cardy density of states q = e πc 3β (7) where c is the central charge and β is the effective temperature which is just a regulator for us. One caveat is that the membrane computation for the random unitary circuits does not have the discontinuities previously mentioned. This discrepancy may either show a difference between the two theories or the analogy may be restored once we account for O(1) contributions in the membrane theory.
d. Holography In Section III C, we identify the geometry dual to the local operator state (3). Because this state lives in two copies of the original Hilbert space, it is natural that the dual geometry has two identical asymptotic boundaries. This is the eternal black hole dual to the thermofield double state with the temperature playing the role of the cutoff. The local operator perturbs the eternal black hole in a similar manner to Refs. [3,25]. It is a massive particle that backreacts on the geometry. We are able to compute the operator mutual information directly from the geometry using the Ryu-Takayanagi formula which precisely matches the CFT calculation.

I. RANDOM UNITARY CIRCUITS AND MEMBRANE THEORY
In this section, we motivate intuition by comparing two effective theories of entanglement dynamics, the quasiparticle picture and the membrane theory, which model integrable and chaotic dynamics respectively.
The quasi-particle picture has been proposed as a universal description of entanglement dynamics in integrable theories [26][27][28]. This posits that when an integrable system is sufficiently excited above its ground state, the entanglement between subsystems may be entirely accounted for by quasi-particle pairs that carry entanglement content that travel at known speeds. These dynamical inputs may be fixed by thermodynamic Bethe ansatz techniques. The entanglement in inherently bipartite by construction because only Bell pair-like correlations are accounted for. This description largely matches our results for free fermions as the local operator state is an excitation above the ground state of the Hamiltonian Severe breakdowns of the quasi-particle picture occur for non-integrable systems because multipartite entanglement becomes increasingly important (see e.g. Refs. [29][30][31][32][33][34][35][36]). Thus, the information about entanglement can no longer be carried by quasi-local objects. Recently, a compelling case has been made that, for quantum chaotic systems, the entanglement dynamics are captured by a codimension-one membrane in spacetime, a manifestly non-local object [37,38]. The dynamical input into this membrane theory is the tension of the membrane which may be explicitly computed in certain cases. A particular instance where this may be computed is for Haar random unitary circuits. In this section, we will study these circuits and adapt the membrane theory to local operator entanglement.

A. Haar random unitary circuits
We begin with a simpler problem of computing just the late-time behavior of local operator entanglement by modeling the random unitary circuit as one big Haar random operator 5 . In the following sections, we will refine these results in order to understand early-time behavior and the membrane theory.
The advantage of modeling chaotic dynamics with Haar random unitary circuits is that analytic results are 5 It has been shown that local random unitary circuits are approximate k-designs at a circuit depth scaling as O(N k), where N is the total number of qudits [39]. In this subsection, we will need at most k = 4, so this quantifies what we mean by "late-time." tractable due to well known results from random matrix theory. In general, we will only need the Weingarten formula which computes the integral of monomials of unitary operators with the Haar measure [40] [ where d is the rank of the unitary. The sum is over elements of the permutation group and Wg is the Weingarten function. We will consider the large system size limit such that the term with στ −1 = e (the identity) is dominant and approximately leading to In Fig. 5, we show the diagrams that compute the average of Trρ 2 (the purity) for the states |U (t) and |O(x, t) . Assuming that the average and the logarithm approximately commute for large system sizes, this computes the average second Rényi entropy. While it is in principle possible to compute all the "average Rényi entropies", which we denoteS and analytically continue to n = 1, we will only consider n = 2 for simplicity. An interesting fact is that while, in general,S n = S n , in the von Neumann limit 6 , they are equal, allowing one to properly find the average entanglement entropy. In fact, such an analytic continuation was computed for operator entanglement of the time evolution operator in Ref. [42]. We warm up with the time-evolution operator and then proceed to local operators. For a q L × q L time-evolution operator, the corresponding state is normalized as Here, the two indices represent the input and output Hilbert spaces respectively as the total state is an element of the Hilbert space of U (q L ). The density matrix is then We can bipartition both the input and output systems arbitrarily The reduced density matrix on AC is where summation over repeated indices is implied. The average purity is then The Wiengarten formula involving only four unitaries needed for the above is simple enough that we may write it out explicitly in terms of Kronecker deltas 6 This leads to Thus S (2) Here, a, b, c, and d are the number of qudits in subsystems A, B, C, and D respectively. Let's look at the tripartite mutual information, taking A to be O(1) and B 1 and B 2 to be semi-infinite (scale as e.g. L/2). Then, where a is the length of the subsystem. Therefore, TOMI tends to −2a log q. If we had taken the output subsystem to be size L − , theñ Thus, the condition for nontrivial mutual information is This means that at late times in a chaotic quantum channel, one needs at least (L − a)/L of the system to recover any information from A. We now progress to local operator entanglement. As shown in Fig. 5, we have twice the number of unitaries to worry about. In the limit of large Hilbert space dimension such that we can make the approximation of (10), this is still tractable by brute force. Our state associated to the local operator is so the density matrix is Again, we bipartition the input and output Hilbert spaces so that we may define the reduced density matrix on AC The Haar averaged purity is We can see from (10) that this integral will involve 24 terms, even after the approximation. After the contraction of many Kronecker delta's, one finds at leading order In general, the second term will be subleading. The immediate consequence is that the answer is operator independent. We then find The mutual information is To find the TOMI, we take c, d = L/2 The second term is subleading in the scaling limit. Because this saturates the bound on tripartite mutual information, we know the state must be approximately maximally entangled, so all of the Rényi entropies are approximately equivalent. Thus, we find This is the lower bound on I 3 allowed by quantum mechanics. Given that this is an important point that will continue come up in this work, we now show the simple derivation. I 3 is composed of three individual mutual informations. The mutual information is positive semidefinite, so The mutual information between subregion A and the entire output B 1 ∪ B 2 is a time independent quantity in finite-dimensional systems. The input and outputs are maximally entangled by construction, so the mutual information is twice the logarithm of the Hilbert space dimension This is only well-defined for subsystems in discrete models. When we move on to quantum field theories, we must regulate the Hilbert space with (6) as the analog. This will lead to strange effects such as the regularized dimension of the Hilbert space associated to a finite region being time-dependent.

B. Random Clifford circuits
The random unitary calculations above effectively capture the late-time behavior of chaotic channels. However, to study interesting early-time behavior, we have a couple of options. One option for modeling stronglyinteracting dynamics for large systems sizes is random Clifford circuits. These are random unitary circuits that are composed of the Clifford group: phase, Hadamard, and CNOT gates. Using the stabilizer formalism, measuring entanglement in these circuits is tractable with computational times scaling polynomially. For entanglement of the unitary evolution operator, it was shown that these circuits maximally scramble and behave extremely similarly to holographic quantum channels [22]. However, it is also known that these circuits have pathological OTOC [38]. Even though the late-time average OTOC is zero, the variance is order 1. This is due to Clifford gates being unitary 3-designs. As is clear from Fig. 5, one must take the fourth moment of the unitary group in order to computeS (2) for local operators. Because Cliffords are 3-designs, a priori, they may have distinct behavior from the Haar random unitaries and chaotic channels in general.
Indeed, this is what we find. The late-time value of the TOMI is a constant, independent on the size of subregion or total system size. However, it does depend on which operator we are evolving. Moreover, the operator has left and right-moving components, so given a configuration where the input subregion spatially overlaps with the partition of the output Hilbert space, twice the information will be scrambled compared to if the operator is initially outside of the interval overlaps. For simplicity, we have used the three local unitary operators that generate the Clifford group as our local operators. This behavior is reminiscent of integrable theories, however we find that there are no recurrences, even for the finite system, a feature of the stochastic time evolution.
In Fig. 6, we show the time evolution of the operator mutual information and operator logarithmic negativity 7 for symmetric intervals. At early times, when the operator has not yet reached the intervals, the correlations are maximal, proportional to the area of the intervals. However, once the operator has time to reach the intervals, the correlations decrease because some information is being scattered as in Fig. 2. We observe the following interesting features that distinguish this quantum channel from chaotic channels, particularly the Haar random unitaries that we have studied. (1) The amount of information scattered is independent of the size of the subregions. (2) The amount of information scattered is operator dependent. In particular, the CNOT gate scatters more than the Hadamard gate. (3) The quantum and classical information delocalize differently. For certain operators, the saturation value of I 3 is equivalent to E 3 , while for others the saturation value of I 3 has greater magnitude than E 3 (see Fig. 7). The latter indicates that some purely classical information has been scrambled.

C. Membrane theory
While the quasi-particle picture is an effective description of entanglement propagation for all integrable systems, it fails to capture the qualitative features of entanglement production in chaotic systems. It is highly desirable to obtain an analogous universal description of entanglement dynamics for chaotic quantum systems. Recently, it has been proposed that these chaotic theories have effective hydrodynamical descriptions where the von Neumann and Rényi entropies may be computed by the area of a spacetime codimension-one brane, M, which is characterized by its tension, T , [17,37,38,44]  . Notably, the CNOT gate scrambles both classical and quantum information while the Hadamard gate only scrambles quantum information. This can be seen by the fact that E3 is larger than I3 for the CNOT gate, but identical to I3 for the Hadamard. Again, the operator is inserted 5 sites to the left of the symmetric intervals.
where x is the position of the membrane and v is the space-time velocity of the membrane (dx/dt). The membrane M A is the extremal surface with respect to the integrand of (35) that is homologous to subregion A. Though derived from finite-dimensional quantum circuits, there are strong parallels of this construction to the holographic description of von Neumann and Rényi entropies in the AdS/CF T correspondence [45][46][47][48]. In essence, both prescriptions require finding the area of an extremal surface that is homologous to the subregions of interest.
With motivations from the holographic description of reflected entropy [49], it was later proposed that this mixed state entanglement measures may also be computed by the area of a different codimension-one membrane [35] S (n) We denote this membrane E W because in the language of AdS/CFT, this surface is the entanglement wedge cross section, a natural geometric object in the bulk that generalizes the Ryu-Takayanagi surface [50]. In the membrane theory, the entanglement wedge of A ∪ B is the codimension-one spacetime region whose boundary is is then defined as the extremal surface separating subregions A and B within the entanglement wedge. We note that for random unitary circuits with large local bond dimension q, the spectrum is effectively flat and the line tension is thus equal for all Rényi's. While these membrane descriptions were well motivated by analysis of random unitary circuits for states after global quenches and for operator states of the unitary evolution operator, for local operator states, only a highly symmetric case has been analyzed [17]; we show that it is straightforward to generalize to generic configurations, giving an intuitive explanation for our late-time result (32) from the previous section. The line-tension for the local operator entanglement is dependent on space and time, not just velocity. The line tension is the same as it was for unitary operator entanglement within the light cone of the local operator while outside the light cone This may be quickly seen by considering the minimal cut through the quantum circuit displayed in Fig. 8. A more sophisticated derivation is explained in Appendix A by mapping the random unitary circuit to a classical spin model. In the large q limit, the number of bonds cut is asymptotically equal to the Rényi entropy. This becomes (37) & (38) in the scaling limit. It is instructive to work out a couple examples. We show three time steps for the entanglement entropy of symmetric intervals of length l in Fig. 9. Initially, the two intervals A and B in the input and output Hilbert spaces, respectively, are maximally entangled with one another, so their total entropy is zero. Using (38), the corresponding minimal membrane is of zero area (left) because it has v = 0. Once the light cone of the operator reaches outside of the intervals, it can break the entanglement between them and entangle them with the rest of the system. This is seen in the intermediate time step where the area of the minimal membrane grows linearly. At sufficiently late times, the entropy saturates to its maximum value, 2l log q, which is described by the disconnected regime on the right in Fig. 9. In summary, we find Because the individual entropies of the intervals are constant in time, we find the mutual information is We can also compute the full time dependence of the tripartite operator entanglement shown in Fig. 10. We take A = (0, l), B 1 = (−∞, 0), B 2 = (0, ∞) for simplicity, but the late-time value will be universal. We find which leads to a tripartite mutual information of The saturation value is identical to the late-time result of the previous section (32) and is of maximum magnitude. We can play the same game for reflected entropy. However, the relevant membrane is now given by the extremal cross section of the codimension-zero region bounded by M and the spacetime boundary. For the symmetric case shown in Fig. 9, we find step function behavior This is dramatically different than the mutual information. Similarly extreme differences were found between the mutual information and reflected entropy for irrational CFTs and random unitary circuits following a global quantum quench [35] and in operator entanglement of the reduced density matrix [51]. Interestingly, this discrepancy has never been observed for integrable theories. We would like to better understand this physically because its information theoretic implications are somewhat puzzling as discussed in Ref. [35]. Multipartite entanglement must play a significant role, but the problem certainly deserves further attention.
For the semi-infinite configuration shown in Fig. 10, we have Because the other terms in the tripartite quantities are constant in time for the given configuration, we find Some of the CFT techniques that we use in subsequent sections are specific to von Neumann entropy and will not apply to the negativity and reflected entropy, so we do not evaluate these in CFT. While these calculations seem tractable, we leave this to future work and assume that the random unitary circuit analysis precisely describes the CFT computations once identifying the bond dimension with the Cardy density of states (7), q = e πc 3β .

II. FREE FERMION SYSTEM
In this section, we will compute the local operator entanglement for a (1+1)-dimensional lattice free fermion system described by a quadratic Hamiltonian, H = x,y c † x H xy c y , where c x /c † x are the real space fermion annihilation/creation operators at site x on the lattice. Specifically, we will take the nearest-neighbor tight-binding Hamiltonian to be our free fermion Hamil-tonian.

This Hamiltonian is diagonalised by a Fourier transform
where L is the total length of the system. The tight-binding dispersion relation is E k = −2t cos k. As the computation of local operator entanglement for generic RCFTs is rather involved, we will study the free fermion system numerically instead.
In order to compare (match) free fermion numerics with field theory results, UV regulators will have to be introduced in the operator state and taken to be much larger than the lattice spacing in order to suppress lattice effects [22,23]. However, the introduction of UV regulators in the local operator state will greatly complicate the expressions so we refrain from doing so. We thus consider the local operator state with no regulators where H = H B − H A , O A = O ⊗ I and |Ω is the infinite temperature thermofield double state. The numerical results are not expected to agree precisely with field theory calculations although they should capture the overall qualitative behaviour. Here, the maximally entangled state can be written in terms of real space fermions as As an local operator, we choose to work with the single site fermion parity operator at some arbitrary site z, This operator is the exponential of a quadratic fermion operator, so it is a Gaussian operator. Hence, we can utilize the correlator method to compute operator entanglement entropies. Since the parity operator squares to one, the state is already normalized. The initial local operator state is then given by the maximally entangled state with a sign flip at site z Noting that the time-evolution of the fermion operators under H B − H A is given by where (−1) A = −1 and (−1) B = 1, the correlation matrices are given by The relevant operator entanglement entropies can then be computed by diagonalizing subblocks of the correlation matrix.
We plot the results in Fig. 11 and find very different behavior than the random unitary circuits. In particular, the BOMI decreases from its initial value once the operator has time to enter the subregions. Then, it relaxes back once it has left the subregion. A similar analysis is made for the TOMI. It is presently unclear whether the values relax all the way back to their initial value because the lattice model has slow quasi-particle modes that take a very long time to travel through the intervals. This, however, is a moot point because we clearly see that the operator scrambles very little (if any) information i.e. the information in the free fermion channel is robust to perturbations, a quality we expect to be generic for integrable systems.

A. Setup
In this section, we compute the local operator entanglement for 2d conformal field theories at large central charge. These represent candidates theories possessing bulk gravitational duals well described by semi-classical physics. We now set up the path integral representation of operator entanglement. We take a local operator situated at position X in the Heisenberg picture and expand in the energy eigen basis as We then perform the state-operator map to create the local operator state in a doubled Hilbert space where N is a normalization constant that ensures that the state has unit norm. More specifically, the normalization squared is Crucially, we have included regulators 1 and 2 in order to smear the operator and cut off the high-energy modes. Consider the density matrix corresponding to the state (54) in Euclidean signature 12 FIG. 11. We show the evolution of the operator entanglement for the fermion parity operator at a single site (−1) nz . Left: the change in the operator mutual information is shown for symmetric intervals of length 50 with the operator inserted 25 sites away from the intervals. After a time corresponding to the distance from the operator to the intervals, the BOMI begins to drop. Once it passes through the intervals, it relaxes but asymptotes back to the value at which it started according to a power law with exponent between −1/2 and −1. This indicates that little information has been scattered. Right: the TOMI is shown for lA = 50 and the partition of B lying at the center of A. The local operator is again located 25 sites away from A. We see minor delocalization of information when the operator has support in region A, but relaxes at late times back to zero.
After performing our computations of BOMI and TOMI in Euclidean space, we will perform the analytic continuation We bipartition the Hilbert spaces and write the density matrix elements in terms of the field configurations on these bipartitions where the states corresponding to complex conjugated fields are defined by Φ a ,Φ b |n = n|Φ a , Φ b , n|Φ a ,Φ b = Φ a , Φ b |n . Consider two subsystems A and B in the first and second Hilbert spaces respectively. The reduced density matrix for the union of these two regions is In order to compute the entropy and thence the mutual information, we must perform the replica trick where we cyclically glue the path integrals defining the above state. This replica manifold is shown in Fig. 12. Equivalently, we can consider a replicated theory on a single cylinder and introduce Z n twist operators that implement the cyclic gluing, where n is the number of replicas. The circumference of the cylinder is In the replicated theory, we consider the operator O n which is the tensor product of the operators from each Contents 1 OTOC in different limit 1 The operator entanglement entropy for two disjoint regions A and B is then computed by the following correlation functions where we introduce the following coordinates on the 13 cylinder: These coordinates are to be analytically continued at the end of the calculations: Similarly, the local operator entanglement entropy for the individual intervals can be obtained by tracing out the other interval, Performing the standard cylinder to plane map z = e 2π β w , the two-point function in the normalization is simply given by

B. Bipartite and tripartite information
The local operator Rényi entropy for each individual set-up must be computed separately as the monodromies of the conformal blocks are highly dependent on the spacetime locations of the operators in the correlators. The computations are thus rather repetitive so we leave them in Appendix B. In computing entropies, we use the four-point functions given by the HHLL vacuum conformal block [52]: Here, in our case, the light operators are the twist operators, h L = h n , h p = 0 for the vacuum conformal block, and where we are considering scalar operators with h O =h O .
The six-point function in (62) can be approximated by two four-point functions using the OPE σ n (1) × σ n (x,x) ≈ I + O((1 − x) s ) where s ∈ Z. Each fourpoint function will contain the local operators O, O † as well as the twist operators σ,σ. If both twist operators in each four-point function correspond to the end-points of the same interval, we say that the six-point function is computed in the disconnected channel. On the other hand, if both twist operator in each four-point function belongs to separate intervals, we say that the six-point function is computed in the connected channel. Holographically, the disconnected channel correponds to bulk geodesics starting and ending on the endpoints of the same interval, while the connected channel corresponds to geodesics beginning on an endpoint of one interval and ending on an endpoint of the other interval.
Combining the various local operator entanglement entropies listed in Appendix B, we obtain the bipartite local operator mutual information. In the following, we list the mutual information (for the connected channel) for various subsystem configurations. Below, we set 1 = 2 .

14
(ii): Partially overlapping intervals I X < X 2 < Y 2 < X 1 < Y 1 : Two comments are in order. (i) The bipartite local operator mutual information is constant until both intervals are within the light cone of the local operator. (ii) We note that in these results, the time-independent part of the bipartite local operator mutual information in the β → 0 limit is simply given by 2πc 3β l A∩B , where l A∩B is the length of the overlap of the two intervals A and B.
With the various bipartite local operator mutual information at hand, we proceed at last to the tripartite local operator mutual information. As a specific setup, we con- . and insert the local operator to the left of subregion A so that Y 3 < X < X 2 < Y 2 < X 1 < Y 1 . We also send Y 1 → ∞ and Y 3 → −∞ so that B 1 and B 2 form a bipartition of the output. The tripartite mutual information is obtained by taking the appropriate limits of (69), (71) and (72). Those expressions are for the connected channel, and the actual bipartite local operator mutual information is given by I AR = Max (I con.
We plot I AB1 , I AB2 , I AB and I(A, B 1 , B 2 ) in Fig. 13 both for light and heavy local operators. The bipartite local operator mutual information for each semi-infinite interval I AB1 and I AB2 vanish after a certain time while the bipartite local operator mutual information I AB for A and the entire output B remains constant even at late times. The tripartite local operator mutual information thus converges to −I AB , where S reg.
is the regulated entanglement entropy for A [23]. This saturates the lower bound for I 3 as fast as is allowed by causality just like the random unitary circuits in Section I and the tripartite unitary operator mutual information of holographic CFTs and random unitary circuits in Refs. [22,23]. Small discontinuities in the local operator mutual informations depend on the weight of the local operator. For a light operator, these discontinuities are small. Here, the dependence of the bipartite and tripartite local operator mutual information on the local operator comes from their conformal dimension, δ =δ = 1 − 24 c h O . When the local operator is light, i.e. h O < c 24 , 0 < δ < 1, as it enters the expressions of bipartite local operator mutual information in the form of Trig. Function δ , its logarithm is much smaller than the kinematical factors that enter the epxressions in terms of exponentials. On the other hand, when the local operator is heavy, h O > c 24 , δ is purely imaginary, and hence the trigonometric functions become hyperbolic functions, and the piecewise constant operator dependent terms give rise to more noticeable discontinuities. (Here, keep in mind that the bar does not refer to complex conjugation but instead refers to the anti-holomorphic conformal dimension.) Since I con.
AB1 and I con. AB2 decay to zero, and the late time value of I con.
AB and S reg.
A are operator independent, the tripartite local operator mutual information for heavy operators still satisfy the equality (73).
In Section I, we found precisely the same results as the large-c calculations for light operators i.e. no discontinuities. This can be quantitatively verified once using the identification (7). It is interesting to consider if and how the discontinuities created by heavy operators can arise in the membrane theory. In Section I, we had neglected O(1) contributions that can arise from the initial state. It is reasonable that by carefully accounting for these O(1) contributions, one can find the discontinuities that depend on the specific operator.
Before concluding the section, let us comment on one of the key differences between the unitary operator entanglement and the local operator entanglement. During the analytic continuation, the cross-ratios for the case of unitary operator entanglement are real and do not follow any non-trivial trajectories [23]. On the other hand, the cross-ratios for local operator entanglement are complex and can encircle the branch point at the origin during analytic continuation as shown in figure 14. The nontrivial time-dependent behaviour of the cross-ratios is a direct consequence of the insertion of local operators. As a result, the dominant conformal block can acquire a monodromy, which contributes to the late-time behavior. This behavior is essentially the same as the cross-ratios that appear in the computation of OTOCs [4,53,54]. In some sense, one can think of the four-point functions in the computation of S AB as an OTOC with the operators V and W from (1) being the twist field and local operator O respectively. The exponential decay of this OTOC (coming from the monodromy of the vacuum conformal block) manifests itself as the linear decrease in bipartite local operator mutual information.

C. Holographic description
We conclude this section by introducing the geometry that is holographically dual to (54). It is the two-sided black hole with a massive object discussed in Refs. [6,25,55].
In the holographic CFT, we study the time evolution of BOMI and TOMI on R 1,1 . Therefore, the gravity dual which we compute the operator entanglement entropies is the geometry in the AdS-Schwarzschild patch where √ M = 2π/β. In order to consider the gravitational dual to (54) (with 1 = 2 = ), we take the period to be β = 2( 1 + 2 ) = 4 . The Kruskal coordinates are related to AdS-Schwarzschild patch of two wedges of AdS 3 as: The resuling metric is where U and V are defined in the region −1 < U V < 1. The conformal boundaries, where the two copies of the CFT live, is at U V = −1, the horizons are at U V = 0, and the singularities are at U V = 1. The regions which correspond to the left and right CFTs are defined by We now place a massive object located at in the coordinate of (74). The geometry back-reacted by the massive object can be constructed by first considering the metric of AdS 3 black hole in the global coordinate, where the black hole (of mass m = µ/(8G N R 2 )) is located at the center of the cylinder. The parameter µ is related to the conformal dimension of local operator O: This metric can then be mapped by the following boost and coordinate transformation so that the resulting metric describes the massive object at the origin of the coordinates in the AdS-Schwarzschild patch where (U, V, ψ) is the Kruskal coordinate. In terms of (U, V ) coordinate, r is given by Since (U, V, ψ) are transformed to (z, t L,R , x) as in (75), the boost parameters Λ 1 and Λ 2 can be determined by requiring the massive object at (z, t, x) = (α, 0, 0) in AdS-Schwarzschild patch corresponds to r = 0 in global coordinate: The above transformation gives the metric in Kruskal coordinates that takes into account the back reaction from the massive object. For our purpose of computing entanglement entropies, let us write the global coordinate for the two wedges (r L , τ, ψ L ) and (r R , τ, ψ R ) in terms of AdS-Schwarzschild coordinate: for the left (p = L) and right boundaries (p = R), respectively. The left and right regions in Kruskal coordinate corresponds to the geometries, which are back-reacted by the massive object and their asymptotic regions are AdS-Schwarzschild. Now, in terms of global coordinates the holographic entanglement entropy S A is given by [29] S A = c 6 log where the boundaries of subsystem A are at (r 2 , τ 2 , ψ 2 ) and (r 1 , τ 1 , ψ 1 ). Since the holographic entanglement entropy is diffeomorphism invariant, we are able to compute holographic entanglement entropy in the two-sided black hole with the massive object by using (88). As the final step, let us determine the parameter α. The expectation value of energy momentum tensor for the state in (4) is equal to the three point function which is universal in 2d CFTs. This can be done by requiring the energy density in the gravity side should be equal to that for the state in (4) with 1 = 2 = : ( This is related to its holographic counterpart T L,R tt hol as T L,R tt hol = T L,R xx hol = 1 2π · T L,R 00 [56,57]. Using this dictionary, α is determined as Then, the massive object is pinned to the horizon of AdS-Schwarzschild coordinate and the origin of Kruskal coordinate as in Fig. 15. In AdS-Schwarzschild coordinate, the subsystems A and B are defined as A = {z, t L , x|z = a, t L = t, X 2 < x < X 1 } and B = {z, t R , x|z = a, t R = t, Y 2 < x < Y 1 }, respectively, where a 1 is the inverse UV cutoff. For M α 2 = 1, the small a expansions of the variables in (88) for A, B and A ∪ B are given to leading order by where i = 1, 2. Here, we assume |X 2 | < |X 1 | and |Y 2 | < |Y 1 |, and we introduced for p, q = L, R. The above equations can be used to compute the operator entanglement entropies in terms of Poincare coordinate. By choosing the minimum one of (88) to be the operator entanglement entropy, we verify that the time evolution of holographic BOMI and TOMI match precisely with those in Section III B.

IV. DISCUSSION
In this work, we have studied a strong version of the butterfly effect in quantum many-body systems from an information theoretic perspective. We have found that local operators in chaotic theories entirely delocalize information, regardless of the details of the operator. In certain "large-N " theories such as holographic CFTs and random unitary circuits with large local Hilbert space dimension, we have found this delocalization process to occur at the fastest possible rate allowed by causality. In contrast, we have found that integrable theories are robust against these perturbations.
There are several interesting avenues for further study of local operator entanglement. These include higherdimensional calculations which may be made possible through the holographic membrane theory [44,58]. Holographically, it should also be tractable to compute the entanglement wedge cross-section in the massive-particle geometry of Section III C. It will be important to understand if the reflected entropy remains parametrically larger than the mutual information as this is a novel phenomenon never seen for simpler quantum systems and is hinting at the fundamental role of multipartite entanglement. Moving beyond holography and maximally chaotic systems, it would be fascinating to understand this notion of the butterfly effect in more generic quantum sys-tems. In particular, it is important to understand the universal features in generic interacting RCFTs and irrational CFTs as has been previously done for OTOC, local quenches, global quenches, and unitary operator entanglement [24,30,32,35,36,53,54,[59][60][61]. Similarly, it is desirable to understand non-conformal theories that are not maximally scrambling such as spin chains and random unitary circuits with finite onsite Hilbert space dimension [17][18][19][20][62][63][64][65][66][67] and interacting integrable systems that exhibit diffusion, a tractable example being the Rule 54 chain [18,68]. of Appendix A will be published elsewhere [70]. In this appendix, we make progress in the derivation of the membrane theory for logarithmic negativity and reflected entropy conjectured in Refs. [22,35] and used in the main text 9 . We do this using the formalism developed in Refs. [38,65] for Rényi entropies. We will work in the q → ∞ limit where significant simplifications may be made. It will be interesting to understand finite-q effects in the future. The q → ∞ limit is relevant to irrational CFTs where the effective bond dimension is determined by the Cardy density of states (7), q = e πc 3β . Progress can be made because of analytic formulas known for averaging over an arbitrary numbers of Haar random unitary matrices (A1) where σ and τ are elements of the permutation group S n and n is the number of unitaries (and duals) averaged over. Once applying this independent averaging on every unitary matrix, we end up with an effective hexagonal lattice of classical S n spins. We are thus instructed to compute the partition function on this lattice. It has been shown that the partition function simplifies by summing over the τ variables and we end up with a triangular lattice with positive three-spin interactions involving the Weingarten function [66] . (A2) These interactions are still quite complicated but they simplify greatly in the q → ∞ limit where they equal q |σ −1 B •σc|−n . | · | for a permutation element denotes the total number of cycles in that permutation.
Rather than considering spin configurations in the partition function, it is more convenient to consider S n domain wall configurations. With these simplifications, we can reduce the problem of computing the partition function to finding the dominant "bulk saddle" completely analogous to the story for conformal field theories in the large-N limit [71], though, so far, the discussion has been quite general and we have not specified to the entanglement entropy. In order to apply this general framework to the specific quantities that we are interested in studying, we must apply appropriate boundary conditions.

Mutual information
For the Rényi entropies, one must impose Z n permutations on the boundaries within the regions of interest while imposing identity elements everywhere else (see Fig. 16). This imposes the correct trace structure. The Z n permutation is represented as where we are using cycle notation to label the elements of S n . To leading order in q, the replica partition function is Here, γ con and γ dis represent the areas of the extremal surfaces in the circuit of different topologies. In the large q limit, when γ con > γ dis , we need to maximize |τ |. This is achieved when τ is the identity element, e, because |e| = n. Analogously, when γ dis > γ con , we need to maximize |σ −1 n • τ | which occurs when τ = σ n . Thus, to leading order, the Rényi entropies are all equal S n (A) = min[γ con , γ dis ] log q.
(A5) This may be described by a membrane theory because it is equivalent to finding the minimal membrane in the circuit with membrane tension log q.

Negativity
For the logarithmic negativity, we must compute even powers of the partial transposed density matrix For this partition function, we must modify the boundary conditions of the S ne spin model to incorporate the simultaneous cyclic and anticyclic gluing of the replica manifold. In cycle notation, the anti-cyclic permutation is We will now perform a 0 th order analysis of this double sum and show various subtleties that arise that are not present in the case of Rényi entropies. The simplest regime is when A and B are sufficiently distant. In this case, we must maximize both |τ 1 | and |τ 2 | which means that we take both of them to be the identity. This contributes to the partition function as In the opposite regime where we expect correlations to be present between A and B, ignoring the term with E W , we must maximize both |σ −1 n • τ 1 | and |σ n • τ 2 | which means that we take τ 1 = σ n and τ 2 = σ −1 n . This contributes to the partition function as Z (P T ) n Z n 1 ⊃ q E W (2−n)+(1−n)(γ A,con +γ B,con ) .
While these two choices for permutations seem most natural, due to the term proportional to E W , there exists another permutation (with degeneracy) that can become important. These permutation elements can be thought of as being halfway between σ n and σ −1 n , leading to the following contribution Z (P T ) n Z n 1 ⊃ q (1−n/2)(γ A,dis +γ B,dis )−(n/2)(γ A,con +γ B,con ) . (A11) The three contributions described above are the most clear leading contributions to the partition function. However, we have not proven that other permutations are not also important or even that it is valid to take a single dominant contribution. If we make the assumption that these are the only important saddles and that it is sufficient to throw away the rest, we find as q → ∞ E (n) ≡ log Z (P T ) n Z n 1 = − min (n − 1)γ con + (n − 2) E W , (n − 1)γ dis , n 2 γ con + n 2 − 1 γ dis ) log q, where γ con ≡ γ con,A + γ con,B and γ dis ≡ γ dis,A + γ dis,B . Here, we have approximated the logarithm of the sum of the three contributions as a "min" function due to the q → ∞ limit. In reality, we are summing the many terms of the form q # . For the disconnected regime (γ con ≥ γ dis ), the middle term will dominate for n > 1, leading to Naively taking the replica limit, this vanishes, implying that the logarithmic negativity is zero. This makes sense because it corresponds to the regime where the intervals are either sufficiently distant or we are at sufficiently late times when thermalization has occurred. The connected regime (γ dis ≥ γ con ) is more subtle. First note that the first and third terms are identical when n = 2. Being linear in n, this means that one of the terms is minimal for all n > 2 and the other is minimal for n < 2 i.e. there is a replica transition. This is a novel phenomenon that we did not see for Rényi entropies. To determine which term is dominant in which regime, take n = 1 for simplicity. In this case, we are comparing the size of E W versus (γ dis − γ con )/2. It is a simple geometric exercise to show that E W is always the greater of these two, thus the first term would appear to be dominant for n < 2 while the third term is dominant for n > 2. If we are to take the analytic continuation seriously, we would determine that the logarithmic negativity is given by This is precisely the membrane theory used in Refs. [22,35]. Moreover, it shows how the entanglement wedge crosssection in the context of negativity can emerge outside of holographic conformal field theories. However, this analytic continuation was too naive. We warn the reader that the above derivation was not rigorous both due to dropping all subleading terms in the double sum and in the analytic continuation to one. We have not proven that other terms are not important. We leave a more rigorous proof (or disproof) that involves computing the full negativity spectrum to future work [70]. One peculiarity in this result that must be resolved is the following: Recently, the logarithmic negativity has been analytically computed for arbitrary tripartitions of Haar random states [72]. Haar random states are expected to accurately describe the late time states after random local quantum circuit evolution, an intuitive idea that has been made precise in Ref. [39]. However, if one compares the results from Ref. [72] with those predicted by (A14) at late-times for a finite size system, one finds disagreements. Rather than E W , Ref. [72] found an answer that is more reminiscent of a Rényi mutual information. Resolving this tension is an important future direction 10 .

Reflected entropy
We are able to run through similar analysis for the reflected entropy by imposing yet another boundary condition on the effective spin system. The partition function is indexed by two replica numbers S (n) R = 1 n − 1 log Z n,m (Z 1,m ) n . (A15)

22
The computation of these partition functions has an associated replica trick (see Ref. [49] for details). For our purposes, we simply need to recall the cycles defined for the relevant twist operators so that we may set appropriate boundary conditions on the S mn spin model. In region B, the permutation element is σ g B = n k=1 (k, k + n, . . . , k + n(m − 1)). (A16) Each factor consists of m − 1 elementary swaps, so in total, the domain wall between this element and the identity is composed of n(m − 1) elementary domain walls. Similarly, the permutation elements acting on region A give n(m − 1) elementary domain walls, but the cycles are of a different form (k, k + n, . . . , k + n(m/2 − 1), k + 1 + nm/2, . . . , k + 1 + n(m − 1)).
Thus, the reflected entropy is zero in the disconnected regime. In the connected regime when γ con < γ dis , we have log Z n,m Z n 1,m = max [2(1 − n)E W , n(1 − m)(γ dis − γ con )] .
The second term is undesirable so we wish to take n → 1 before m → 1. This was previously noted in Refs. [60,61] as necessary for picking the correct "entanglement wedge." However, this order of limits subtlety is only present because we have taken q → ∞. At finite q, there should be no ambiguity. A better understanding of the necessity of this order of limits deserves further attention. Taking the proper order of limits, the reflected entropy becomes as advertised. Again, we stress that we have made serious assumptions about the dominating terms in the sums over the permutation group. A more thorough analysis of this issue is being pursued [73].