Many-body chaos in the antiferromagnetic quantum critical metal

TL;DR

This paper calculates chaos and scrambling rates at an antiferromagnetic quantum critical point, revealing a hierarchy of Lyapunov exponents for bosons and fermions, with implications for understanding quantum chaos in strongly coupled systems.

Contribution

It introduces a detailed analysis of many-body chaos at a quantum critical point, highlighting the different scrambling behaviors of bosonic and fermionic degrees of freedom.

Findings

Bosonic Lyapunov exponent scales as √w, larger than fermionic exponent.

Fermionic Lyapunov exponent scales as w^2, smaller than bosonic.

Chaos front exhibits infinite butterfly velocity due to non-local bosonic propagator.

Abstract

We compute the scrambling rate at the antiferromagnetic (AFM) quantum critical point, using the fixed point theory of Phys. Rev. X $\boldsymbol{7}$, 021010 (2017). At this strongly coupled fixed point, there is an emergent control parameter $w \ll 1$ that is a ratio of natural parameters of the theory. The strong coupling is unequally felt by the two degrees of freedom: the bosonic AFM collective mode is heavily dressed by interactions with the electrons, while the electron is only marginally renormalized. We find that the scrambling rates act as a measure of the "degree of integrability" of each sector of the theory: the Lyapunov exponent for the boson $\lambda_L^{(B)} \sim \mathcal O(\sqrt{w}) \,k_B T/\hbar$ is significantly larger than the fermion one $\lambda_L^{(F)} \sim \mathcal O(w^2) \,k_B T/\hbar$, where $T$ is the temperature. Although the interaction strength in the theory is of order unity, the larger Lyapunov exponent is still parametrically smaller than the universal upper bound of $\lambda_L=2\pi k_B T/\hbar$. We also compute the spatial spread of chaos by the boson operator, whose low-energy propagator is highly non-local. We find that this non-locality leads to a scrambled region that grows exponentially fast, giving an infinite "butterfly velocity" of the chaos front, a result that has also been found in lattice models with long-range interactions.

Equations

$$f_{A,B}(t)=\mathrm{Tr}\left[e^{-\beta H/2}\left[A(t),B\right]e^{-\beta H/2}{\left[A(t),B\right]}^{†}\right].$$

$$\begin{split}\mathcal{S}=&\sum_{n=1}^{4}\sum_{m=\pm}\sum_{\sigma=\uparrow,\downarrow}T\sum_{\omega_{k}}\int\frac{d\vec{k}}{(2\pi)^{2}}\,{\psi}^{(m)*}_{n,\sigma}(k)\left[i\omega_{k}+e^{m}_{n}(\vec{k};v)\right]\psi^{(m)}_{n,\sigma}(k)\\ &+\frac{1}{2}T\sum_{\omega_{q}}\int\frac{d\vec{q}}{(2\pi)^{2}}\,\Bigl{[}\left|{\omega_{q}}\right|+c(v)(\left|{q_{x}}\right|+\left|{q_{y}}\right|)+M(T,\Lambda,v)\Bigr{]}\vec{\phi}(q)\cdot\vec{\phi}(-q)\\ &+\sqrt{\frac{\pi v}{2}}\sum_{n=1}^{4}\sum_{\sigma,\sigma^{\prime}=\uparrow,\downarrow}T\sum_{\omega_{k}}\int\frac{d\vec{k}}{(2\pi)^{2}}\,T\sum_{\omega_{q}}\int\frac{d\vec{q}}{(2\pi)^{2}}\Bigl{[}\vec{\phi}(q)\cdot{\psi}^{(+)*}_{n,\sigma}(k+q)\vec{\tau}_{\sigma,\sigma^{\prime}}\psi^{(-)}_{n,\sigma^{\prime}}(k)+c.c.\Bigr{]}.\end{split}$$

$$c(v)=\frac{1}{4}\sqrt{v\log (1/v)},$$

$$w(v)\equiv \frac{v}{c(v)}\ll 1.$$

$$\begin{array}{rl}{f}_F(t,\vec{x})& =\theta (t)\frac{1}{2^2}\sum _{\sigma ,{\sigma }^{\prime }=\uparrow ,\downarrow }\mathrm{Tr}[e^{-\beta H/2}\{{\psi }_{\sigma }(\vec{x},t),{\psi }_{{\sigma }^{\prime }}^{\ast }(0)\}\\ & e^{-\beta H/2}\{{\psi }_{\sigma }(\vec{x},t),{\psi }_{{\sigma }^{\prime }}^{\ast }(0){\}}^{†}],\end{array}$$

$$\begin{array}{rl}{f}_B(t,\vec{x})& =\theta (t)\frac{1}{3^2}\sum _{i,j=1}^3\mathrm{Tr}[e^{-\beta H/2}[{\varphi }_i(\vec{x},t),{\varphi }_j(0)]\\ & e^{-\beta H/2}[{\varphi }_i(\vec{x},t),{\varphi }_j(0){]}^{†}].\end{array}$$

$$G_{(n,m)}^W(k,T)=\frac{A_{n,m}(k)}{2\cosh \frac{\beta k_0}{2}}=\frac{\pi \delta (k_0+e_n^m(\vec{k}))}{\cosh \frac{\beta k_0}{2}},$$

$$D^W(q,T)=\frac{B(q)}{2\sinh \frac{\beta q_0}{2}}$$

$$=\frac{q_0}{q_0^2+{\left(c(v)(|q_x|+|q_y|)+M(T,\Lambda ,v)\right)}^2}\frac{1}{\sinh \frac{\beta q_0}{2}},$$

$$\begin{array}{rl}& G_{(n,m)}^R(k,T)=[k_0+e_n^m(\vec{k})\\ & +i\frac{3}{4}w(v)T(2\frac{-i{k}_0-\pi T}{2\pi T}\log \left(\frac{\Lambda }{2\pi T}\right)\\ & {-2\log \left[\Gamma \left(1+\frac{-i{k}_0-\pi T}{2\pi T}\right)\right]+\log \frac{\Lambda }{M(T,\Lambda ,v)})]}^{-1},\end{array}$$

$$\displaystyle D^{R}(q,T)=\left[-iq_{0}+c(v)\Big{[}|q_{x}|+|q_{y}|\Big{]}+M(T,\Lambda,v)\right]^{-1},$$

$$f_F(\omega )=\frac{1}{2}\sum _{n=1}^4\sum _{m=\pm }\int dk{f}_F^{(n,m)}(\omega ,k),$$

$$\begin{array}{rl}{f}_F^{(n,m)}(\omega ,k)& =G_{(n,m)}^R(k)G_{(n,m)}^A(k-\omega )\\ & \times \left[1+\int d{k}^{\prime }{K}_F(k,k^{\prime },\omega )f_F^{(n,\bar{m})}(\omega ,k^{\prime })\right],\end{array}$$

$$\int dk{f}_F^{(n,m)}(\omega ,k)=\int dk{G}_{(n,m)}^R(k)G_{(n,m)}^A(k-\omega )$$

$$=i\int \frac{d{k}_x}{2\pi }\int \frac{d{k}_0}{2\pi }[\omega +i\frac{3}{4}wT(-\frac{i\omega +2\pi T}{\pi T}\log \left(\frac{\Lambda }{2\pi T}\right)$$

$$-2\log \left[\Gamma \left(1+\frac{-i{k}_0-\pi T}{2\pi T}\right)\Gamma \left(1+\frac{i(k_0-\omega )-\pi T}{2\pi T}\right)\right]$$

$${+2\log \frac{\Lambda }{M})]}^{-1}.$$

$$\begin{array}{rl}{f}_F^{(n,m)}(\omega ,k)& =G_{(n,m)}^R(k)G_{(n,m)}^A(k-\omega )\\ & \times \int d{k}^{\prime }{K}_F(k,k^{\prime },\omega )f_F^{(n,\bar{m})}(\omega ,k^{\prime }).\end{array}$$

$$K_F(k,k^{\prime },\omega )=\frac{3\pi }{2}v{D}^W(k-k^{\prime })$$

$$=\frac{3\pi }{2}v\frac{k_0-k_0^{\prime }}{(k_0-k_0^{\prime }{)}^2+{\left(c(v)(|k_x-k_x^{\prime }|+|k_y-k_y^{\prime }|)+M\right)}^2}$$

$$\times \frac{1}{\sinh \frac{\beta (k_0-k_0^{\prime })}{2}}.$$

$$\begin{array}{rl}{f}_F^{(\pm )}(\omega ,k)& =G_{\pm }^R(k)G_{\pm }^A(k-\omega )\\ & \times \frac{3\pi }{2}v\int d{k}^{\prime }{D}^W(k-k^{\prime })f_F^{(\mp )}(\omega ,k^{\prime }).\end{array}$$

$$\begin{array}{rl}{f}_F^{(\pm )}(\omega ,k_0,k_y)& =\frac{1}{k_0\pm k_y+i\frac{3}{4}wT\left(2\frac{-i{k}_0-\pi T}{2\pi T}\log \left(\frac{\Lambda }{2\pi T}\right)-2\log \left[\Gamma \left(1+\frac{-i{k}_0-\pi T}{2\pi T}\right)\right]+\log \frac{\Lambda }{M}\right)}\\ & \frac{1}{k_0-\omega \pm k_y-i\frac{3}{4}wT\left(2\frac{i(k_0-\omega )-\pi T}{2\pi T}\log \left(\frac{\Lambda }{2\pi T}\right)-2\log \left[\Gamma \left(1+\frac{i(k_0-\omega )-\pi T}{2\pi T}\right)\right]+\log \frac{\Lambda }{M}\right)}\\ & \times \frac{3}{4}v\int \frac{d{k}_0^{\prime }d{k}_y^{\prime }}{(2\pi {)}^2}\frac{\pi -2\arctan \left(\frac{c(v)|k_y-k_y^{\prime }|+M}{|k_0-k_0^{\prime }|}\right)}{c(v)|k_0-k_0^{\prime }|}\frac{k_0-k_0^{\prime }}{\sinh \frac{\beta (k_0-k_0^{\prime })}{2}}{f}_F^{(\mp )}(\omega ,k_0^{\prime },k_y^{\prime }).\end{array}$$

$$f_F(\omega ,k_0)$$

$$\times i\frac{3}{4}v\int \frac{d{k}_0^{\prime }}{2\pi }\frac{\pi -2\arctan \left(\frac{M}{|k_0-k_0^{\prime }|}\right)}{c(v)|k_0-k_0^{\prime }|}\frac{k_0-k_0^{\prime }}{\sinh \frac{\beta (k_0-k_0^{\prime })}{2}}{f}_F(\omega ,k_0^{\prime }).$$

$$\displaystyle\int\frac{dk_{0}^{\prime}}{2\pi}\Bigg{(}i\frac{3}{4}w(v)\frac{\pi-2\arctan\left(\frac{\bar{M}}{|k_{0}-k_{0}^{\prime}|}\right)}{|k_{0}-k_{0}^{\prime}|}\frac{k_{0}-k_{0}^{\prime}}{\sinh\frac{k_{0}-k_{0}^{\prime}}{2}}-2\pi\delta(k_{0}^{\prime}-k_{0})$$

$$\displaystyle\times\Bigg{[}\omega-i\frac{3}{2}w(v)\left(\left(1+\frac{i\omega}{2\pi}\right)\log\left(\frac{\bar{\Lambda}}{2\pi}\right)+\log\left[\Gamma\left(\frac{1}{2}-\frac{ik_{0}}{2\pi}\right)\Gamma\left(\frac{1}{2}+\frac{i(k_{0}-\omega)}{2\pi}\right)\right]-\log\left(\frac{\bar{\Lambda}}{\bar{M}}\right)\right)\Bigg{]}\Bigg{)}f_{F}\left(\omega\,T,k_{0}^{\prime}\,T\right)=0,$$

$${\lambda }_L^{(F)}=0+\mathcal{O}(w^2)T.$$

$$\begin{array}{rl}{f}_B(\omega ,k)& =D^R(k)D^A(k-\omega )\\ & \times \left[1+6\int d{k}^{\prime }{K}_B(k,k^{\prime },\omega )f_B(\omega ,k^{\prime })\right].\end{array}$$

$$\begin{array}{rl}{K}_B(k,k^{\prime },\omega )& =\frac{i}{\omega }\frac{{\pi }^{2}v}{16}\frac{k_0-k_0^{\prime }}{\sinh \frac{\beta (k_0-k_0^{\prime })}{2}}\\ & \times \sum _{n=1}^4\sum _{m=\pm }\delta (k_0^{\prime }-k_0+e_n^{(m)}({\vec{k}}^{\prime }-\vec{k})).\end{array}$$

$$\begin{array}{rl}{f}_B(\omega ,k)& =\frac{1}{-i{k}_0+c(v)\left(|k_x|+|k_y|\right)+M}\frac{1}{i(k_0-\omega )+c(v)\left(|k_x|+|k_y|\right)+M}\\ & \times \frac{i}{\omega }\frac{6{\pi }^{2}v}{16}\int d{k}^{\prime }\frac{k_0-k_0^{\prime }}{\sinh \frac{\beta (k_0-k_0^{\prime })}{2}}\sum _{n=1}^4\sum _{m=\pm }\delta \left(k_0^{\prime }-k_0+e_n^{(m)}({\vec{k}}^{\prime }-\vec{k})\right)f_B(\omega ,k^{\prime }).\end{array}$$

$$\displaystyle\int dk_{0}^{\prime}\,dk_{x}^{\prime}\,dk_{y}^{\prime}\Bigg{[}\frac{i}{\omega}\,\frac{3\,v}{2^{5}\pi}\frac{k_{0}-k_{0}^{\prime}}{\sinh\frac{\beta(k_{0}-k_{0}^{\prime})}{2}}\sum_{\pm}\Bigg{(}\delta(k^{\prime}_{0}-k_{0}\pm(k_{y}^{\prime}-k_{y}))+\delta(k^{\prime}_{0}-k_{0}\pm(k_{x}^{\prime}-k_{x}))\Bigg{)}$$

$$\displaystyle-\delta(k^{\prime}-k)\left[-ik_{0}+c(v)\left(|k_{x}|+|k_{y}|\right)+M\right]\left[i(k_{0}-\omega)+c(v)\left(|k_{x}|+|k_{y}|\right)+M\right]\Bigg{]}f_{B}(\omega,k^{\prime})=0,$$

$$\displaystyle\mathcal{M}_{3\text{D}}(\omega)=\Bigg{[}\Delta k^{\prime}_{x}\Delta k^{\prime}_{y}\,\frac{i}{\omega}\frac{3\,v}{2^{5}\pi}\frac{k_{0}-k_{0}^{\prime}}{\sinh\frac{k_{0}-k_{0}^{\prime}}{2}}\Bigg{(}\delta_{k^{\prime}_{0},k_{0}+k_{y}^{\prime}-k_{y}}+\delta_{k^{\prime}_{0},k_{0}-k_{y}^{\prime}+k_{y}}+\delta_{k^{\prime}_{0},k_{0}+k_{x}^{\prime}-k_{x}}+\delta_{k^{\prime}_{0},k_{0}-k_{x}^{\prime}+k_{x}}\Bigg{)}$$

$$\displaystyle\qquad\qquad-\delta_{k^{\prime}_{0},k_{0}}\,\delta_{k^{\prime}_{x},k_{x}}\,\delta_{k^{\prime}_{y},k_{y}}\left[-ik_{0}+c(v)\left(|k_{x}|+|k_{y}|\right)+\bar{M}\right]\left[i(k_{0}-\omega)+c(v)\left(|k_{x}|+|k_{y}|\right)+\bar{M}\right]\Bigg{]}.$$

$$\mathcal{M}_{2\text{D}}(\omega)=\Bigg{[}\Delta k_{x}\,\Delta k_{y}\frac{i}{\omega}\frac{3\,v}{2^{4}\,\pi}\Bigg{(}\frac{k_{y}-k_{y}^{\prime}}{\sinh\frac{k_{y}-k_{y}^{\prime}}{2}}+\frac{k_{x}-k_{x}^{\prime}}{\sinh\frac{k_{x}-k_{x}^{\prime}}{2}}\Bigg{)}-\delta_{k^{\prime}_{x},k_{x}}\,\delta_{k^{\prime}_{y},k_{y}}\left[c(v)\left(|k_{x}|+|k_{y}|\right)+\bar{M}\right]\left[-i\omega+c(v)\left(|k_{x}|+|k_{y}|\right)+\bar{M}\right]\Bigg{]}.$$

$${\lambda }_L^{(B)}(v,\Lambda ,T)=Th(\bar{\Lambda })\sqrt{w(v)}{\left(\log \frac{1}{w(v)}\right)}^{1/4},$$

Full text

I Introduction

Thermalization in closed quantum systems Deutsch (1991); Srednicki (1994); Tasaki (1998); **Rigol **et al. (2008) is an important topic in quantum dynamics that has recently received much attention **Langen **et al. (2013); **Jurcevic **et al. (2014); **Kaufman **et al. (2016); **Schreiber **et al. (2015); **Choi **et al. (2016); **Kondov **et al. (2015); **Meldgin **et al. (2016); **Smith **et al. (2016). Thermalizing quantum many-body systems are “non-integrable", i.e. they do not have a number of conserved quantities that is of the order of the large number of degrees of freedom in the system.

A candidate quantitative measure of physical processes that destroy integrability and lead to thermalization is the degree of “many-body chaos" or “scrambling" as computed from certain special correlators Larkin and Ovchinnikov (1969); Shenker and Stanford (2014); Kitaev (2014). These correlators measure the growth of the commutator of two local operators in time. Given operators $A,B$ of a quantum many-body system with a Hamiltonian $H$ at inverse temperature $\beta$, the correlator at time $t$ is defined by

[TABLE]

As opposed to the usual thermal trace, the thermal density matrix has been split up into two parts and inserted between the two commutators in order to avoid short-distance divergences of the continuum field theory, thereby “regularizing" the correlator. The operators $A$ and $B$ are chosen such that $f_{A,B}(t)$ starts out very small at $t=0$. At late times, $f_{A,B}(t)$ will saturate to an ${\mathcal{O}}(1)$ value. At intermediate times, for non-integrable systems, the out-of-time-ordered (OTO) terms in $f_{A,B}(t)$ are expected to lead to a rapid growth $f_{A,B}(t)\sim\epsilon\,e^{\lambda_{L}^{(A,B)}\,t}$ Kitaev (2014). The exponent $\lambda_{L}^{(A,B)}$ is called a “Lyapunov exponent" and generically quantifies the non-integrability of the quantum many-body system 111There have been recent works showing that there are some exceptions to this Gopalakrishnan (2018); **Gopalakrishnan **et al. (2018). The small parameter $\epsilon\ll 1$ enables the exponential growth to be defined as $f_{A,B}(t)$ grows from an $\mathcal{O}(\epsilon)$ value at $t=0$ to an $\mathcal{O}(1)$ value at later times. Certain many-body quantum systems with a large number, $N$, of local degrees of freedom have $\epsilon\sim\mathcal{O}(1/N)$ even if $A,B$ are spatially close to each other. However, in systems with an ${\mathcal{O}}(1)$ number of local degrees of freedom, a large spatial separation between local operators $A,B$ is required to have a small $\epsilon$ Kitaev (2014); Stanford (2016); Patel and Sachdev (2017); Chowdhury and Swingle (2017); **Werman **et al. (2017); **Klug **et al. (2018); Jian and Yao (2018); **Patel **et al. (2017); **Aleiner **et al. (2016); Banerjee and Altman (2017); **Khemani **et al. (2018); Xu and Swingle (2018).

The largest Lyapunov exponents are expected to occur in many-body quantum systems without quasiparticles, since the complete lack of particle coherence in any basis of operators implies the non-existence of a thermodynamic number of conserved quantities and strong non-integrability. The upper bound of $\lambda_{L}=2\pi k_{B}T/\hbar$ (from here on we set $k_{B}=\hbar=1$) **Maldacena **et al. (2016) is saturated by black holes Shenker and Stanford (2014) and certain Sachdev-Ye-Kitaev models of strange metals Sachdev and Ye (1993); Kitaev (2014). On the other hand, Fermi liquids, for example, have $\lambda_{L}\sim T^{2}$ **Aleiner **et al. (2016); Banerjee and Altman (2017). Several works so far have computed $\lambda_{L}$ for various weakly and strongly coupled systems Kitaev (2014); Stanford (2016); Patel and Sachdev (2017); Chowdhury and Swingle (2017); **Werman **et al. (2017); **Klug **et al. (2018); Jian and Yao (2018); **Patel **et al. (2017); **Aleiner **et al. (2016); Banerjee and Altman (2017). However, $\lambda_{L}$ has not yet been computed analytically in a strongly coupled system without a large number of flavors.

In this work we study many-body quantum chaos at the antiferromagnetic (AFM) quantum critical point (QCP) of a two-dimensional metal, believed to exist in a wide range of layered compounds **Helm **et al. (2010); **Hashimoto **et al. (2012); **Park **et al. (2006). It is described by the spin-fermion model Abanov and Chubukov (2000); **Abanov **et al. (2003); Abanov and Chubukov (2004); Metlitski and Sachdev (2010); Abrahams and Wolfle (2012); Sur and Lee (2015); Maier and Strack (2016); **Berg **et al. (2012); **Schattner **et al. (2016); **Wang **et al. (2016); **Lunts **et al. (2017); **Schlief **et al. (2018, 2017), in which the fluctuations of the AFM collective bosonic mode scatter electrons between pairs of “hot spots”. The stable low energy fixed point of this model was found in Ref. **Schlief **et al. (2017). This strongly coupled fixed point occurs near perfect nesting of the paired hot spots, and deviations from it are perturbative in an emergent dynamical control parameter, $w$, which is a ratio of natural parameters of the theory. The electrons are perturbed away from free particles only by marginal corrections of strength $w\,\log(\Lambda/T)$ (where $\Lambda$ is a UV cutoff), whereas the bosonic collective excitation is made incoherent and non-local by strong renormalization from the electrons.

We compute two Lyapunov exponents at this fixed point, for the scrambling of both the fermion and boson operators. We do this perturbatively in $w$, without a large number of flavors. Since the interaction strength is ${\mathcal{O}}(1)$, one might expect the two scrambling rates to be the same and that both scale as $\lambda_{L}\sim{\mathcal{O}}(1)\,T$. However, the strong coupling is felt disproportionately by the boson, as the dressing of the fermion is of the order of the small parameter $w$. We find that the fermion exponent $\lambda_{L}^{(F)}\sim{\mathcal{O}}(w^{2})\,T$ is much smaller than the boson exponent $\lambda_{L}^{(B)}\sim{\mathcal{O}}(\sqrt{w})\,T$ (both relations are up to logarithms in $w$ and $T$). Their large difference is a reflection of the fact that the degrees of renormalization of the two sectors are wildly different, even though at low energies they interact strongly. In this sense, the Lyapunov exponents here act as measures of the “degrees of integrability” of the two sectors of the theory. However, contrary to the intuition about systems with strong interactions, the larger exponent $\lambda_{L}^{(B)}$ is parametrically smaller than ${\mathcal{O}}(1)\,T$. This is therefore a critical point where the effective coupling is strong (${\mathcal{O}}(1)$), and the boson is heavily renormalized at low energies, but the degree of non-integrability is determined by a parameter other than the effective coupling, and the Lyapunov exponents are not close to the upper bound.

Since the boson scrambles much more strongly, we additionally compute the spatial structure of its scrambling, by separating out the two boson operators in Eq. (1) at large distances. We find that because of the non-local nature of the boson propagator, boson operators spread exponentially fast through the system. In particular, the front of this “chaotic wave” is not a linear light-cone, and a finite “butterfly velocity” cannot be defined. This is similar to systems with non-local interactions, where Lieb-Robinson bounds Lieb and Robinson (1972) break down Hastings and Koma (2006); **Gong **et al. (2014); **Foss-Feig **et al. (2015).

II Model

The spin-fermion model Abanov and Chubukov (2000); **Abanov **et al. (2003); Abanov and Chubukov (2004); Metlitski and Sachdev (2010); Abrahams and Wolfle (2012); Sur and Lee (2015); Maier and Strack (2016); **Berg **et al. (2012); **Schattner **et al. (2016); **Wang **et al. (2016); **Lunts **et al. (2017); **Schlief **et al. (2018, 2017) describes a system of two-dimensional electrons interacting with a collective spin excitation at a finite ordering wavevector. This interaction is strongest at certain “hot spots” on the Fermi surface that are connected by the ordering wavevector. For simplicity, we take the case of a square lattice with one band, and the wavevector to be $(\pi,\pi)$ (this is most similar to the situation in the cuprates). The Gaussian tree level action has been extensively studied, so instead we start by writing the effective Euclidian action for the low-energy fixed point of this model found in Ref. **Schlief **et al. (2017) at finite temperature,

[TABLE]

Here, $k=(\omega_{k},\vec{k})$ denotes the fermionic Matsubara frequency and the two-dimensional momentum $\vec{k}=(k_{x},k_{y})$. $\psi_{n,\sigma}^{(m)}$ are the fermion fields that carry spin $\sigma=\uparrow,\downarrow$ at the hot spots labeled by $n=1,2,3,4,\leavevmode\nobreak\ m=\pm$.

With this choice of axis the ordering wave vector is $\vec{Q}_{AFM}=\pm\sqrt{2}\pi\hat{k}_{x},\pm\sqrt{2}\pi\hat{k}_{y}$ up to the reciprocal lattice vectors $\sqrt{2}\pi(\hat{k}_{x}\pm\hat{k}_{y})$. See Fig. 1 for details. The fermion dispersions are given up to linear order by $e^{\pm}_{1}(\vec{k};v)=-e^{\pm}_{3}(\vec{k};v)=vk_{x}\pm k_{y}$, $e^{\pm}_{2}(\vec{k};v)=-e^{\pm}_{4}(\vec{k};v)=\mp k_{x}+vk_{y}$, where $\vec{k}$ is the momentum deviation from each hot spot. The curvature of the Fermi surface can be ignored, since the patches of Fermi surface connected by the ordering vector are not parallel to each other for $v\neq 0$, i.e. the problem is fully two-dimensional. The component of the Fermi velocity along the ordering vector has been set to unity by rescaling $\vec{k}$. $v$ is the component of Fermi velocity that is perpendicular to $\vec{Q}_{AFM}$. It controls the degree of nesting between coupled hot spots. A necessary criterion for the validity of Eq. (2) is that $v\ll 1$, i.e. the fermions are close to perfect nesting. This is because, as explained below, a power of $v$ acts as a control parameter for the theory. $\vec{\phi}(q)$ is the boson field with three components which describes the AFM collective mode with frequency $\omega_{q}$ and momentum $\vec{Q}_{AFM}+\vec{q}$. $\vec{\tau}$ represents the three generators of the $SU(2)$ group. Due to the irrelevance of all local (in space and time) boson terms, there is a freedom to re-scale the boson field (the coefficients of the non-local terms are generated from the fermions and are not independent parameters). This freedom is used to set the Yukawa coupling between the collective mode and the electrons to $\sqrt{\pi v/2}$. The non-local dynamics of the boson is generated from the leading order contributions to the Schwinger-Dyson equation, which are shown in Fig. 2.

The ‘velocity’ $c(v)$ of the strongly damped boson is given by

[TABLE]

to the leading order in $v$. $M(T,\Lambda,v)$ is the thermal mass with $\Lambda$ being the momentum cutoff in Eq. (2). The leading order in $v$ non-zero contribution to $M(T,\Lambda,v)$ is computed in Appendix A.

Although the Yukawa coupling, originally denoted as $g$, is scaled to be $\propto\sqrt{v}\ll 1$, the theory is strongly coupled. The actual strength of interactions, or effective coupling, is given by $g^{2}/v$ Abanov and Chubukov (2000), which after our rescaling is ${\mathcal{O}}(1)$. The strength of interactions is what determines the prefactor of the leading kinetic term for the boson, which comes from the RPA correction and is given by $\left|{\omega_{q}}\right|$. By choosing this scaling, the boson kinetic term is of the same order as the fermion one, $i\omega_{k}$. Despite this, the leading order momentum dependence of the boson comes with a small prefactor of $c(v)\sim\sqrt{v}$. This dependence states that at this fixed point the boson motion is entirely due to the motion of fermions in the direction perpendicular to $\vec{Q}_{AFM}$. Of course, the generated momentum dependence is still larger than that of the bare kinetic term, $c_{0}^{2}\left|{\vec{q}}\right|^{2}/\tilde{\Lambda}$, below momenta of order $q\sim c(v)\tilde{\Lambda}/c_{0}$, where $c_{0}$ is the bare velocity and $\tilde{\Lambda}$ is a large UV scale.

The action in Eq. (2) obeys the $z=1$ critical scaling. All operators are marginal (this is called the interaction-driven scaling Sur and Lee (2014); **Schlief **et al. (2018)). It turns out that all quantum corrections to Eq. (2) are controlled by

[TABLE]

A renormalization group analysis shows that $v$ flows to zero with decreasing energy scale 222A similar ratio of velocities was found to act as a control parameter for the low energy fixed point of the theory of nematic ordering in $d$-wave superconductors Huh and Sachdev (2008). Therefore, there exists a basin of attraction around $v=0$ of finite size where this low energy fixed point is stable. In our work we assume a bare value of $v$ within this basin of attraction and small enough to satisfy Eq. (4).

There are several points to make here regarding the energy scales present. First, since the interaction in Eq. (2) is marginal, the perturbative in $w(v)$ renormalization of the fermion propagator leads to marginal Fermi liquid behavior at energies/temperatures low enough for quantum corrections to become important. Therefore, the renormalizations of the collective mode and fermion are highly asymmetric: the fermion renormalizes the collective mode to a highly incoherent and non-local form at relatively high energies, i.e. below the cutoff for the theory in Eq. (2), while the feedback onto the fermion remains weak down to much lower energies. In between these two scales is a superconducting transition temperature. Finally, the flow of $v(\mu)\rightarrow 0$ with decreasing energy scale $\mu$ will introduce an additional energy/temperature dependence into the propagators, beyond the usual logarithmic renormalization corrections. However, this flow is logarithmic in $\mu$, and therefore in order to see $v$ change significantly from its bare value $v_{0}$ the energy scale $\mu$ must be extremely small, in particular much less than the scale at which the fermions lose coherence. These various energy scales are shown in Fig. 3.

We are interested in the large energy window above the superconducting transition temperature and below the cutoff for the fixed point theory. Therefore, in particular, we can ignore the flow of $v(\mu)$ and treat it as $v=v_{0}$.

III Computing many-body chaos

We are interested in computing the many-body chaos (or “scrambling") in both the fermion and boson sectors of the theory. This requires the computation of two correlators of the type in Eq. (1): one with fermion and one with boson creation/annihilation operators. They are given by

[TABLE]

These correlation functions are generated from an action defined on a complex-time contour that has two real-time folds separated by $i\beta/2$. We compute $f(t,\vec{x})$ perturbatively in $w(v)$ by solving a Bethe-Salpeter equation for the Fourier transform $f(\omega,\vec{p})$ that gives a resummation of an infinite number of ladder diagrams of different kinds, and then Fourier transforming back Stanford (2016); Patel and Sachdev (2017); **Patel **et al. (2017); Chowdhury and Swingle (2017). Because of the unique time contour, the rungs in the Bethe-Salpeter equations are correlation functions between different time folds and are called Wightman functions, given by

[TABLE]

where $A_{n,m}(k)$ and $B(q)$ are the fermion and boson spectral functions, respectively. The rails in the Bethe-Salpeter equation are the standard retarded Green’s functions, with leading order in $w(v)$ self-energy corrections,

[TABLE]

and the advanced Green’s functions which are simply the complex conjugate of the retarded ones, $G^{A}_{(n,m)}=G^{R\,*}_{(n,m)},D^{A}=D^{R\,*}$. The fermion self-energy at finite $T$ is computed in Appendix B, and we note that to leading order it is independent of $\vec{k}$. To get Eq. (10) we simply analytically continue the propagator from Eq. (2) to real frequency, since the boson self-energy to leading order in $w$ is already given in Eq. (2). There are no other types of propagators to consider, since the interaction vertices in the expansion are only placed on the real time folds and not on the thermal circle, as it is believed those corrections will not change the spectrum of growth exponents Stanford (2016).

IV Chaos of the fermions

We start by computing the spatially averaged correlator for the fermion, $f_{F}(t)\equiv\int d^{2}x\,f_{F}(t,\vec{x})$. We first decompose the Fourier transform $f_{F}(\omega)$ as

[TABLE]

where $dk\equiv\frac{d^{3}k}{(2\pi)^{3}}$. We study the Bethe-Salpeter equation for each component, and then integrate over $k$ to get back $f_{F}(\omega)$. The Bethe-Salpeter equation for $f_{F}^{(n,m)}(\omega,k)$ to leading order can be expressed as

[TABLE]

or diagrammatically in Fig. 4.

The first diagram in the series is the zero-rung diagram, which includes only self-energy corrections. The second diagram is the leading order rung contribution.

We examine the zero-rung diagram separately. It is given by the sum over hot spots of products of a retarded and advanced Green’s functions. For a given hot-spot, to leading order, we can set $v=0$ in $e_{n}^{m}(\vec{k})$. The divergent integral over the momentum with velocity component $v$ contributes to the normalization of $f_{F}(t)$, but does not affect it’s growth rate. We can then do the integral over the other momentum component, noting that the quantity $\log(\Lambda/M)\gg 1$ dominates the location of the pole. This gives

[TABLE]

The Fourier transform of this contribution decays exponentially in time.

The exponential growth of $f_{F}(t)$ must come from the rung diagram and we can therefore neglect the zero-rung diagram, which gives

[TABLE]

This equation is equivalent to the statement that we can add a rung to $f_{F}(\omega,k)$ without changing the long-time behavior of $f_{F}(t)$ Stanford (2016). The kernel $K_{F}$ is given by

[TABLE]

Without loss of generality, we focus on a single hot spot pair, $(1,\pm)$, for which the coupled set of equations is (for ease of notation we omit the ’$1$’ from the superscripts and subscripts)

[TABLE]

As for the first diagram, we can approximate $e_{1}^{\pm}\approx\pm k_{y}$ and integrate over $k_{x}$. The equation for $f_{F}^{(\pm)}(\omega,k_{0},k_{y})\equiv\int\frac{dk_{x}}{2\pi}f_{F}^{(\pm)}(\omega,k)$ is given by

[TABLE]

Since we expect that $f_{F}$ depends on $k_{y}^{\prime}$ with no small prefactor, to leading order we can set $c(v)|k_{y}-k_{y}^{\prime}|\rightarrow 0$, which enables us to integrate over $k_{y}$ in the same way as we did for the first diagram. The equation for $f_{F}^{\pm}(\omega,k_{0})\equiv\int\frac{dk_{y}}{2\pi}f_{F}^{\pm}(\omega,k_{0},k_{y})$ is entirely independent of the hot spot index and so we remove it,

[TABLE]

We can see from the scaling that $\lambda_{L}^{(F)}\propto T$ up to logarithmic corrections, and we can scale out the temperature. We convert Eq. (18) to a matrix equation of the form $\mathcal{M}(\omega)f_{F}(\omega)=0$,

[TABLE]

where $\bar{\Lambda}\equiv\Lambda/T$ and $\bar{M}\equiv M/T$. To find the fastest exponential growth of $f_{F}(t)$ we look for the eigenvectors $f_{F}(\omega,k_{0})$ of $\mathcal{M}(\omega)$ with the largest eigenvalue, for $\omega$ on the positive imaginary axis. This is done by discretizing $k_{0}$ in Eq. (19) and diagonalizing the resulting finite matrix numerically. The details of the numerical solution are outlined in Appendix C. We find that to the leading order, which is ${\mathcal{O}}(w)$, the Lyapunov exponent for the fermion is zero,

[TABLE]

Here, the ${\mathcal{O}}(w^{2})$ term is understood to be up to logarithms in $w$ and $T$. Computing $\lambda_{L}^{(F)}$ to ${\mathcal{O}}(w^{2})$ involves the computation of higher-order self-energy diagrams **Schlief **et al. (2017) and the higher order rung diagrams shown in Appendix D, and is beyond the scope of this work. The reason for this null result is that the self energy contribution exactly cancels that coming from the kernel $K_{F}$. Although we don’t understand the reason for this completely, it seems to be a product of the structure of the theory, the linear fermion dispersion, and the momentum-independence of the self-energy at leading order Patel and Sachdev (2017). Because the Lyapunov exponent vanishes to the order we are working at, we cannot compute the spatial dependence of $f_{F}(t,\vec{x})$ in a controlled way.

V Chaos of the bosons

For the boson, we are able to compute the full spatiotemporal dependence of $f_{B}(t,\vec{x})$ to the leading order in $w(v)$. We start by computing $f_{B}(t)=\int d^{2}x\,f_{B}(t,\vec{x})=\int\frac{d\omega}{2\pi}e^{-i\omega t}\frac{1}{3}\int dk\,f_{B}(\omega,k)$. The Bethe-Salpeter equation for $f_{B}(\omega,k)$ to leading order is shown diagrammatically in Fig. 5,

and is given by

[TABLE]

The zero-rung contribution is trivial to compute and gives an exponential decay of $f_{B}(t)\propto e^{-2\,M\,t}/(c(v)^{2}t^{2})$. Therefore, as for $f_{F}$, it can be ignored when looking for the exponential growth of $f_{B}(t)$. The rung contribution has no minus sign coming from the fermion loop because two of the legs are Wightman functions. We note that there is no crossed diagram because of the time ordering properties of the expansion. The kernel $K_{B}$ is computed in Appendix E and is given by

[TABLE]

Ignoring the zero-rung term we have

[TABLE]

We convert Eq. (23) to a matrix equation $\mathcal{M}(\omega)f_{B}(\omega)=0$,

[TABLE]

where we have set $v=0$ in the fermion dispersions, since it is a subleading contribution. As for the fermion, the fastest exponential growth of $f_{B}(t)$ will be given by the largest eigenvalue of $\mathcal{M}(\omega)$ for $\omega$ on the positive imaginary axis. We have to find this eigenvalue numerically, by discretizing the integration variables. We scale out temperature and work with dimensionless units of frequency and momentum. The matrix $\mathcal{M}_{3\text{D}}(\omega)$ to diagonalize is given by

[TABLE]

By solving this equation for small system sizes we can see that the eigenvector $f_{B}(i\lambda_{L}^{(B)})$ of the largest eigenvalue is almost a delta function peak at $k_{0}=0$. We can therefore use the ansatz $f_{B}(\omega,k)=g_{B}(\omega,\vec{k})\delta(k_{0})$. Then, we can integrate both sides of the equation over $k_{0}$, which uses the delta functions to fix $k_{0}-k_{0}^{\prime}$. This gives an integral equation for $g_{B}(\omega,\vec{k})$, for which the new matrix is

[TABLE]

$\mathcal{M}_{2\text{D}}(\omega)$ is in one lower dimension, and is therefore computationally more manageable. The details of the numerical computation are in Appendix F. The scaling dictates that $\lambda_{L}^{(B)}\propto T$ up to logarithmic corrections. From analyzing the numerically obtained solution $\lambda_{L}^{(B)}(v,\bar{\Lambda})/T$ ($\bar{\Lambda}\equiv\Lambda/T$), we can see that for a fixed $\bar{\Lambda}$ the dependence on $v$ is well fit by

[TABLE]

with $h(\bar{\Lambda})$ plotted in Fig. 6.

We can fit the small and large $\bar{\Lambda}$ regions to an inverse logarithmic form, $\frac{\alpha}{\log(\bar{\Lambda})+\beta}-\gamma$. From extrapolating to the infinite cutoff (here the cutoff of the integral is taken to be different from $\bar{\Lambda}$, see Appendix F for details) and thermodynamic limits we find the coefficients to be

[TABLE]

The result of Eq.(27) is larger than the expected ${\mathcal{O}}(w)$ by a factor of $w^{-1/2}$. This enhancement comes from the thermal mass, which plays an important role: without it $\lambda_{L}^{(B)}$ would be infinite. This is qualitatively similar to Ref. Stanford (2016), where at criticality the thermal mass changed both the scaling with coupling and with temperature of the Lyapunov exponent. However, here the scaling of $\lambda_{L}^{(B)}$ with small $M$ is not as straightforward. Despite this enhancement, $\lambda_{L}^{(B)}$ still vanishes in the small $v$ limit, which indicates that as the effective velocity of the boson $c(v)\rightarrow 0$, the ${\mathcal{O}}(1)$ Landau damping is not enough to induce chaos.

V-1  Propagation of chaos in space

To compute the spatial dependence of $f_{B}(t,\vec{x})$ we need to solve the same Bethe-Salpeter equation as in Fig. 5, but with an external momentum $\vec{p}$ injected into the correlation function, $f_{B}(\omega,\vec{p},k)$. By the same reasoning as for $\vec{p}=0$ (see Appendix F) we arrive at the resulting two-dimensional matrix,

[TABLE]

Finding the largest eigenvalue of $\mathcal{M}_{2\text{D}}(\omega,\vec{p})$ for small $\vec{p}$ gives us the deviation

[TABLE]

where $\vec{p}=(p\cos(\theta),p\sin(\theta))$. We note that $\lambda_{L}^{(B)}(p,\theta)$ doesn’t develop an imaginary part. The exponent $\alpha$ varies between $1<\alpha<2$, starting out close to one for larger $v$ and monotonically approaching two as $v$ approaches zero. The coefficient $a$ grows monotonically as $v\rightarrow 0$, eventually scaling as $a(v)\sim v^{-1/4}$ when $v$ is small enough (“small enough” depends on the value of $\bar{\Lambda}$). The dependence of $\alpha$,$a$ on $\theta$ is weak. The dependence of $\alpha$,$a$ on $\bar{\Lambda}$ is logarithmic, and in the range we study $10<\bar{\Lambda}<10^{5}$ there is not much of a change. The plots of $a(v,\bar{\Lambda},\theta)$ and $\alpha(v,\bar{\Lambda},\theta)$ as functions of $v$ for various values of $\bar{\Lambda}$ are shown in Fig. 7 for $\theta=0$. More plots detailing the $\theta$ dependence of the parameters are shown in Appendix F.

The full form of the leading time dependence of $f_{B}(t,\vec{p})$ is $f_{B}(t,\vec{p})\sim e^{(\lambda_{L}^{(B)}-a\,p^{\alpha})\,t}\tilde{f}_{B}(t,\vec{p})$, where $\tilde{f}_{B}(t,\vec{p})$ is the eigenvector of the dominant eigenvalue, and we have again set $T=1$. The $\vec{p}$ dependence of the eigenvector in certain situations is known to have singularities for complex $\vec{p}$ that modify the real space structure of scrambling at large $|\vec{x}|$ and $t$ **Gu **et al. (2017); Gu and Kitaev (2019); **Guo **et al. (2019). Namely, there could exist a complex $\vec{p}$ for which $\lambda_{L}^{(B)}(\vec{p})$ is at its maximal value of $2\pi$ and $\tilde{f}_{B}$ is singular. However, Ref. **Guo **et al. (2019) found that (i) such a singularity in $\tilde{f}_{B}$ does not occur for the OTO correlation functions of the form in Eq. (6) that we compute (“retarded" OTO correlators in the language of Ref. **Guo **et al. (2019)), and, independently, (ii) if $\lambda_{L}^{(B)}(\vec{p}=0)\ll 2\pi$, as is the case in this paper, the singularity in the amplitude $\tilde{f}_{B}$, if it occurs, would occur at large values of complex $\vec{p}$ that would render the contribution to the Fourier transform of $f_{B}$ from those $\vec{p}$ severely suppressed by the exponential factor for most values of $\vec{x},t$. Therefore, it is safe for us to disregard the $\vec{p}$ dependence of $\tilde{f}_{B}$.

Computing the Fourier transform of the exponential factor then gives

[TABLE]

where $\theta_{x}$ is the angle of $\vec{x}$. Since the dependence of $a$ and ${\alpha}$ on $\theta$ is weak, we can ignore it, as the qualitative $\left|{\vec{x}}\right|$ dependence will not be affected. Integrating over $\theta$ in this approximation gives,

[TABLE]

where $J_{0}$ is the zeroth Bessel function of the first kind, and we have scaled out $a(v,\bar{\Lambda})\,t$ from the $p$ integral. This last integral is hard to compute analytically, but from numerics we can see that for $\left|{\vec{x}}\right|\gtrsim 10\,(a\,t)^{1/\alpha}$ the integration gives an inverse power-law scaling of

[TABLE]

The proportionality constant is fairly flat until $\alpha$ is extremely close to $2$, and since we are studying small but finite $v$ we treat the proportionality constant as effectively $\alpha$-independent.

From Eq. (33) we can define a typical “operator radius” for the boson operator, which is defined as the $R(t)$ for which $f_{B}(t,R(t))\sim 1$ (here we forget about the angular dependence as it is weak). This is the radius within which the initially local operator at $\vec{x}=0$ has spread, and is given by

[TABLE]

The form of Eq. (34) tells us how fast the scrambling of the boson operators spreads through the system. For large $R(t)$, which we are interested in, the “scrambled” region grows exponentially with time, as illustrated in Fig. 8.

Since the boundary of this region cannot be linearized, we cannot define a finite butterfly velocity. The reason for this strange behavior is due to the highly non-local nature of the boson propagator: $D^{R,-1}(t,\vec{x})=\pi^{2}e^{M\,t}(1+(\frac{x_{1}}{c(v)\,t})^{2})(1+(\frac{x_{2}}{c(v)\,t})^{2})$. The algebraic decay of $D^{R}(t,\vec{x})$ with $\left|{\vec{x}}\right|$ for any non-zero $t$ means that the boson is highly delocalized, and therefore an initial perturbation at the origin propagates with “infinite” velocity. As $c(v)\sim\sqrt{v}\rightarrow 0$, the non-locality in the propagator (slowly) disappears and instead the boson becomes completely localized in space. In the limit $v\rightarrow 0$, $\alpha(v)\rightarrow 2$ and Eq. (34) is modified to $R(t)\sim 2\sqrt{a\,\lambda_{L}^{(B)}}\,t$. In this case, after restoring factors of $v_{F}$, the butterfly velocity would be $v_{B}^{(B)}=2\sqrt{a\,\lambda_{L}^{(B)}}\,v_{F}$, which is independent of $v$ since the $v$-dependence from $a$ and $\lambda_{L}^{(B)}$ would cancel. This again makes sense, since the only way a completely localized boson can propagate is by fermionic particle-hole production, which travel at a speed of $v_{F}$. We note, though, that the complete localization of the boson for $c(v)\rightarrow 0$ is not real, since at extremely small $c(v)$ the bare term $c_{0}^{2}\left|{\vec{q}}\right|^{2}/\tilde{\Lambda}$ would control the boson dynamics.

VI Discussion

We have computed the Lyapunov exponents that describe the growth rate of both fermion and boson OTO correlators. We worked at the leading perturbative order in the control parameter $w$, and naively it seemed that both exponents might be $\sim{\mathcal{O}}(w)$. Instead, it turns out that the boson scrambles significantly faster than the fermion, $\lambda_{L}^{(B)}\sim{\mathcal{O}}(\sqrt{w})\,T\gg\lambda_{L}^{(F)}\sim{\mathcal{O}}(w^{2})\,T$. This large discrepancy stems from their different degrees of renormalization. The reason for $\lambda_{L}^{(F)}$ being smaller than ${\mathcal{O}}(w)$ seems to be accidental, and not a generic feature of Yukawa-type field theories. The enhancement of $\lambda_{L}^{(B)}$, however, is not accidental and is due to the irrelevance of the bare terms in the bosonic propagator. However, even though the Landau-damped frequency dependence of the boson kinetic term is ${\mathcal{O}}(1)$, the degree of chaos still vanishes in the $v\rightarrow 0$ limit. This indicates that strong Landau damping alone is not enough to rapidly scramble the system. Even though this fixed point has an effective coupling of order unity that controls the low-energy boson dynamics, the degree of non-integrability of the boson is still proportional to a positive power of $w$, and its scrambling rate is not close to the maximal one.

We also find that the boson scrambling spreads in space via a logarithmic light cone and an infinite butterfly velocity. Initially, one would expect that this is at odds with the Lieb-Robinson bound Lieb and Robinson (1972) on the spread of information in quantum systems. However, these features have been seen before in systems with long-range interactions Hastings and Koma (2006); **Gong **et al. (2014); **Foss-Feig **et al. (2015), which violate the assumptions set forth in the Lieb-Robinson bound. Here, the effective action of Eq. (2) violates those same assumptions because of the non-local boson propagator, and our results are consistent with recent findings.

Acknowledgments

We thank Igor Aleiner, Antoine Georges, Martin Claassen, Vadim Oganesyan and Subir Sachdev for useful discussions. The Flatiron Institute is a division of the Simons Foundation. AAP was supported by the US Department of Energy under Grant No. DE-SC0019030, and by a Harvard-GSAS Merit Fellowship.

Appendix A Thermal mass

Here we compute the thermal mass $M(T,\Lambda,v)$ to the leading non-vanishing order in $w(v)$. At the leading order of ${\mathcal{O}}(w\log(w))$ it was computed in Ref. **Schlief **et al. (2017) and found to vanish. However, we expect is to have a non-zero value at $O(w)$. The only contribution at this order comes from the diagram in Fig. 2. Is it given by

[TABLE]

where

[TABLE]

with $\omega_{p},\omega_{k}$ being bosonic and fermionic Matsubara frequencies, respectively, and the zero temperature polarization is the straightforward analogue. We first change variables to $k_{+}=e_{n}^{(m)}(\vec{k}),k_{-}=e_{n}^{(\bar{m})}(\vec{k})$ (which has a Jacobian of $1/2v$). The integrations over $k_{\pm}$ are done via poles. The Matsubara sum over $\omega_{k}$ is straightforward and is equal to the integral at zero temperature. This bring us to

[TABLE]

where $\varepsilon(\vec{p},M)=c(\left|{p_{x}}\right|+\left|{p_{y}}\right|)+M$ and $\varepsilon(\vec{p})=c(\left|{p_{x}}\right|+\left|{p_{y}}\right|)$. The Matsubara summation and the frequency integration can be done directly, which gives

[TABLE]

where $\psi(z)=\Gamma^{\prime}(z)/\Gamma(z)$ is the digamma function. The thermal mass is proportional to temperature (up to logarithms), as can be seen from scaling $T$ out of Eq. (A4). Since the only two scales are $T$ and $\Lambda$, the ratio $M(T,\Lambda,v)/T$ is a function of $v$ and $\bar{\Lambda}=\Lambda/T$ only. We can therefore set $T=1$ and compute $\bar{M}(v,\bar{\Lambda})\equiv M(v,\bar{\Lambda})/T$.

We solve Eq. (A4) numerically for $\bar{\Lambda}\geq 10$ and $v\leq 10^{-2}$. For a fixed $\bar{\Lambda}$ the scaling with $v$ is hard to fit systematically. We therefore fix $v$ and fit the resulting curves as functions of $\bar{\Lambda}$. These curves are shown in Fig. A1.

From fitting the various regions and values at their junctions we find that

[TABLE]

The functions $a_{i}(v)$ and $b_{i}(v)$ are well estimated by

[TABLE]

For a fixed $\bar{\Lambda}$, as we take $v$ to be smaller, $\bar{M}(v,\bar{\Lambda})$ changes between the three different regimes indicated in Eq.(A5), and therefore the scaling with $v$ is not uniform for all $v$.

Appendix B Fermion self-energy at finite temperature

The leading order contribution to the fermion self-energy is independent of momentum. This is because the momentum dependence carries an extra factor of $c$ Sur and Lee (2015); **Lunts **et al. (2017); **Schlief **et al. (2017, 2018). The temporal part of the fermion self-energy after scaling out $1/c$ is given by

[TABLE]

where $\omega_{n}=2\pi Tn$ with $n\in\mathbb{Z}$, and $k_{0}=2\pi T(m+\frac{1}{2})$ with $m\in\mathbb{Z}$. Since we do not expect any logarithmic divergences in $v$, we can set $w=c=0$ inside the integral, and do the $p_{y}$ and $p_{x}$ integrals

[TABLE]

where $\Lambda$ is the momentum cutoff. The infinite sum can be simplified

[TABLE]

where we have used the notation $\bar{k}_{0}\equiv|\frac{k_{0}}{2\pi T}|-\frac{1}{2}$. Up to now we could compute everything exactly. Now we take the limit of large momentum, $\Lambda\gg\left|{k_{0}}\right|+M$. We then have

[TABLE]

Since $\frac{M}{2\pi T}\rightarrow 0$ as a power of $v$ (c.f. Appendix A), we can neglect it to leading order,

[TABLE]

We convert back to the frequency $k_{0}$ and get

[TABLE]

Appendix C Numerical calculation of $\lambda_{L}^{(F)}$

Here we give some details on the numerical calculation of $f_{F}$ from Eq. (19). We convert the integral into a discrete sum and introduce a cutoff $\bar{\Lambda}_{0}$ for the frequency summation. The equation becomes of the form $\mathcal{M}(\omega)f_{F}(\omega)=0$, where $\mathcal{M}$ is a finite matrix given by

[TABLE]

where $k_{0},k_{0}^{\prime}\in[-\bar{\Lambda}_{0},\bar{\Lambda}_{0}]$. We sweep values of $\omega$ on the positive imaginary axis, and plot the eigenvalue with the smallest magnitude, $\left|{E_{0}}\right|$. The Lyapunov exponent is the largest value of $\lambda_{L}^{(F)}=-i\omega$ for which $\left|{E_{0}}\right|=0$. For a large enough $\bar{\Lambda}_{0}$ the value of the integral doesn’t change any more. We take $\bar{\Lambda}_{0}=15$ and for each other parameter we decrease $\Delta k_{0}$ until the integral converges. We find that $\lambda_{L}^{(F)}\sim v\sim w^{2}$ up to logarithms. This means that the $O(w)$ terms in Eq. (C10) cancel. The nonzero result is due to the thermal mass $M(T,\Lambda,v)$. However, our calculation is not controlled up to $O(w^{2})$, since we have not computed higher order rung corrections to the Bethe-Salpeter equation and higher order self-energy terms. Therefore, our result is null at the order up to which we have control.

Appendix D Higher-order graphs for $f_{F}$

Here we consider the higher order contributions to $f_{F}$ shown in Fig. D2 and argue that they are at least of order ${\mathcal{O}}(w^{2})$ up to logarithms.

The three kernels are given by

[TABLE]

where $\tilde{K}_{F}$, $\bar{K}_{F}$ and $\hat{K}_{F}$ refer to diagrams (a), (b) and (c) in Fig. D2, respectively, and we have excluded numerical factors, since we are interested in the power of $v$ that the kernels scale with. Note that $\tilde{K}_{F}$ connects Green’s functions at hot spot $(n,m)$ to $f_{F}^{(n^{\prime},m^{\prime})}$ at all others hot spots $(n^{\prime},m^{\prime})$, while the other two only connect the same hot spots. It is easy to check that once $1/c$ is scaled out of the $\vec{q}$ variables the integrals are UV finite (the exponential decay with large $q_{0}$ from the Wightman functions makes them even more UV safe than the equilibrium perturbative corrections). This implies the kernels come with a factor of $\sim v^{2}/c$. The kernels $\tilde{K}_{F}^{(n,m),(n^{\prime},m^{\prime})}$ that involve fermions from nearly perpendicular hot spots are finite even without scaling out $1/c$ and are therefore suppressed even further **Schlief **et al. (2017); **Lunts **et al. (2017). After the $q$ integrations, the $\vec{k}^{\prime}$ integration will give another factor of $1/c$, as in the lowest order rung diagram in Fig. 4. Since all matrix elements of $\tilde{K}_{F},\bar{K}_{F},\hat{K}_{F}$ are then parametrically smaller than those of $K_{F}$, the resummation of the insertions in Fig. D2 will only make a contribution to $\lambda_{L}^{(F)}$ of $\sim O(v^{2}/c^{2})\sim O(w^{2})$.

Appendix E Calculation of $K_{B}$

The kernel $K_{B}$ contains a fermion loop. Summing over the contributions from all the hot-spot pairs we get

[TABLE]

Here, the Green’s functions and the Wightman functions are of the free fermions, since self-energy corrections would be higher order than the order we are working at. We start with $K_{B}^{(1,+)}(k,k^{\prime},\omega)$,

[TABLE]

We change variables to $q_{+}=q_{x}+q_{y},\,q_{-}=q_{x}-q_{y}$ and do the integration over $q_{-}$ by using the first delta function,

[TABLE]

The integration over $q_{+}$ can be done via poles, and then the $q_{0}$ integration is trivial as well,

[TABLE]

The result is trivially extended to all other $K_{B}^{(n,m)}(k,k^{\prime},\omega)$, and we can sum all contributions to get

[TABLE]

Appendix F Numerical calculation of $\lambda_{L}^{(B)}$

Here we provide some more details about the numerical calculation of $\lambda_{L}^{(B)}$. We first note the transition from $\mathcal{M}_{3\text{D}}$ to $\mathcal{M}_{2\text{D}}$ is additionally justified by the fact that for a fixed $\Delta k_{x},\Delta k_{y}$ the largest $\lambda_{L}^{(B)}$ computed from $\mathcal{M}_{2\text{D}}$ is always larger than that of $\mathcal{M}_{3\text{D}}$. Since we are ultimately interested the $k$-integral of the fastest growing eigenvector, $\int dkf(i\lambda_{L}^{(B)},k)$, this implies that the ansatz we make is safe.

The ratio of momentum cutoff over temperature $\bar{\Lambda}=\Lambda/T$ appears in two different places in Eq. (26): as the largest momentum index, and as a parameter in $\bar{M}(\bar{\Lambda},v)$. The actual physical $\bar{\Lambda}$ can be quite large, but we don’t have the computational power to go to the necessary resolution except for modest $\bar{\Lambda}$. However, the dependence of the integral (or sum) on the cutoff is weak, i.e. the result is a series in the inverse cutoff, since the integral is convergent (we confirm this). On the other hand, the dependence on $\bar{\Lambda}$ from $\bar{M}(\bar{\Lambda},v)$ is logarithmic. Therefore, in order to find the Lyapunov exponent for large values of $\bar{\Lambda}$, we can treat the integration cutoff and $\bar{\Lambda}$ as separate. For each given $10\leq\bar{\Lambda}\leq 3\times 10^{5}$, we find the Lyapunov exponent for several (mostly) much smaller cutoffs and then extrapolate the answer to the infinite cutoff limit. Our exact extrapolation procedure is as follows: we first compute $\lambda_{L}^{(B)}$ for some fixed cutoff and momentum spacing $\Delta k$, and from Eq. (27) we extract the fit for $h(\bar{\Lambda})$ of the form in Eq. (28). Then holding the cutoff fixed we recompute the numerical coefficients of $h(\bar{\Lambda})$ for several decreasing values of $\Delta k$, and extrapolate these results to the $\Delta k\rightarrow 0$ limit using a linear fit. We then assume that the slope of that fit is independent of $\bar{\Lambda}$, which holds for the values we have checked, and perform the infinite cutoff extrapolation on the $\Delta k\rightarrow 0$ values of the numerical coefficients in $h(\bar{\Lambda})$, giving us the values in Eq. (28).

For a finite external momentum $\vec{p}$, Eq. (24) becomes (before setting $v=0$ in the dispersions)

[TABLE]

As before, we set $v=0$ in the dispersions and assume an ansatz of the form $f_{B}(\omega,\vec{p},k)=g_{B}(\omega,\vec{p},\vec{k})\,\delta(k_{0})$. The resulting two dimensional equation is given in Eq. (29). Solving it gives the deviation from the $\vec{p}=0$ Lyapunov exponent in Eq.(30). In Fig. F3 we include some plots of the weak $\theta$ dependence of $a(v,\Lambda,\theta)$ and ${\alpha}(v,\Lambda,\theta)$.