Quantum butterfly effect in weakly interacting diffusive metals

We study scrambling, an avatar of chaos, in a weakly interacting metal in the presence of random potential disorder. It is well known that charge and heat spread via diffusion in such an interacting disordered metal. In contrast, we show within perturbation theory that chaos spreads in a ballistic fashion. The squared anticommutator of the electron field operators inherits a light-cone like growth, arising from an interplay of a growth (Lyapunov) exponent that scales as the inelastic electron scattering rate and a diffusive piece due to the presence of disorder. In two spatial dimensions, the Lyapunov exponent is universally related at weak coupling to the sheet resistivity. We are able to define an effective temperature-dependent butterfly velocity, a speed limit for the propagation of quantum information, that is much slower than microscopic velocities such as the Fermi velocity and that is qualitatively similar to that of a quantum critical system with a dynamical critical exponent $z>1$.


I. INTRODUCTION
Elucidating the physics of thermalization in isolated quantum systems [1][2][3][4] represents an ongoing challenge in quantum many-body physics, and great progress has been made in recent years due to advances in both theory and experiments [5][6][7][8][9][10][11][12]. In this work we are interested in the process of thermalization in interacting disordered metals, specifically in the physics of quantum information scrambling. Starting from a local perturbation, scrambling describes the spreading of quantum entanglement and information across all of the degrees of freedom in a system [13][14][15][16], leading to a loss of memory of the initial state. The onset of scrambling is associated with the growth of chaos and is an intermediate step in the eventual global thermalization at late times of an isolated quantum many-body system.
It has become clear recently that certain special correlation functions can probe the onset of scrambling [17,18]. While such correlators first appeared in the literature many decades ago [19], there has been a revival in their interest, partly due to their relevance in studying information scrambling in black holes [17,18,20]. For two local operators X and Y in a system described by a Hamiltonian H, these correlation functions are defined as where ρ ∝ e −H/T is the density matrix of an equilibrium state at temperature T and X(t) = e iHt Xe −iHt . The intuition for considering this object is that local operators must grow in time if information is to spread across a system and the commutator measures this growth. Furthermore, in order to access generic matrix elements of the commutator, one considers the average of the square of the commutator, f (t), which is non-negative and avoids phase cancellations. In contrast, the average of the commutator is a response function, and these tend to decay to zero at late times in a chaotic system.
A few comments about f (t) may be helpful. When expanding out f (t) in terms of 4-point functions one finds that it contains both time-ordered and out-of-timeorder (OTO) pieces. When dealing with fermionic operators, it is more convenient to study instead the squared anti-commutator. For non-interacting fermions the anticommutator is proportional to the single particle propagator and encodes causality. More generally, one can relate commutators of composite bosonic operators, e.g. fermion bilinears, to the basic fermion anti-commutator. For a field-theory defined in the continuum, we use the 'regulated' version of the correlator above, where two of the operators have been moved halfway along the thermal circle to deal with spurious divergences. In a chaotic system with a local Hamiltonian, one expects f (t) to start out small when X and Y are spatially separated, and to grow exponentially in time, f (t) ∼ e λ L t , where is a small parameter that may depend on time and the distance between X and Y . By considering an appropriate analytic continuation of f (t), one can show that there is a fundamental upper bound on λ L (≤ 2πk B T / ) [21]; black-holes and certain random fermion models [18,22] saturate the bound. On the other hand in glassy systems, or in systems that simply fail to thermalize (but are not fully integrable), f (t) may have a power-law form [23]. While a measurement of such correlation functions is highly non-trivial, naively requiring a 'time-machine' in the laboratory, a few novel protocols have been proposed [24][25][26] and three preliminary experiments [27][28][29] have already been carried out within the last couple of months.
In this paper, we study scrambling in (weakly) interacting diffusive metals [30]. We consider the case of Coulomb interactions as well as short-range interactions in two and three spatial dimensions. Based on the gen-eral intuition that disorder slows the spread of charge and heat, one might also expect that operators spread more slowly in space in a disordered metal relative to a clean metal. Relatedly, we expect the effects of interactions to be enhanced relative to the clean metal since diffusive electrons move slowly compared to ballistic electrons and the effects of interactions can build up. We compute the growth exponent to lowest order in the strength of the interaction while carrying out an infinite resummation over disorder and indeed find that λ L is larger at low T than the corresponding result for a clean Fermi liquid. We also find that chaos grows in a ballistic fashion, with a velocity that is parametrically smaller than the Fermi-velocity at low temperatures.
These computations confirm a recent argument by two of us [23] that even though the transport of charge and energy is diffusive in such metals, generic operators grow ballistically (see also Refs. 31-33 for a related observation in one-dimensional systems). This is not too surprising, since there is no reason for the motion of charge and energy to be tied to the growth of chaos in interacting systems; an extreme example being that of a many-body localized (MBL) phase [34,35] where partial scrambling occurs even in the absence of any transport of charge and heat [23,[36][37][38][39][40]. When considering long range Coulomb interactions, it is particularly interesting that we find ballistic growth of operators since the microscopic model does not have a Lieb-Robinson bound [41]. There are variants of the Lieb-Robinson bound for systems with power law interactions [42,43], but these bounds allow exponential growth of operators with time while we find only linear growth.
On general scaling grounds, the butterfly velocity can be estimated to be v B ∼ √ Dγ in , where D ∼ l 2 τ −1 is the diffusion constant (l ≡ mean-free path, τ −1 ≈ elastic scattering rate) and γ in is a small interaction-induced inelastic scattering rate [23]. In the presence of weak interactions, we expect the exponential growth reflects the onset of chaos in an ergodic system as discussed above while the latter contribution is a result of diffusion. Solving for f (t, R(t)) ∼ 1, where R(t) is a typical 'operator-radius'-which, given an initial perturbation, defines the region in space over which information has spread over time t-leads to R 2 ∼ 4Dλ L t 2 ( Figure 1). One therefore obtains a light-cone like growth of f with a butterfly velocity We show in this paper, by carrying out a perturbative 'ladder' computation [44], that the disordered metal does obey Eq. (1.2), and the growth exponent, λ L is indeed mostly given by the inelastic scattering rate with a singular temperature dependence. Note that the unitarity of quantum mechanics prevents f (t, x) from growing to values 1 and thus it saturates at very long times. We make an analogy to the spread of an epidemic [31,44]: An 'infected' electron inserted into the center of the figure at t = 0 diffuses outwards (fuzzy red circle). As it encounters other diffusing electrons, it infects them. These newly infected electrons further infect other electrons and so on (fuzzy green circles). The flight paths of the butterflies track the spread of the infection. The radius of the region containing infected electrons (bounded by the dashed red circle) grows ballistically as vBt. Although not shown in the figure, the electrons also have a finite lifespan, given by the inverse of the quasiparticle decay rate. This needs to be taken into account when considering the population of infected electrons as a function of time. The function f (t, x) is roughly equivalent to the local fraction of infected electrons at a point x. (b) The behavior of f (t, x) for one operator placed at the center of the figure (red dot) and the other at a position x shown as a function of x at a given time t. f (t, x) displays a light-cone (a time slice of which is bounded by the dashed red circle; this region exclusively contains infected particles) within which it has saturated and no longer grows. The radius of this region grows as vBt.
Eq. (1.2) and the ladder computation are valid only for the pre-saturation growth of f . The rest of this paper is organized as follows: in Section II, we define our model of interacting electrons in the presence of static disorder and set up the basic elements required for carrying out perturbation theory to leading order in the coupling strength. Section III deals with the perturbative computation of the important terms contributing to λ L and v B for the case of Coulomb interactions in three spatial dimensions. In Section IV, we consider some additional effects in perturbation theory, as well as the case of short-range interactions, and show that our main results are unchanged by these modifications. Finally, in Section V, we study the two-dimensional version of the problem, and point out a subtle difference between λ L and the inelastic scattering rate. Unless explicitly mentioned, = k B = 1 in the rest of this paper.

II. PRELIMINARIES
We consider a model of N species of electrons in d ≥ 2 spatial dimensions subject to random potential disorder and weak interactions. We do not take any kind of large-N limit; N is a finite number (N = 2 for the case of spinful electrons). For most of this work, we shall focus on the physically relevant case of long-range Coulomb interactions in a metal; we also analyze the case of shortrange interactions in Section IV. From now on, we focus on the three-dimensional problem with d = 3 unless otherwise stated, but will analyze the case of two spatial dimensions with d = 2 in Section V.
The Hamiltonian of interest is, ) represent fermionic creation (annihilation) operators satisfying the usual anticommutation algebra, µ is the chemical potential and m is the effective mass of the electrons. The disorder potential U (x) breaks translational invariance and we assume where ... denotes averaging over disorder realizations and U 0 denotes the strength of disorder. We shall treat the interaction, V b (|x − x |), perturbatively, but will allow for strong disorder via the resummation of various classes of Feynman-diagrams with disorder lines. For Coulomb interactions in any number of dimensions V b (|x − x |) = e 2 /|x − x |, where e 2 will be the small parameter in our perturbative treatment.
Let us now review the key features of the above theory before setting up the computation for the correlation functions describing chaos in Section III. The remainder of this section closely follows the discussion in standard references (see e.g. Ref. 30).
The bare electron imaginary time Green's function after including the impurity self-energy (Figure 2a) is is the elastic electron scattering rate due to disorder (g(0) is the density of states at the Fermi level; we use the convention g(0) = 2π d 3 p (2π) 3 δ( p 2 2m − µ)). The real time Green's functions are defined as (ψ(0) ≡ ψ(0, 0)) As is well known in the theory of non-interacting disordered metals, the disorder averaged product of Green's functions in the particle-hole polarization bubble (density-density correlator) gives rise to the 'diffuson' mode at low frequencies and momenta (|ω|, v F q τ −1 ), . The noninteracting compressibility is dn/dµ = N g(0)/(2π). In the presence of interactions, the above diffuson mode introduces large vertex corrections ( Figure 2b) to the electron-interaction vertices where ω m = 2πmT and n = π(2n + 1)T are Matsubara frequencies at a temperature T and effectively screens ( Figure 2c) the long-range Coulomb interaction to, . (2.7) In the above expression, K 2 = 4πe 2 dn/dµ is proportional to the charge compressibility. Note that despite the factor of e 2 , we still treat K 2 as an O(1) quantity while doing perturbation theory in e 2 , since dn/dµ ∝ k F in d = 3 (k F ≡ Fermi-momentum) is large. Equivalently, this amounts to doing perturbation theory in Let us also review the computation of the disorderaveraged electron lifetime, which provides the inelastic scattering rate [30]. The process in Figure 2d, which includes the effects of the dynamically screened interaction and the vertex corrections, gives, via Fermi's Golden Rule, the following expression for the out-relaxation rate or the 'inelastic scattering rate' γ in ( ) for particles with energy of a given flavor i [45][46][47]: (2.8) Here n F (...) is the Fermi-Dirac distribution function.
Here, the incoming particles are on-shell while the outgoing particles are allowed to be off-shell due to the dynamical interaction V R (Ω, q). At the Fermi level ( = 0), this simplifies to where we made the reasonable assumptions q K and Ω ∼ T DK 2 . We also used the non-interacting result N g(0) ≈ 2πdn/dµ, as corrections due to interactions will only correct γ in (0) at higher orders in e 2 . At a finite energy away from the Fermi level,

III. MANY-BODY QUANTUM CHAOS
To study the onset of quantum chaos for the model introduced in Eq. (2.1), we compute the flavor-averaged squared anticommutator of electron field operators perturbatively to leading non-trivial order in the coupling e 2 , The prefactor of 1/N is inserted so that the bare contribution to f (t, x) is free of factors of N . The splitting of e −βH into two factors of e −βH/2 ensures that all operator insertions occur at distinct complex time points, thus avoiding short-distance divergences. The strict positivity of f (t, x) also guarantees exponential growth at a rate equal to that of the correlator where e −βH is not split [44,48]. These "regularized" correlators have also been shown to obey fluctuation-dissipation-like relations [49]. Computing f (t, x) involves defining the action on a complex-time contour with real time folds separated by iβ/2 [44,48,50,51]. We must then solve a Bethe-Salpeter equation arising from the resummation of different classes of ladder diagrams to determine f (ω, q), which after a Fourier transform yields information about the spatial and temporal structure of growth of chaos. An outline of the derivation of the Feynman rules for Eq. (3.1) required to set up the following diagrammatic calculation is presented in Appendix A.
Let us first quote the results for the non-interacting case, where V = 0. Here we do not expect chaotic growth of entanglement because the many-body state can be written as a Slater determinant of exact eigenstates of the one-body Hamiltonian. Summing the simplest class of where f 0 is a rapidly decaying function of time with a rate set by τ −1 and is dominated by the second term, which grows diffusively but then decays as a power law at long times. The diffusive behavior is expected as we have merely computed the particle-hole polarization bubble in real-time. As expected, there is no exponential growth.
We note here an important point, namely that we are actually computing f (t, x) averaged over different realizations of disorder. In a disordered metal for which the localization length of the eigenstates is far larger than the typical length scale over which the disordered potential varies, the disorder self-averages, and it hence makes sense to consider the disorder average of f (t, x) within a single copy of a system.
We now consider the effects of interactions, using a diagrammatic formalism which sums all the singular terms associated with diffuson and 'Cooperon' modes perturbatively in the intereaction strength [52][53][54][55]. Our perturbative computation sums all singular disorder corrections while working at O(e 2 ) in the interaction, and is formally identical to the theory of Altshuler and Aronov [47]. We will examine two effects: (i) dissipative 'self-energy' corrections ( Figure 4) that lead to decay, and, (ii) 'ladder' corrections ( Figure 5) that lead to an exponential growth of the squared anticommutator [44]. In order to obtain a non-trivial chaotic growth, the effect of the latter has to overwhelm the former. Let us discuss them now one by one.

A. Self-energy corrections
The non-interacting L(ω, q) (Figure 3a) is given by .
FIG. 4: The dominant Fock-type self-energy corrections to L(ω, q), as described in Refs. [52,53]. Each diagram has a partner diagram generated by reflecting about the horizontal axis. Also not shown are the Hartree-type contributions, which are suppressed for sufficiently long-range interactions.
The dissipative self-energy corrections to the above quantity were considered by Castellani et. al. [52,53]. These renormalize L(ω, q) at small ω, q to .

(3.4)
For T = 0, Σ L (0, 0) = 0. The field renormalization Z and the renormalization of D →D are not of particular concern to us as they will provide a correction to the growth exponent at O(e 4 ); from now on we take Z = 1 andD = D. The finite-temperature lifetime is important, and corrects the growth exponent downwards.
To compute Σ R L (0, 0) for the correlator spread across the two time folds, we note the Fock-type diagrams in Figure 4 (and their partners obtained by reflection about the horizontal axis). We ignore the corresponding Hartree-type diagrams, which are relatively suppressed by a factor of K 2 /k 2 F [47]. In the Fock diagrams, the time folds are connected only by static disorder lines and not the dynamical interaction, and hence there is no distinction between the two time fold correlator and the realtime retarded particle-hole correlator. Thus, in the end we only need to focus on the contribution arising from Figure 4(d) and its partner, as was noted in Refs. [52,53], to get Σ L (0, 0). We have We do this sum by contour integration, noting the branch cut in V (Ω m , k) as iΩ m crosses the real axis and that n is a fermionic Matsubara frequency. The non-vanishing contribution upon analytically continuing n , ω l → 0 is [53] Σ R L (0, 0) = 2g(0)τ 2 (3.6) We have Note that −Σ R L (0, 0) is also the decay rate of the Cooperon at zero external pair momentum and frequency [53,56], which has been interpreted as the decay rate of electrons in exact eigenstates near the Fermi level [56][57][58].

B. Ladder diagrams
In the ladder diagrams with interaction rungs (Figure 5), the disorder correction to the interaction vertices occurs on a single time fold. Therefore the second term of Eq. (2.6) does not apply, as it would correspond to the bare interaction vertex connecting Green's functions on opposite time folds before being corrected by disorder, a possibility that is ruled out by the locality of the bare interaction vertex in time. Since the dynamic interaction (which can be interpreted to be mediated by a dynamically fluctuating boson) rung connects two time folds on opposite sides of the thermal circle, its propagator is given by a bosonic Wightman function [44,48,51] Note that only the dynamical part of the interaction V (ω m , q) − 4πe 2 /q 2 (which behaves like a Landaudamped boson) contributes to the Wightman function.
Direct Insertion.-We first consider the simplest summation of the ladder diagrams with alternating interaction and 'diffuson' rungs, L(ω, q), given by Figure 5a. By explicitly considering the series of diagrams, we see that the resulting unit, F , depends only upon the frequencies passing through it, but not the momenta. The Bethe-Salpeter equation for F reads The overall sign of the rung term is +1, coming from i 2 (−i) 2 , where the factors of i are generated by the Hubbard-Stratonovich transformation of the Coulomb interaction to a fermion-boson interaction and the two real-time fermionboson interaction vertices. After some manipulations, and assuming k 1 , k 2 are close to the Fermi surface, this becomes where we also set ω, q = 0 in the internal Fermion Green's functions, because the leading dependence on ω, q for small ω, q comes from the 1/(−iω + Dq 2 ) multiplying the integral. We can rewrite this for small k 0 , k 0 τ −1 µ as (3.12) As a matrix equation where the elements of A 0 are given by the integral kernel of the previous equation (3.12). Note that the translationally invariant structure of A 0 implies plane wave eigenstates. The growing part of f (ω, q) is obtained by appending external lines to F , capping off the ladder sum and integrating over momenta (which just provides two factors of g(0)τ = for ω τ −1 ) and frequencies ( Figure 5c): f (ω, q) = (g(0)τ ) 2 dk0 2π dk 0 2π F (ω, q, k 0 , k 0 ). Therefore F and A 0 have the same eigenvectors and the largest positive eigenvalue of A 0 (= πT 2 ) provides the growth exponent (3.14) Thus the growth exponent produced by the simplest 'direct' ladder insertion considered above is insufficient to overwhelm the T 3/2 decay rate from the self-energy corrections. We need to thus consider other ladder insertions at O(e 2 ) and check to see if they generate an exponent that successfully competes with the decay rate. Henceforth, we ignore the contribution of A 0 to the ladder sum. Exchange insertion.-As discussed above, we need to consider additional ladder insertions at the same order in perturbation theory which at least compete with the previously computed decay rate. At O(e 2 ), these come from Figure 5b. The sum of the two insertions gives the following integral equation: The overall sign of this contribution is +1 for the same reasons as above. Moreover, the two contributions are equal to each other. As before, we ignore the small ω, q contribution coming from within the integrand, and throw out the short-wavelength/high-frequency parts of the interaction. Since the interaction is long-ranged, the largest contribution to the integrals comes when the momentum k 2 − k 3 appearing in the internal interaction and in the 'diffuson' rungs is small compared to the momenta flowing through the internal fermion lines, which are O(k F ). We thus shift k 3 → k 3 + k 2 and then ignore k 3 everywhere except in the interaction and 'diffuson' rungs, which are singular at small k 3 . Then we have, This gives the matrix equation where the matrix elements of A 1 are given by the integral kernel in the last line of the above equation. As was the case with A 0 , the largest positive eigenvalue of A 1 comes from an eigenvector with constant entries. We thus obtain the net growth exponent after taking into account the dissipative self-energy: This returns Eq. (1.2) after a Fourier transform.

IV. ADDITIONAL CONSIDERATIONS
In the previous section, we computed the squared anticommutator and the leading O(e 2 ) correction to the growth exponent by doing an infinite resummation of the disorder lines. It is natural to ask the following questions: (i) Do ladder diagrams with a different skeleton structure of the disorder lines affect the exponent? (ii) What is the contribution of the other diagrams at O(e 2 ) that have been ignored in Figure 5 above? (iii) How sensitive are the above results to the specific form of the (Coulomb) interaction, V (|r − r |)?
We address all of these concerns one by one in this section.

A. Crossed disorder rungs
Instead of using the 'diffuson' rung, L(ω, q), considered thus far, we can sum diagrams with 'maximally-crossed' disorder rungs ( Figure 6). As is well known, this gives where Q is the total momentum of the incoming or outgoing particle-particle pairs. As with L(ω, q), ω is still the net lateral frequency transfer above as the disorder rungs cannot transfer frequency. At the non-interacting level, this gives It is easy to see that this expression does not have a pole at small ω, and hence we do not need to consider contributions with L c as the base unit (In two spatial dimensions, there is a logarithmic singularity at small ω that is still weaker than the pole in the contribution with L). We can also insert L c as an internal rung in the series with L as the base unit, such as by replacing in the integrand of Eq. (3.15). However, in this case, the same small momentum then does not appear in both the interaction and L c rungs, and the resulting contribution is thus less singular than the one in Eq. (3.15), scaling as subleading powers of T starting at T 2 . Similarly, we can consider insertions such as those in Figure 5, but with additional internal L rungs. These are also less singular than the ones shown for the same reason. At O(e 2 ) we have to also consider the diagrams shown in Figure 7. In the diagrams given by Figure 7(a), (b), the internal interaction line carries only the external frequency and momentum. We assume that the Coulomb interaction actually has a long static screening length ξ l, where l = v F τ is the disorder mean free path, and that we are probing scrambling at length scales x ξ.
The interaction line in Figure 7a is given by 3) The insertion in Figure 7(a) in the limit of small external frequency and momentum (ω, q) is then given by 1 cosh(βk 0 /2) cosh(βk 0 /2) The factor of i comes from (−i) 3 from the three advanced Green's functions, and the additional factor of N arises because the flavor indices on the left and the right sides of the diagram are decoupled. The partner insertion obtained by reflection about the horizontal axis is the complex conjugate of this, so their sum vanishes. For the insertion in Figure 7(b), the internal Wightman line is given by so this diagram is not important. In Figure 7 (c), the internal Wightman line carries the external frequency ω λ L T , so it can be approximated by 4πe 2 DK 2 p 2 T , where p is an internal momentum. However, in this case, once again, the same small momentum p does not appear in both the interaction and the internal L or L c , so this diagram ends up being less singular and scales as subleading powers of T starting at T 2 . For Figure 7 (d), the internal interaction line is just −i 4πe 2 K 2 and we get for the insertion, after appropriately shifting momenta, for both the internal L and internal L c cases Reflecting this insertion about the horizontal axis produces its complex conjugate, and reflection about the vertical axis effectively interchanges k 0 , k 0 . The four contributions then sum to zero.

C. Short-range interactions
Based on the analysis of Section III, we see that the Lyapunov exponent is simply given by as Im[V R (k 0 , k)] and Im[L(k 0 , k)] are both odd functions of k 0 for the interactions we consider. Since |1/ sinh(βk 0 /2)| > |2/ sinh(βk 0 )|, the first term of the above (coming from the ladder sum of Figure 5) always dominates the second (coming from the self-energy corrections), and the exponent is thus always positive if sgn(Im[V R (k 0 , k)]) = −sgn(k 0 ). For a short-range interaction that does not vanish as q → 0 (we take a contact interaction for which V R bs (q) = V 0 ), screening by the diffuson produces Inserting this into Eq. (4.7), we see that all the integrals converge, and that λ L ∼ +V 2 0 T 3/2 for d = 3. Thus, shortrange interactions behave qualitatively in the same way as Coulomb interactions from the point of view of scrambling, consistent with previous work on the inelastic scattering rate [30].

V. TWO SPATIAL DIMENSIONS
In two spatial dimensions, the diffuson-screened Coulomb interaction is [ We probe scrambling at length scales x much larger than the mean free path l and the screening length K −1 smaller than the eventual localization [59] length l e k F l of the electron wavefunctions [30] (The light-cone like growth of f (t, x) will be arrested beyond this localization length, i.e. the operator-radius R(t) is bounded by this length). Then, the same approximations and lines of reasoning we used in three dimensions also work in two dimensions, and the Lyapunov exponent is still given by Eq. (4.7) with d = 2. Inserting this dynamically screened Coulomb interaction, we obtain the leading contribution where R = 1/(e 2 D(dn/dµ)) is the sheet resistivity [47] and we restored factors of k B and . This cannot saturate the universal bound λ L ≤ 2πk B T / unless the effective coupling e 2 R /h becomes large, which also determines crossover or transition to an insulating state. According to experimental results reported in Ref. 60 and theory discussed in Ref. 61, the density-tuned metal-insulator crossover/transition occurs at around R ≈ 3h/e 2 , which is smaller than the value required to saturate the bound by about a factor of 3. This indicates that the metallic state has a Lyapunov exponent numerically, but not parametrically, smaller than the bound. From Eq. (5.2) above, we see that it contains the difference of two terms. The term being subtracted is the decay rate of electrons in exact eigenstates of the disorder potential [56,57], whereas the term being added gives the rate at which chaos spreads, i.e. how electrons would be infected within an epidemic picture (See Figure 1) if there were no electron 'deaths'. Both these terms individually contain a logarithmic infrared divergence, which cancels when their difference is taken. The logarithmic divergence in the exact eigenstate decay rate was removed in a self-consistent computation [58], by using the rate itself as an infrared energy cutoff, but this is not required here. For the exact eigenstate decay rate, the self-consistent computation provides instead a regularized logarithmic factor of ln(πD(dn/dµ)) [30,47,58], which doesn't appear in the Lyapunov exponent.
Let us comment now on why the logarithmic divergence cancels out in the expression for the Lyapunov exponent but appears in the exact eigenstate decay rate. It arises from an infrared divergence in the collision integral in Eq. (2.8) when the energy transfer in a collision approaches zero. At zero energy transfer, the interaction of the electrons with another particle-hole excitation (or equivalently the boson representing the Coulomb interaction) is like the electrons scattering off a random static potential. Each instance of such a scattering event can be described by a quadratic integrable Hamiltonian, and is hence incapable of producing chaos. However, this process still leads to decoherence of the individual electron wavepackets and hence contributes to the decay rate. A similar cancellation between singular pieces of self-energy and ladder contributions coming from zero energy transfer collisions was first pointed out by two of us in the computation of the Lyapunov exponent of a Fermi surface coupled to a gapless fluctuating gauge field in Ref. 51. For the short-range interactions considered in the previous section, the logarithmic factor still cancels in the Lyapunov exponent, and we obtain λ

VI. DISCUSSION
We have studied the spread of many-body quantum chaos due to electron-electron interactions in diffusive metals. We find that chaos spreads ballistically, even though quasiparticles are transported diffusively. This is because the spread of chaos is linked only to the propagation of quantum information about inelastic collisions of quasiparticles, which does not require the transport of quasiparticles themselves. In three dimensions, we found that the Lyapunov exponent scales as the inelastic scattering rate of quasiparticles, whereas in two dimensions the inelastic scattering rate is larger than the Lyapunov exponent by a logarithmic factor arising from 'classical' collisions that do not involve quantum fluctuations. In d spatial dimensions, we find λ L ∼ T d/2 , which leads to v B ∼ T d/4 . Comparing the form of the butterfly velocity to a scaling form v B ∼ T 1−1/z , where z is the dynamical exponent, we find that our result is qualitatively similar to that of a critical system with z > 1. While our computations in d = 2 and 3 were carried out with the 1/r Coulomb interaction, we expect similar results to hold in d = 2 for the ln r Coulomb interaction.
Remarkably, we find the above ballistic growth of operators even though the Coulomb interaction is long ranged and no microscopic Lieb-Robinson bound exists. This result is a particularly striking example of the idea that the butterfly velocity can function like a low energy Lieb-Robinson velocity [62]. It raises the question of what other long range models might be harboring an emergent ballistic growth of operators at low energy.
We note the recent experimental measurement by Kapitulnik et. al. of local thermal diffusivity using an optical method [63]. It would be interesting to measure the local heat diffusion constant in an interacting diffusive metal using this method. The heat diffusion constant is given by the ratio of thermal conductivity and specific heat; at low enough temperatures, in a regime where both of these quantities are dominated by the electronic contribution, it would be interesting to compare the measured diffusion constant to the known quasiparticle diffusion constant D that appears to be relevant to quantum chaos. While in the non-interacting case one expects the thermal diffusivity to be equal to D, significant deviations may arise due to interactions, especially in two dimensions [64].
In this work we only focused on disorder averaged cor-relation functions in the diffusive, ergodic phase. However, one could also ask about rare-region effects [65,66]. For example, can rare 'localized' regions in the ergodic phase impede the spread of chaos? How does the different inelastic scattering rate in these regions [67] affect the Lyapunov exponent? Alternatively could there be rare-regions, with very little disorder, that lead to an even faster butterfly velocity? In dimensions greater than one, the effect of such rare-regions are expected to be significantly suppressed, but we leave a detailed study for future work. Finally, it would also be interesting to study the growth of entanglement in an interacting diffusive metal, and compare it to the spread of chaos. We also leave this question for future study.
For H clean free , Eq. (3.1) simply factorizes by Wick's theorem into a product of a retarded Green's function and an advanced Green's function. When disorder is included, we have where T denotes time-ordering. The exponentials containing H dis free may now be expanded, this produces corrections to Eq. (3.1) with H = H clean free that can be contracted by Wick's theorem and the disorder average Eq. (2.2). Since the disorder is time-independent, this produces to lowest order the disorder self-energy corrections to the Green's functions (Figure 2a), and also the disorder ladder corrections in Figure 3a. These corrections can then be resummed to obtain the non-interacting f (t, x) as shown in Appendix B.
With the inclusion of interactions, we use It is helpful to consider for the purposes of this illustration the decoupling the four-Fermion interactions using a bosonic field ϕ(ω, k) with a propagator given by the unscreened Coulomb interaction V b (k). The perturbative expansion now generates the corrections shown in Figure 4 that involve the usual correction to the Green's functions due to interactions, along with a new set of corrections that involve the contraction of the boson field across the e −βH/2 thermal factors of Eq. The expression for this Wightman propagator V W is provided in Eq. (3.9), and its relation to the spectral function is derived in detail in Refs. [48,51]. These new corrections generate the diagrams shown in Figure 5. Note that interaction corrections to the e −βH/2 thermal factors in Eq. (3.1) correspond to the dressing of the Wightman propagators, which we take into account since we use the dynamically screened Coulomb interaction for V W in Eq. (3.9).
where D = v 2 F τ /d, and in the intermediate steps, we expanded in small q assuming that the largest contributions to the integrals come from the regions with k ∼ k F = mv F q, and that µ τ −1 |ω|. We assumed that L(ω, q) does not depend on any other combinations of momenta and frequencies passing through it apart from (ω, q), which turns out to be self-consistent. Each disorder rung is multipled by a factor of −i 2 = 1, where the i's come from the real-time electron-disorder vertices. We thus see that f (t, x) ∼ f 0 (t, x) + f 1 (t, x)e −x 2 /(4Dt) , where f 0 decays rapidly in time at a rate given by τ −1 and f 1 is a slowly varying function of space and time. Henceforth we ignore f 0 as we are interested in long times t τ and set f 1 to 1. Since there is no exponential growth in f (t, x) we conclude that the non-interacting disordered metal does not have many body quantum chaos.