Onset of many-body chaos in the $O(N)$ model

TL;DR

This paper calculates how chaos emerges in the 2+1 dimensional $O(N)$ model at finite temperature, showing exponential growth of operator commutators with a temperature-dependent rate and ballistic operator spreading.

Contribution

It provides the first leading-order calculation of chaos growth rate and butterfly velocity in the $O(N)$ model using $1/N$ expansion at finite temperature.

Findings

Chaos growth rate $oxed{ ext{~} rac{3.2 T}{N}}$

Butterfly velocity $v_B/c oxed{ ext{~} 1}$

Growth of chaos is slow at large $N$ due to weak interactions.

Abstract

The growth of commutators of initially commuting local operators diagnoses the onset of chaos in quantum many-body systems. We compute such commutators of local field operators with $N$ components in the $(2+1)$-dimensional $O(N)$ nonlinear sigma model to leading order in $1/N$. The system is taken to be in thermal equilibrium at a temperature $T$ above the zero temperature quantum critical point separating the symmetry broken and unbroken phases. The commutator grows exponentially in time with a rate denoted $\lambda_L$. At large $N$ the growth of chaos as measured by $\lambda_L$ is slow because the model is weakly interacting, and we find $\lambda_L \approx 3.2 T/N$. The scaling with temperature is dictated by conformal invariance of the underlying quantum critical point. We also show that operators grow ballistically in space with a "butterfly velocity" given by $v_B/c \approx 1$ where $c$ is the Lorentz-invariant speed of particle excitations in the system. We briefly comment on the behavior of $\lambda_L$ and $v_B$ in the neighboring symmetry broken and unbroken phases.

Figures

Equations

$$C(t)=\mathrm{Tr}\left\{\rho [W(t),V{]}^{†}[W(t),V]\right\}=-\mathrm{Tr}\left\{\rho [W(t),V{]}^2\right\}$$

$$C(t)\sim ϵe^{{\lambda }_{L}t}$$

$$\mathcal{C}(t)=\mathrm{Tr}\left\{{\rho }^{1/2}[W(t),V{]}^{†}{\rho }^{1/2}[W(t),V]\right\}=-\mathrm{Tr}\left\{{\rho }^{1/2}[W(t),V]{\rho }^{1/2}[W(t),V]\right\}.$$

$$\mathcal{C}(t)$$

$$G(t)$$

$$F(t)$$

$$\displaystyle{\cal{L}}=\bigg{[}\frac{1}{2}(\partial\varphi_{a})^{2}-\frac{v}{2N}\bigg{(}\varphi_{a}^{2}-\frac{N}{g}\bigg{)}^{2}\bigg{]}.$$

$$\displaystyle\tilde{Z}=\int{\cal{D}}\varphi_{a}~{}\textnormal{exp}\bigg{(}i\int_{x}{\cal{L}}\bigg{)},$$

$$\mathcal{C}(t,\mathbf{x})=-\frac{1}{N^2}\sum _{a,b}\mathrm{Tr}\left\{{\rho }^{1/2}[{\phi }_a(\mathbf{x},t),{\phi }_b]{\rho }^{1/2}[{\phi }_a(\mathbf{x},t),{\phi }_b]\right\},$$

$${\lambda }_L\approx 3.2\frac{T}{N}.$$

$$\lambda (\mathbf{k})={\lambda }_L-D_L|\mathbf{k}{|}^2+...$$

$$\mathcal{C}(t,\mathbf{x})\sim {\int }_{\mathbf{k}}{\chi }_{\mathbf{k}}{e}^{i\mathbf{k}\cdot \mathbf{x}+{\lambda }_L(\mathbf{k})t}$$

$$\mathcal{C}(t,\mathbf{x})\sim \exp \left({\lambda }_{L}t-\frac{|\mathbf{x}{|}^2}{4D_{L}t}\right)$$

$$ϵ\sim \exp \left(-\frac{|\mathbf{x}{|}^2}{4D_{L}t}\right),$$

$${\partial }_t\mathcal{C}(t,\mathbf{x})={\lambda }_L\mathcal{C}(t,\mathbf{x})+D_L{\nabla }^2\mathcal{C}(t,\mathbf{x})+...$$

$$v_B=\sqrt{4D_L{\lambda }_L}.$$

$$\frac{v_B}{c}\approx 1.$$

$$\frac{1}{\mathcal{C}}\frac{d\mathcal{C}}{dt}={\lambda }_L+\frac{|\mathbf{x}{|}^2}{4D_L{t}^2},$$

$$\displaystyle\mathcal{C}(t,\boldsymbol{x})\big{|}_{|t-|\boldsymbol{x}|/v_{B}|\ll t}\sim\exp\left(2\lambda_{L}[t-|\boldsymbol{x}|/v_{B}]\right).$$

$$\displaystyle{\cal{L}}\rightarrow\overline{{\cal{L}}}=\frac{1}{2}\bigg{[}(\partial\varphi_{a})^{2}+\frac{\lambda}{\sqrt{N}}\bigg{(}\varphi_{a}^{2}-\frac{N}{g}\bigg{)}+\frac{\lambda^{2}}{4v}\bigg{]},$$

$$\tilde{Z}=\int \mathcal{D}{\phi }_a\mathcal{D}\lambda \text{exp}\left(i{\int }_x\overline{\mathcal{L}}\right).$$

$$\displaystyle{\cal{L}}_{E}=\frac{1}{2}\bigg{[}(\partial\varphi_{a})^{2}-\frac{\lambda}{\sqrt{N}}\bigg{(}\varphi_{a}^{2}-\frac{N}{g}\bigg{)}-\frac{\lambda^{2}}{4v}\bigg{]}$$

$$Z=\int \mathcal{D}{\phi }_a\mathcal{D}\lambda \text{exp}\left(-\int d\tau d^2\mathbf{x}{\mathcal{L}}_E\right).$$

$$\mathcal{G}(\tau ,\mathbf{x}){\delta }_{a,b}=⟨{\phi }_a(\mathbf{x},\tau ){\phi }_b(0,0)⟩.$$

$$\mathcal{G}(i{\omega }_n,\mathbf{k})$$

$${ϵ}_{\mathbf{k}}^2$$

$${\mathcal{G}}_R(t,\mathbf{x}){\delta }_{a,b}=-i⟨[{\phi }_a(\mathbf{x},t),{\phi }_b(0,0)]⟩\theta (t)$$

$${\mathcal{G}}_R(\omega ,\mathbf{k})=-\mathcal{G}(i{\omega }_n\to \omega +i{0}^{+},\mathbf{k})=\frac{1}{(\omega +i{0}^{+}{)}^2-{ϵ}_{\mathbf{k}}^2},$$

$$A(\omega ,\mathbf{k})=\frac{\pi }{{ϵ}_{\mathbf{k}}}[\delta (\omega -{ϵ}_{\mathbf{k}})-\delta (\omega +{ϵ}_{\mathbf{k}})].$$

$${\mathcal{G}}_R(\omega ,\mathbf{k})=\int \frac{d\nu }{2\pi }\frac{A(\nu ,\mathbf{k})}{\omega +i{0}^{+}-\nu }.$$

$${\mathcal{C}}_0(t,\mathbf{x})={\left[{\mathcal{G}}_R(t,\mathbf{x})\right]}^2.$$

$${\mathcal{G}}_W(t,\mathbf{x}){\delta }_{a,b}=\text{Tr}\left\{{\rho }^{1/2}{\phi }_a(\mathbf{x},t){\rho }^{1/2}{\phi }_b(0)\right\}.$$

Full text

1 Introduction

Understanding thermalization in closed quantum systems is one of the great challenges of quantum many-body physics [1, 2, 3, 4]. Issues of thermalization have recently received renewed attention, due in part to numerous experiments probing quantum dynamics in approximately closed systems [5] [6] [7, 8] [9, 10] [11] [12] as well as new theoretical tools, like AdS/CFT duality, which are able to address dynamics in isolated strongly interaction systems (see Refs. [13],[14] for a review). In general, there are numerous time scales associated with thermalization, which proceeds through distinct stages. In the simplest scenario, one can identify three distinct stages of thermalization.

At early times there is a process of relaxation, which describes the initial decay of local perturbations. The onset of relaxation can be diagnosed from the time dependence of simple auto-correlation functions of local observables [15]. The time-scale for relaxation is typically independent of the system size. Relaxation is followed by a process called scrambling at intermediate times [16], which describes the spreading of quantum entanglement and information across all of the degrees of freedom of the system. Accessing scrambled information necessarily requires making non-local measurements; the memory of the initial out-of-equilibrium state is effectively lost as far as local probes are concerned [17]. For local Hamiltonians the time-scale to reach a scrambled state necessarily depends on the total number of degrees of freedom. Relaxation and scrambling are distinct phenomena, though relaxation processes and time-scales will influence how fast the system scrambles. After the system is effectively scrambled, the quantum state continues to move in Hilbert space and can continue to grow in complexity up to extremely long times [18].

Scrambling is fundamentally an information theoretic notion [19][20], but the growth of scrambling can also be probed using special correlation functions as we discuss below [21, 22]. The relevant correlation functions can be related to information theoretic measures of scrambling, but are easier to access, both theoretically and experimentally. Furthermore, in systems with a semi-classical limit, the growth of scrambling, as measured by these special correlation functions, can be heuristically related to the growth of chaos in the corresponding classical model as measured by classical Lyapunov exponents [22]. However, there are some important subtleties with this connection and it remains incompletely understood [23]. Because the growth of scrambling describes the effective loss of memory of the initial state and because of its relation to growth of classical chaos, the onset of scrambling can be regarded as a fully quantum avatar of the growth of chaos.

To define the physical correlation functions of interest, consider two local Hermitian operators $W$ and $V$ and a local Hamiltonian $H$. From these objects we can form the time evolution operator $U(t)=e^{-iHt}$ and the Heisenberg operator $W(t)=U(-t)WU(t)$. We take the state of the system to be in thermal equilibrium at temperature $T=1/\beta$: $\rho\propto e^{-H/T}$. The “unregulated” squared commutator in state $\rho$ is

[TABLE]

where in the second equality we used $[W(t),V]^{\dagger}=-[W(t),V]$ valid for Hermitian $W$ and $V$. This commutator is called “unregulated” because some of the operators entering $C(t)$ are inserted at the same spacetime point, a situation which typically leads to additional short-distance divergences in a quantum field theory. The expectation in a quantum chaotic system is that $C(t)$ begins small if $W$ and $V$ are well separated and that it subsequently grows exponentially in time,

[TABLE]

for some small parameter $\epsilon$ with an exponent denoted $\lambda_{L}$. The small parameter $\epsilon$, which in general may also depend on time in some way, is typically related to some measure of the number of relevant degrees of freedom, e.g. the number of components of a field or the distance between $W$ and $V$.

For the calculations in this paper, it is more convenient to consider a “regulated” squared commutator defined as

[TABLE]

This object has the virtue that all operator insertions occur at distinct spacetime points thus removing short-distance divergences from coincident operators. The regulation of divergences is achieved by moving two of the operators half way around the thermal circle which is equivalent to an analytic continuation off the real time axis.111To be clear, consider the four point function ${\rm Tr}\left\{\rho O_{1}(t_{1})O_{2}(t_{2})O_{3}(t_{3})O_{4}(t_{4})\right\}$. We can move the $j$th operator along the thermal circle by an amount $\tau$ by shifting $t_{j}\rightarrow t_{j}-i\tau$ or $O_{j}(t_{j})\rightarrow\rho^{-\tau/\beta}O_{j}(t_{j})\rho^{\tau/\beta}$. Since the thermal circle is periodic with period $\beta=1/T$, $\tau=\beta/2$ corresponds to half way around the thermal circle. Shifting two operators, say $1$ and $2$, half way around amounts to considering ${\rm Tr}\left\{\rho^{1/2}O_{1}(t_{1})O_{2}(t_{2})\rho^{1/2}O_{3}(t_{3})O_{4}(t_{4})\right\}$. We furthermore expect $\mathcal{C}(t)$ to grow exponentially in time at the same rate as $C(t)$. This is because $\mathcal{C}(t)$ is strictly positive, $\mathcal{C}={\rm Tr}\left\{A^{\dagger}A\right\}$ where $A=\rho^{1/4}[W(t),V]\rho^{1/4}$, so analytically continuing $t$ in $e^{\lambda_{L}t}$ still yields exponential growth given the guarantee that $\mathcal{C}$ is positive.

It is also conventional to expand the commutators appearing in $\mathcal{C}(t)$ or $C(t)$ to give two terms:

[TABLE]

Note that $G$ is positive and time-ordered. On the other hand, $F$ is called an out-of-time-order (OTO) correlator because of the unusual time order of its operator insertions. The placement of the operators in $F$ is shown in Figure 1. $F$ is approximately equal to $G$ at early times, when $W(t)$ and $V$ still approximately commute, but differs from $G$ at later times. If $G(t)$ reaches its equilibrium value after some local relaxation time $t_{r}$, then the subsequent dynamics of $\mathcal{C}(t)$ will be controlled by $F(t)$.

OTO correlators showed up decades ago in the study of semi-classical methods in superconductivity [24]. More recently, they have been studied in much greater detail in the context of the AdS/CFT correspondence where they have been shown to diagnose quantum chaos in black hole physics [21] [25][22]. In particular, black holes are supposed to be the fastest scramblers in nature [26]. At the other extreme, disordered quantum systems typically scramble much more slowly than their clean counterparts [27] [28] [29] [30] [31] [32]. Several experimental proposals have also appeared recently [33] [34] [35] that enable one to measure OTO correlators and scrambling and three preliminary experiments have already been carried out [36] [37] [38]. There has been a flurry of recent calculations of out-of-time-order correlators in a variety of models [39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51]. There have also been some recent works exploring other formal aspects of these correlators, including fluctuation-dissipation-like theorems [52][53][54][55].

In this paper we study scrambling in a class of quantum critical field theories, namely vector models associated with an $N$ component real scalar $\varphi_{a}$ field in (2+1)-dimensions with $O(N)$ symmetry. These theories describe a zero-temperature quantum phase transition from a disordered phase where $\langle\varphi_{a}\rangle=0$ to an ordered symmetry breaking phase where $\langle\varphi_{a}\rangle\neq 0$. A large class of problems at long wavelengths are described by this model, e.g. the theory with $N=2$ describes the transition between a superfluid and a bosonic Mott insulator [56] and the theory with $N=3$ describes the transition from a paramagnet to a Heisenberg antiferromagnet [57] [58] [59]. While the system can be described by well-defined quasiparticles on either side of the transition, the critical point itself is distinguished by the absence of any quasiparticle-like excitations. Moreover, the quantum many-body system at a non-zero temperature ‘above’ the quantum critical point (QCP) [60] is described by only one relevant scale, namely the temperature itself. It is then natural to ask how fast this broad class of theories scrambles at a given temperature.

Before stating our results, we will explicitly describe the model. We study the theory of an $O(N)$ symmetric real vector field, $\varphi_{a}$, $a=1,..,N$, governed by the real time Lagrangian

[TABLE]

The real time generating functional is

[TABLE]

where $\int_{x}\equiv\int dt~{}d^{2}\boldsymbol{x}$ is the integral over $(2+1)-$ dimensional space-time222We reserve $\tau~{}(i\omega_{n})$ for Euclidean-imaginary time (Matsubara frequency) and $t~{}(\omega)$ for real time (frequency).. The parameter $g$ is used to tune across the QCP (which occurs at $g=g_{c}$) and $v>0$ is a self-interaction coupling constant. We have set the “speed of light” $c=1$.

We study the index-averaged commutator of the fields, $\varphi_{a}(\boldsymbol{x},t)$,

[TABLE]

where $\varphi_{b}=\varphi_{b}({\boldsymbol{0}},0)$ and $\rho=e^{-\beta H}/Z$ is the thermal density-matrix. This corresponds to the thermally regulated commutator discussed above with $W=\varphi_{a}(\boldsymbol{x})$ and $V=\varphi_{b}({\boldsymbol{0}})$. The normalization in Eq. (1.9) amounts to averaging the square commutator over $a,b=1,...,N$.

Our main results are summarized below. All results are to leading non-trivial order in a $1/N$ expansion at the quantum critical point at temperature $T$. We first compute the chaos exponent $\lambda_{L}$ in the spatially averaged case, i.e. the growth rate of $\int d^{2}\boldsymbol{x}~{}\mathcal{C}(t,\boldsymbol{x})$, and find

[TABLE]

As discussed below, our calculation of the prefactor is approximate; we roughly estimate the accuracy to be about $10$%.

We next analyze the commutator, $\mathcal{C}(t,\boldsymbol{k})=\int d^{2}\boldsymbol{x}~{}e^{-i{\boldsymbol{k}}\cdot\boldsymbol{x}}\mathcal{C}(t,\boldsymbol{x})$, at a fixed external momentum, $\boldsymbol{k}$. This grows exponentially in time as well with a $\boldsymbol{k}$ dependent growth exponent given by

[TABLE]

The spacetime structure of the commutator is determined from the Fourier transform

[TABLE]

which we evaluate in a saddle point approximation (assuming no singular structure in $\chi_{\boldsymbol{k}}$). In holographic and related studies [61, 62], there is an additional pole-like singularity in $\chi_{\boldsymbol{k}}$ which, when closer to the real-axis, determines the spatial decay. As we discuss later, we see no evidence for such a pole in our numerical computations and speculate that it may be related to the presence of gravitational modes.

Within the above approximation, the space-time dependence of $\mathcal{C}(t,\boldsymbol{x})$ then takes the form

[TABLE]

where $\lambda_{L}$ is the uniform growth exponent and $D_{L}$ has units of a diffusion constant. Comparing Eq. (1.13) to Eq. (1.2) shows that the role of $\epsilon$ is played by

[TABLE]

which for large $|\boldsymbol{x}|$ indeed represents a suppression factor in the square commutator.

From Eq. (1.13) it follows that $\mathcal{C}(t,\boldsymbol{x})$ obeys a local equation of motion of the form

[TABLE]

where $...$ denotes higher order terms in a gradient expansion and non-linear terms. In practice, we directly compute the momentum dependent chaos exponent and only obtain the spacetime structure indirectly.

The form in Eq. (1.13) is consistent with a ballistic growth of operators with a “butterfly velocity” [63]

[TABLE]

We find that $v_{B}$ is independent of $T$ and $N$ at large $N$; expressed as a ratio with the speed of light, our approximate numerical calculation gives

[TABLE]

As discussed below, our calculation of the ratio is approximate; we roughly estimate the accuracy to be about $10$%. Of course, causality constrains $v_{B}/c$ to be less than one. The butterfly velocity has been argued [42] to play the role of a low energy state dependent Lieb-Robinson velocity [64].

Interestingly, the derivative of the diffusive-plus-growth form in Eq. (1.13) is

[TABLE]

which for sufficiently large $|\boldsymbol{x}|/t$ is greater than the chaos bound of $2\pi T$ [26]. Of course, the chaos bound does not directly address the dynamics of $\mathcal{C}$, but the bound does suggest that the diffusive form must break down for sufficiently large $|\boldsymbol{x}|/t$. For example, higher order terms in the ${\boldsymbol{k}}$ expansion of $\lambda({\boldsymbol{k}})$ will become important as the saddle point momentum increases away from the real axis. Additional singularities, e.g. poles and branch cuts, and the exact vanishing of the commutator outside the lightcone may also become relevant before the chaos bound is violated.

There is another feature of note concerning the growth of chaos in space. Expanding Eq. (1.13) near the wavefront $|\boldsymbol{x}|=v_{B}t$ gives

[TABLE]

Note the peculiar factor of $2$ in front of $\lambda_{L}$ which arises from the expansion. The full behavior of the squared commutator near the front is not fully understood at this point and is beyond the scope of this work. Moreover, the nature of the crossover of the squared commutator outside of the causal region will be the topic of future work.

We also estimate $\lambda_{L}$ and $v_{B}$ in the proximate phases. In the symmetry unbroken phase at low temperature we find $\lambda_{L}\sim e^{-2\mu/T}$ ($\mu$ is the mass of the $\varphi$ particles) and $v_{B}\sim\sqrt{T/\mu}$. In the symmetry broken phase at low temperature we find $\lambda_{L}\sim T^{3}/\rho_{s}^{2}$ and $v_{B}/c\sim 1$. A word of caution is appropriate here: in two spatial dimensions at finite $N$ and non-zero temperature the symmetry broken phase is disordered and the symmetry is restored. Our estimates are based on the existence of a Goldstone mode which will not persist at the longest distances at finite temperature. We discuss this subtlety more carefully in Section 6. Our results for the quantum critical region and our estimates for the proximate phases are summarized in Figure 2; in general, there is a complex crossover function which connects the three regimes discussed above.

The rest of this paper is organized as follows: in Section 2, we set the stage for carrying out a $1/N$ expansion for the $O(N)$ theory described by Eq. (1.7). Section 3 develops the perturbative expansion used to compute the regulated commutator $\mathcal{C}(t,\boldsymbol{x})$. Sections 4 and 5 discuss the numerical solution of the ladder series derived in Section 3 to obtain the growth exponent and the butterfly velocity. We briefly discuss scrambling in the symmetry broken and unbroken phases in Section 6, leaving a detailed computation for future work. We end with a discussion where we highlight some interesting future directions. The appendices contain a variety of technical details related to our calculations.

Notation guide: In the remainder of the paper, $N$ denotes the number of components of the field $\varphi_{a}$ and $\lambda$ denotes the Lagrange multiplier in the $O(N)$ action. The imaginary time $\varphi$ correlator is ${\cal G}$, the retarded $\varphi$ correlator is ${\cal G}_{R}$, and the Wightman $\varphi$ correlator is ${\cal G}_{W}$. Similarly, the imaginary time $\lambda$ correlator is ${\cal G}_{\lambda}$, the retarded $\lambda$ correlator is ${\cal G}_{R,\lambda}$, and the Wightman $\lambda$ correlator is ${\cal G}_{W,\lambda}$. The thermal mass of the $\varphi$ field is denoted $\mu$ and $\Lambda$ is a momentum cutoff where needed. $\lambda_{L}$ is the chaos exponent and $v_{B}$ is the butterfly velocity. We distinguish position and time and momentum and frequency with spatial vectors denoted by bold face, e.g. $\boldsymbol{x}$ and $\boldsymbol{k}$. Integrals over momentum are denoted $\int_{\boldsymbol{k}}$ and include the $(2\pi)^{-2}$ phase space factor in $d=2$ spatial dimensions. The thermal state is $\rho\propto e^{-\beta H}$ where $H$ is the Hamiltonian and $\beta=1/T$ is the inverse temperature. We set the speed of light $c=1$, unless explicitly stated.

2 Preliminaries

In this and the next section, we set up the formalism required to compute $\mathcal{C}(t)$ in a $1/N$ expansion in the theory Eq. (1.7). We are particularly interested in the limit of strong coupling $v\rightarrow\infty$ at the quantum critical point where Eq. (1.7) describes a $2+1$-dimensional conformal field theory [60].

Decoupling the quartic term in Eq. (1.7) using a Hubbard-Stratonovich field, $\lambda(\boldsymbol{x},t)$, leads to a new Lagrangian

[TABLE]

such that the generating function becomes

[TABLE]

Although we are ultimately interested in real time evolution, it is also quite useful to consider the imaginary time formulation of the above theory. As usual we set $it\rightarrow\tau$ and construct the Euclidean Lagrangian

[TABLE]

and the partition function

[TABLE]

The derivative term now uses the Euclidean $+++$ metric instead of the $+--$ metric in real time.333Note that for comparison with some of the literature, the field $\lambda$ is sometimes traded for the field $\tilde{\lambda}=i\lambda$. We will work with $\lambda$ consistently throughout because it is naturally regarded as Hermitian.

To leading order in $N$, the effect of $\lambda$ is simply to generate an effective mass for the particles described by the $\varphi$ field: $-\langle\lambda\rangle/\sqrt{N}=\mu^{2}$. The imaginary time $\varphi_{a}$ propagator is

[TABLE]

In momentum space ${\cal G}$ has the simple form,

[TABLE]

It is important to note that at the critical point $g=g_{c}$, the mass term $\mu$ for $\varphi_{a}$ is proportional to the temperature since there are no other scales. The mass, $\mu=\Theta T$, where $\Theta=2\log\frac{1+\sqrt{5}}{2}$; the standard computation is reviewed in Appendix A.

The retarded propagator is defined in the usual way as

[TABLE]

where $\theta(t)$ is the Heaviside step function. ${\cal G}_{R}$ is related to ${\cal G}$ by analytic continuation,

[TABLE]

and the spectral function $A(\omega,{\boldsymbol{k}})=-2\,\text{Im}[{\cal G}_{R}(\omega,{\boldsymbol{k}})]$ is

[TABLE]

To fix notation, recall that the inverse formula is

[TABLE]

At infinite $N$, this is the whole story. The $\varphi$ particles are massive but non-interacting and chaos does not develop. In particular, the squared commutator is

[TABLE]

This implies that the chaos exponent is of order $1/N$ at large $N$. Hence we must study $1/N$ corrections in the theory Eq. (2.1).

In order to evaluate the squared commutator to higher orders in $1/N$, we need additional real-time propagators. Thus we also define the symmetrized Wightman function

[TABLE]

As reviewed in Appendix C, this function may be written in terms of the corresponding $\varphi$ spectral function as

[TABLE]

For our subsequent computations involving $1/N$ corrections, we also need the propagator for the $\lambda$ field. The bare Euclidean propagator is ${\cal G}^{(0)}_{\lambda}(i\omega_{n},{\boldsymbol{k}})=-4v$. At infinite $v$, we may re-sum and dress the $\lambda$ propagator as shown in Figure 3. Then,

[TABLE]

where $\Pi(i\omega_{n},{\boldsymbol{k}})$ is the one-loop $\varphi_{a}$ bubble,

[TABLE]

We evaluate the above diagram in Appendix B. We define the retarded bubble $\Pi_{R}(\omega,{\boldsymbol{q}})=\Pi(i\omega_{n}\rightarrow\omega+i0^{+},{\boldsymbol{q}})$. The real and imaginary parts of the retarded bubble $\Pi_{R}(\omega,q)$ are evaluated numerically as discussed in Appendix G and serve as inputs to the later calculation of $\lambda_{L}$ and $v_{B}$.

From the basic Euclidean $\lambda$ propagator we extract the real time correlators of the $\lambda$ field. First, the retarded $\lambda$ correlator is

[TABLE]

The $\lambda$ spectral function is then

[TABLE]

The imaginary part of the retarded bubble is positive for positive frequency, so the $\lambda$ spectral function is also positive for positive frequency. From this spectral function we immediately obtain the symmetrized Wightman function for $\lambda$,

[TABLE]

The $\lambda$ Wightman function is also positive.

Finally, in order to carry out the computation consistently to leading order in $1/N$, we need to include the leading $1/N$ correction to the self-energy for the $\varphi_{a}$ propagator (Figure 4). This is given by,

[TABLE]

The two contributions in the first line above correspond to the two diagrams in Figure 4. The self-energy modifies the $\varphi$ propagator,

[TABLE]

As with the bubble, we define the retarded self-energy by $\Sigma_{R}(\omega,{\boldsymbol{q}})=\Sigma(i\omega_{n}\rightarrow\omega+i0^{+},{\boldsymbol{q}})$. Converting the self-energy modified Euclidean propagator to the corresponding retarded correlator yields

[TABLE]

For our purposes, the important physical effect of $\Sigma_{R}$ is to introduce a finite lifetime to the $\varphi$ particles. The poles of ${\cal G}_{R}$ are shifted by $\Sigma_{R}$, so working to first order in $\Sigma_{R}$ gives poles at

[TABLE]

Since the real part of $\Sigma_{R}$ is symmetric and the imaginary part is anti-symmetric, we find the poles

[TABLE]

This suggests defining the inverse lifetime

[TABLE]

so that $\Gamma_{\boldsymbol{q}}>0$ corresponds to decay and ${\cal G}_{R}$ remains analytic in the $\omega$ upper half-plane.

We have evaluated the imaginary part of the self-energy in Appendix D; the final expression is given by,

[TABLE]

As a check on the computation, notice that positivity of the $\lambda$ Wightman function implies positivity of the decay rate $\Gamma_{\boldsymbol{q}}$.

The inverse lifetime is

[TABLE]

where the function ${\cal{R}}_{1}({\boldsymbol{k}},{\boldsymbol{q}})$ is given by,

[TABLE]

The significance of ${\cal{R}}_{1}$ will become clear shortly.

We can define an inverse “(on-shell) phase-coherence time”,

[TABLE]

$\tau_{\varphi}$ is associated with the loss of phase coherence and can be deduced from a two-point function. It is the analog of the relaxation time discussed in the introduction. We find

[TABLE]

This is a high precision value good to three digits444We thank W. Witczak-Krempa and S. Sachdev for helping benchmark the numerical computation of $\tau_{\phi}$ [65]. which serves as a check on our numerical calculation of $\Pi$ and as a benchmark for the more complex calculations of $\lambda_{L}$ and $v_{B}$ since these calculations also output approximate values for $1/\tau_{\varphi}$.

3 Diagrammatic expansion for the squared commutator

We now assemble the pieces discussed in the previous sections to compute $\mathcal{C}(t)$ in theory Eq. (2.1) to leading non-trivial order in $1/N$. Recall that as shown in Figure 1 the squared commutator is naturally associated with a complex time contour which includes the imaginary time thermal circle and two “real time folds”. The perturbation expansion is naturally organized in terms of retarded correlation functions within the time folds and Wightman correlation functions between the two time folds.

To generate the perturbation series, we expand $\mathcal{C}$ in powers of the interaction vertex on both time folds. As far as possible, we keep only diagrams which affect the result at order $1/N$. However, we re-sum a number of contributions, so we are not performing a strict $1/N$ expansion since diagrams with arbitrarily high powers of $1/N$ are included. Some additional classes of diagrams which are not $1/N$ suppressed are nevertheless ignored, as discussed below, because they do not affect the leading $1/N$ growth rate.

The rules of our calculation are discussed in detail in Appendix E. They are summarized as follows:

Vertex insertions are restricted to lie along the real time folds. Each real time vertex insertion is associated with a factor of $\frac{i}{2\sqrt{N}}$. Vertex insertions along the thermal circle dress the thermal state but do not directly lead to real time growth; the relevant effects are included in the various finite temperature correlation functions obtained via analytic continuation from imaginary time. Mixed insertions involving contractions between operators on the thermal circle and operators on a time fold also do not grow in time. In fact, they decay in time as the operator on the time fold moves away from the thermal circle. 2. 2.

Horizontal lines correspond to retarded propagators, $i{\cal G}_{R}$ or $i{\cal G}_{R,\lambda}$, and vertical lines correspond to Wightman propagators, ${\cal G}_{W}$ or ${\cal G}_{W,\lambda}$. The direction of the lines is not meaningful in Euclidean diagrams. The Euclidean $\lambda$ propagator is given by $-1/\Pi$ to leading order. 3. 3.

The fastest growing diagrams correspond to a set of ladder diagrams with two types of rungs as shown in Figure 5. The type-I rung (dashed line) corresponds to ${\cal G}_{W,\lambda}$ while the type-II rung (wavy line) corresponds to the box insertion which is schematically $\sim{\cal G}_{R,\lambda}^{2}{\cal G}_{W}^{2}$. When considering the averaged squared commutator, the sum over $a$ and $b$ effectively caps off the ladder as shown in Figure 5. 4. 4.

Accounting for the appropriate closed loops of $\varphi$ which appear when considering the type-II rung, both rungs are associated with a factor of $1/N$. This leads to a factor of $1/N^{\ell}$ for a diagram with $\ell$ rungs. Hence $\lambda_{L}\sim 1/N$. 5. 5.

In addition to the explicit rungs, the $\varphi$ sides of the ladder are dressed with self-energy corrections. The self-energy is computed to order $1/N$. 6. 6.

Vertex corrections and other higher order diagrams are ignored. Crossed ladder diagrams are of the same order in $1/N$ as the ladder diagrams we keep, but we argue in Appendix F that these diagrams do not contribute to the leading order growth rate. Roughly speaking, this is because they are far from being on-shell.

A few comments are in order. The focus on ladder diagrams is the standard perturbative approach to the calculation of chaos exponents [39, 40, 43, 41]. Our calculation extends a recent perturbative calculation in a $\varphi^{4}$ matrix-model in which the fields $\varphi$ are $N\times N$ matrices [43]. There are some important differences between the vector and matrix models. The growth exponent is $1/N$ suppressed in the vector model but not the matrix model due to differences in the $N$ counting. The crossed ladder diagrams are $1/N$ suppressed in the matrix model because they are non-planar. The type-II rungs are also suppressed due to the weak-coupling treatment of the matrix model. There are additional re-summations in the vector model, e.g. to compute the renormalized $\lambda$ propagator, and we are able to access a conformal fixed point in the vector model.

The explicit $N$ counting is straightforward. A ladder with $\ell_{1}$ rungs of type-I and $\ell_{2}$ rungs of type-II, denoted $D_{\ell_{1},\ell_{2}}$, is associated with a factor of

[TABLE]

Hence the $N$ suppression of a ladder depends only on the total number of rungs of either type, $\ell=\ell_{1}+\ell_{2}$. Summing over all ladder diagrams then gives a square commutator of order

[TABLE]

We also see that the growth rate will have an explicit factor of $1/N$. The scrambling time, $t_{*}$, defined as $\mathcal{C}(t_{*})\sim 1$ will thus be of order

[TABLE]

This is consistent with earlier arguments [66].

We are now ready to evaluate the ladder sum for $\mathcal{C}$. We start by writing down a Bethe-Saltpeter equation for $\mathcal{C}(t,\boldsymbol{x})$ (shown in Figure 6). The first diagram on the right (without any rungs) is just given by,

[TABLE]

Chaotic growth comes from the effects of interactions which are captured by the other two diagrams in Figure 6 which correspond to the two types of rungs discussed above.

We first review some features associated with the type-I rung diagram from the top row in Figure 6:

[TABLE]

The overall prefactor is positive since the two factors of $i$ from the vertex cancel the overall minus sign in the definition and the four factors of $i$ associated with the retarded correlators give unity. The product of the retarded propagators, ${\cal G}_{R}(\nu-\omega,-\boldsymbol{p})~{}{\cal G}_{R}(\omega,\boldsymbol{p})$, starting with their free-particle forms, is given by

[TABLE]

Evaluating the $\omega$ integral using the method of residues gives terms proportional to: $(\nu-2\epsilon_{\boldsymbol{p}})^{-1},~{}\nu^{-1},~{}(\nu+2\epsilon_{\boldsymbol{p}})^{-1}$.

In order to extract the leading contribution to the long-time form of ${\cal{C}}(t)$, we focus on the most singular $\nu^{-1}$ term. Then we can approximate,

[TABLE]

The overall positive sign arises from performing the $\omega$ integral by closing it with a clockwise contour in the lower half-plane combined with an explicit minus sign in the residue. At this point, restoring the lifetime from the self-energy correction, we obtain,

[TABLE]

This combination can be recast as

[TABLE]

Let us now proceed by summing the ladder series on the first line of Figure 6 to all orders (i.e. ignoring the diagrams on the second line). We start by defining a function,

[TABLE]

Then, the Bethe-Saltpeter equation becomes,

[TABLE]

To estimate the behavior at large time, we drop the homogenous term; assuming that $f(t)$ grows exponentially, the ladder sum should remain invariant under the addition of an extra rung to the series. Therefore, substituting the form for the propagators,

[TABLE]

It is natural to postulate the following on-shell form of the function $f$:

[TABLE]

The integral equation for the ladder sum then becomes,

[TABLE]

Finally, plugging in the explicit form of $\Gamma_{\boldsymbol{p}}$ from Eq. (2.28), we obtain the following integral equation,

[TABLE]

Having understood the general structure of the ladder sum, let us now include the type-II rung. Just as in the previous computation, we begin by considering the one-block diagram:

[TABLE]

Note that the overall sign is again positive. It is easiest to proceed by defining,

[TABLE]

so that the one-block diagram becomes,

[TABLE]

The one-block and one-rung diagrams can be grouped into one contribution by shifting ${\cal G}_{W,\lambda}(\omega^{\prime}-\omega,\boldsymbol{p}^{\prime}-\boldsymbol{p})\rightarrow{\cal G}_{W,\lambda}(\omega^{\prime}-\omega,\boldsymbol{p}^{\prime}-\boldsymbol{p})+{\cal G}_{\textnormal{eff}}(\omega^{\prime},\omega,\boldsymbol{p}^{\prime},\boldsymbol{p})$; therefore we can simply proceed by setting up the ladder sum as before after shifting the Wightman function as above. Note that we do not assume ${\cal G}_{\textnormal{eff}}$ to be a function of only the relative frequency or momentum.

The next step is to evaluate ${\cal G}_{\text{eff}}$. Recall Eq. (2.14) which is ${\cal G}_{W}(\omega,{\boldsymbol{k}})={\cal{Q}}(\omega)~{}A(\omega,{\boldsymbol{k}})$, where ${\cal{Q}}(\omega)=[2\sinh(\beta\omega/2)]^{-1}$. In this formula it is sufficient to use the bare $\varphi$ spectral function since the rung already contains an explicit factor of $1/N$. Inserting the spectral function gives

[TABLE]

Denote the relative momentum, $\boldsymbol{p}^{\prime}-\boldsymbol{p}=\bar{\boldsymbol{p}}$, and the centre-of-mass momentum, ${\boldsymbol{P}}=(\boldsymbol{p}+\boldsymbol{p}^{\prime})/2$. Also define $\bar{\omega}=\omega^{\prime}-\omega$ and shift $\boldsymbol{p}^{\prime\prime}\rightarrow\boldsymbol{p}^{\prime\prime}+{\boldsymbol{P}}$. We can now use one of the delta functions above to integrate over $\omega^{\prime\prime}$ in Eq. (3.19). This leads to,

[TABLE]

The next step is to set $\nu=0$ in the above integrand; one can ignore the $\nu$ dependence to leading order since $\nu$ is the external frequency and it will ultimately be of order $1/N$. Denoting $\epsilon_{\frac{\bar{\boldsymbol{p}}}{2}\pm\boldsymbol{p}^{\prime\prime}}=\epsilon_{\pm}$, the above can be expressed as,

[TABLE]

The delta function constraints in the above equation are of the form,

[TABLE]

Squaring both sides, we have,

[TABLE]

leading to

[TABLE]

where we have chosen the physical positive root and $\theta$ is the angle between $\boldsymbol{p}^{\prime\prime}$ and $\overline{\boldsymbol{p}}$. We can carry out the integral over the magnitude of $\boldsymbol{p}^{\prime\prime}$ using the delta function. The final integral over the orientation, i.e. the angle $\theta$ between $\boldsymbol{p}^{\prime\prime}$ and $\bar{\boldsymbol{p}}$ must be done numerically. The use of the delta function and the subsequent angular integral parallel the computation of the bubble in Appendix B.

The fully re-summed ladder with both rung types is obtained from the solution of

[TABLE]

Postulating the on-shell form in Eq. (3.13) for $f$, the integral equation for the ladder sum becomes,

[TABLE]

Finally, together with the explicit form of $\Gamma_{\boldsymbol{p}}$ from Eq. (2.28), we obtain the following integral equation,

[TABLE]

where we have factored out the explicit $1/N$ dependence of $\Gamma_{\boldsymbol{p}}$ above.

Some comments are in order. First, note that both ${\cal{R}}_{1}$ and ${\cal{R}}_{2}$ are positive and tend to cause chaos to grow while $\Gamma_{\boldsymbol{p}}$ suppresses the growth. Besides the explicit momenta, the only energy scale appearing in ${\cal{R}}_{1}$ and ${\cal{R}}_{2}$ is the thermal mass $\mu(\propto T)$ and the temperature itself, so in fact the only scale is the temperature. By rescaling all momenta and frequencies by a factor of $T$, it immediately follows that the self consistent ladder equation is

[TABLE]

where $\widehat{\mathcal{M}}$ is a scaled integral kernel in which the temperature has been set to one. Hence without any further analysis it follows that if a growing mode exists, then its growth exponent is proportional to $T/N$.

4 Chaos exponent

By substituting $-i\nu\rightarrow\partial_{t}$, we now view Eq. (3.31) as a matrix differential equation

[TABLE]

where $\widehat{\cal{M}}$ represents the contributions from both rungs and the damping. By further rescaling the frequency, we obtain a dimensionless integral equation which is suitable for numerical analysis. The growth of $f$ and hence $\cal{C}$ will then be determined by the largest positive eigenvalue of $\widehat{\cal{M}}$. All of the individual ingredients in $\widehat{\cal{M}}$, i.e. the rung functions, can be computed numerically starting from the polarization bubble.

We now summarize the numerical procedure used to solve the eigenvalue equation

[TABLE]

Our main assumption is that the eigenvector $f_{1}$ associated with the largest eigenvalue $\lambda_{1}\equiv\lambda_{L}$ is rotationally invariant, $f_{1}(\boldsymbol{p})=f_{1}(|\boldsymbol{p}|)$. Then the eigenvalue equation can be simplified by projecting both sides onto states with angular momentum $l=0$, i.e. ‘$s$-wave’ states. The resulting eigenvalue equation involves an integral transform of $f$ which is one dimensional, involving only the magnitude of the momentum. The assumption of rotational invariance essentially permits a more precise numerical computation of the spectrum of $\cal{M}$ within the rotationally invariant subspace.

The strategy is then three-fold. First, we compute the spectral function $A_{\lambda}$ numerically on a fine grid in $(\omega,|{\boldsymbol{q}}|)$ space. Second, we use those data to construct the rung functions and to perform an approximate angular integral in the eigenvalue equation. Third, we solve the eigenvalue equation on a discretized grid in the radial momentum. Since the rung functions effectively vanish when the norms of their arguments greatly exceed $T$, we use a simple linearly spaced grid with a hard momentum cutoff of order $10~{}T$. As discussed above, the temperature is the only scale in the problem, so the entire numerical computation is set up in terms of the scaled variables $\mathfrak{q}=q/T$ and $\mathfrak{w}=\omega/T$.

The data for $1/\tau_{\varphi}$ and $\lambda_{L}$ as a function of grid size are shown in Figures 7 and 8. We roughly estimate, based on a comparison with the high precision value of $1/\tau_{\varphi}$, that our results are accurate to roughly $10$%. However, we caution the reader that $\tau_{\varphi}$ is not sensitive to the type-II rung, so this is only a partial benchmark.

We emphasize that a positive growth exponent proportional to $T/N$ is a robust prediction of our calculation. To better explain the limitations of our calculation and the sources of error, we provide additional details of the calculation in Appendix G.

5 Butterfly velocity

In this section we go back to studying $\mathcal{C}(t,\boldsymbol{x})$ when $\boldsymbol{x}$ is not integrated over, which is given by

[TABLE]

This amounts to evaluating the ladder-sum as earlier, but at a finite external momentum. For the purpose of illustration, we start again by considering only the type-I rung (first line of Figure 6); the generalization due to the type-II rung is straightforward.

At finite external momentum, the product of the $\varphi$ propagators is

[TABLE]

Evaluating the integral over $\omega$ by residue again gives rise to a variety of terms; the long time behavior is well approximated by retaining just two terms,

[TABLE]

The overall sign is again positive due to the clockwise contour and the explicit minus sign in the residue.

These two terms generalize the $1/\nu$ terms from the ${\boldsymbol{k}}=0$ calculation; they are the most singular terms at small $\nu$ and ${\boldsymbol{k}}$ and hence are expected to lead to the fastest growth. Restoring the self-energy correction but ignoring the correction due to a finite ${\boldsymbol{k}}$ on the damping, i.e. $\Gamma_{{\boldsymbol{k}}-\boldsymbol{p}}\approx\Gamma_{\boldsymbol{p}}$, we have,

[TABLE]

Let us now proceed to evaluate the ladder-sum. We first define (as was done previously),

[TABLE]

such that the Bethe-Saltpeter equation for the first line in Figure 6 for $f$ becomes,

[TABLE]

We ignore the homogenous term as before and make the following ansatz for $f(\nu,{\boldsymbol{k}};\omega,\boldsymbol{p})$:

[TABLE]

Roughly speaking, the quantities $f_{+}$ and $f_{-}$ are analogous to the density of particles and holes in a Boltzmann-like equation, to be derived shortly. The “on-shell” condition now accounts for propagation with two opposite velocities, as will become clear soon. Note that the above reduces to the previous ansatz in the limit of ${\boldsymbol{k}}\rightarrow 0$, where $f_{+}(\nu,\boldsymbol{0},\boldsymbol{p})=f_{-}(\nu,\boldsymbol{0},\boldsymbol{p})=f(\nu,\boldsymbol{p})$.

The ladder sum in Eq. (5.5) becomes,

[TABLE]

where ${\cal{R}}_{1,\pm}(\boldsymbol{p}^{\prime},\boldsymbol{p})$ are given by Eq. (3.16). We have made the following simplifications in arriving at the above equations: $\epsilon_{\boldsymbol{p}^{\prime}}\epsilon_{{\boldsymbol{k}}-\boldsymbol{p}}\approx\epsilon_{\boldsymbol{p}}\epsilon_{\boldsymbol{p}^{\prime}}$ and $\epsilon_{\boldsymbol{p}}-\epsilon_{{\boldsymbol{k}}-\boldsymbol{p}}\equiv\delta\epsilon_{{\boldsymbol{k}}}\approx{\boldsymbol{k}}\cdot v_{\boldsymbol{p}}$, where $v_{\boldsymbol{p}}=\nabla_{\boldsymbol{p}}\epsilon_{\boldsymbol{p}}$ is the group velocity.

The analogous equations with the second rung included are obtained by replacing ${\cal{R}}_{1,\pm}$ by ${\cal{R}}_{1,\pm}+{\cal{R}}_{2,\pm}$. From now on, the second rung will be included in the calculation. It is convenient to define ${\cal{R}}_{\pm}={\cal{R}}_{1,\pm}+{\cal{R}}_{2,\pm}$.

The above two equations can be combined into a compact expression by introducing the vector $\Psi^{T}(\nu,{\boldsymbol{k}};\boldsymbol{p})=(f_{+}(\nu,{\boldsymbol{k}};\boldsymbol{p})~{}~{}f_{-}(\nu,{\boldsymbol{k}};\boldsymbol{p}))$, and the matrix,

[TABLE]

The equation for $\Psi$ then takes a form that is similar to a kinetic theory Boltzmann equation,

[TABLE]

Finally, recall that $\Gamma_{\boldsymbol{p}}~{}(\sim 1/N)$ can itself be expressed in terms of ${\cal{R}}_{1}$.

The matrix structure can be further simplified with the ansatz that the “particle” density at momentum $\boldsymbol{p}$ is equal to the “hole” density at momentum $-\boldsymbol{p}$, i.e. $f_{+}(\nu,{\boldsymbol{k}};-\boldsymbol{p})=f_{-}(\nu,{\boldsymbol{k}};\boldsymbol{p})$. One can then decouple the equations for $f_{\pm}$ in Eq. (5.8) as,

[TABLE]

The final form of the finite ${\boldsymbol{k}}$ equation, let us say for $f_{+}$, is thus

[TABLE]

where we use the shorthand notation $f_{+}(\boldsymbol{p})\equiv f_{+}(t,{\boldsymbol{k}};\boldsymbol{p})$ and the expression for $\widehat{\mathcal{K}}(\boldsymbol{p}^{\prime},\boldsymbol{p})$ can be read off from Eq. (LABEL:decoupled),

[TABLE]

Note that we have factored out the $1/N$ dependence in $\Gamma_{\boldsymbol{p}}$ above.

We analyze Eq. (5.15) by thinking of the ${\boldsymbol{k}}$ term as a perturbation. To do this consistently, we should first render ${\boldsymbol{k}}$ dimensionless. Let $\lambda_{0}\sim T/N$ denote the largest eigenvalue of $\widehat{\mathcal{K}}$ times the temperature at ${\boldsymbol{k}}=0$. Then $1/\lambda_{0}$ defines a length-scale (recall that $c=1$) associated with the inelastic physics of growing chaos. Just like the dephasing-time or length-scale, $1/\Gamma_{0}=\tau_{\varphi}$, the length-scale $1/\lambda_{0}$ is proportional to $N/T$ and hence grows at large $N$. For simplicity we thus introduce the dimensionless variable

[TABLE]

and consider perturbation theory in small ${\boldsymbol{u}}$.

Treating ${\boldsymbol{u}}$ as small amounts to considering distances that are long compared to $1/\lambda_{0}$ and operators that are smeared on the scale of $1/\lambda_{0}$. This is a physical requirement because there is an important order of limits issue. For distances smaller than $1/\lambda_{0}$, the physics is that of free particles, so our analysis of the chaotic growth of operators only applies for times and distances longer than $1/\lambda_{0}$.

We now consider the leading eigenvalue $\lambda_{L}({\boldsymbol{u}})$ for small ${\boldsymbol{u}}$. The leading result is of course just $\lambda_{0}$ and the first order correction vanishes because of the spherical symmetry of the leading eigenvector at ${\boldsymbol{k}}=0$. Equivalently, we are assuming that the largest eigenvalue is a function of only $u=|{\boldsymbol{u}}|$. The first correction occurs at second order in $u$,

[TABLE]

In this formula, $\lambda_{0}$ and $\lambda_{2}$ are positive and proportional to $T/N$.

Returning to $\mathcal{C}(t,{{\boldsymbol{k}}})$, we find that it grows at the rate $\lambda_{L}({\boldsymbol{u}})$. $\mathcal{C}(t,\boldsymbol{x})$ is obtained by taking the Fourier transform,

[TABLE]

As discussed in Eq. 1.12, we assume that the physically relevant dependence on ${\boldsymbol{k}}$ comes from the exponential growth factor (i.e. is dominated by the saddle-point $\lambda_{L}(\boldsymbol{k})$, and not by the eigenvector). We then estimate

[TABLE]

The result of the ${\boldsymbol{k}}$ integral is

[TABLE]

where

[TABLE]

and

[TABLE]

Recall that $\lambda_{0},\lambda_{2}\sim T/N$, which implies that $D_{L}\sim N/T$.

Balancing the two terms in the exponent of $\mathcal{C}(t,\boldsymbol{x})$ in Eq. (5.21) gives a ballistic condition

[TABLE]

which translates to a butterfly velocity

[TABLE]

Since $\lambda_{0}$ and $\lambda_{2}$ both scale like $T/N$, it follows that $v_{B}$ is independent of $T$ and $N$. This is consistent with $v_{B}\sim v_{0}T^{1-1/z}$ obtained from simple scaling arguments with $[x]=-1$ and $[t]=-z=-1$, where $z$ is the dynamical exponent ($v_{0}$ is some natural velocity scale in the problem, appropriately dimensionalized) and has been pointed out in earlier holographic studies of quantum chaos [67] [42]. There is a subtle order of limits issue here as discussed above; the result for $v_{B}$ is only valid on scales longer than $1/\lambda_{0}$, which diverges as $N\rightarrow\infty$, i.e. as we approach the free particle limit.

We may also take the exponent of $\mathcal{C}(t,\boldsymbol{x})$, write $t=|\boldsymbol{x}|/v_{B}+(t-|\boldsymbol{x}|/v_{B})$, and expand in $t-|\boldsymbol{x}|/v_{B}$ (valid near the wavefront) to yield

[TABLE]

Our numerical computation of $v_{B}$ proceeds as in Section 4 except that we no longer assume rotational invariance. This is because the leading eigenvector is not rotationally invariant at finite external momentum ${\boldsymbol{k}}$. We approximately compute the largest eigenvalue of the integral kernel $\hat{\mathcal{K}}$ as a function of momentum ${\boldsymbol{k}}$ and find excellent agreement with the form postulated in Eq. (5.18). From these data we extract $\lambda_{L}$ and $v_{B}$; we also extract $1/\tau_{\varphi}$. The data are shown in Figures 9, 10, and 11.

Curiously, the extrapolated value of $v_{B}/c$ is very slightly greater than one, but of course we must have $v_{B}/c\leq 1$ by causality. However, the largest value obtained in the calculation for the largest grid size considered by us is $v_{B}/c\approx.999$. In fact, if the form of $\mathcal{C}$ in Eq. (5.21) is valid, we claim that $v_{B}/c$ should be strictly less than one. This is because Eq. (5.21) does not exactly vanish outside the “butterfly cone” $|\boldsymbol{x}|=v_{B}t$. Since the commutator must exactly vanish outside the light cone $|\boldsymbol{x}|=ct$, it follows that either Eq. (5.21) is incorrect or $v_{B}<c$. All our numerical results are consistent with the basic form Eq. (5.21), so we argue that $v_{B}<c$.

6 Onset of chaos in proximate phases

Although we leave the full quantitative calculations to future work, the framework developed above enables a qualitative discussion of the behavior of $\lambda_{L}$ and $v_{B}$ in the symmetry broken and unbroken phases separated by the quantum critical point. We focus only on the low temperature limit inside each phase; more generally there is a complicated crossover function which describes the entire phase diagram but is beyond the scope of this work.

First consider the symmetry unbroken phase. Here the $\varphi$ particles have a finite mass, so at the lowest temperatures the system will correspond to a dilute gas of weakly interacting well-defined $\varphi$ quasiparticles. The scattering rate (or inverse dephasing time) will be proportional to the square of the particle density since two particles must approach close to each other in order to scatter. Assuming the chaos exponent is comparable to the scattering rate (as we found above), we presumably have $\lambda_{L}\sim e^{-2\mu/T}$ where $\mu\sim\sqrt{|g-g_{c}|}$ is the mass of $\varphi$. Such a dilute gas of $\varphi$ particles has an interparticle spacing much larger than the thermal wavelength, hence the gas is effectively classical. Expanding the dispersion $\epsilon_{{\boldsymbol{k}}}$ near zero momentum gives $\epsilon_{{\boldsymbol{k}}}=\mu+|{\boldsymbol{k}}|^{2}/2\mu$, so equipartition of energy gives a typical speed of the particles of order $\sqrt{T/\mu}$. Keeping in mind the order of limits issue discussed in Section 5, we expect that at the longest times there is a butterfly velocity going like $v_{B}\sim\sqrt{T/\mu}$. This is also consistent with the scaling form of $v_{B}$ with $z=2$.

Now consider the symmetry broken phase. Our analysis here is schematic since the symmetry broken phase does not survive to non-zero temperature in two spatial dimensions. The symmetry broken phase can be stabilized at non-zero temperature by considering a three dimensional system, e.g. weakly coupled layers of the two dimensional system.

The dominant low energy degrees of freedom are the massless Goldstone modes denoted $\pi_{a}$, $a=1,...,N-2$ (the amplitude mode is gapped and by the discussion in the previous paragraph, leads to slower growth of chaos at low temperatures). However, in addition to being massless, they are also derivatively coupled which will suppress the growth of chaos at the lowest temperatures. For $N>2$ (the non-linear sigma model) the leading correction to the free Goldstone action is determined just by symmetry,

[TABLE]

where $\rho_{s}$ is the stiffness. The dimension of $\rho_{s}$ is $[\rho_{s}]_{2}=1$ in two spatial dimensions.

Introducing the canonically normalized field $\chi_{a}=\sqrt{\rho_{s}}\pi_{a}$, we see that the second term is an irrelevant interaction with a universal coefficient determined by $\rho_{s}$. Assuming the basic rung contains two factors of the interaction vertex, dimensional analysis then suggests that in two dimensions the chaos exponent goes like $\lambda_{L,2}\sim T^{3}/\rho_{s}^{2}$. The Goldstone bosons are massless and hence move at the speed of light, however, by analogy with the Boltzmann-like equation derived above, the butterfly velocity may in general be smaller than but of the same order as the speed of light.

The case $N=2$ is special since the leading symmetry controlled irrelevant operator vanishes. The leading correction to free Goldstone action, which is $\rho_{s}(\partial\theta)^{2}$, is now determined by non-universal physics via the operator $g(\partial\theta)^{4}$. In two spatial dimensions $g$ has energy dimension $[g]=-1$. As before, introducing $\chi=\sqrt{\rho_{s}}\theta$ gives a canonically normalized free scalar with an irrelevant interaction $(g/\rho_{s}^{2})(\partial\chi)^{4}$. Dimensional analysis now suggests that the chaos exponent is of order $g^{2}T^{7}/\rho_{s}^{4}$. However, we still expect a finite butterfly velocity of the order of the speed of light.

However, this analysis of the $N=2$ case is complicated by the presence of vortex excitations in two dimensions. Unlike for $N>2$, the physics of the symmetry broken phase is not completely destroyed at non-zero temperature. Phase fluctuations of the complex $N=2$ order parameter do destroy the long-range order, but there is a distinct finite temperature phase in which vortices have not proliferated.

For completeness, we record our expectations for the three dimensional case where the symmetry broken state is stable at small non-zero temperature. The dimension of $\rho_{s}$ is now $[\rho_{s}]_{3}=2$ in three spatial dimensions. Again introducing a canonically normalized scalar field and assuming that the basic rung has two copies of the interaction vertex, dimensional analysis gives the estimate $\lambda_{L,3}\sim T^{5}/\rho_{s}^{2}$ and $v_{B}\sim c$.

7 Discussion

In this work we studied the onset of chaos in the $O(N)$ nonlinear sigma model in a controlled $1/N$ expansion. As the model is weakly interacting at infinite $N$, we found a chaos exponent $\lambda_{L}=3.2~{}T/N$ which vanished in the limit of $N\rightarrow\infty$ and trivially obeys the chaos bound. The calculation predicts that models of smaller $N$ scramble faster and, making an uncontrolled extrapolation all the way to $N=1$, suggests that the Ising quantum critical point scrambles with a growth exponent $\lambda_{L}^{(N=1)}\approx 3.2~{}T$. We also studied the butterfly velocity and found that at weak coupling the butterfly velocity (in units of the speed of light) is a universal number independent of $T$ and $N$ at large $N$. At finite $N$ we expect $v_{B}$ to remain independent of $T$ but to potentially depend on $N$. The results for the $O(N)$ model can be contrasted with a different large $N$ problem of a Fermi surface of $N$ component fermions coupled to a $U(1)$ gauge field, where $\lambda_{L}$ is also proportional to $T$ but not suppressed by $1/N$ (although the precise numerical prefactor is difficult to reliably compute due to the breakdown of the $1/N$ expansion which removes the $1/N$ suppression in $\lambda_{L}$) and $v_{B}\sim T^{1/3}$ [45].

Our calculation is uncontrolled at small $N$, which includes the most physically relevant cases, so in this case it would be quite interesting to carry out experimental measurements in suitable model systems. It would also be interesting to carry out a complementary $\epsilon=3-d$ expansion at a fixed $N$ in order to study scrambling at the Wilson-Fisher conformal fixed point [60].

Recent work has also focused on the possible relation between diffusion, the chaos exponent, and the butterfly velocity [27] [67] [62]. In particular, a generalized version [62] of the Sachdev-Ye-Kitaev (SYK) [68] [22] finds an exact relationship between the energy diffusion constant, $D_{E}=v_{B}^{2}/\lambda_{L}$, while no such relationship exists for the charge diffusion constant, $D_{c}$ [69]. The energy diffusion constant is undefined at the $O(N)$ quantum critical point due to an absence of momentum relaxation. If the stress tensor is given by $T^{\mu\nu}$, then the energy current $T^{0i}$ is also the momentum density. Hence the total energy current cannot relax unless momentum is also not conserved. However, if we choose a “charge” $U(1)$ subgroup of the $O(N)$ symmetry, then the charge (for $N=2$), or, spin (for $N=3$) diffusion constant is well defined and known at large $N$.

The Einstein relation is $D_{c}=\sigma_{\textnormal{dc}}/\chi$ where the charge susceptibility is $\chi=\sqrt{5}\Theta T/2\pi$ [60] and the dc conductivity is $\sigma_{\textnormal{dc}}=0.085N$ [60] [70]. The Einstein relation then gives $D_{c}=0.249N/T$. By comparison, we find that $v_{B}^{2}/\lambda_{L}\approx 0.31N/T$ which is reasonably close to the value for $D_{c}$. Note, however, that in two spatial dimensions so-called long-time tails ultimately lead to divergences in transport coefficients [71] [72], although this effect is suppressed at large $N$. Hence in the scaling limit the true diffusion constant is infinite while $v_{B}^{2}/\lambda_{L}$ remains finite. It continues to be an interesting problem to explore the extent to which charge and energy diffusion constants in generic strongly correlated systems without a well-defined quasiparticle description are related to quantum chaos. If any relationship exists, this would be interesting for studying transport in various ‘strange’ metallic phases that are ubiquitous in correlated materials, such as the cuprates, the pnictides and other heavy-fermion materials [73].

Given that the calculation focused on the elementary field operators, one may also wonder about the extent to which commutators of other operators grow at the same rate. That commutators of other operators cannot grow more slowly is easily seen. To illustrate the physics, consider the example of the composite operator $\varphi^{2}$; in the nonlinear sigma model the field $\lambda$ plays the role $\varphi^{2}$. Consider then

[TABLE]

The commutator satisfies

[TABLE]

where $...$ denotes three other terms with different placements of the operators. In perturbation theory we can then contract the extra factors of $\varphi$ to leave the original expression for $\mathcal{C}$, i.e.

[TABLE]

where $...$ denotes contributions from additional diagrams that don’t factorize. Note also that the first term will be exponentially small at finite temperature if $|\boldsymbol{x}|\gg\beta$, so only the second term contributes, but this term amounts to $\mathcal{C}$ up to a local renormalization equal to $\langle\varphi^{2}\rangle^{2}$. Assuming the other $...$ contributions do not magically cancel this term, we have shown that $\mathcal{C}_{\varphi^{2}}$ grows as fast as $\mathcal{C}$ (note, however, that the time for $\mathcal{C}_{\varphi^{2}}$ to become large may differ owing to the renormalization).

In diagrammatic terms, given any set of composite operators, we can always dress a copy of the elementary $\varphi$ ladder with some local renormalizations (contractions between operators at the same spacetime point) to obtain a set of diagrams that yield an exponentially growing composite operator square commutator. Thinking more broadly about weakly coupled theories, this general argument suggests that commutators of composite operators should grow as fast as commuators of elementary operators (fields that give a good description of the unperturbed theory). What this general argument does not show is that commutators of composite operators must grow at the same rate as commutators of elementary operators. One could worry that they could potentially grow faster, e.g. due to some special collective modes that were not included in the “elementary ladder”. However, this is unlikely on physical grounds. Diagrammatically, we can imagine inserting the “composite ladder” into the diagrammatic expansion for the commutator of elementary fields, thus showing that they actually grow at the same rate. Physically, products of local composite operators generate sums over local operators, so generic sets of operators are expected to all have commutators which grow at the same rate. In short, chaos in one local operator easily spreads to all other local operators. Possible exceptions to this rule include non-local operators like defect operators or line operators.

Finally, the computations presented here can be extended to a variety of other models, including to fermionic QED3 in (2+1)-dimensions [74] and to the weak coupling Banks-Zaks fixed point in (3+1) dimensions [75]. It will be interesting to understand, for example, if a more highly entangled critical point like the one described by QED3 scrambles faster than the corresponding $O(N)$ critical point [76].

8 Acknowledgements

We thank J. McGreevy, A.A. Patel, S. Sachdev, D. Stanford, J. Steinberg, W. Witczak-Krempa and L. Zou for valuable disussions. We especially thank W. Witczak-Krempa and S. Sachdev for helping benchmark some of our numerical computations of the phase-coherence time and L. Zou for a careful reading of our manuscript. DC is supported by a postdoctoral fellowship from the Gordon and Betty Moore Foundation, under the EPiQS initiative, Grant GBMF-4303, at MIT. BGS is supported by the Simons Foundation as part of the It From Qubit collaboration and through a Simons Investigator Award to Senthil Todadri and by a MURI grant W911NF-14-1-0003 from ARO.

Appendix A Thermal mass

Here we review the standard computation for the thermal contribution to the mass, $\mu$, at the critical point $g=g_{c}$. In Eq. (2.1), and at infinite $v$, the $\lambda$ field exactly implements the hard constraint

[TABLE]

At large $N$, the constraint is implemented on average and so we may replace this hard constraint by its expectation value,

[TABLE]

It is convenient to adopt a regulator scheme where we keep all frequencies but implement a hard momentum cutoff $\Lambda$. The quantum critical point is then determined by setting $\mu^{2}=0$:

[TABLE]

For $d=2$ the above implies,

[TABLE]

The thermal mass is then determined by

[TABLE]

The analysis is easiest if we add and subtract the zero temperature integral but with $\mu^{2}$ restored:

[TABLE]

Both the LHS and the RHS are manifestly UV convergent; the result is a non-linear equation for $\mu=\Theta T$.

The equation for the thermal mass $\mu=\Theta T$ is (the integral is done by changing variables to $\epsilon/T$)

[TABLE]

$\Theta$ thus obeys

[TABLE]

with the usual physical solution

[TABLE]

Appendix B One-loop bubble

Here we evaluate the $\varphi_{a}$ bubble:

[TABLE]

where $\int^{\Lambda}_{\boldsymbol{k}}\equiv\int_{|{\boldsymbol{k}}|<\Lambda}\frac{d^{d}{\boldsymbol{k}}}{(2\pi)^{d}}$. We first rewrite the bubble as,

[TABLE]

where $b(z)=(e^{\beta z}-1)^{-1}$. Evaluating the residues and simplifying the final expressions, we obtain,

[TABLE]

The overall minus sign arises because the non-Matsubara poles are encircled clockwise instead of counterclockwise.

The imaginary part of the retarded bubble after analytically continuing to real frequencies $(i\nu_{n}\rightarrow\nu+i0^{+})$ is,

[TABLE]

We note that the above is an odd function of $\nu$, as should be the case. It is positive for positive frequency $\nu$.

As a simple check, let us evaluate the bubble in the zero momentum limit, ${\boldsymbol{q}}\rightarrow 0$. This limit yields

[TABLE]

In two spatial dimensions, $d=2$, we find

[TABLE]

[TABLE]

This matches nicely with what we would expect from the Cutkosky cutting rules. This result can be compared with the zero temperature result,

[TABLE]

Continuing to real frequencies, $i\nu_{n}\rightarrow\nu+i0^{+}$, the imaginary part of $\Pi_{T=0}$ is non-zero for $\nu>q$. We have the standard result that

[TABLE]

In the limit that $\nu\gg T$, the zero temperature and finite temperature results agree, e.g. as can be seen from the finite temperature formula at ${\boldsymbol{q}}=0$.

Returning to the general case, since $\textnormal{Im}[\Pi_{R}(\nu+i0^{+},q)]$ is odd in $\nu$, we only have to evaluate it for $\nu>0$. In that case, we can drop the last term in the expression Eq. (B) for $\textnormal{Im}[\Pi_{R}(\nu+i0^{+},q)]$ since the argument in the delta function is positive. Moreover, rotational invariance implies that the first term is equal to the second term, so we find

[TABLE]

In order to evaluate the above quantity numerically, it is useful to recast it as follows. Shift the ${\boldsymbol{k}}$ integral so that the energies are

[TABLE]

The first delta function requires

[TABLE]

while the second requires

[TABLE]

These may be put in the common form

[TABLE]

where $s=\pm 1$. Squaring both sides gives

[TABLE]

In radial coordinates $\{|{\boldsymbol{k}}|,\theta\}$, the energies are

[TABLE]

The above condition then translates to

[TABLE]

Note that the integral is only over positive $|{\boldsymbol{k}}|$, so there is only one root to consider. Furthermore, in order to have $|{\boldsymbol{k}}|_{0}$ real, we must have $\nu^{2}<|{\boldsymbol{q}}|^{2}\cos^{2}\theta$ or $\nu^{2}>|{\boldsymbol{q}}|^{2}+4\mu^{2}$. Finally, to actually satisfy the delta function, we must check that $\epsilon_{+}+s\nu=s\epsilon_{-}$. For $\nu<q$ the only solution is for $s=1$ in the range $\theta\in(\pi-\cos^{-1}\frac{\nu}{q},\pi+\cos^{-1}\frac{\nu}{q})$. For $\nu>q^{2}+4\mu^{2}$ the only solution is for $s=-1$ in the range $\theta\in(0,2\pi)$.

Appendix C Wightman function

In order to express the symmetric Wightman function, ${\cal G}_{W}(\omega,{\boldsymbol{k}})$ in terms of the spectral function, we first insert a complete set of (many-body) eigenstates, $|n\rangle$ (with energy $E_{n}$ and momentum ${\boldsymbol{P}}_{n}$)

[TABLE]

where $Z$ is the partition-function and $|\varphi_{mn}|^{2}=|\langle m|\varphi|n\rangle^{2}$.

Now recall that for ${\cal G}_{R}(\omega,{\boldsymbol{k}})$ a similar procedure leads to

[TABLE]

and so we can read off the spectral function as,

[TABLE]

The above can be simplified further by using the delta-function constraint,

[TABLE]

Therefore, the Wightman function can be expressed as,

[TABLE]

which leads to the useful relation

[TABLE]

Appendix D One-loop self-energy

Here we evaluate the self-energy of $\varphi_{a}$:

[TABLE]

Let us focus on the contribution from the first term above, which we denote $\tilde{\Sigma}_{1}(i\omega_{n},{\boldsymbol{q}})$,

[TABLE]

We can re-express the above as

[TABLE]

Evaluating this leads to,

[TABLE]

The minus signs in front of the second two terms again arise from the clockwise encircling of the poles.

Let us shift ${\boldsymbol{k}}\rightarrow{\boldsymbol{k}}-{\boldsymbol{q}}$ and further simplify the above expressions,

[TABLE]

We have used the formula Eq. (2.17) relating the Euclidean $\lambda$ correlator with the retarded correlator (note the minus sign). The analytic continuation is $i\omega_{n}\rightarrow\omega+i0^{+}$ and the imaginary part is

[TABLE]

The plus signs in the last two terms arise because $\text{Im}[{\cal G}_{\lambda}(x-i0^{+})]=-\text{Im}[{\cal G}_{\lambda}(x+i0^{+})]$. The overall positive sign comes from the relative sign relating ${\cal G}_{\lambda}$ and ${\cal G}_{R,\lambda}$.

Further simplification yields

[TABLE]

Finally,

[TABLE]

It should be clear from the structure of the above expression that $\textnormal{Im}[\tilde{\Sigma}]=\textnormal{Im}[\tilde{\Sigma}_{1}]$, i.e. the second piece in Eq. (D.1) doesn’t contribute to the imaginary part.

It is useful to go one step further and express the self-energy in terms of the Wightman function. Eq.  (D.8) then becomes,

[TABLE]

Recall that the spectral function,

[TABLE]

and using the relationship between the Wightman function and the spectral function in Eq. (C.7), we have

[TABLE]

Tracking through the minus signs, we see that the imaginary part of the self-energy is positive for positive frequency.

Appendix E Feynman rules for ladder diagrams

Here we show how the diagrammatic rules for performing the ladder sum are derived. We focus on the “bare” ladder; the dressing of propagators due to self-energy effects is standard and will not be reviewed here. The basic object of interest is again the squared commutator ($\varphi(0)\equiv\varphi(\boldsymbol{0},0)$):

[TABLE]

Recall that at infinite $N$ the fields $\varphi_{a}$ are free, so $\mathcal{C}$ reduces to a free particle result. There are two sources of $1/N$ corrections to the $N\rightarrow\infty$ form of $\mathcal{C}$, one from the thermal state $\rho$ and one from the time-evolution operator $e^{-iHt}$ defining the Heisenberg operators in Eq. (E.1). As discussed in the main text, we may ignore the corrections to the thermal state when computing the chaos exponent to leading order in $1/N$. The nontrivial contributions to scrambling thus arise from the time-evolution operator.

It is convenient to go to the interaction picture with respect to the non-interacting $N\rightarrow\infty$ Hamiltonian. Then the non-trivial part of the time evolution operator is given by a time-ordered exponential of the interaction term in the interaction picture:

[TABLE]

(the overall minus sign in the Hamiltonian arises because the interaction has a positive sign in the Lagrangian). The subscript ‘[math]’ indicates that these fields evolve under the non-interacting Hamiltonian. To be precise, the interaction picture time-evolution operator is

[TABLE]

For clarity of presentation, we henceforth drop the factors of $N$ and the index structure. These can be easily restored using the rules discussed in the main text. Our purpose in this appendix is to understand the overall sign of diagrams, the various combinatorial factors, and the analytic structure i.e. which Greens functions appear.

We may thus rewrite $\mathcal{C}$ as (after dropping factors of $N$ and various indices)

[TABLE]

This formula is the basis of the perturbative expansion. The perturbative expansion for $U_{I}$ is

[TABLE]

where we have been careful with the time ordering. The perturbative expansion for $U_{I}^{\dagger}$ is

[TABLE]

where we have been careful to reverse the terms in the second order expression. We also used $\lambda_{0}^{\dagger}=\lambda_{0}$ and $\varphi_{0}^{\dagger}=\varphi_{0}$ and have dropped the spatial arguments in the expressions above.

The leading contribution to $\mathcal{C}$ is obtained by setting $U_{I}=1$. Then since

[TABLE]

is proportional to the identity operator, we obtain

[TABLE]

The first order contribution is obtained by expanding $U_{I}$ to first order,

[TABLE]

Since we need two powers of $\lambda_{0}$ to get a non-zero thermal expectation value (the expectation value giving the mass has been subtracted out), we need to either expand $U_{I}$ on a single time fold to higher order or expand $U_{I}$ to first order on both time folds. The first option leads to the standard self-energy correction along a single time fold and will not be reviewed here. The second option gives the first non-trivial ladder contribution. We have

[TABLE]

Using the free particle results $[\lambda_{0},\varphi_{0}]=0$, $[\varphi_{0}(\boldsymbol{x},t),\varphi_{0}(0)]=i{\cal G}_{R}(\boldsymbol{x},t)$, and $[\varphi_{0}(\boldsymbol{x},t),\varphi_{0}^{2}(\boldsymbol{y},s)]=[i{\cal G}_{R}(\boldsymbol{x}-\boldsymbol{y},t-s)]2\varphi_{0}(\boldsymbol{y},s)$ (valid provided $t>s$) gives,

[TABLE]

or

[TABLE]

Up to factors of $1/N$ this is precisely the claimed one rung contribution of the type-I rung. Note the absence of any extra combinatorial prefactors.

The two rung contribution to the ladder is obtained by expanding $U_{I}$ to second order on both time folds. Suppressing the spatial labels for ease of notation, the second order expansion is

[TABLE]

Extracting an overall factor of $i^{2}/4$ and breaking the $s_{1}$ integral in the second term into two pieces ($s_{2}<s_{1}$ and $s_{2}>s_{1}$) and relabelling time integration variables, we recognize the formula

[TABLE]

There are now two types of terms that arise from the nested commutators. The first type generates the two-rung contributions from the type-I rung and the second type is responsible for the one-rung contribution from the type-II rung (see Figure 5). The expansion of the nested commutator from above is,

[TABLE]

which gives

[TABLE]

As indicated above, the first line of Eq. (E.16) gives rise to a two-rung contribution with type-I rungs as well as a self-energy correction for $\varphi$ (depending on whether the $\lambda$ fields are contracted within a time fold or between time folds). The second line of Eq. (E.16) gives rise to a one-rung contribution with a type-II rung and a self-energy correction for $\lambda$ (coming from the two terms in $[\varphi_{0}(s_{1})\varphi_{0}^{2}(s_{2}),\varphi_{0}(0)]$).

Notice that for the two-rung contribution, the factors of $4$ from the vertex and nested commutators cancel. The type-II rung contribution arises from the case where $\varphi_{0}^{2}(s_{2})$ commutes with $\varphi_{0}(0)$ to give $2\varphi_{0}(s_{2})[\varphi_{0}(s_{2}),\varphi_{0}(0)]$ and the remaining $\varphi_{0}(s_{1})\varphi_{0}(s_{2})$ factor is contracted with a corresponding pair on the other time fold. For this type two rung, the factors of $4$ also cancel. It remains to establish the overall sign. The expression is

[TABLE]

or

[TABLE]

Again, this is the claimed one rung contribution of type two including the overall positive sign and combinatorial factor of unity.

Following the manipulations for the first and second order terms, the general term starts from the expansion

[TABLE]

The various multi-rung ladder diagrams all arise from these nested commutators. The following claims can be straightforwardly verified from the general expansion. First, the combinatorial factors always cancel to give unity in ladder diagrams with either rung type. Second, the time ordering guarantees that every commutator can be replaced by a retarded Green’s function; furthermore, including the explicit factor of $\theta(t-s)$ in $G_{R}(t-s)$ allows us to extend all time integrals to run over the entire line $(-\infty,\infty)$. Third, the overall sign is positive when written as an integral over products of ${\cal G}_{R}$ and ${\cal G}_{W}$. This follows because in an $\ell$ rung ladder, there are $2\ell+2$ commutators which each give a factor of $i$, there are $2\ell$ factors of $i$ from the interaction vertices, and there is an overall minus sign, which gives a total sign of $-i^{2\ell+2}i^{2\ell}=+1$.

Appendix F Ladder versus crossed diagrams

In this appendix, we discuss an important issue, namely whether it is justified to retain only the ladder series (see Figure 12a) as discussed earlier or if one needs to include the crossed diagrams as shown in Figure 12b. In order to study this aspect in more detail, consider the two-rung ladder and two-rung crossed diagrams respectively, and focus only on the internal “loop”.

The internal loop in the two-rung ladder is given by,

[TABLE]

Similarly, the loop in the two-rung crossed diagram is given by,

[TABLE]

Let us now compare and contrast the pole structures of the above diagrams. For the ladder diagrams, the retarded propagators always occur in pairs of the form, ${\cal G}_{R}(\nu-\omega,-\boldsymbol{p})~{}{\cal G}_{R}(\omega,\boldsymbol{p})$, with an integral over $\boldsymbol{p}$, which upon doing the $\omega$ integral leads to a pole structure dominated by the $\nu^{-1}$ term, as far as the long time form of ${\cal{C}}(t)$ is concerned. Moreover, the contribution to $\chi_{\textnormal{ladder}}$ above from all three pieces then goes as $\nu^{-3}$, since the different pairs of the retarded propagators depend only on one frequency that is being integrated over (e.g. $\varepsilon,\omega,\omega^{\prime}$) 555The frequencies are, of course, not entirely decoupled. They are coupled via the rungs of the ladder, ${\cal G}_{W,\lambda}$..

On the other hand, notice that for the crossed diagram, the combination of the retarded propagators coming from the “cross” (or the flipped; see fig. 12b) region is,

[TABLE]

Note that the momenta in the two legs are not directed opposite to each other. Carrying out the $\omega$ integral leads to pole structures of the form: $(\nu-\omega^{\prime}-\varepsilon\pm\epsilon_{\boldsymbol{p}}\pm\epsilon_{\boldsymbol{p}-\boldsymbol{p}^{\prime}-\boldsymbol{\ell}})^{-1}$. This can be combined with the contributions from the other pieces in $\chi_{\textnormal{crossed}}$. It is not hard to see that upon doing the $\omega^{\prime}$ integral leads to a pole structure: $\left[\nu(\nu-\varepsilon\pm\epsilon_{\boldsymbol{p}^{\prime}}\pm\epsilon_{\boldsymbol{p}}\pm\epsilon_{\boldsymbol{p}-\boldsymbol{p}^{\prime}-\boldsymbol{\ell}})\right]^{-1}$. Finally, doing the $\varepsilon$ integral leads to: $\left[\nu^{2}(\nu\pm\epsilon_{\boldsymbol{\ell}}\pm\epsilon_{\boldsymbol{p}^{\prime}}\pm\epsilon_{\boldsymbol{p}}\pm\epsilon_{\boldsymbol{p}-\boldsymbol{p}^{\prime}-\boldsymbol{\ell}})\right]^{-1}$.

There is an equivalent real time picture of the suppression of the crossed ladder relative to the uncrossed ladder. Taking for simplicity two type-I rungs, the uncrossed ladder is schematically

[TABLE]

Roughly speaking, the effect of the Wightman propagators is to set $s_{i}=s^{\prime}_{i}$. Said differently, the Wightman functions constrain the relative variables $s_{i}-s^{\prime}_{i}\approx 0$ to a precision of about $\beta$ but the center of mass variables $\frac{s_{i}+s^{\prime}_{i}}{2}$ remain unconstrained. Then the ladder contribution is schematically

[TABLE]

On the other hand, the crossed ladder is schematically

[TABLE]

The Wightman functions now constrain $s_{1}-s^{\prime}_{2}\approx 0$ and $s_{2}-s^{\prime}_{1}\approx 0$. However, since the integrand is zero unless $s_{1}>s_{2}$ and $s^{\prime}_{1}>s^{\prime}_{2}$, the effective range of integration is severely restricted by these constraints. We may write the causality constraints as

[TABLE]

so we also have $s_{1}-s_{2}\approx 0$ to a precision of about $\beta$. This further restriction of the time arguments in the integral implies that $\chi_{\text{crossed}}$ is parametrically smaller than $\chi_{\text{ladder}}$.

As a simple model of this suppression, consider a free particle with position $X$ and decay rate $\gamma$. The equation of motion is $\partial_{t}^{2}X=-\gamma\partial_{t}X$. The retarded correlator is

[TABLE]

For $\gamma t\gg 1$ and $\gamma\beta\ll 1$ the ratio of the ladder to the crossed ladder is then of order

[TABLE]

Appendix G Details of the numerical computations

Here we discuss the details of our numerical computation of $\lambda_{L}$ and $v_{B}$. The first step is to compute the analytic continuation of the bubble. The imaginary part is computed by numerical integration of Eq. (B) using a standard MATLAB routine. The real part is then obtained by a stable version of Kramers-Kronig,

[TABLE]

Here we have subtracted a multiple of the zero integral $\mathbb{P}\int d\nu\frac{1}{\omega^{2}-\nu^{2}}=0$ so that the resulting integrand is non-singular at $\omega=\nu$ and the principal part can be removed.

This integral is computed in three parts because the imaginary part is non-analytic at $\omega=|{\boldsymbol{q}}|$ and $\omega=\sqrt{|{\boldsymbol{q}}|^{2}+4\mu^{2}}$. The numerical integral defining the real part is then carried out by breaking it into three pieces, $\nu\in(0,|{\boldsymbol{q}}|)$, $\nu\in(|{\boldsymbol{q}}|,\sqrt{{\boldsymbol{q}}^{2}+4\mu^{2}})$, and $\nu\in(\sqrt{{\boldsymbol{q}}^{2}+4\mu^{2}},\infty)$. The upper limit is taken to be a large number of order $10^{4}$ so that the result is converged to several digits. We have also checked the result by treating the discontinuity in $\text{Im}[\Pi_{R}]$ at $\omega=\sqrt{{\boldsymbol{q}}^{2}+4\mu^{2}}$ exactly and performing Kramers-Kronig on the continuous subtracted imaginary part. The methods agree to several digits.

The values of the real and imaginary parts are precomputed on a relatively fine grid in frequency and momentum for later use (this precomputing is necessary to avoid prohibitive slowdowns in later steps). For the results quoted in the main text, the precomputed grid is $400\times 400$ with $\omega/T,|{\boldsymbol{q}}|/T\in[0,12]$.

Given the real and imaginary parts of the bubble, we are then able to compute the $\lambda$ spectral function $A_{\lambda}$, the $\lambda$ retarded correlation function ${\cal G}_{\lambda,R}$, and the $\lambda$ Wightman function ${\cal G}_{\lambda,W}$. If values of $\Pi_{R}(|\boldsymbol{q}|,\omega)$ not in the precomputed grid are needed, then we use two dimensional spline interpolation to determine approximate values. If $|\boldsymbol{q}|$ or $\omega$ lie outside the grid, then we set the value of $\Pi_{R}$ to zero. These are additional approximations.

From here the calculation follows two routes depending on whether we are interested in the $\boldsymbol{k}=0$ chaos exponent or the $\boldsymbol{k}\neq 0$ chaos exponent.

G.1 $\boldsymbol{k}=0$

The basic assumption underlying this calculation is that the growing mode is spherically symmetric and depends only on $|\boldsymbol{p}|$. Hence we use radial coordinates so that the angular variable drops out of the final eigenvalue equation.

Using the computed values of the propagators, we first obtain the dephasing rate $\Gamma_{{\boldsymbol{q}}}$ by numerically integrating the rung function ${\cal{R}}_{1}$ over angle and then over radial momentum. The angular grid is taken be uniform with $n_{\theta}=60$ points, while the radial grid is also uniform and has up to $n_{\text{1d}}=200$ points.

We then average ${\cal{R}}_{1}(\boldsymbol{\ell},\boldsymbol{p})$ and ${\cal{R}}_{2}(\boldsymbol{\ell},\boldsymbol{p})$ over the angle between $\boldsymbol{\ell}$ and $\boldsymbol{p}$ for $|\boldsymbol{\ell}|$ and $|\boldsymbol{p}|$ on the radial grid. Note that there is an extra complication here, since the calculation of ${\cal{R}}_{2}(\boldsymbol{\ell},\boldsymbol{p})$ already involves an angular integral coming from the integral over the loop momentum. This angular integral is also done approximately on a grid with $n_{\theta}=60$ points.

Finally, using the dephasing rate computed before and the angular averaged versions of ${\cal{R}}_{1}$ and ${\cal{R}}_{2}$, we numerically solve a discretized version of the integral equation Eq. (3.31). We generically find one or a few real positive eigenvalues and then many more complex eigenvalues with negative real parts. By studying the system size dependence of the results, we attempt to approximately extrapolate to a converged limit.

G.2 $\boldsymbol{k}\neq 0$

In this version of the calculation we cannot assume that the growing mode is spherically symmetric. Hence we work on a general two-dimensional momentum grid. Another possibility is to break the equation down into partial waves and solve a set of coupled radial integral equations, but we did not pursue it. We define a square grid of vectors $(q_{x},q_{y})$ with $|q_{x}|,|q_{y}|<k_{\max}$ and $n_{\text{2d}}^{2}$ points. We consider up to $n_{\text{2d}}=61$ in our calculation, corresponding to matrices of size $n_{\text{2d}}^{2}\times n^{2}_{\text{2d}}=3721\times 3721$. Since the matrix is not sparse, memory sets a limit to how fine a grid we can consider. However, the time needed to set all the values of the matrices is by far the limiting step.

The dephasing rate is again computed by numerically integrating the formula for $\Gamma_{{\boldsymbol{q}}}$ over the two-dimensional grid. The value of ${\cal{R}}_{1}$ and ${\cal{R}}_{2}$ are also computed on the grid and the additional angular integral entering the definition of ${\cal{R}}_{2}$ is done approximately using $n_{\theta}=60$. We then solve the eigenvalue equation as a function of external momentum $\boldsymbol{k}$. To be precise, as discussed in the main test, by working with the scaled variable $N\boldsymbol{k}/T$, the entire eigenvalue equation is rendered proportional to $T/N$. We find that the leading eigenvalue is positive and real and obeys

[TABLE]

for $N\boldsymbol{k}/T\ll 1$. This is shown in Figure 13. There are also subleading eigenvalues including complex eigenvalues with negative real parts. Again, by studying the system size dependence of the results, we attempt to approximately extrapolate to a converged limit.

G.3 Run times

With the $400\times 400$ $\omega,|{\boldsymbol{q}}|$ grid, the full precomputation step including the calculation of the type-II rung takes a few days in MATLAB running on a laptop with two $2.4$ GHz cores and $8$ GB of memory. The calculation of the rungs in the spherically symmetric case is comparatively quick, with the matrix element setup and subsequent diagonalization taking about an hour for the largest grid sizes considered ($n_{\text{1d}}=200$). By contrast, the matrices involved in the calculation of $v_{B}$ are much larger and take much longer to initialize. The matrix element setup and subsequent diagonalization took about a day and a half for the largest grid sizes considered ($n_{\text{2d}}=61$). All told, the setup and data collection took about a week for the data presented in the main text.

As indicated above, the most expensive step by far is the initialization of the matrices. The largest matrices considered were less than $4000\times 4000$, so memory was not the primary limitation in our calculation. Given more time and further optimizations of the code and possibly parallelization, the calculation could be further improved as needed.

G.4 Extrapolation in to infinite $\omega,|{\boldsymbol{q}}|$ grid size

In the main text we reported results for the extrapolation of $1/\tau_{\phi}$, $\lambda_{L}$, and $v_{B}$ to infinite grid size in the solution of the integral equation. However, there are two additional grid sizes in the problem, the initial fine grid in $\omega,|{\boldsymbol{q}}|$ and the angular grid used to approximately compute the type-II rung. The results did not appear to depend too sensitively on the angular grid size, but because of limited time we did not pursue a systematic study of this dependence.

We did study the results as a function of the precomputed grid size. We repeated the calculations in the main text for a precomputed grid size of $200\times 200$. We then took the extrapolated values for $400$ and $200$ and performed a further linear fit to one over the precomputed grid size. Of course this is only two data points and we have not carefully studied the dependence on the angular grid, so our conclusions in this appendix are tentative.

The estimated infinite grid size values for the ${\boldsymbol{k}}=0$ calculation and the ${\boldsymbol{k}}\neq 0$ calculation are

[TABLE]

The agreement between the two methods is reasonable. Comparison with the high precision value of $1/\tau_{\phi}$ suggests errors at the level of at least $5$%. Obviously the butterfly velocity must obey $v_{B}\leq 1$, although the extrapolation gives a value slightly larger than one. Within the accuracy of the calculation we simply say that $v_{B}\approx 1$.