Four-point function in the IOP matrix model

The IOP model is a quantum mechanical system of a large-$N$ matrix oscillator and a fundamental oscillator, coupled through a quartic interaction. It was introduced previously as a toy model of the gauge dual of an AdS black hole, and captures a key property that at infinite $N$ the two-point function decays to zero on long time scales. Motivated by recent work on quantum chaos, we sum all planar Feynman diagrams contributing to the four-point function. We find that the IOP model does not satisfy the more refined criteria of exponential growth of the out-of-time-order four-point function.


Introduction
Matrix models are useful toy models of gauge theories and holography. Strongly coupled quantum field theories are difficult to understand directly, having a prohibitively large set of Feynman diagrams that must be summed. A good model should have a sufficiently small and well-organized set of diagrams, allowing for the computation of the full planar correlation functions. The diagrammatic structure should, however, be sufficiently nontrivial so as to capture the essential features of the bulk.
The IP model [1] is a simple large-N system of a harmonic oscillator in the U (N ) adjoint representation plus a harmonic oscillator in the U (N ) fundamental representation, coupled through a trilinear interaction. It has the same graphical structure as the 't Hooft model of two-dimensional QCD [2]. The IOP model [3] is a more tractable variant of the IP model. It possesses the same degrees of freedom, but the trilinear interaction is replaced by one that is quartic in the oscillators but quadratic in the U (N ) charges. Building on ideas of [4], the IP and IOP models were introduced in [1,3] as toy models of the gauge theory dual of an AdS black hole. These models capture a key property of black holes: the long time decay of the two-point function at infinite N , but not at finite N [5].
In this paper we compute the thermal four-point function in the IOP model in the planar limit. The motivation for studying the four-point function comes from recent work in quantum chaos and holography [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22]. A signature of quantum chaos in a large-N theory is the exponential growth in time of the connected out-of-time-order four-point function [23]. The growth rate is identified as a Lyapunov exponent. A black hole has a Lyapunov exponent of 2πT [7,9], which is the maximal possible Lyapunov exponent [8]. The significance of the out-of-time-order four-point function as a diagnostic for the viability of a model of holography was recognized in [6].
In Sec. 2.1 we review the role of the two-point function as a diagnostic of thermalization. In Sec. 2.2 we review the role of the out-of-time-order four-point function as a diagnostic of chaos. In Sec. 2.3 we briefly mention the Sachdev-Ye-Kitaev model [7,24], which was recently recognized to be maximally chaotic [7]. We point out that the random coupling can, to leading order in 1/N , be replaced by a quantum variable.
In Sec. 3.1 we review the calculation of the planar two-point function in the IP model. In Sec. 3.2 we compute the planar four-point function. This involves summing ladder diagrams, which can only be done analytically in the limit of small adjoint mass, to which we restrict ourselves.
In Sec. 4.1 we review the planar two-point function in the IOP model. In Sec. 4.2 we compute the planar four-point function. Diagrammatically, the IOP model is more involved than the IP model. However, it has the advantage of allowing analytic computations for any adjoint mass. For both the IP and IOP models, we work in the limit that the mass of the fundamental is heavy, as compared to the temperature.
In the regimes considered, we find that the IP and IOP models are not chaotic. Some speculations on why this is so, and possible modifications of the models, are mentioned in Sec. 5.
2 Thermalization, chaos, and large N

Thermalization
Holography has provided useful insights into both strongly coupled field theories, as well as their gravity duals. A well-studied property of a black hole is its approach 2) studied in [4].
to equilibrium after a perturbation. A two-point function computed in a black hole background exhibits late time decay of the form [25,26], where c is an order-one constant and β is the inverse temperature. The late time decay of the two-point function has a clear interpretation in the bulk: matter falls into the black hole, but classically nothing escapes. Computing subleading corrections in G N to (2.1) does not prevent the late time decay. As recognized in [5], the late time decay to zero of a two-point function is inconsistent with the properties of a finite entropy quantum mechanical system. On the field theory side, one thus has the statement that, to all orders in 1/N , the two-point function decays to zero at late times, even though this property does not hold nonperturbatively in 1/N . The two-point function φ(t)φ(0) β can be regarded as the overlap between the states φ(0)|β and φ(t)|β ; its decay is a probe of thermalization. Therefore, the large N limit acts like a thermodynamic limit [4].
This late time breakdown of perturbation theory was studied in the context of matrix quantum mechanics in [4]. Reducing Yang-Mills on a sphere in terms of spherical harmonics, one obtains a Hamiltonian whose essential features can be captured by considering just two interacting matrices. For instance, [4] considers two large N matrices M 1 , M 2 with a Hamiltonian, The relevant diagrams for the decay of the two-point function are the sunset diagrams shown in Fig. 1. The model (2.2) has the drawback of still being too complicated to allow the summation of all planar Feynman diagrams. The goal of [1] was to find a matrix model i i j Figure 2. The basic graphical unit of the IP model (2.3) studied in [1]. It is like the diagram in Fig. 1, but cut in half. A single line is a fundamental, a double line is an adjoint.
that is more tractable, while still exhibiting the late time decay of the planar two-point function. The IP model [1] is given by the Hamiltonian, where a i is the annihilation operator for a harmonic oscillator in the fundamental of U (N ), while A ij is the annihilation operator for an oscillator in the adjoint, and As we review in Sec. 3, the planar two-point function can be found if one takes the mass of the fundamental to be large compared to the temperature, M T . For a general mass m for the adjoint, the planar Schwinger-Dyson equation for the two-point function can be solved numerically, exhibiting the desired late time exponential decay. In the limit of small mass for the adjoint, m → 0, the two-point function can be found analytically, giving late time power law decay.
A variant of the IP model, the IOP model, was introduced in [3], This model has the feature that analytic computations are possible for any mass m. It again exhibits power law decay of the two-point function at long times.

Chaos
Chaos in deterministic systems is understood as aperiodic long-term behavior that exhibits sensitive dependence on initial conditions. Two points in phase space, characterized by a separation δx(0), will initially diverge at a rate, where κ is the Lyapunov exponent. For a number of reasons [27], there is no straightforward extension of chaos to quantum systems. In the semiclassical regime, [23] gave an intuitive definition of chaos.
Replacing the variation in (2.5) by a derivative, and noting that this is given by a Poisson bracket, the generalization to quantum systems consists of replacing the Poisson bracket by a commutator. The commutator is generally an operator, so seeing exponential growth requires taking an expectation value. The expectation value of the commutator in the thermal state will vanish, as a result of phase cancelations. A simple way to cure this is to consider the square of the commutator [23], 2 Alternatively, one can consider the thermal expectation value of the commutator times the anticommutator; this will scale as . Either of these consist of sums of out-oftime-order four-point functions. The important point is that a chaotic system has an out-of-time-order four-point function that exhibits exponential growth. The exponential growth persists until a time of order −κ −1 log , at which point the commutator saturates at an order one value. For a large N field theory, 1/N plays the role of , and the classical limit is the infinite N limit. For matrix models, such as the IP and IOP models, leading order in 1/N corresponds to keeping the planar Feynman diagrams. The criteria of chaos for evaluating the viability of a model is a powerful one, that was recognized in [6]. A good model of a strongly coupled gauge theory should having an exponentially growing out-of-time-order four-point function. Moreover, if it is to be dual to a black hole, the Lyapunov exponent must match that of a black hole [7,9].

Thermalization and chaos
There is generally an intimate connection between thermalization and chaos. In the context of classical systems, there is a precise version of this statement [28], which we now review.
Letting A and B be regions of phase space, occupying phase space volumes µ(A) and µ(B), respectively, and letting φ t denote time evolution, a dynamical system is said to be mixing if µ [φ t A ∩ B] → µ(A) · µ(B) as t → ∞, for all sets A and B. In the notation of quantum mechanics, this is the statement that a system is mixing if the (connected) two-point function of any two operators decays to zero at late time. A system is defined to be ergodic if for every function f , the time mean of f (x) is equal to the space mean of f (x). It is shown in [28] that mixing implies ergodicity, but ergodicity does not necessarily imply mixing.
It is important to note that for a system to be mixing, the two-point function of all operators must decay. In fact, the IP and IOP models do not satisfy this criteria, as it is only the two-point functions of the fundamentals that exhibit late time decay. 3 The adjoints have a two-point function of a free harmonic oscillator; they have no self-interaction, and the interaction generated via the fundamentals is 1/N suppressed. Thus, exponential growth of the out-of-time-order four-point function for the fundamentals is a more refined criteria than the decay of the two-point function of the fundamentals at long times.

SYK model
Kitaev has proposed a variant of the Sachdev-Ye model [24] as a model of holography [7]. The SYK model consists of N 1 Majorana fermions χ i with a quartic interaction with random coupling J jklm , where couplings are drawn from the distribution, giving a disorder average of, Remarkably, one can analytically compute the disorder averaged large-N correlation functions in the SYK model at finite temperature and strong coupling, βJ 1. The two-point function exhibits exponential late time decay, see [7,24,29,30]. The out-oftime-order four-point function exhibits exponential growth [7], For studies of the four-point function, see [7,21,31,32]. An important aspect of the SYK model is the quenched disorder: if the coupling J jklm where instead a fixed constant, there would be additional Feynman diagrams that would contribute at leading order in 1/N . Here we simply point out that the disorder J jklm can be replaced by a quantum variable, as the quantum corrections are 1/N 3 suppressed.
Recall that the disorder averaged expectation value of an operator O composed of the fields χ i is, The interpretation of (2.12) is that one first computes the expectation value O with some coupling J jklm drawn from the distribution (2.9), and then averages over the J jklm . If one were to instead treat J jklm as a static quantum variable, then the expectation value of O would be given by, In terms of Feynman diagrammatics, if J jklm is a classical Gaussian-random parameter, then it has a two-point that is exactly 3!J 2 /N 3 . If instead J jklm is a quantum variable, then its leading two-point function can be chosen to be 3!J 2 /N 3 , however this will receive quantum corrections, as shown in Fig. 4. Thus, generally (2.12) and (2.13) are different. However, for the SYK model, the first quantum correction is suppressed by a factor of 1/N 3 : the loop diagram in Fig. 4 has two J jklm propagators, giving a factor of (3!J 2 /N 3 ) 2 . So, at leading order in 1/N , (2.12) and (2.13) are the same. Equivalently, the effective action for J jklm is at leading order in 1/N . Note that the structure of the vacuum is different depending on if J jklm is quenched disorder or a quantum field: the vacuum loop scales like N , and receives a correction of the same order from interactions with χ i , as there is now a summation over the indices. This, however, is irrelevant for the purposes of connected correlation functions. The variable J jklm is still not yet a standard quantum variable, due to the constraint that it be static. There are a few somewhat artificial ways to achieve this. One could add to the action a termJ jklm φ, where φ is some Lagrange multiplier field. A better option is to regard J jklm as the momenta of harmonic oscillators for which the frequency is taken to zero. Consider a harmonic oscillator with the standard Lagrangian, (mẋ 2 − mω 2 x 2 )/2. The Euclidean two-point function for the momentum is p(t)p(0) = mωe −ωt /2. Now take the limit of ω → 0, so as to remove the time dependance. Letting mω = 12J 2 N −3 , the momenta have the same two-point function as (2.10).

IP model
The IP model [1] is a quantum mechanical system, with a harmonic oscillator in the adjoint of U (N ) and a harmonic oscillator in the fundamental of U (N ), coupled through a trilinear interaction. The Hamiltonian for the IP model is given by (2.3). The twopoint function is found by summing rainbow diagrams (see Fig. 5) and is reviewed in Sec. 3.1. The four-point function is given by a sum of ladder diagrams (see Fig. 6), which we evaluate in Sec. 3.2.

Two-point function
The bare zero temperature propagator for the fundamental is defined as, Trivially, one has that, It will be assumed that fundamental has a large mass, M T , where T is the temperature. In this case, the bare finite temperature two-point function is the same as the bare zero temperature two-point function.
The adjoints have no self-interaction, and the backreaction from interactions with the fundamental is suppressed by 1/N . Thus, the propagator for the adjoint is that of a free oscillator in a thermal bath, where we have defined y = e −m/T . It will be useful for later to note that in the limit that the adjoints become massless, m → 0 (y → 1), their two-point function becomes, The planar two-point function for the fundamental is found by summing rainbow diagrams. The Schwinger-Dyson equation for the two-point function is given by (see Fig. 5): where the 't Hooft coupling is λ = g 2 N . In general, such an integral equation is difficult. However, the assumption that M T implies that G(t) = 0 for t < 0. As a result, G(ω) has no poles in the upper half plane, allowing us to close the integration contour in (3.5) in the upper half plane ω plane. Picking up the residues at ω = ω ± m, the Schwinger-Dyson equation turns into a difference equation, This can be solved numerically [1], however to proceed analytically we take the limit of small adjoint mass and small 't Hooft coupling, In this limit one finds [1], Here the ω should really be an ω + i ; we will generally suppress the i , remembering that all the poles are in the lower half complex plane. The Fourier transform of the two-point function is a Bessel function, We will later encounter integrals of a similar form, so we show (3.9) in some detail. For positive times, the ω contour in (3.9) wraps around the branch cut stretching from − √ 2ν < ω < √ 2ν. 4 Using (3.8) and moving the square root to the numerator, we rewrite (3.9) as, (3.10) The integral of the first term vanishes, while the second is twice a line integral, which gives (3.9). Now let us redo the calculation for the Fourier transform (3.9) slightly differently. Taking G(ω) in the form (3.8) and changing variable to, Our ω integral was from −∞ < ω < ∞. For positive time, we close the contour in the lower half plane. The branch cut is slightly below the real axis, and so is inside the contour. We can shrink the contour so that it hugs the branch cut. For negative times, the ω integral is closed in the upper half plane, and so gives zero. Also, our choice of location for the branch cut corresponds, for instance, to writing √ ω 2 − 2ν 2 = ω exp 1 2 log(1 − 2ν 2 /ω 2 ) and taking the principal branch for the logarithm. gives, The ω contour in (3.9) that hugs the branch cut maps into an x contour that is a circle of radius √ 2ν and centered around the essential singularity at the origin. Using the integral representation of the Bessel function, where the contour circles clockwise around the origin, we have, which is equal to (3.9). At late time, νt 1, the propagator decays as G(t) ∼ t −3/2 .

Four-point function
We now turn to the connected four-point function. In the planar limit, it consists of a sum of ladder diagrams, as shown in Fig. 6. The ingoing momenta are ω 1 , ω 2 , while the outgoing momenta are ω 3 , ω 4 . 5 As in the case of the two-point function, to proceed analytically we must work in the limit specified in (3.7). In this limit, the propagator for the adjoint is given by (3.4).
Consider the ladder diagram that consists of a single rung. It is given by, into (3.16), evaluating the integral, and then taking the m → 0 limit, yields for (3.16), We now sum all the ladder diagrams. As a result of the limit (3.7), all the pieces appearing in the Feynman diagrams are on-shell. Defining G 4 (ω 1 , ω 2 , ω 3 , ω 4 ) = δ(ω 1 − ω 3 )δ(ω 2 − ω 4 )G 4 (ω 1 , ω 2 ), and letting n denote the number of rungs, we have , (3.19) 5 The ingoing momenta are drawn in Fig. 6 as coming from the upper left and lower right in order for the diagram to look planar.  where ν was defined in (3.7).
The Fourier transform of (3.19) gives the position space four-point function, where we have defined t 31 ≡ t 3 − t 1 , t 42 ≡ t 4 − t 2 . In addition, G(ω) really denotes G(ω + i ); we suppress the i , remembering that, if we are using G in a time-ordered correlator, all the poles are in the lower-half complex plane.

Free propagator
The propagator entering the four-point function (3.20) is given by (3.8). As a warmup, it is useful to study (3.20) with the free propagator (3.2), rather than the dressed one.
In this case we have, Performing the ω 2 integral, and closing the contour in the lower half plane, we pick up poles at ω 2 = 0 and ω 2 = ν 2 /2ω 1 , Using the integral representation of the Bessel function (3.14), we get, Eq. 3.23 is the time-ordered four-point function, as evidenced by the theta functions.
We can obtain the out-of-time-order four-point function by dropping the theta functions. In particular, setting t 31 = −t 42 = t gives, In the limit νt 1, By summing only a subset of the Feynman diagrams: the ladder diagrams with undressed propagators, we have found exponential growth in the out-of-time-order fourpoint function. While intriguing, using the free propagator is certainly not legitimate, as it violates unitarity; classically it would be equivalent to violating Liouville's theorem. However, before evaluating (3.20) with the dressed propagator, it will be instructive to study (3.21) a bit further.
Returning to (3.22), and taking the limit of large t 31 , t 42 , we approximate the integral via the method of steepest descent (see Appendix A). This involves deforming the contour of integration in order for it to pass through the saddle point, at an angle so as to maintain constant phase. The saddle point of the exponent, occurs atω 1 = ±ν t 42 /2t 31 . As we continue from a time-ordered four-point function, to an out-of-time-order four-point function, t 42 → −t 42 , the saddle moves off of the real axis and onto the imaginary axis. At t 31 = −t 42 = t, the saddle is atω 1 = ±iν/ √ 2. The leading exponent in the integral in (3.22) can therefore be approximated by, reproducing (3.25). Let us also reproduce (3.23) by returning to (3.19) and computing the Fourier transform of each term before taking the sum. From (3.19) and (3.2) we have, (3.28) The Fourier transform gives, where we have made use of the series definition of the Bessel function. The expression (3.29) is easy to see directly in time-space. Since the free two-point function for the fundamental is simply θ(t) (3.2), a ladder diagram with n rungs will have n + 1 propagators for the fundamentals on each of the two sides. For one such side we have a product of theta functions, with the time insertions of the rungs integrated over. For the top side, and similarly a factor of t n 42 /n! from the bottom side. Accounting for the coupling at each vertex, −ig, as well as the sum over indices, and the factor of m −1 (1−y) −1 coming from the adjoint propagator, we recover the sum in (3.29).
If we wish to form an out-of-time-order four-point function, for instance with t 42 < 0, then on the bottom edge of the ladder diagrams, time runs backwards: we must use a two-point function that is θ(−t) rather than θ(t). In addition, since time is running backwards on the bottom edge, the interactions come with a factor of ig, instead of −ig. This results in the elimination of the minus sign in the sum in (3.29), and correspondingly gives exponential growth.

Dressed propagator
We now return to the frequency-space four-point function (3.19), and evaluate the Fourier transform (3.20), this time using the full dressed propagator. Inserting the propagator G(ω) (3.8) into (3.20) gives, where we have split off a G(ω 1 )G(ω 2 ) from (3.19), giving the first term in (3.31). Changing integration variables to (3.32) Our goal is to see if (3.32) exhibits exponential growth; if this does occur, it will be in the out-of-time-order regime, such as t 31 = −t 42 = t. We consider the late time limit, 6 and approximate (3.32) via the saddle point method (Appendix A): we seek to deform the contours of integration of x 1 , x 2 such that they pass through a saddle, at an angle such that the phase is constant. If we are away from the poles of the integrand, the saddle point occurs at x i = ± √ 2ν, which clearly only gives oscillatory behavior. Now consider the regions at the poles of the integrand, at x 1 x 2 = 2ν 2 . This a peculiar region, as is invariant under x → 2ν 2 /x. Inserting this x 2 = 2ν 2 /x 1 into the exponent in (3.32) , the exponent becomes,

IOP model
We now turn to the IOP model [3]. Like the IP model, this is a quantum mechanical system, with a harmonic oscillator in the adjoint of U (N ) and a harmonic oscillator in the fundamental of U (N ). However, the interaction is now quartic in the oscillators (2.3), and quadratic in the U (N ) charges. The latter property makes the IOP model more analytically tractable than the IP model, although diagrammatically it is more involved. As in the IP model, we consider the limit in which the fundamental is heavy, M T . However, we can now obtain analytic results at any mass m for the adjoint. We review the two-point function in Sec. 4.1, and compute the four-point function in Sec. 4.2.

Two-point function
The bare propagator for the fundamental is the same as in the IP model (3.2). The propagator for the adjoint is that of free harmonic oscillator in a thermal bath, defined by L(t)δ il δ jk = T A ij (t)A † kl (0) , and giving, The Schwinger-Dyson equation for the planar two-point function for the fundamental is (see Fig. 7), As G only has poles in the lower-half plane, we can close the ω i integrals in the lowerhalf plane and pick up residues only from L. This leads to an algebraic equation for G, with the solution [3], where the 't Hooft coupling is λ = hN . The propagator has a branch cut from ω − to ω + , leading to late-time power law decay, t −3/2 .

Four-point function
We now turn to the four-point function in the planar limit. The connected four-point function is found by summing ladder-like diagrams, shown in Fig. 10, where each "rung" of the ladder is found by summing an infinite number of diagrams, like the ones shown in Fig. 8. We warm up by computing the four-point function in the limit of small adjoint mass m, before doing the calculation for arbitrary m.

Small adjoint mass
We start with the limit m → 0. In particular, where κ is held constant. In this limit, the two-point functions for the adjoint (4.1) and the fundamental (4.4) become, To compute the four-point function, we first sum the diagrams shown in Fig. 8, to get where ω ij ≡ ω i − ω j and, , (4.9) where the index n/m labels the number of intermediate fundamental propagators on the top/bottom edge. As in the IP model, as a result of the m → 0, all intermediate propagators are on-shell. The four-point function is given by the ladder-like sum of the Γ (see Fig. 10), (4.10) Inserting (4.9) into (4.10) gives the frequency-space four-point function G 4 (ω 1 , ω 2 , ω 3 , ω 4 ) = (2π) 2 δ(ω 13 )δ(ω 24 )G 4 (ω 1 , ω 2 ) where, . (4.11) Like in the IP model, we find exponential growth in the out-of-time-order fourpoint function if we only sum the diagrams containing the free propagator: (4.11) with (3.2) and t 31 = −t 42 = t gives a four-point function ∼ N −1 exp(2κt) for large t.
Now consider (4.11) with (4.7). The position-space four-point function is thus, (4.12) Changing integration variables to x i = ω i + ω(ω − 4κ), (4.12) becomes (4.13) We approximate the integral by taking the limit of large time separations, and looking for saddle points which could give rise to exponential growth. Picking up the pole at x 1 x 2 = 2κ(x 1 + x 2 ), the exponent becomes, Like in the IP model, there is no exponential growth.

Arbitrary adjoint mass
We now compute the four-point function, with the adjoints taking arbitrary mass m. The Feynman diagrams contributing to "rung" Γ are shown in Fig. 8. A term in this sum, having n fundamental propagators on the upper edge and m on the lower, is given Figure 9. One of the diagrams entering Γ in Fig. 8, given by n = 2, m = 1 in (4.15). by, (4.15) where the ingoing frequencies are ω 1 , ω 2 , the outgoing frequencies are ω 3 , ω 4 , and we have defined r 2 = r 1 + ω 1 − ω 3 , and suppressed an overall factor of N −1 . In Fig. 9 the n = 2, m = 1 diagram from Fig. 8 is shown in more detail. Since G(ω 1 − p i ) has poles in the upper half p i plane, we close the contour in the lower half plane. Similarly for the q i integral. This gives for (4.15), Evaluating the integral over r 1 by closing the contour in the upper half-plane, (4.16) becomes, 7 To sum over all the diagrams contributing to Γ (see Fig. 8), we must sum (4.17) over n, m from 0 to infinity. This gives Γ = yΓ where,  where we have defined, , (4.19) and have simplified notation to denote ω j by j, and recall that κ ≡ λ/(1 − y). One can also rewriteΓ in (4.18) as, 20) which, recalling that ω 1 +ω 2 = ω 3 +ω 4 , is manifestly symmetric under ω 1 ↔ ω 2 , ω 3 ↔ ω 4 .
Attaching external propagators to (4.18) gives the first term in the sum for the four-point function shown in Fig. 10. The second term requires gluing two of theΓ together, where ωā = ω a + ω 4 − ω 1 . Performing the integral in (4.21) by closing the ω a contour in the upper-half plane gives, where both theΓ ×Γ on the left, and theΓ on the right, are functions of the external ω i .
(4.26) The four-point function is given by S, with external propagators attached.
Thus, the connected four-point function for the IOP model in the planar limit is, where (1, 4) , (4.29) where j denotes the frequency ω j , the propagator G(i) for the fundamental is given by (4.4), the constant y is the Boltzmann factor y = e −m/T where m is the mass of the adjoint and T is the temperature, and z(j, l) was defined in (4.19) and is a function of G(j), G(l), and κ = λ/(1 − y), where λ is the 't Hooft coupling. In the limit of small adjoint mass m (y → 1), the first term in (4.27) survives and reproduces the earlier result (4.11). The out-of-time-order four-point function does not exhibit exponential growth with time, for reasons similar to those seen in the small adjoint mass limit (4.13, 4.14); see Appendix B.

Discussion
The absence of exponential growth in the out-of-time-order four-point function implies that the IOP model is not chaotic. In fact, there is a heuristic way to understand the absence of chaos in the IOP model. The interacting part of the Hamiltonian (2.4) can be written as, As a result of the absence of self-interactions for the adjoints, combined with the assumption of large fundamental mass M T , the number of fundamentals is timeindependent and, a i (t) = e −ihQ il t a l (0) .

(5.2)
Since Q is a hermitian matrix, it has real eigenvalues, and so the norm of the a i operators does not grow. If we relax the assumption that M T , the above argument is no longer applicable, though this may not be sufficient to make the model chaotic. Heuristically, chaos is associated with rapid growth. As we evolve a fundamental, it is emitting and absorbing adjoints. Since the adjoints have no self-interaction, and conversion of an adjoint into two fundamentals is suppressed by 1/N , the only way for the adjoints to continue evolving in between emissions and absorptions is if they interact with fundamentals in the thermal bath.
It may be useful to modify the IOP model, so as to have several flavors of fundamentals. Also, the interaction (5.1) can written in terms of the quadratic Casimirs, −hq · Q = 1 2 hTr(q 2 + Q 2 − (q + Q) 2 ), allowing a computation of the two-point function at finite N through a sum over Young tableaux [3]. One could study the four-point function in this way as well.

A Steepest descent
In this appendix, we review some aspects of evaluating integrals by the method of steepest descent, see e.g. [33]. Consider an integral of the form, where the integral is evaluated along some contour. For now, let g(z), f (z) be smooth functions; we will discuss later how to relax this assumption. Since t 1, the integrand generically undergoes rapid oscillations which cancel out. The idea will be to deform the contour of integration so as to follow a path for which the phase remains constant. As long as we do not cross any singularities, we are free deform the contour. Splitting f (z) into a real and imaginary part, we need to deform the contour to follow a path of constant u(z). The most relevant region of the integrand is one in which the real part is maximized. Letting z = a + ib, As a result of the Cauchy-Riemann equations, this amounts to finding the saddle points, f (z) = 0. Therefore, the prescription for approximating (A.1) is to focus on the vicinity of the dominant saddle point, and choose a direction for the contour that moves away from the saddle point so as to maintain constant phase u(z).
As an example, consider the integral representation of the Bessel function, This has a branch cut, x ∈ (−i∞, −i) ∪ (i, i∞). We perform a change of variables, Extermizing f (u) = sinh u, the saddle points are at u = ±πi/2. The line of constant phase passing through the saddle points is one that runs along the imaginary axis. We deform the contour so that it runs along −∞ < u < −iπ/2. Moving downward from u 0 = −iπ/2 is a direction of steepest descent. In the vicinity of the saddle, Defining a new variable z as u = u 0 − iz, (A.5) becomes, which is the correct large t expansion of K 0 (t).
We have so far discussed approximating (A.1) by the behavior near the saddle point. There are several contexts in which other regions may be relevant. If the contour has endpoints, then one must analyze the behavior near the endpoints. Additionally, if g(z) has singularities, then one must analyze the integrand near those regions as well. In particular, it may happen that there is no way to deform the contour into the relevant steepest descent contour, without passing through singularities. If the singularity of g(z) is a simple pole, then we may simply deform through it, picking up the contribution of the pole. If, instead, g(z) has a branch cut or an essential singularity, we must analyze the integrand in the vicinity of these regions.
For instance, consider again approximating (A.4), but without changing variables. In this case, g(x) = (1 + x 2 ) −1/2 and f (x) = x. The exponential has no saddle points, so we focus on the regions where g(x) is large: near x = ±i. We integrate along a direction running parallel to the imaginary axis, as we still need to maintain constant phase for the exponent. Letting x = −i − ρi, with ρ 1 so that √ 1 + x 2 ≈ √ 2ρ, (A.4) is approximated by, where we have extended the range of integration to infinite ρ, as its contribution is negligible. Evaluating (A.8) gives (A.7).

B Four-point integral
The four-point function for the IOP model is, where ω 4 = ω 1 + ω 2 − ω 3 and G 4 (ω 1 , ω 2 , ω 3 , ω 4 ) is given by (4.27). Our eventual interest is the out-of-time-order four-point with time separations t 41 = 0, t 34 = −t 42 = t and large t. At large t, the exponent in (B.1) undergoes rapid oscillations as ω 2 , ω 3 are varied. Since the exponent clearly has no saddle point, the only regions which could lead the four-point function to grow exponentially are those in which G 4 (ω 1 , ω 2 , ω 3 , ω 4 ) is singular. We thus hold ω 1 fixed, and scan over ω 2 , ω 3 , looking for regions in which the frequency-space four-point function is divergent. The relation between ω 2 and ω 3 where this occurs then determines the form of the exponent in (B.1), which can then be written just as a function of ω 2 . This function may have saddles, which will either lead to an oscillatory exponent or a growing one.