Scrambling with conservation law

In this article we discuss the impact of conservation laws, specifically $U(1)$ charge conservation and energy conservation, on scrambling dynamics, especially on the approach to the late time fully scrambled state. As a model, we consider a $d+1$ dimensional ($d\geq 2$) holographic conformal field theory with Einstein gravity dual. Using the holographic dictionary, we calculate out-of-time-order-correlators (OTOCs) that involve the conserved $U(1)$ current operator or energy-momentum tensor. We show that these OTOCs approach their late time value as a power law in time, with a universal exponent $\frac{d}{2}$. We also generalize the result to compute OTOCs between general operators which have overlap with the conserved charges.


Motivation
Quantum information scrambling [1,2,3] is a fundamental phenomenon in chaotic many-body systems that has been under wide discussion in recent years. On the theoretical side, this activity focused on the study of out-of-time-order correlators (OTOCs) [4,5,6,7,8,9,10,11,12,13,14,15,16]. For a chaotic system with large N number of degrees of freedom per unit volume, the OTOC displays an exponentially increasing deviation from its initial value which is characterized by a quantum Lyapunov exponent. This period of growth occurs after local equilibrium is achieved but before the scrambling time, when the system approaches global equilibrium. Schematically, given simple Hermitian operators W and V , one has For a large N conformal field theory (CFT) holographically described by Einstein gravity, earlier works [5] [6][17] [18] [19] calculated OTOCs through geometric methods, by studying shockwave geometries. The Lyapunov exponent in such a theory is equal to 2π β , saturating the conjectured chaos bound [20]. For times much larger than scrambling time (the late time regime), the OTOC decays to zero exponentially fast, with a different but related exponent. By contrast, it was shown in [21] that in a random circuit model with local charge conservation law, the OTOC between the charge density operator and a non-conserved operator displays a power law tail at late time. While such power law tails are expected to be generic, they have not yet been seen in holographic systems. This missing piece of physics motivated the present study. Other studies of the interplay between conservation laws, hydrodynamics, and OTOCs include [22,23,24,25,26,27,28,29,30,31,32,33,34,35].
The bulk of this work focuses on the U (1) case: we compute OTOCs between the conserved charge density and a (non-conserved) scalar operator and show that holographic systems also exhibit power law tails consistent with the random circuit result. We also argue that power law tails will be induced by energy conservation. This has not been shown before, but it is important since Hamiltonian systems generically have energy conservation but not charge conservation. In the remainder of the introduction, we review existing holographic calculations, focusing on the scattering approach. Then, in the rest of the article, we show how these calculations are modified due to U (1) charge conservation and we interpret the result physically. We also discuss the extension to energy conservation. Overall, we view this work as a study of the interplay between slow modes, as in hydrodynamics, and the fast dynamics of scrambling.

Review of the holographic calculation
The holographic calculation of OTOCs for scalar operators is discussed in many works. The approach that we follow here is based on the scattering approach discussed in [6]. In this approach, the boundary OTOC is written as an inner product of in and out asymptotic states.
W (t 1 , x 1 )V (t 2 , x 2 )W (t 1 , x 1 )V (t 2 , x 2 ) = out|in |in = W (t 1 , x 1 )V (t 2 , x 2 )|T F D |out = V † (t 2 , x 2 )W † (t 1 , x 1 )|T F D . From the bulk perspective, these in and out states may be written in terms of bulk wavefunctions as |in = dp u dx dp Here u and v are null coordinates in the black hole geometry dual to the thermofield double state |T F D . Note that the momentum states are only well-defined near the black hole horizon. We think of the scattering as taking place in the approximately flat near horizon region as described by Kruskal coordinates. This corresponds to the case of large time separation between t 1 and t 2 , as shown in Fig 1. The wave functions are given by bulk-to-boundary propagators. So These formulae have a direct interpretation as bulk scattering states sourced by boundary operators. The relevant energy scale is determined by the boundary time separation through the Mandelstam variable . For scrambling physics, we are interested in time scales that are larger than the relaxation time. The dominate contribution in this regime comes from summing ladder diagrams with graviton exchanges [36]. To leading order in s, the S-matrix approaches a pure phase, obtained from the eikonal approximation, The physical interpretation of this phase factor is interaction of particle 1 with a gravitational shockwave sourced by particle 2 [37] [38]. The shockwave metric is While passing through the shockwave located near u ∼ 0, the scalar wave function receives a jump in the v coordinate, The phase δ(s, b) is then identified with the displacement factor h( For scalar operators, the OTOC is evaluated in earlier works. In the limit ∆ W ≫ ∆ V ≫ 1, the heavy particle's wavefunction is not much affected by the shockwave sourced by the light particle. So the amplitude is simply an inner product between V particle wavefunctions with and without the shockwave, Γ is a constant depending on the regularization of the correlator. At early time, when the second term in denominator is much smaller than the first term, the OTOC is decreasing as shown in Eq. (1). Note that e 2π β t is roughly the colliding energy, and it enters through p dependence of h(x). The exponent is 2π β , independent of the operator. The correlator has decayed significantly when the second term in the denominator is O(1). After this time, the correlator experiences an expontial decay with an operator-dependent exponent. We refer to this as the late time regime of OTOC.

Results and outline
The situation is changed when we consider an OTOC that involves conserved charges. For example, the R-charge in N = 4 SYM theory, and more generally the energy-momentum tensor. We will see that due to the hydrodynamical property of the conserved current, the particle sourced by these operators in the bulk spreads over a large region of space-time. As a result, the collision responsible for scrambling happens at a wide range of space-time points in the classical picture (see Fig 2). When the collision occurs near the horizon, the center of mass energy is large, but when the collision occurs further away from the horizon, the center of mass energy is smaller. This leads to a smearing of the exponential factor in Eq. (8), effectively replacing the OTOC formula with where c and α are some constants. One can then show that at late time, the OTOC becomes ∼ t − d 2 . The article is organized as follows. For simplicity, we start with the conserved charge density of a U (1) symmetry. In section 2 we give an expression for the dual photon's wave function, to lowest order in transverse momentum and frequency. Then we analyze the inner product and equation of motion for the photon field. Using these ingredients and some additional approximations, we calculate the OTOC. In section 3 we interpret these calculations in terms of the classical picture in Fig. 2. In section 4 we discuss the case of the energy-momentum tensor. Finally, we generalize the result in section 5 to the case of generic operators with overlap with the conserved currents.

Solution to equation of motion
In an AdS-Schwarzschild black hole background, the metric is where R = R + is the horizon and boundary is at R = 0. The inverse temperature is β = 4πR+ d+1 and f (R) = 1 − ( R R+ ) d+1 . It's more convenient to use the tortoise coordinate, r, defined by O(t 1 ) The domain of r is r ∈ (−∞, 0] with −∞ corresponding to the horizon and 0 corresponding to the boundary. In this coordinate, the metric components are g rr = −g tt = L 2 f (R) R 2 . Here we use x to denote the coordinates of the transverse directions. We also consider a Maxwell field propagating in this geometry whose dynamics determines the behavior of a U (1) current operator J 0 on the boundary. In this note, we will focus on the four-point correlation function with two insertions of J 0 and two insertions of a scalar operator O.
The Maxwell-Einstein equation is Consider the ansatz A i = 0 for all the transverse directions. The Fourier mode decomposition is Using the trick in [39], we pick a transverse coordinate frame for each k, such that the x-axis is parallel to k direction, with the other axes perpendicular to it. Then ∂ x can be replaced by ik when acting on that mode, and derivatives in other transverse direction are replaced by 0. We also assume that the A field is spherically symmetric and excited by a point source sitting at x ′ = 0, so the momentum space A only depends on norm of k. Written in component form, the equations are From the second equation, one can deduce Then one can use the first equation to find a second order differential equation for A t , This equation has two independent solutions. As r → ∞, they behave like e ±iωr , corresponding to out-going and in-falling boundary condition at the horizon. We focus on one of them, since the other one is obtained by complex conjugation. The solution can be found explicitly in the small frequency and small momentum limit by setting ω → λω, k → λk, and taking λ ≪ 1. Up to first order in λ, we have The normalization constant C is fixed by requiring lim r→0 A t (r, ω, k) → 1. Then we find that The factor 1 d−1 R + is identified as the diffusion constant D. H(r) is the indeterminate integral of the second term in Eq. (18). If we ignore the ωH(0) term in the denominator, then this function has a pole at ω = −iDk 2 . With this ω ∼ k 2 scaling, the ωH(0) is indeed subleading compared to the ik 2 D/ω term at small ω and k. Moreover, we can anticipate that ω 2 is of order k 4 after using the residue theorem. Hence, to the leading order in small k and ω, we can neglect ω 2 and higher order terms. Another approximation is to set ( R R+ ) ∼ 1, corresponding to the near horizon region. Since the scrambling time is large, the relevant physics indeed happens within this region. In summary, we have The real space form is where t ′ and x ′ label the position of the boundary source. We expect this expression to hold near the horizon when the transverse separation from the source is large. Now we use Eq. (15) to obtain the other components, Since g tt → 0 near horizon, F tr is also very small there. (Note that the physical electric field E r = g tt √ g rr F tr is still finite).

Inner product of gauge field
The gauge invariant inner product between two gauge field configurations A 1 and A 2 is Where the integration is over a Cauchy surface. Following [6], we choose to integrate on constant u ∼ 0 slice. Then this inner product becomes with u and v related to t and r by The various components of gauge field are related to the old ones by Following Eq. (19) and Eq. (21), Plug this into the inner product expression Eq. (23), we get where φ(v, x) is the diffusion kernel given in Eq. (26). We will see that the inner product written in this way is very convenient when discussing the shockwave.

Scattering states
The original definition of OTOC is complicated by a UV divergence arising from coincident operator insertions. To avoid this, one can consider a regularized version obtained by inserting operators at different values of imaginary time, To simplify the discussion, we choose a symmetric regularization. In the following, we will consider the correlator It has a representation as an inner product of in and out states where |S − AdS denotes the Schwarzschild-AdS black hole thermal double state, Figure 3: in and out sates in a symmetric regularization scheme The superscript L and R on each operator means that the excitation is created on the left or right side, respectively. For example, J R 0 (t 2 ) creates a photon on the right. We denote the corresponding gauge field by A R µ . In general it is a linear superposition of in-falling and out-going solutions, depending on the state we are constructing. Following [40][41], we choose the coefficients such that the field has positive Kruskal frequency for in-falling mode and negative Kruskal frequency for out-going mode. For this purpose, we use the following combination n(ω) is the Boltzmann factor. For the in-falling part, we use the ansatz: This is just rewritten from Eq. (26), except that we have inserted an extra 1 in the denominator. This change doesn't modify the long time behavior of the wave function, but it does provide some convenience in the analysis. On the other hand, the out-going part in A v is proportional to u. To evaluate the OTOC, we choose to calculate the inner product on the surface u ∼ 0. Hence, the out-going mode's contribution can be neglected. The thermal factor (1 + n(ω)) is proportional to , evaluated at the diffusive pole 1 . It suggests the following ansatz for A nice property of the above expression is that the real space form of φ R and φ R in−f alling are related by analytical continuation, if we treat it as a complete solution of equation of motion. To see this, note that Similarly, on the left side boundary, the operator J L 0 sources a wave function A L v . By symmetry, Then A L v can be written as The scalar operator O R and O L create scalar mode in the bulk. For simplicity, we assume that the scalar operator O has a large conformal dimension, corresponding to a particle in bulk with a large mass. As a consequence, we can treat the scalar particle semi-classically. As it moves deep into the bulk, the scalar mode carries the shockwave along with it, which modifies the photon's wave function as we will analyze in the next section. In contrast, we will neglect the back-reaction from the photon on this scalar mode.

Shockwave geometry
Since we are interested in the large time limit of the OTOC, we will take t 1 ≫ β and t 2 ≪ −β. As a result, the geodesic of the scalar particle is approximated by u = ǫ ∼ 0. The shockwave geometry is described by a metric that contains a singularity near the horizon, where where the i's run from 1 to d and label the transverse directions. The effect of the shockwave is encoded in the g vv component. γ is the metric determinant in the transverse directions, so −g = g 2 uv γ. As before, we simplify the equation by approximating R R+ ∼ 1 such that we can set γ = ( L R ) 2d and In the region u < ǫ and u > ǫ, the first two terms are finite. At u = ǫ, we are searching for a solution where A v jumps and A u may contain a delta function δ(u). On the other hand, from the first equation of Eq. (38), F uv must be finite. Therefore the delta functions in ∂ u A v and ∂ v A u must cancel each other at u = ǫ. Applying this condition, Eq. (39) can be written as Integrating over u we obtain Note that this shift in v doesn't commute with derivatives in the transverse directions, so this simple shift rule only applies to ∇A instead of A itself. The first equation in Eq. (38) gives the same relation between F uv and A as in Eq. (21). Comparing with Eq. (26), we find that this simple shift rule also applies to which is the reason to write the inner product in the form of Eq. (27). In summary, after scattering with the shockwave, the wave function sourced by J R 0 is modified according to the rule Eq. (42). If we choose to evaluate the inner product, Eq. (27), just after the scattering, along u = ǫ + 0 + , then φ 1 and φ 2 (in Eq. (27)) should be

Calculation of OTOC
In this section, we evaluate the inner product in Eq. (27). It turns out that the relevant integral in is easier to do in momentum space. Define so that Note that if we only keep the leading order terms in k 2 , then the term sin( βDk 2 2 ) will cancel with Γ( β 2π Dk 2 ). As long as we are considering large transverse coordinate separation, this approximation should be qualitatively correct. Plugging these into Eq. (27) and approximating R R+ ∼ 1, we obtain To save space, we have suppressed the factor β 2π setting the units of time. After some changes of variable and approximations shown in Appendix B, we obtain the final expression for the OTOC as At early time, this expression admits a large N expansion in which the leading term still grows exponentially with time. Moreover, one can identify the same Lyapunov exponent and butterfly velocity as in the non-conserved OTOC in Eq. (8). In the late time limit, the ln of the exponentially growing part in the denominator gives rise to a power law time decay behavior. Hence, we find a significant difference from the OTOC of non-conserved operators in the late time regime.

Physical interpretation
In this section we try to understand the result in Eq. (47) in a more intuitive way. In the last section, we calculated the OTOC by integrating first over the radial momentum. Although this makes the calculation easier, the physical reason why we expect a power law decay at late time is somewhat obscured. As an alternative approach, we can perform the spatial momentum integral at the beginning and rewrite the integral in Eq. (46) as Here s is defined as s = ln( e −t 2 p ) and we have cut-off the large momentum contribution for p > e −t2 . This formula has a direct physical interpretation. The squared term can be viewed as the photon's wave function. From the solution Eq. (20), we see that the photon's wave function is extended both in the radial direction and in the transverse directions. Although we obtain Eq. (48) in momentum space, it's more inspiring to think of s as the radial coordinate. Small values of s correspond to regions close to the horizon, while larger values of s correspond to regions further from the horizon. The last term is a phase induced by the gravitational scattering, and e t12−s measures the relative scattering energy of the two particles. Therefore, it is natural to think of the scattering between the photon's large wave-packet and the scalar particle's localized wave-packet as taking place over a large range of radial depths. Then as the scalar particle scans through the photon's wave function, the colliding energy effectively becomes smaller and smaller.
The term in the middle can be thought of as a regulator for large momentum (small s). Since we didn't obtain the complete solution of the photon's wave function in the black hole geometry, the regulator might be replaced by a more general function in the full solution. However, if we're only interested in the late time behavior, meaning t 12 − µ| x 12 | ≫ ln(N ), then the integral in s receives its dominate contribution from s > t 12 due to the fast oscillation of the phase term for small s. The regulator is therefore not important, and we can directly see that the result is proportional to 1 t12 d 2 .

Correlator with stress-energy tensor
The low energy hydrodynamics of the stress-energy tensor shares some similarity with that of U (1) charge. In this section, we show that the OTOC has the same late time behavior. As shown in [42][43] [44], in hydrodynamic limit, small perturbation of stress-energy tensor splits into sound and shear modes. The corresponding graviton wave functions in the bulk have a sound pole and a diffusive pole, respectively. The shear mode is relatively easier to analyze, as it involves less components of the metric perturbation, but we will see that they have the same qualitative effect on the OTOC in the late time regime.
As above, we take the k direction as the x-axis. The shear modes consist of T ty , T xy . Here y can be any direction perpendicular to x, Sound modes are more involved, including T tt , T tx , T xx and T yy . The shear and sound mode operators on the boundary excite bulk metric fluctuation corresponding to vector and scalar perturbations, respectively. Taking the operator T ty as an example, it sources the bulk metric perturbations h ty and h yr . Other choices of non-zero metric components are related to this by a gauge transformation. The spherically symmetric choice of scalar perturbation sourced by T tt involves h tt , h tr , h xx = h yy , and h rr . For simplicity, we mainly discuss the shear mode and comment on the sound mode in the end.

Wave function of graviton
For a given momentum k, choose the coordinate system such that the x-axis is parallel to k. The shear modes involve the T ty (ω, k) and T xy (ω, k) components, and are described by which together imply In AdS/CFT, the dynamics of these modes can be found by solving the linearized Einstein equation, where d is the spacial dimension of boundary theory. h µν is the metric perturbation, δg µν = h µν . The equations are simplified if we set to zero all the components except h ty and h ry . Then there are only two independent equations, where R is the radial coordinate in the metric Eq. (10) and r is the tortoise coordinate defined in Eq. (11). The two first order equations then lead to a second order differential equation for h y t , This equation is the same as equation in U (1) charge case except with d−2 replaced by d (see Eq. (115)). Thus we have the same solution, except the diffusion constant D T = 1 d+1 R + is different. Also, while the photon wave function excited by J 0 is spherically symmetric in the boundary spatial plane, in this case, since the shear mode operator T ty contains a spatial index, it breaks the spherical symmetry. So we expect that the wave function given by Eq. (20) only captures the dependence on directions x satisfying x⊥y. In the y direction, the mode propagates as sound. To avoid the complexity of mixing sound and shear modes, we restrict to consider the string operator living in d − 1 spatial dimensions, T s (t, x) := +∞ −∞ dyT ty (t, y, x). Note that the choice of line operators over point operators can change late-time exponents by shifting the effective dimension of space.

Interaction with the shockwave
We continue to focus on the OTOC between energy-momentum tensor and a scalar operator with large conformal dimension. As above, we approximate the shockwave as sourced by the heavy scalar without backreaction from the graviton. Then it remains to consider the evolution of the graviton wave function in the geometry. In the following, to distinguish the metric perturbation from the shockwave, we will use f ( x) as the displacement in the metric, We solve the linearized equation This is a very complicated equation in general. Some other components have to be generated even if we start with only the shear mode perturbation. However, if we restrict the shockwave term g uu to only depend on directions that are perpendicular to y, the equation becomes easier to deal with. So we require S to be a function of x i with x i ⊥y. Then the equations are simplified to where 'finite terms' denotes terms that don't contain a delta function. Here we are searching for a solution such that ∂ u h y v and h y u are proportional to δ(u). Requiring the cancellation of all δ(u)s, we find that Hence, we again have the simple shift rule, after the graviton passing the shockwave. Finally, in order for f to only depend on x⊥y, we have to also consider a string operator built from the scalar O s (t, x) = dyO(t, x, y).

Inner product
As before, we construct the gauge invariant inner product using the symplectic form, We find that Knowing that the trace of h 1 and h 2 are zero, we only keep the first two terms. Then, after some cancellations, a much simpler form remains We have chosen n to be ∂ ∂v . This expression is very similar to the photon case (see Eq. (22)), except for second term, which is proportional to the connection. Then we observe that Γ y yv is of order u near the horizon, so if we choose to evaluate the inner product along the hyper-surface u ∼ 0, the second term can be neglected.

Shear mode
Just like the U (1) case, the OTOC, can be written as an inner product between the graviton wave functions, before and after passing the shockwave. Plugging the solution Eq. (20) and Eq. (58) into Eq. (61), we get exactly the same expression as in the U (1) case, and the subsequent calculation is completely parallel. Note that these string operators live in an effective d − 1 spacial dimension. The exponent of power law tail is also modified to d−1 2 .

Sound mode
We can also consider a local insertion of a stress-tensor operator that creates spherically symmetric wave function propagating in the bulk. At low energy, these modes have the dispersion relation of a sound wave. For instance, we can choose the insertion T tt and d i=1 T ii . The pole in the Green's function is at where D T is the diffusion constant (same as that of the shear mode). The quadratic term has the same effect as the shear mode, broadening the wave-packet. Hence, we expect the OTOC to have the same power law tail as in the photon case. Another interesting regime for the sound mode is at early time, where we have the shockwave that propagates at the butterfly velocity v B as well as the hydrodynamical mode that propagates at the sound speed v s . For the gravity model we considered in this work, v B > v s for physically sensible spatial dimension. However, in other models the sound speed might be larger. 2 In this case we find that the information carried by sound mode can scramble faster, the spreading speed of which is determined by v s instead. We solve the sound mode gravitaional perturbation equation and provide a detailed analysis of the OTOC in Appendix C.

Higher order corrections
We have calculated the photon wave function by just keeping the leading order in ω and k. This gave an OTOC with a power law dependence on t at large time. We may include higher order terms in the wave function, for example, consider The pole is now at ω = −iDk 2 (1 + γk).
After Fourier transformation with respect to ω, the wave function becomes φ(v, k) = dωφ(ω, k)e −iωv , (66) Following the steps above, we should multiply the integrand of Eq. (46) by a factor (1 + αk + βk 2 ). Also we should modify the exponent of [2 + f ( x)e t12 ] in second equation of Eq. (78) to −Dk 2 (1 + γk). After performing integration with respect to k, these modifications only contribute factors of higher order in t −1 12 , and do not change the leading long time behavior. One can also include loop corrections due to a graviton, as we discuss in Appendix D. This correction creates a branch cut in the wave function but doesn't modify the long time behavior of the OTOC.

Hydrodynamics and OTOCs of non-conserved operator
In this section, we explore the question of how hydrodynamic modes may affect scrambling of non-conserved operators. In the process shown in Fig 4 (with the photon line representing either a photon mode or a sound mode), imagine that a massive scalar particle (created by a non-conserved operator) emits a hydrodynamic mode through a coupling of order O(g c ) in the near boundary region. This mode grows in size as it falls into the black hole. The in state before the shockwave scattering would contain a term where f s is the sound mode wave function and f m is the massive scalar wave function. The state with momentum q is the second massive mode. The collidision induces a phase factor e iqph(x,y)+iqp ′ h(x ′ ,y) . So we expect the scattering amplitude to be the one that involves hydrodynamic modes (like in Eq. (29)) multiplied with the one that involves only massive scalars. As discussed in the previous sections, the hydrodynamic OTOC has a power law tail at late time. However, the non-conserved OTOC is multiplied with it and the combination decays faster. Hence, the complete OTOC at late time is still controlled by massive particle scattering. On the other hand, the early time behavior can be modified, if we consider the non-conserved mode coupling with a sound mode, becuase the scattering amplitude between the sound mode and the second massive particle starts to decay earlier if v s > v B (see Appendix C for details). Since the hydrodynamic mode's wave function is created with amplitude O(g c ), its influence on the OTOC is of order O(g 2 c ). So even if v s > v B , the fast propagating wave-front cannot grow to exceed the same order, Since all modes couple to gravity, they also couple to the sound mode of the gravitational perturbation. In this case, g 2 c ∼ 1 N , and the last term can grow with time up to order O( 1 N ). A similar situation has been discussed in [22], where they gave a bound on the squared commutator of the form where the function a(x, t) is non-zero when |x| < v s t, and is bounded by an O(1) quantity. However, the picture we considered here (Fig 4) is slightly different.

Summary
In conclusion, we have explored OTOCs between hydrodynamic operators and generic operators. In the late time regime, these OTOCs obey a power law scaling, while at early time the deviation still grows exponentially. In some models where the sound speed v s is large, the information near the wave-front scrambles with a velocity that depends on v s . Finally, when generalizing to OTOCs of generic operators that couple to hydrodynamic modes, we found that the late time power law decay is absent, while there can be a small amount of information scrambling faster than v B (when v s > v B ) near the wave-front. We can also understand the conclusion intuitively in the boundary picture. Starting with a nonconserved operator O. A part of it evolves into g c JO, where J represents a hydrodynamic operator. Then a small part of the information is carried by the hydrodynamic mode. A pure hydrodynamic operator J may spread fast (in sound mode case) but release its information slowly (which causes the late time power law tail). Therefore, although the order g 2 c amount of information may propagate fast and lead to a rapidly moving wave-front in the OTOC, the scalar operator accompanied with J releases most of the information and breaks the power law tail at late time. Due to conservation law constraints, the dynamics prevents O from turning into pure J's. In bulk theory, this is a constraint from gauge symmetry.

Acknowledgements
This work is supported in part by the Simons Foundation through the It From Qubit Collaboration.

A Solution to the Maxwell-Einstein differential equation
In our metric, −g tt = g rr = R 2 f and √ −g = R −(d+2) f . The Maxwell equation simplifies to There is a singular point at f → 0, where the equation becomes with solution A t (r, ω, k) ∼ e −iωr . We have picked the in-falling mode on this boundary. The ansatz yields an equation for F , where (ω, k) = (λω, λk) and λ ≪ 1. The solution F should have an expansion in the form F = F 0 + λF 1 + · · · . The leading F 0 satisfies F 0 goes like C 1ω 2 r as r → −∞, so a regular solution should have C 1 = 0. Similarly, the first order term F 1 should satisfy The integration constant can be fixed by requiring regularity at horizon, then one obtains To first order in λ, we find

B Details of the OTOC calculation
Start from the inner product, As in the main text, we will just keep the leading order dependence in k for the gamma function. We change the integration variables to K = k+k ′ 2 and κ = k−k ′ 2 . This gives where g and E are defined as Since both of the two terms contain e − | x| 2 2D ln g , we expect the integral to receive its dominant contribution from |x| ∼ 0. Integrating x over this saddle point gives a factor of (2D ln g) d 2 . Therefore, we finally obtain In fact, the saddle point approximation in last step is only valid for very small diffusion constant D and for not too large values of the function g(x). Thus it is necessary to discuss the late and early time limits separately from the above treatment. Looking back at Eq. (78), in the large t 12 limit, t 12 ≫ LnN + | x12| vB . We can expand the function ln g ∼ t 12 (1 − | x| vB t12 ), for |x| < cv B t, where c is some finite constant smaller than 1. Then the integral over x can be evaluated as 2D ln g , For the first part, the change of variable to x → The second part is bounded by From these results, the late time amplitude indeed decays in a power law manner, and is consistent with the result Eq. (80).
On the other hand, when t 12 − | x21| vB ≪ ln N , we expand ln g ∼ ln 2 , as well as The leading order deviation is . We can see that the velocity of the wave-front is still given by v B . For completeness, we give the exact result of this integral in 1d, (87)

C OTOCs of sound mode operators
In this section, we evaluate the OTOC between a sound mode operator (eg. T tt , T xx + T yy ) and a scalar operator with large conformal dimension ∆. For this purpose, we first solve the sound mode equation in a symmetric gauge to second order of ω and k. By doing this, we obtain the dispersion relation as in [46] and the near horizon form of the wave function. Unlike in [46], we apply a gauge fixing condition that preserves boundary spatial rotational symmetry. If the boundary source respects rotational symmetry (for example, consider a boundary insertion of T tt or d i=1 T ii ), we expect a bulk configuration satisfying h ti = h ri = 0 and h ij proportional to identity matrix, for i, j ∈ {1, · · · , d}. Without loss of generality, we will take d = 2 in the following. Then the non-zero components are h tt ,h tr , h rr and h xx = h yy . They satisfy a set of differential equations. Using the tortoise coordinate r = − dR f (R) , and the Fourier decomposition h MN (r, t, x) = dωdkh MN (r, ω, k)e −iωt+ikx , these equations are written as − These seven equations are not independent. They reduce to three independent first order differential equations together with an algebraic equation, Eq. (90). Then the last three equations, Eq. (92)-(94), give an algebraic constraint, (95) There are two independent solutions, with out-going and in-falling conditions near the horizon. To see this, make the substitution lim r→−∞ h µν (r) = e νr F µν to the equations and take the limit R → 1. The equations become There are three eigenvalues, ( 3 2 , iω, −iω), corresponding to a spurious solution, out-going, and in-falling solutions, respectively. The eigenvalue 3 2 is discarded, since the corresponding eigenvector is not compatible with the constraint Eq. (95). The eigenvector of the in-falling solution tells us that One can check that the gauge condition fixes the gauge completely. As a result, the in-falling solution is unique. To second order in momentum, we find the solution for h xx The integration constants c 1 x and c 2 x are fixed by requiring that the r → −∞ limit of h xx matches with Eq. (97).
The solution for h tt is where h (2) tt is a complicated function. Again, the integration constants can be chosen according to Eq. (97). Finally, solution for h tr is Close to the boundary (r → 0), h tt ∼ O( 1 r ), and h xx has the asymptotic form Therefore, for the boundary souce h 0 xx (ω, k) = h 0 yy (ω, k) = 1, we obtain the overall constant C(ω, k) = −3iω (1+ iω 3 Ln (3))(2ω 2 −k 2 )+ 4iω 3 (k 2 −ω 2 ) . The near horizon limit of all the components is also clear. Due to rotational symmetry, the x-axis can be any direction. So we have the following ansatz for the wave function, with the correct sound pole, Using this, we can estimate the late time tail of the sound mode OTOC, given that the inner product contains a term where is a solution to the linearized Einstein equation on the shockwave background. The calculation is slightly different from the previous one due to the linear in k term in the sound pole dispersion. In the end, it becomes the following integral where we have introduced an unknown cutoff factor G(p) to regulate the large momentum behavior. Similar to the method in Section 3, one defines s = − ln(p) and performs the integral over the transverse momentum first. Then the inner product can be written as In the final, we performed the integral over | x| assuming that it receives its dominant contribution around | x| ∼ v s s. The factor G(s) regulates small s region. When t 12 is large, the integral in s becomes significant only when s > t 12 . So we conclude that the late time OTOC has a power law tail. It's also interesting to discuss the early time regime, where we can expand the integrand in 1 N , and obtain We need to analyze the integral over s. Denote the angle betweenn and x 12 by θ. The exponential reaches its maximum at s * = x12 vs (cos θ −sin θ tan ϕ), with ϕ = arcsin vB vs . So for v B > v s , or cos θ < vB vs , then s * = 0, and we conclude that the information still scrambles at the speed v B . On the other hand, for cos θ > vB vs , we expand s = s * + l and approximate the exponential factor as vs (cos θ+sin θ cot ϕ) .
From this expression, we learn that when v s > v B , a portion of the information, corresponding to solid angle θ 0 , can scramble at a speed vs cos θ0+sin θ0 cot ϕ > v B , with cos θ 0 < sin ϕ = vB vs . For θ 0 sufficiently small, the information spreading speed is close to speed of sound.

D Non-linear corrections
In this section, we consider the leading non-linear correction to the photon's propogator. The correction is from graviton dressing, as shown in Fig. 5. Expanding the Einstein-Maxwell action, one obtains a gauge-gauge-graviton vertex of the following form, Following [47], we consider the gravitational dressing in the near horizon region, and only include the interaction between the photon and the diffusive mode of the graviton. In other words, we only keep h tx and h rx non-zero. Many terms in Eq. (111) are suppressed by small ω and k, so the leading order term involving A t is 1 t t F xr h tx Figure 5: gravitational dressing to photon's wave function To be consistent with the main text, we still use the gauge A x = 0. To find the bulk correlators, we have to solve the Maxwell-Einstein equation, Eq. (12), with a source term such that where J t and J r are 4-vectors that satisfy To obtain G tt , we require J r = 0, and J t (r, ω, k) = δ(r − r ′ ). Then J x (r, ω, k) = ω k δ(r − r ′ ). One can thenderive the following equation, Applying the in-falling condition at the horizon and a Neumann boundary condition on the boundary, we obtain an approximate solution by expanding in small ω and k. Moreover, if we focus on the near horizon region, the solution takes a simple form, G tt (ω, k, r, r ′ ) ∼ ω 2 k 2 iω − Dk 2 (1 + O(ωr)), G tx,tx (ω, k, r, r ′ ) ∼ Similarly, to solve for G rr , we use a source such that J t = 0, and J r (r, ω, k) = δ(r − r ′ ). So we have J x (r, ω, k) = 1 −ik δ(r − r ′ ), as well as the equation whereÃ r contains a contact term and is defined as 1 R d−2 A r + 1 k 2 δ(r − r ′ ). In the near horizon region, one can also find that G rr (r, r ′ , ω.k) ∼ ω 2 k 2 iω − Dk 2 (1 + O(ωr)). (118) To simplify the calculation, we replace ω by its value at the pole and set ω 2 k 2 to Dk 2 . Then we evaluate the bubble diagram by summing over imaginary frequencies and continuing to real frequency in the end, where D m = DDT D+DT is a mixed diffusion constant. The loop correction creates a branch cut along −iω < −D m k 2 , as expected by considering the on-shell condition of the two internal states. Then we suppose that the wave function in the near horizon region also receives a correction in the denominator, of the form The branch cut structure is similar to that considered in [48]. Close to the line −iω < −D m k 2 , the pole is split into two poles, at ω = −iDk 2 ± α N |k| d+4 . So the modification to real time functions is sub-leading in 1/|t|. When going to the real time, we need to add a integral along a contour that circulates the branch cut. This integral is proportional to The integrand is like a Lorentzian function, so we can estimate it by multiplying its peak height and width, which is again e −Dk 2 t (1 + O(k 2 )). So we conclude that the loop correction doesn't change the conclusion in the main text.