Many-body chaos at weak coupling

The strength of chaos in large $N$ quantum systems can be quantified using $\lambda_L$, the rate of growth of certain out-of-time-order four point functions. We calculate $\lambda_L$ to leading order in a weakly coupled matrix $\Phi^4$ theory by numerically diagonalizing a ladder kernel. The computation reduces to an essentially classical problem.


Introduction
Non-time-ordered four point functions can be used to diagnose chaos in many-body quantum systems [1][2][3][4][5][6][7][8]. For example, we can understand the butterfly effect as the statement that for rather general operators W, V , the thermal expectation value of the square of a commutator should become large at late time, of order 2 V 2 β W 2 β . The way this function grows can be interesting. At least for simple operators in large N systems, one expects a long period of exponential growth [6], plus higher orders in N −1 that make the function eventually saturate. The rate defined by the exponential, λ L , is a measure of the strength of chaos. It satisfies λ L ≤ 2π β [9]. The purpose of this paper is to evaluate λ L in a weakly coupled large N quantum field theory at finite temperature. Ref. [7] suggested an approach to the weak coupling calculation, based on analogy to the BFKL [10,11] analysis of high energy scattering in --+ t β Φa'b' Φab Φab Φa'b' Figure 1: The squared commutator (4) can be expanded to four terms represented by the path integral contours shown. The vertical segment (ends should be identified) represents the imaginary-time circle, and the horizontal folds implement the real time evolution to produce Φ ab (t). The two folds are separated by half of the thermal circle.
gauge theories. In this paper we set up and carry out this calculation for a simple model system. Specifically, we consider the theory of a single Hermitian matrix field Φ ab in four spacetime dimensions, with the Lagrange density 1 The 't Hooft coupling is defined as λ = g 2 N . The goal is to compute λ L to leading order in λ. We will take a very direct approach, evaluating a subset of thermal Feynman diagrams for an index-averaged and spatially-averaged version of the squared commutator: Here ρ = ρ(β) is the thermal density matrix at inverse temperature β. Splitting ρ into the two √ ρ factors amounts to putting the two commutators on opposite sides of the thermal circle. It is a choice that does not affect λ L , but that makes some of the equations below a little simpler.
To generate the perturbation theory for C(t), it is helpful to think about expanding out the two commutators to give four terms. Each term could be computed by a particular analytic continuation of the Euclidean correlator. Equivalently, we can represent the four terms by path integral contours in complex time, where we append some real-time folds to the Euclidean thermal circle. This is illustrated in Fig. 1. In principle, for each term we should follow the usual procedure of expanding down powers of the interaction vertex, integrating them along the contour and connecting fields by contour-ordered propagators.
The quantity λ L is defined by the asymptotic rate of growth of the N −2 term in C(t). In order to compute this, we can restrict to planar diagrams. We can also restrict the region of integration for the interaction vertices to the real-time folds. The integral over the thermal circle implements corrections to the thermal state; such corrections would be important for getting the exact C(t), but we believe that they do not affect the spectrum of growth exponents. These two simplifications would be valid at any value of the coupling, but they are not enough to make the problem tractable. To compute λ L to leading order in the coupling, we can make another simplification, which is to keep only the fastest-growing function of time at each order. Roughly, we will sum all powers of λ 2 t, but ignore terms proportional to e.g. λ 3 t or λ. This simplification restricts the class of diagrams that we need to consider, and it also allows us to get by with simplified versions of the diagrams that we keep. The structure is very similar to the leading-log approximation in high energy scattering, but with t playing the role of log s [7].
The diagrams that must be summed consist of a set of dressed ladder diagrams. The rate of growth of the sum of ladders can be determined by finding the largest eigenvalue of a one-dimensional integral equation, which we diagonalize numerically. In the case where the bare mass m is nonzero but small compared to the temperature, we find The fact that λ L is proportional to 1/m indicates that the important degrees of freedom are the highly populated, frequently colliding low energy quanta with E ∼ m, not the thermal scale quanta that one might have expected. Naively, this result diverges for a massless field, but if the tree level mass vanishes we must include the one-loop thermal mass m 2 th = 2λ/3β 2 , giving This is still small, but it is parametrically enhanced relative to the naive λ 2 scaling. In addition to the parallels to BFKL, our calculation shares much in common with the analysis by Jeon [14] of the shear viscosity in weakly coupled φ 4 theory. Both the work of Jeon and the review of BFKL by Forshaw and Ross [15] were very useful guides.
Kitaev [16] has computed λ L in a strongly-coupled fermion quantum mechanics similar to the Sachdev-Ye model [17][18][19]. The diagrams are structurally similar to the ones that we study in this paper, but in that context they give the exact O(1/N ) answer as a function of the coupling, saturating the bound 2π/β as the coupling goes to infinity.

The ladder diagrams
In this main section of the paper, we will study a set of ladder diagrams for C(t) and derive an eigenvalue equation that determines the growth rate λ L at leading order in the coupling λ. This equation arises from the fact that an infinite ladder is not changed if we add one extra rung. There are two slightly subtle points in the analysis. The first is related to the fact that we are doing perturbation theory on a pair of folded time contours, and the interaction vertices should be integrated over both sides of each fold. The sum over the two sides turns the side rails in the ladder diagrams into retarded propagators, while the rungs remain Wightman correlators. The second subtlety is that we have to include self-energy corrections for the retarded propagators. In order to explain both of these points, we will go through the first couple orders of perturbation theory explicitly, in §2.2, §2.3 and §2.4. We will then analyze the ladder diagrams in §2.5, discussing corrections in §2.6.

Free propagators
For the computations below, we will need two types of correlation functions: the retarded propagator G R and a Wightman functionG with the operators separated by half of the thermal circle. The functions are defined by Here and elsewhere in the paper, a field operator without spacetime coordinates is assumed to be at the origin Φ ≡ Φ(0, 0). In momentum space, it is often useful to express these propagators in terms of the spectral function ρ(ω, |k|). (We will distinguish the spectral function ρ(k) = ρ(k 0 , |k|) from the thermal density matrix ρ by making the arguments explicit.) For a free field, we have ρ(ω, , and therefore We will now use these propagators to discuss the first couple of orders of perturbation theory for C(t).

Order λ 0
In the free theory, we simply sum over the three ways of contracting the four operators on the contours of Fig. 1. Two of these cancel when we take the sum of the terms to form the square of the commutator. The only nonvanishing contraction is the one in which the two fields on the bottom fold are contracted with each other, and likewise for the top fold. For this pattern of contractions, the sum over the four terms gives where the minus sign comes from the † reversing the order of the operators inside the commutator in Eq. (4).

Order λ
The first correction to the free theory comes from a term where we integrate one copy of the interaction vertex over both of the real-time folds. In Fig. 2  self energy diagram that results. Similar diagrams appear at higher orders in perturbation theory as well, decorating all propagators. The effect is a one-loop temperature-dependent correction to the mass of the field. In principle, this can be absorbed by using the "thermal mass" m 2 th (β) in all propagators. However, for nonzero m, λ L depends smoothly on the mass, so to leading order in λ we can just use the tree-level mass and ignore the one-loop self energy altogether. The exception is the case where the tree-level mass is zero; there we must include the thermal mass. We work out the value in appendix A.1:

Order λ 2
At order λ 2 , we integrate two interaction vertices over the contour. When these vertices are on the same Lorentzian fold, as in Fig. 3(a), we get the second term in the geometric series of one-loop self energies, as well as the first term of the two-loop self energy. We can account for the two-loop diagram (as well as similar diagrams dressing propagators at higher orders in perturbation theory) by including a λ 2 self energy correction in the propagators. The real part of the two-loop self energy is a momentum-dependent shift in the mass that can be ignored relative to the tree-level or one-loop mass discussed above. The imaginary part leads to exponential decay of correlation functions due to the finite Figure 4: The sum over indices in C(t) is equivalent to contracting the external operators in the four point function with semicircle caps. The one rung diagram then has three index structures that each contribute 16N 4 , where 16 = 4 · 4 is a combinatoric factor that arises from the possibility of "rotating" each of the vertices in the plane of the diagram. The total factor is 48N 4 ; dividing by N 4 to turn the sum into an average, we get 48. lifetime of a single particle state; at finite temperature, a particle can be knocked into a different momentum mode by a collision with a thermal excitation. We will see below that this process will contribute to the leading order λ L . We can include the relevant effect by modifying the retarded propagators to The computation of the two-loop width Γ k is standard [20,14], and we review it in Appendix A.2. We will also have a configuration in which one interaction vertex is attached to each of the Lorentzian folds. This leads to a qualitatively new diagram, the "one rung" diagram, illustrated in Fig. 3b. To analyze this diagram we would like to study the Fourier transform C(ω). This is slightly delicate, because we anticipate that the result of our computation will be an exponentially growing function C(t), for which the Fourier transform does not exist. Instead, we will actually study a Laplace transform C(ω) = ∞ 0 e iωt C(t) where the integral is over positive t only. We can recover C(t) at positive time by taking C(t) = dω 2π e −iωt C(ω) along a contour that runs above all singularities in the complex ω plane.
Including the combinatoric factor explained in Fig. 4, the one rung diagram is equal to where the rung function R contains the loop integral and the product of Wightman correlators This function R will be important in what follows. We evaluate the integral explicitly in Appendix B. The overall sign in (16) is −i 2 = 1, where the two factors of i come from the real-time interaction vertex, and the minus sign comes from the † in (4). We would like to understand how to extract the leading large-time behavior of (16). It is helpful to begin by understanding the expression with the free propagators. In this case, the leading large-time behavior of C one rung (t) is linear in t. In frequency space, this comes from a double pole in ω. Let us see how we get this behavior. The above expression contains two pairs of retarded propagators. With the free propagators, the first The integral of the complete expression over p 0 will have a contribution from taking the residues of these poles. We get terms proportional to ω −1 , to (ω + 2E p ) −1 and to (ω − 2E p ) −1 . The first term is the important one, because when we combine it with a similar term from the second pair of retarded propagators, we will get the desired double pole ω −2 .
The lesson is that we can get the correct large-time behavior by making the replacement Now we add back in the self energy correction, replacing the i factors in (18) with iΓ p . The shift in the location of the poles will affect the on-shell delta functions, but only by a small amount that we can ignore at leading order. The only effect that we need to include is to modify the replacement as A convenient feature is that the delta functions in this substitution, together with the delta functions in the Wightman functionsG, force each of the four momenta p, p and (p − p )/2 ± to be on shell. This will be important in comparing to a classical calculation below. Notice also that we are supposed to sum over positive and negative energy for each of the momenta.

Higher orders
At higher orders in λ the leading diagrams consist of dressed ladders. To sum these diagrams, it is useful to define the function f through The sum of ladders illustrated in Fig. 5. satisfies a simple integral equation This is an inhomogeneous equation. At large time we expect f (t, p) to be growing, whereas the homogeneous term will be decaying. We can therefore get the correct large-time rate of growth by dropping the homogeneous term. We also substitute in the approximate form for G R discussed above. The integral equation then becomes Because of the on-shell delta functions in all pairs of retarded propagators, f (ω, p) will be entirely supported for on-shell p. We therefore write where f (ω, p) and f (ω, p) are distinguished by their second argument. Substituting into (23), we get where the kernel m(k, p) is given in terms of the rung function R as For on-shell k and p, one can show that R(k + ) comes only from the first term in equations (46) and (49) while R(k − ) comes only from the second term. As shown in (43), we can actually write the decay rate Γ p in terms of the rung matrix, so the whole equation is In real time, this equation is d dt f = M f , where M represents the integral operator on the RHS. The eigenvalues of M determine the spectrum of growth exponents; the largest positive exponent is λ L . We do not know how to solve this integral equation analytically, but we can solve it numerically. Appendix C and Fig. 7 give some details on this. The results we find were already reported in the Introduction: for small mβ we find λ L ≈ 0.025λ 2 /mβ 2 . For larger values of mβ we find that λ L is exponentially small in βm.

Corrections
There are two qualitatively different types of corrections to the analysis above. The first class consists of corrections suppressed by the 't Hooft coupling λ. These come from other planar diagrams that have more than two powers of λ for each power of t. An example of such a correction is shown in Fig. 6(b). By studying decorated ladder diagrams of this type, we expect that one could compute successive corrections to λ L in powers of λ.
The second class is 1/N corrections. These may also correct λ L , but they have a second qualitatively different effect, which is to give N −4 e 2λ L t corrections to C(t). In fact, we expect a full power series in N −2 e λ L t that will sum to a function that saturates for large t. An example of a diagram that will give this type of correction is shown in Fig. 6(c). The two ladder sub-diagrams will give a contribution ∝ e 2λ L t , and the overall factor is 1/N 4 . These corrections are analogous to the multi-Pomeron exchange diagrams in high energy scattering, or multi-graviton exchanges for holographic theories. Actually, the N counting here is a little subtle; we haven't shown the double line diagram, but there is a special subset of index contractions that give a 1/N 2 piece in the double-ladder diagram. These correspond to taking just one of the three index structures in Fig. 4 for each rung in the ladders. The expected contribution from this subset of contractions would be proportional to N −2 e 2yt , where y is the leading eigenvector of an equation where we divide the first term in (27) by three. We believe that y = 0, both based on numerics and intuition to be described in point 2 of the Discussion; this contribution is then N −2 uniform in time, and therefore ignorable relative to the N −2 e λ L t piece that we calculated.
We can make this more precise by using (46) to write (25) as where we defined We would like to interpret f + (p, t) as the expected number of infected particles of type p at time t and f − (p, t) as the expected number of infected holes. This includes the sum over all indices, so the expected number of a given index type would be f ± /N 2 . The first line in the evolution equation (28) represents the loss of infected particles of type p due to collisions with other particles. The third line represents the fact that if an infected particle of type k collides with a thermal particle, we gain two new infected particles. The factor 96λ 2 is 2|T | 2 for the theory we are considering. In such a collision, we also gain one infected hole, for the state k that is no longer occupied. Similarly, the scattering evolution of an infected hole results in two infected holes and one infected particle. The second line represents this process. The only surprise in the collision integrals of (28) is the following: we might have expected to have final state bose factors 1+n(E p ) = (1 − e −βEp ) −1 for the new infected particle in the collision integrals of equation (28). We do not have such factors.
2. From this perspective, we can understand the relative factor of three between the two terms in (27): each collision involving an infected particle results in the loss of one and the creation of three (two particles and one hole). It is because of this factor of three that λ L is positive and C(t) grows.
3. A surprising feature of our result for λ L is the 1/m behavior, a type of IR divergence. We can understand where this comes from by looking at (28), and putting in the massless form E k = |k|. If we go to very low momentum we can expand the exponentials in the denominator. We then have six powers of momentum in the denominator and five in the numerator. So if we scale a candidate eigenvector f (p) toward smaller momentum, we increase the eigenvalue.
This IR effect also shows up in a way in the imaginary part of the self energy. For a massless field, the two-loop decay rate Γ p is proportional to 1/|p| at small p [20,14]. This means that low-momentum modes are experiencing lots of collisions. This does not lead to IR sensitivity for the shear viscosity, because these modes don't carry much momentum, but it does affect λ L .
4. The rate of growth of the epidemic is a measure of the many-body butterfly effect, but it is not exactly a Lyapunov exponent for the gas of particles. That would be defined as the growth of an infinitesimal perturbation. In the case at hand, the perturbation is only small in the sense that it initially affects a small fraction of the total number of degrees of freedom N 2 .
It would be nice to know if there is some other effective classical description where λ L approaches a Lyapunov exponent. Presumably the holographic dual has this property. Another possibility, in case of very small mass where the relevant modes are highly occupied, would be to study the nonlinear classical field equations in the spirit of recent work on one dimensional matrix systems of [24]. Previous work [25][26][27] relating the quantum Boltzmann equation (at high occupation number) to nonlinear classical field equations may be useful.
5. As pointed out in [7], there is another weakly coupled field theory where λ L was computed long ago. This is four dimensional gauge theory on hyperbolic space at temperature β = 2πR. Up to beta functions, which are higher order in the coupling, this problem is equivalent to a high energy scattering problem in the Minkowski vacuum. The BFKL analysis of this problem determines λ L as This is proportional to λ at weak coupling, rather than λ 2 or λ 3/2 as in this paper. The key difference is that the BFKL ladders are cubic ladders because of the cubic Yang-Mills coupling. The vertices come with √ λ, so one rung costs a single factor of λ. One might wonder whether weakly coupled gauge theory at finite temperature in flat space would have λ L ∝ λ for the same reason. We believe that it would not, since in flat space one cannot make on-shell ladders out of cubic vertices.
6. It would be nice to compute λ L in other weakly coupled theories, and to explore the behavior of C(t, x) where we do not integrate over the spatial separation of the operators. This might give some insight into the diffusive vs. ballistic behavior in the growth of operators discussed in [5,7].
Subtracting off the zero temperature integrand and then doing the integral, we find the thermal mass A.2 The imaginary part of the two-loop self energy To lighten the notation, we will suppress the spatial momenta, restoring them below. Then the two-loop self energy is We can use the spectral representation to write the correlators as The integral over τ is now simple. Using e iωnβ = 1, we get We now do the following: (i) we continue iω n → p 0 + i , (ii) we use Im[ 1 x+i ] = −πδ(x) to find the imaginary part, (iii) we restore the integral over spatial momenta and (iv) we use the energy conservation delta function to replace the (1 + e −βE ) factors by 2 sinh βE 2 . The result is This can be reduced further, to a one-dimensional integral [20,14]. However, we are interested in deriving an expression for Γ p = −Im[Π(E p + i , p)]/2E p in terms of the rung function. To get this, we substitute in for ρ in terms ofG: where the p in the argument of R can be either of (±E p , p). We can use the free expression forG (12) to write this in a form useful in the main text:

B The rung function
We will evaluate the rung function A very similar integral was provided in [14]. Plugging in forG using (12) and doing the p 0 integral with one of the delta functions, we find One of the four possible configurations of s,s gives zero. The sum of the other three is At this point, we use rotation invariance to set k = (|k|, 0, 0), and we decompose p = (p 1 , p ⊥ ). The integrand depends on p ⊥ only through p 2 ⊥ , so we can replace d 3 p → πdp 1 d(p 2 ⊥ ). Specializing to the dispersion relation E p = m 2 + p 2 , we can write p 1 , p 2 ⊥ in terms of E ± . Scanning over all values of p 1 , p ⊥ , we cover all positive values of E ± satisfying the constraint One also finds the Jacobian for the change of variables In these variables the integral can be done straightforwardly, with the upper and lower limits of integration determined by the constraint (47) and the remaining delta function. One finds that the first term in (46) is zero unless (k 0 ) 2 ≥ k 2 + 4m 2 , and the second term is zero unless k 2 ≥ (k 0 ) 2 . All together, we find where k 2 = −(k 0 ) 2 + k 2 and Note that the variable x ± is real whenever one of the θ functions is nonzero.

C Some details on finding λ L
We start by converting the three-dimensional integral equation to a one-dimensional equation, for the function f (ω, |p|). In this section we will use the notation P = |p| and K = |k|.
To get the one-dimensional equation, we have to integrate over the angle between k and p. It is simplest if we introduce a variable y = |k − p| so that The eigenvalue equation we want to solve is then where the angle-averaged kernel is written in terms of In this expression, E ± = |E K ± E P | and For nonzero mass, the integral in (53) must be done numerically. The matrix m 1 (P, K) is not symmetric, because of the factor of K/P . We will fix this below. We now change variables from P, K to u, v in order to map the half-line into the unit interval, which we will then discretize uniformly. A convenient choice is where P 0 determines the scale of momenta that receive the most attention in the discretization. For small mass we found P 0 ≈ 3m to be good, and for large mass P 0 ≈ m. Changing variables in the integral, and defining (the below contained a typo in v1) the eigenvalue problem becomes This is now a symmetric integral equation that can discretized as matrix multiplication. One finds a "continuum" of negative eigenvalues corresponding to approximately localized eigenvectors, and a discrete set of positive eigenvalues corresponding to smooth delocalized eigenvectors. λ L is defined as the largest positive eigenvalue. This can be determined by exact diagonalization or by iteration. For small mass, the largest negative eigenvalue is larger in magnitude than λ L , so a simple iteration of (59) will not work, but one can iterate a procedure where you update f to be a weighted average of the previous f and the result of applying the matrix to the previous f . For an appropriate choice of weighting this converges to the leading positive eigenvector. See Fig. 7 for some plots of λ L and the leading eigenvector.  Figure 7: Top: the combination (β 2 m/λ 2 )λ L is shown as a function of βm. This combination approaches a constant ≈ 0.025 at small mass. Although this is no obvious from the log scale plot here, it decreases exponentially in βm for βm > 1. We used a uniform discretization of the u interval into 640 lattice points. Bottom: the eigenvector f 2 (u) corresponding to λ L for small mass and large mass.