Pole skipping away from maximal chaos

Pole skipping is a recently discovered subtle effect in the thermal energy density retarded two point function at a special point in the complex $(\omega,p)$ planes. We propose that pole skipping is determined by the stress tensor contribution to many-body chaos, and the special point is at $(\omega,p)_\text{p.s.}= i \lambda^{(T)}(1,1/u_B^{(T)})$, where $\lambda^{(T)}=2\pi/\beta$ and $u_B^{(T)}$ are the stress tensor contributions to the Lyapunov exponent and the butterfly velocity respectively. While this proposal is consistent with previous studies conducted for maximally chaotic theories, where the stress tensor dominates chaos, it clarifies that one cannot use pole skipping to extract the Lyapunov exponent of a theory, which obeys $\lambda\leq \lambda^{(T)}$. On the other hand, in a large class of strongly coupled but non-maximally chaotic theories $u_B^{(T)}$ is the true butterfly velocity and we conjecture that $u_B\leq u_B^{(T)}$ is a universal bound. While it remains a challenge to explain pole skipping in a general framework, we provide a stringent test of our proposal in the large-$q$ limit of the SYK chain, where we determine $\lambda,\, u_B,$ and the energy density two point function in closed form for all values of the coupling, interpolating between the free and maximally chaotic limits. Since such an explicit expression for a thermal correlator is one of a kind, we take the opportunity to analyze many of its properties: the coupling dependence of the diffusion constant, the dispersion relations of poles, and the convergence properties of all order hydrodynamics.


Introduction and summary of results
There is a rich variety of phenomena and signatures associated to quantum chaos. Sharpening our understanding of one of these, the growth of operator complexity and its probe, the outof-time order correlation function (OTOC) has lead to great advances in the understanding of many-body chaos and its relation to gravity through AdS/CFT in recent years.
Recently, a novel signature of chaos was proposed, the pole skipping phenomenon [1,2]. While this effect was very convincingly demonstrated in [1][2][3][4][5][6] for theories saturating the chaos bound of [7], it was not understood what the fate of this effect was away from maximal chaos. In this paper, we propose a generalization of pole skipping away from maximal chaos. We provide a stringent test of the proposal by exact computations in a solvable chaotic theory, the large-q SYK chain [8,9].
We spend the rest of the introduction reviewing the behavior of OTOCs in large-N theories, introducing the pole skipping phenomenon, formulating it away from maximal chaos, and providing evidence for the proposal. Since our investigation led us to a precious simple closed form thermal correlator in an interacting chaotic theory, we close the introduction with this formula and a brief discussion of its properties.

Pole skipping at maximal chaos
In maximally chaotic theories in the regime β t, |x| t scr = O(log N ) the OTOC behaves as where the λ L = 2π/β is the maximal Lyapunov exponent and u B is the butterfly velocity that determines the boundary of the butterfly cone, the region in which the OTOC grows.
Recently the pole skipping phenomenon in the thermal energy density retarded two point function G R εε (ω, p) was discovered [1,2] and tested in maximally chaotic theories [1][2][3][4][5][6]. The energy density retarded two point function has a family of hydrodynamic poles on the complex ω plane as we vary p: ω pole (p) = ±c s p + . . . where the dispersion relation depends on whether the system conserves momentum or not. The discovery of [1] was that this family of poles goes through a special point determined by λ L and u B : (ω, p) p.s. = iλ L 1, 1 u B , maximal chaos (1.3) and at this precise location its residue vanishes. The proposed explanation of this pole skipping is that the growth of the OTOC comes from the same hydrodynamic mode that is responsible for energy transport [2]. In the OTOC, this pole gives rise to the exponential growth, while in the energy-energy correlator, the growth must be absent, and therefore the pole must be skipped. This prediction of pole skipping in the thermal energy density two point function for maximally chaotic system was verified in a wide class of holographic systems dual to Einstein gravity with a general matter content [3].
Hence it is natural to expect pole skipping to take place at the point (1.3) also away from maximal chaos. This would make pole skipping extremely exciting, as it would suggest that one can extract some data about chaos from a two point function, which is a much simpler observable than the OTO four point function. However, the naive conjecture (1.3) cannot be true. The most elementary way to see this, is to realize that while not all 2d CFTs are maximally chaotic, their retarded energy density two point function G R εε is universal [4], 1 where > 0 is infinitesimal. This expression exhibits pole skipping at the point (ω, p) p.s. = i (2π/β)(1, 1), which is unrelated to λ L in general. Note that u B = 1 in any 2d CFT [11]. We explain below how to modify (1.3) so that it holds away from maximal chaos.

Pole skipping away from maximal chaos
Away from maximal chaos (1.1) is replaced by where we introduced the velocity dependent Lyapunov exponent (VDLE) λ (u) [11][12][13]. The ordinary Lyapunov exponent λ L is obtained by setting u = 0: λ L = λ (0), and u B is defined by the edge of the growing region, λ (u B ) = 0. Note that while in the maximally chaotic case chaos is characterized by two numbers, λ L and u B , in the non-maximal case there is a whole function λ (u) worth of freedom. In terms of this function, the chaos bound translates to the pointwise bound λ(u) ≤ 2πβ −1 (1 − |u|/u B ) [11]. Importantly, in all examples we know λ (u) is a result of a competition between two terms, one coming from the "leading Regge trajectory" and another describing the contribution of the stress tensor. 2 For small |u| the Regge contribution dominates. For larger |u| ≤ u B there are two possible scenarios depending on the details of λ (u) [11]: In the first case the Regge contribution dominates for all 0 ≤ |u| ≤ u B . In the second case, above some critical velocity u * ≤ |u| ≤ u B the stress tensor dominates, and we have (1.6) The stress tensor contribution to chaos is determined by two numbers λ B . 3 In maximally chaotic theories one has u * = 0, i.e. the stress tensor dominates for all u, and we recover (1.1).
We are now ready to formulate our proposal. We conjecture that the pole skipping point in the retarded energy density two point function is at 4 (ω, p) p.s. = iλ i.e. it encodes the stress tensor contribution to chaos. While the stress tensor determines chaotic dynamics at maximal chaos, its contribution can be decreased in non-maximally chaotic theories or completely cancelled in integrable systems; an example of the latter is provided in [14]. The proposal is consistent with the literature on pole skipping in maximally chaotic theories, since in that case λ L = λ (T ) B . It is also consistent with the 2d CFT result discussed around (1.4).
Let us discuss the existing evidence for (1.7) in the literature. The main evidence comes from conformal field theories on Rindler space. Rindler space is a patch of flat spacetime but it is conformal to S 1 × H d−1 and therefore it is an example of thermal physics that can be studied using vacuum correlators in a conformal field theory. For example, the thermal energy density two point function can be obtained from the vacuum one via a Weyl transformation and therefore is universal to all CFTs. Pole skipping in this context was recently studied in [6] and it was shown that it happens at the values corresponding to maximal chaos. 5 Note that most CFTs on Rindler space do not display maximal chaos, instead the Lyapunov exponent is determined by the resummation of operators on the leading conformal Regge trajectory and is related to the analytic continuation of the spin of this trajectory to nonphysical operator dimensions [7,11]. However, for a large class of theories the stress tensor 3 This inequality follows from the concavity of the contribution of the leading Regge trajectory to the VDLE. In all the examples we know, this is true, and we suspect it holds generally. 4 This formula is valid in the presence of (possibly discrete) translational symmetry. 5 Since the spatial manifold H d−1 is curved, (1.7) must be interpreted so that it applies to the appropriate notion of Fourier transform. In flat space the Fourier mode at the poles skipping point takes the form x . In analogy with this result, the right interpretation in hyperbolic space is that for large geodesic separations, the Fourier mode at the pole skipping point must behave as exp − ρ is the spatial geodesic distance. contribution dominates near the butterfly front and u (T ) B = (d − 1) −1 gives the true butterfly speed; for example this happens in planar N = 4 SYM theory whenever the 't Hooft coupling is greater than 37.7384 [11], while the theory is only maximally chaotic at infinite coupling.
Another piece of non-trivial evidence comes from holographic theories with higher derivative corrections performed using shockwaves [7,[15][16][17]. 6 Such corrections do not affect the Lyapunov exponent but they change the butterfly speed. It was shown in [5] that both for Gauss-Bonnet coupling and the leading α 3 R 4 correction the change in the butterfly speed and the energy-energy correlator are such that pole skipping remains valid. This is consistent with (1.7) since the Lyapunov exponent stays maximal. In fact, when we consider finite coupling corrections to the thermal OTOC in N = 4 SYM 7 we need to include worldsheet effects that correct the Lyapunov exponent [19] in addition to the leading higher derivative term α 3 R 4 that corrects the butterfly speed. 8 The fact that in this case pole skipping still happens at maximal Lyapunov exponent is additional evidence for (1.7) (and contradicts the naive propsal (1.3)).
The main technical result of this paper is the confirmation of (1.7) in the large-q limit of the SYK chain introduced in [8], where we obtain both λ(u) and G R εε in closed form, and can prove (1.7) analytically. The model has two dimensionless coupling constants 0 < v < 1 and 0 < γ ≤ 1. The latter controls the strength of inter-site coupling, while the former controls the interaction strength; as a function of v the model interpolates between a free and a maximally chaotic theory. There exists a critical line in coupling space v * (γ), above which (in the regime v * (γ) ≤ v < 1) chaos is non-maximal, but the butterfly speed is maximal B are nontrivial functions of the couplings. 6 Shockwave computations only capture t-channel graviton exchange, but correctly compute the VDLE λ(v) in Einstein gravity. It was found in [18] that there are contributions to the four graviton S-matrix in higher derivative gravity theories that are not captured by the shockwave computations and that change the value of the Lyapunov exponent. It is an open problem to determine how these terms change the VDLE and vB. We thank Shiraz Minwalla and Douglas Stanford for a discussion on this point. 7 Note that here we put the theory on S 1 × R d−1 as opposed to the Rindler discussion of the previous paragraph. 8 The first correction to the Lyapunov exponent comes at order 1/ √ λ while to the butterfly speed at order 1/λ 3/2 , where λ is the 't Hooft coupling. The two dimensional coupling space of the model. There is a critical line v * (γ) that separates the blue region, where above a critical velocity u * the VDLE is dominated by the pole and is maximal, from the white region, where the saddle always dominates, and the VDLE is nowhere maximal.
We end this section with some comments about (1.7).
• The proposal is somewhat anticlimactic, as it shows that just from determining the pole skipping point, we cannot read off the true λ L and u B characterizing the OTOC. On the other hand, while λ (T ) L ≡ 2π/β does not carry any non-trivial information, u B still depends on the theory. Moreover, for a large set of strongly coupled but not maximally chaotic theories, one has u B = u (T ) B . In these theories, the OTOC butterfly speed can be read from the location of the pole skipping point, however, we caution that there is no way to tell whether u (T ) B is the true butterfly speed just by having access to the energy density two point function.
• The proposal (1.7) refers to the stress tensor contribution to the growth of the OTOC.
This can be defined in models where we can analytically compute the OTOC, but it is not entirely clear how to define it in complete generality. However, when there is a critical velocity u * < u B , we can read this data directly from the behavior of the VDLE near u = u B .
• There is strong evidence that u B ≤ u (T ) B is true generally [11], in which case we can read a bound on the butterfly effect from the pole skipping point using the VDLE bound of [11] (given in the text after (1.5)).
• It would be very interesting to prove (1.7) in full generality. This would likely require the development of an effective field theory framework for the four point functions of nonmaximally chaotic theories generalizing the works [2,20] that apply at maximal chaos. The theory would have to apply for large Lorentzian times for OTOC configurations, and would also have to govern energy transport at finite frequency and momentum.

The energy density two point function in the large-q SYK chain
A byproduct of our investigation is a closed form expression for the energy density two point function in the large-q SYK chain, see Sec. 2 for the definition of the model. It is given by the simple formula: where v is the SYK coupling constant, we chose the inverse temperature β = 2π (without loss of generality), and the function ψ n (θ) is a sum of two hypergeometric functions given explicitly by: where γ is the inter site coupling constant and p is the momentum. We believe this result is precious: this is the only closed form thermal correlator in a chaotic system, which changes nontrivially as we go from weak to strong coupling. The only analogous known expressions is (1.4) in 2d CFT (and similar results for higher dimensional CFTs in hyperbolic space [6]), which however is universal for any CFT integrable or chaotic, and is completely determined by symmetry. Beyond this, there exist results either for small coupling [21][22][23][24][25], or in the holographic setting corresponding to large coupling [21,[26][27][28][29]. We extract (1.8) from a complicated formula for the four point function by taking an OPE limit. The main technical result of the paper is that we find the unique analytic continuation of this Euclidean Green function defined at discrete Matsubara frequencies to the entire complex plane, which allows us to prove (1.7) and to study its analytic structure.
In more detail, is easy to verify that (1.8) exhibits pole skipping at (1.10) With the explicit simple formulas for λ(u), u B , λ The only singularities of the correlator (1.8) are poles. For small p, the closest pole to the real ω axis is a diffusion pole with dispersion relation ω(p) = −iDp 2 +. . . . As we increase p this pole can collide with a non-hydrodynamic pole on the imaginary ω axis, and subsequently the poles can depart the imaginary ω axis to the complex plane. While this scenario is reminiscent of what happens in theories with weakly broken momentum conservation [30], in our model these poles never become propagating quasiparticles (sound modes), as they always have comparable real and imaginary frequencies. We also show that besides a Drude peak (present for small momentum) the spectral function is featureless: it does not show the presence of quasiparticles.
The collision of poles delineates the range of momenta for which all order hydrodynamics converges [31,32]. We determine the radius of convergence by finding pole collisions and the reconnection of pole trajectories for complex momenta. Unlike in more familiar examples, here the radius of convergence of hydrodynamics decreases with increasing coupling v.

Outline
The outline of the paper is as follows. In Sec. 2 we review the SYK chain and its large-q limit, where the quantum many-body problem can be reduced to solving differential equations. In Sec. 3 we use the retarded kernel approach for computing λ(u). This is considerably simpler than computing the full four point function that we undertake in Sec. 4. In Sec. 5 we extract a simple formula for the Euclidean energy density two point function from the complicated formula for the four point function by taking an OPE limit. In Sec. 6 we find all the pole skipping points in the energy density correlator, and verify the proposal (1.7). In Sec. 7, we perform a detailed study of the analytic structure of this correlator, and extract information about the (all order) hydrodynamics of the SYK chain. In Appendix A, we give a new derivation of the coordinate space four point function of the large-q SYK dot first derived in [33,34] using two different methods.

Large-q SYK chain
We will be studying the SYK chain introduced in [8], which is a 1+1 dimensional generalization of the SYK dot [9,[35][36][37]. It is a disordered chain of Majorana fermions with a large N number of on-site degrees of freedom that is chaotic, yet solvable using methods resembling Dynamical Mean Field Theory. The fermions interact in groups of q, and by taking this parameter to be large, we can analytically solve the model for all values of the two remaining dimensionless coupling constants of the theory: v that controls the strength of interactions (and can be thought of as a proxy for temperature measured in units of the dimensionful coupling constant J ) and γ that controls the relative strength between on-site and inter-site interactions and that enters in all results in a simple kinematical way.
The Hamiltonian of the system is given by Here, χ i,x , i = 1, ..., N are Majorana fermions obeying the commutation relation {χ i,x , χ j,x } = δ xy δ ij and periodic boundary condition, χ i,0 = χ i,M . q is an even integer, and the first line of (2.1) is an on-site term, while the second line is an interaction between two groups of q/2 fermions on neighboring sites. J i 1 ...iq,x , J i 1 ...i q/2 j 1 ...j q/2 ,x are independent Gaussian random variables with zero mean. Their variances are given by 2) The SYK model is known to be self-averaging [35,36,38]. This means that while we are interested in the theory with quenched disorder, we can solve the technically much simpler problem of the theory with annealed disorder. That is, we average the generator of disconnected correlation functions instead of correlation functions themselves. After integrating out the disorder and the fermions, one obtains an effective action for bilocal fields This action is further simplified and becomes local (in the two time dimensions) in the large q limit, when N q 2 1. 10 The non-trivial limit is obtained by fixing J 2 0,1 ≡ qJ 2 0,1 /2 q−1 . In this limit, the Schwinger-Dyson equations for G x , Σ x suggest that the system is almost free, with interactions entering at order 1/q. Hence following [40], we proceed by changing variables according to ( expanding the action in 1/q and integrating out σ x . The resulting action at leading order in 1/q is a periodic lattice of coupled Liouville theories The field g x (τ 1 , τ 2 ) is symmetric under the exchange τ 1 ↔ τ 2 and satisfies the boundary condition g x (τ 1 , τ 1 ) = 0.

Linearized action
We are interested in fluctuations around the large N saddle point. Following [8], we look for a translation invariant saddle, We look for a finite temperature solution satisfying g(0) = g(β) = 0. We will set β = 2π for simplicity, this can be easily restored by dimensional analysis. The solution satisfying the boundary conditions is [9] e g s (τ ) = (2.7) We now expand around the saddle g x = g s + δg x so that the second order action becomes (2.8) 10 It would also be interesting to analyse the regime q 2 ∼ N using the methods of [39].
The inverse of the kernel in this bilinear action is related to the Euclidean time-ordered connected four point function F xy : (2.9) We proceed by Fourier transforming δg (from now on we set the lattice spacing to be 1) and taking the limit of an infinite lattice M → ∞ so that the momentum variable p = 2πk/M is continuous. We obtain (2.10) Note that each momentum sector is decoupled at leading order and thus we can evaluate the fermion four point function in the momentum sector p as where we used translational invariance to define (2.12) The boundary condition and symmetries of δg imply just as in the case of the single SYK dot.
The kernel of the quadratic action (2.10) reads as To further simplify the kernel, we introduce 2θ . Then the eigenvalue problem for the g p,n (X) reduces to the following equation (we separate the overall factor of N/q 2 for simplicity) The above kernel admits the following zero modes: (2.17) These eigenfunctions are even/odd under θ → π − θ. For integer n, these modes do not satisfy the boundary condition that the fluctuations must vanish at X = 0, 2π, or θ = π 2 (1 ± v). We can obtain eigenfunctions instead by replacing n by a non-integral m that is defined by requiring these boundary conditions to hold. In this case, the eigenvalue µ p,n in (2.16) is given by (n 2 − m 2 )/4 (this is the route taken in [34] for the SYK dot). This will turn out to be impractical below, instead, we will invert the kernel using the zero modes (2.17) with integer n.

Velocity dependent Lyapunov exponent
We would like to extract the velocity dependent Lyapunov exponent from the ladder kernel (2.14). The simplest way to read the leading growing contribution of a mode with momentum p to the four point function is to follow the retarded kernel method of [9,36]. The rule is to replace the sides of the ladder with retarded two point functions (which are now just Heaviside theta functions) and the rungs with left-right propagators (corresponding to the Lorentzian 11 For generic γ the range of h is smaller, and we only fill the range for γ = 1. continuation X → iχ + π in the two point function, or θ → iξ + π/2 in (2.16)). This gives rise to a Schrödinger problem. The maximal growth exponent in Lorentzian Y admitted by the kernel corresponds to the ground state energy −κ 2 of The ground state wave function of this problem is ψ(ξ) ∼ 1 (cosh ξ) h−1 , and the corresponding ground state energy gives rise to the Lyapunov exponent where the momentum dependence of h is defined via taking the 1 ≤ h ≤ 2 branch of (2.15). The relevant strong coupling limit of this result can be obtained by taking v = 1 − δv, p 2 = δvp 2 and expanding (3.2) in δv, which yields κ(p) = 1− 1 + γp 2 /3 δv+. . . , whose analog for q = 4 was derived in [8].
Here we obtain κ(p) for all values of the coupling and all momenta. The growing piece of the four point function in real lattice space x is then represented by an integral of the form where we have used the ladder identity of [41] relating the prefactor to the exponent. For large t, this integral is either dominated by a saddle point or a pole at κ(p) = 1, depending on x.
We associate the pole with the stress tensor contribution. 12 This mechanism is quite generic [8,11,12,19,41,42], and as explained in [11], leads to the velocity dependent Lyapunov exponent (we use u = x/t for velocity to avoid confusion with SYK coupling v) where the stress tensor butterfly speed u (T ) B and critical velocity u * are defined via solving Explicit expressions for these quantities read as Let us, for completeness, also write down the explicit formula for the Legendre transform part of the VDLE (first line of (3.4)) B , in which case the true butterfly speed (defined via λ(u B ) = 0) is the stress tensor one u B = u (T ) B , and the velocity dependent exponent is maximal for u > u * [11]. On the other hand, when u * > u (T ) B , one only has the first branch of (3.4) and the exponent is nowhere maximal. In this case, the true butterfly speed is determined by solving λ saddle (u B ) = 0. Note that due to concavity of λ saddle (u) B for all couplings. There is a critical line in coupling space marking the boundary between these two types of behaviors, given by equating the two velocities in (3.6). This is shown on Fig. 1.
It is perhaps educational to give explicit formulas in the strong coupling limit [11]: We see that u * goes to zero, while u B , we see that the theory is maximally chaotic, while if we take very small u's, we can detect the deviation from maximal chaos (that itself goes to zero as δv → 0). We close this section by noting that (3.3) only determines the t and x dependence of the OTOC, where t = i(τ 1 + τ 2 − τ 3 − τ 4 )/2 is the separation between the "center of mass" times, while x is the spatial distance between the two operator pairs (see (2.9)). However, the OTOC also depends on the relative times ξ = iv 2 (τ 1 − τ 2 + π), ξ = iv 2 (τ 3 − τ 4 + π). This dependence is encoded in the wave function factors in the retarded kernel eigenvalue equation (3.1) [9,20]. The full position dependence of the growing piece of the OTOC is therefore (3.9) In the chaos limit, the relative positions ξ, ξ remain bounded, so this extra dependence on them does not affect the steepest descent contour and the discussion in this section.

Four point function
In this section we invert the kinetic kernel of (2.10). Instead of starting from its spectral representation, we solve Green's equation with the appropriate delta sources on the right hand side. We switch to Matsubara representation using δ(Y − Y ) = 1 2π n e in(Y −Y ) , and expand so that Green's equation for the kerenel (2.14) reads as This equation is supplemented with the boundary condition We can find G n using the homogeneous solutions Ψ ± n that satisfy Ψ − n (θ v , v) = Ψ + n (π − θ v , v) = 0. We explain how to do this on Fig. 2. The resulting formula is Here the Wronskian is given by W = Ψ − n ∂ θ Ψ + n − Ψ + n ∂ θ Ψ − n and is independent of θ and is inserted to correctly normalize the delta functions on the RHS of (4.2). We can easily construct homogeneous solutions with the required properties using the zero modes (2.17): Since the Wronskian is independent of θ, we may evaluate it at θ = θ v that results in where prime denotes θ derivative. We recognize in the brackets the Wronskian between ψ e n and ψ o n which must be independent of θ v since they solve the same second order homogeneous ODE without first order terms. We may therefore replace θ v → π/2 in the square brackets, which (using (2.17)) then evaluates to one. In summary, we simply have It would be interesting to obtain a position space expression by summing (4.3). We did not manage to do this in general, however, note that for p = 0 (h = 2) the result must agree with the four point function in the case of the SYK dot [33,34]. We confirm that this is the case in Appendix A.

Prescription
Here we would like to compute the Matsubara amplitudes of the connected energy density two point function (recall that we take the lattice spacing to be 1 and set β = 2π) where we define the energy density as a local term in the Hamiltonian (2.1) The idea is to extract this correlator from an OPE limit of the time ordered four point function (TOC) . This method was used in [8] in the conformal limit, with a prescription inspired by the factorization of the p = 0 contribution to the TOC. We cannot rely on this here, away from the conformal limit. Instead we can find the right prescription using the observation of [43], i.e. that the Clifford anticommutation relations imply Note that the energy density is not a uniquely defined operator in a lattice model; our definition in (5.2) was chosen such that the above identity holds. Next, we represent the four point function (2.9) as a trace and plug in the formula χ i,x (t) = e iHt χ i,x (0)e −iHt for Heisenberg operators. We can derive the following identity: where we exploit that time derivatives are implemented by commutators with H, which in the OPE limit τ 1 → τ 2 , τ 3 → τ 4 can be combined to become two copies of the (Heisenberg evolved) energy density according to (5.3). A contact term appears due to the non-triviality of taking the coincident limit τ 2 − τ 4 → 0 first. We can write this formula in Fourier space in terms of (4.3) as In Fourier space, the possible contact terms are polynomials in n and cos(p). We will determine the right contact term in Sec. 5.4 after obtaining the analytically continued expression of the regular part of energy two point function.

Calculation
We proceed to evaluate (5.6). Using (4.3), (4.4), (4.6) we find that up to a possible contact term where we again used that the Wronskian of the two solutions is independent of θ and evaluated it at θ v → π/2. The other combination can be evaluated as (5.9) Therefore we find the Matsubara Green's function for the energy density + contact term n ∈ 2Z + 1 .

(5.10)
Next, we proceed to obtain the analytic continuation in n after which we determine the contact term.

Analytic continuation
In order to obtain retarded correlators we need to determine the analytic continuation of (5.10) from integer n to the complex n plane. This is unique due to Carlson's theorem once we eliminate the exponential growth in the n → ±i∞ directions. We can do this by defining a master function that unifies ψ e/o n of (2.17) for even/odd n and does not grow exponentially in the limit n → ±i∞ at θ = θ v . It turns out that the right master function is 13 ψ n (θ) = c o ψ o n (θ) + c e ψ e n (θ) , (5.12) For integer n, ψ n (θ) ∝ ψ e/o n (θ) depending on the parity of n, up to a θ independent prefactor that cancels in the ratio (5.10). We choose to introduce this prefactor so that c e/o can be entire functions in n; this will soon be useful. With the use of this master function, we may write the retarded correlator as One may confirm that this is analytic in ω in the upper half plane for real momentum (1 ≤ h ≤ 2), as it should be.

Contact term
Now we determine the contact term in the energy density two point function. One necessary condition on the contact term comes from the Ward identity corresponding to energy conservation which gives [44]: This condition however is not sufficient to fix the contact term, since it leaves the freedom of adding a p dependent contact term which vanishes at p = 0. We can fix this freedom by matching to the expected short time behavior of the energy correlator. Since we are analyzing a quantum mechanical model of Majorana fermions, which in the UV become asymptotically free (at leading order in N ), and the operator ε x (τ ) in (5.2) is a polynomial of fermions, but not their time derivatives, we expect that G R εε (t = 0 + , p) = finite . (5.15) 13 One may arrive at this by taking the v → 1 limit in which case ψ e/o in (2.17) are simple functions of θ, but involves certain gamma function prefactors. We obtain the answer for general v by replacing n → n/v inside the resulting gamma functions. This leaves us with an undetermined phase: in the sin πh 2 + πn 2v term in co in (5.12), the phase shift πh 2 is not fixed by this argument (but it must agree with the phase shift in the cos πh 2 + πn 2v term in ce). We guessed the right value based on matching with (rather high order) perturbation theory in δv ≡ 1 − v, where the analytic continuation is straightforward. This is only possible, if the correlator decays fast enough in frequency space. For convenience we take the UV limit in the Euclidean theory, which corresponds to taking ω → i∞ with 0 < v < 1 and the momentum (or equivalently h) fixed in (5.13). We get: The Fourier integral producing (5.15) will only converge, if the contact term cancels the constant in the above equation. This can indeed be done by a contact term, since h(h − 1) is a (first order) polynomial in cos(p), see (2.15). Fourier transforming to real lattice space, this expression produces a combination of δ(τ 1 − τ 2 )δ x,x and δ(τ 1 − τ 2 )δ x,x ±1 , which are contact terms in both time and space. Getting rid of these terms, we get the complete expression of the retarded energy density two point function as announced in (1.8): 6 Pole skipping

Pole skipping points
We are interested in points where zero and pole lines of (5.17) meet, which are called pole skipping points. The necessary condition of pole skipping is that the denominator and the numerator of the first term in the parenthesis of (5.17) goes to zero simultaneously Let us denote this equation as A c = 0. Since det A = W(ψ o n , ψ e n )(θ) = 1, pole skipping can only happen when c e = c o = 0. Inspecting (5.12) we find that there are six classes of solutions (i) cos nπ THus we see that pole skipping points are only possible for real-valued h and so it is natural to divide into the following three cases {p ∈ iR, p ∈ R, p ∈ (2Z + 1)π + iR}. We note that only the first case has pole skipping points on the upper half ω plane, while all three cases have additional pole skipping points on the lower half ω plane.
• Purely imaginary momentum p ∈ iR: h ≥ 2 Here pole skipping happens both on the upper and the lower half of the complex frequency plane. Let us focus on the pole skipping points on the upper half plane n > 0 with h ≥ 2. These points can only arise from (i), (iii) in (6.2) and can be further simplified to The pole skipping point connected to the diffusion pole is the one at n = 1 and h = 1 + 1 v , we show this on Fig. 3. Our modified pole skipping conjecture (1.7) is that pole skipping on this pole line happens at n = 1 and p = i/u (T ) B , which indeed translates to h = 1 + 1/v based on our discussion of the OTOC in Sec. 3. Therefore, the conjecture (1.7) indeed holds exactly in the SYK chain at any coupling.
Furthermore, we note that there are additional pole skipping points on the lower half plane n < 0 corresponding to (ii), (iv), (v), (vi) and the remaining solutions of (i), (iii) in (6.2) satisfying h ≥ 2. The complete collection of pole skipping points for imaginary momentum can be observed in Fig. 4, where we see that pole skipping points on the lower half plane exist for both integer and non-integer values of n, whereas on the upper half plane, all pole skipping points have integer n.
• Real momentum p ∈ R: 1+ For real p the retarded Green's function can only have poles (and hence pole skipping points) in the lower half plane. In this case, the list of pole skipping points coming from (i), (ii), (iii), (iv) of (6.2) simplifies into the following expression which corresponds to integer n Furthermore (v) and (vi) of (6.2) generate pole skipping points at h = 2 (equivalently p = 0) and generically non-integer n given by We plot the pole skipping point for the real momentum case on Fig. 5.
Finally, there is one more class of pole skipping points which happens on the line of complex momentum whose real part lies on the edge of the Brillouin zone Re(p) = (2Z + 1)π. Similarly to the real momentum case, pole skipping only happens on the lower half of the complex frequency plane and (i), (ii), (iii), (iv) of (6.2) correspond to the following unified expression    B . Therefore they do not contribute to the growth of the OTOC, except for m = 0. We show this in Fig. 6.

Pole skipping on the lower half plane
In the discussion above, we have found that the large q SYK chain has pole skipping points on the lower half plane. These points can be divided into two classes: one class involves negative integer Matsubara frequencies, while the others are at non-integer values of n frequencies.
Note that both classes occur both at real and complex momenta p. While pole skipping points on the lower half plane cannot contribute to an exponential growth of the OTOC, we remark that the existence of pole skipping at non-integer multiples of the unit Matsubara frequency is interesting in its own right.
Pole skipping in the thermal energy two point function at negative integer multiples of the unit Matsubara frequency was discovered in [45] in the context of holography. The authors found from the near horizon expansion that at special values of (ω, p), the quasinormal modes of linearized Einstein gravity are not unique, implying that at these frequencies, the holographically dual retarded Green's function [26] is indefinite, which is another indication of pole skipping. In the discussion of [45], the universal structure of the black hole metric and the near horizon expansion always seems to yield such locations at negative Matsubara frequencies ω = −i 2πn β , n ∈ N (this is also robust under higher derivative corrections [46]). It would be interesting to see whether non-integral valued pole skipping points like (6.5), (6.7) van be understood in the context of holography by finding singular quasinormal modes at non-integral frequencies.
7 Hydrodynamics and analytic structure of the two point function The energy density retarded two point function determines the linear response behavior of the system. If we take |ω|, |p| T , we are in the hydrodynamic regime, and from the pole closest to the origin, we can determine the hydrodynamic transport coefficients. We find that for all values of the coupling the transport of energy is diffusive, and is controlled by a diffusion pole.
We find that the two point function is meromorphic: it only has poles, but no branch cuts. We investigate the motion of the hydrodynamic and non-hydrodynamic poles on the complex ω plane as we change p. We already analyzed this problem partially: it is the continuation of the diffusion pole line to p ∼ T that participates in pole skipping (1.7). Here we ask, whether poles collide as we change p; the collision between the diffusion pole and a non-hydrodynamic pole delineates the applicability of the (all order) hydrodynamic expansion. These phenomena have been thoroughly analyzed in the planar four-dimensional N = 4 super Yang-Mills (SYM) theory both at weak and at strong coupling [21, 26-28, 31, 32, 47]. Our system allows us to perform the analysis at all values of the coupling, and we comment on similarities and differences between the SYK model and SYM theory.

Diffusion
We can extract the diffusion constant by examining the ω = in 1, p 1 → h ≈ 2 limit of (5.13). In this limit one has which by using the relation between h and p in (2.15) leads to the diffusion constant: The strong coupling limit of this result is D = γ 6 δv with δv ≡ 1 − v, whose analog for q = 4 was derived already in [8]. (See also the discussion around (3.2).) Note that D is an increasing function of the coupling v, unlike in the more familiar cases of field theories, where the diffusion constant diverges at weak coupling.
Note that one has D ≤ (u B ) 2 when we reinstate the temperature) for all values of the couplings v and γ, with saturation at strong coupling v → 1. This is consistent with (in fact stronger than) the diffusivity bound of [48]. One may also examine the bound in terms if the true butterfly speed when it is not given by u (T ) B . In these cases, u B can be determined numerically by equating (3.7) to zero. It turns out that D ≤ u 2 B can be violated at weak coupling, but we get a correct bound by dividing with the Lyapunov exponent, that is, D ≤ u 2 B /v is always true, again consistently with [48]. One may confirm this analytically in the weak coupling limit, where we can solve for the true butterfly speed where e is the natural number, while D = 1 6 γv + O(v 2 ) in this limit. Note that the pole line giving rise to the diffusion pole is the same as the pole line participating in pole skipping at (1.7).

Movement of poles
Let us concentrate on the movement of poles for real p first. To explore the full range of h ∈ [1, 2] as we move around in the Brillouin zone p ∈ [−π, π], we set γ = 1. All the plots in this section can be easily converted to any value of γ using the relation p γ = arccos cos(p here )−1+γ γ that follows from (5.13). Relatedly, for γ < 1 we only cover part of the possible h range.
From Fig. 5, we can already start building intuition about the movement of pole trajectories (hot lines). In Fig. 7 we plot the dispersion relations of the first few poles for representative values of v as we go from weak to strong coupling. We see that at weak coupling the first few poles stay at pure imaginary ω. At stronger coupling, we encounter collisions of poles, which is followed by a gap in momentum, after which another pair of poles appears. Next we seek to understand the details of pole collisions. It turns out that after the collision the poles move off to the complex plane. We illustrate this on Fig. 8, by examining the location of poles on the complex ω plane for the two values of momenta marked by gridlines on the v = 0.8 figure of Fig. 7. As the dashed lines show on Fig. 8, as we increase p the pair of poles return to the imaginary axis and continue to live there until we reach the boundary of the Brillouin zone.
The fact that poles never wander too far from the imaginary axis implies that the spectral function ρ(ω, p) = 2ImG R εε (ω, p) is rather featureless: there are no quasiparticle peaks emerging even at weak coupling, see Fig. 9. Besides the usual sharp Drude peak for small fixed momentum (at small frequency), which is the signature of the hydrodynamic diffusion pole close to the real ω axis, a noteworthy feature is that for larger momentum and sufficiently strong coupling, the peak can shift away from zero frequency. This can be observed on the right panel of Fig. 9. This happens when the retarded correlator has a zero on the negative imaginary ω axis that is closer to the real axis than the first pole. Note that this is very different behavior than what one finds in SYM theory slightly away from infinite coupling [49].  Next, following [31,32], we ask about the convergence radius of all order hydrodynamics. As explained in [31], we can use the analytic implicit function theorem to find the radius of convergence. This amounts to finding a point on the diffusion pole line n(p) such that ∂ n ψ n (θ v )| n(pc),pc = 0 . (7.4) We plot the solution p c to this equation on Fig. 10. We find that the radius of convergence |p c | increases as we decrease the coupling. 14,15 For strong coupling v c < v < 1, where v c = 0.6875 ± 0.0005, the critical momentum p c is real, while for v < v c it develops an imaginary part. This corresponds to the fact that we have a dispersion relation for the diffusion mode that has finite p support for v > v c , while it reaches the edge of the Brillouin zone for v < v c , see Figs. 7 and 5. There is another distinguished value of the coupling which we denote by v c = 0.241 ± 0.005. For v < v c we have Re p c = π, therefore the region of convergence of all order hydrodynamics contains the entire Brillouin zone. Note that the convergence radius shrinks to zero size as we go to maximal coupling v → 1. To gain more intuition about how this happens, we sketch the analog of Fig. 7 for very strong coupling in Fig. 13. We see that the diffusive hydrodynamic pole with D diverging as v → 1 is intersecting with a non-hydrodynamic pole that is sitting at iω = O(1) that has mild momentum dependence. This clearly leads to a decreasing p c as we increase the coupling. coupling, the reconnection always happens at real p, and we can read off the convergence radius of all order hydrodynamics from Fig. 7. Figure 11: Motion of the hydrodynamic and a non-hydrodynamic pole on the complex ω plane for v = 0.8 at fixed p m and varying phase of the momentum φ. Left: For p m = 1.297 the poles sit on the imaginary axis of the complex ω plane (horizontal plane) close to each other for φ = 0, π. As we increase φ (vertical axis) from 0 to 2π the projection of their trajectory onto the ω plane goes around a circular curve twice, since the function only depends on p 2 . Note that the vertical axis represents a periodic direction, so the upper and lower end planes are identified. Right: For p = 1.35 the two curves reconnect, and we wind around the projection of the curve to the ω plane once as we change the phase from 0 to 2π. There is also another curve whose plot coincides with this one, but has phase shifted by π. Insets: The insets zooms onto the point where the reconnection takes place, and are marked with red cuboids on the large figures.
At weaker coupling, we find a pole collision of a different flavor. Instead of the reconnection of pole trajectories, they are unlinked at small p, and become linked at a complex p c plotted on Fig. 10, without the reorganization of each individual curve. We show this on Fig. 12.

The strong coupling limit
For fixed h, n we can take the strong coupling limit corresponding to small δv ≡ 1 − v, and obtain the following simple formula valid for 2 > Re(h) > 1/2: + O δv 2 , δv 2h−1 .

(7.5)
Note that the terms in the first line are analytic in ω; they play the important role of cancelling the singular contribution of the second line at h = 3/2, where the powers δv and δv 2h−2 coincide. We can simply verify pole skipping on this function, and we plot the motion of the poles for real momentum on Fig. 13. This plots is helpful in that it explains where some features we saw on Fig. 7 originate from. Note however that this formula misses the important diffusion pole, as that happens for 2 − h = O(δv). If we take another limit, we can recover the results of [8] including the diffusion pole, but we miss the higher poles discussed above. To analyze the hydrodynamic limit (p, ω → 0) of the retarded correlator at strong coupling, we rescale the momentum as p 2 =p 2 δv and expand (5.13) in terms of δv (see also the discussion around (3.2)). Together  Figure 13: The dispersion relation of the poles of (7.5) is drawn with blue lines for real momenta. We also include a sketch of the diffusion pole with an orange line, which is not captured by the formula (since the pole degenerates to the vertical axis). Compare with the last plot of Fig. 7.
with a first order approximation h = 2 − γ 6 p 2 and the diffusion constant D = γ 6 δv we get This matches the result obtained in [8] up to the contact terms that they did not keep track of. We note that we do not see the poles that were plotted in Fig. 13, as they have subleading residues compared to the diffusion pole. This leads to the curious conclusion that the convergence radius of hydrodynamics in this scaling limit is infinite, as opposed to the result of Fig. 10 that took into account the collision with the aforementioned subleading poles. Another consequence is that the dispersion relation of the pole, ω(p) = −iDp 2 is an entire function of p, and according to the result of [51] this implies the relation D = (u (T ) B ) 2 that we found in the strong coupling limit Sec. 7.1. 16

Analytic structure at weak coupling
Here we wish to analyze (1.8) in the weak coupling limit v → 0. In analogy with the strong coupling limit, we will consider two different scalings of the parameters of the Green function.
Let us first consider fixed real momentum, h ∈ [1, 2]. Already from Fig. 7 it is visible that the poles on the negative imaginary axis become denser as we decrease v, and in the limit v → 0 they proliferate. We did not find a nice limiting function, in particular the accumulation of poles does not seem to form a branch cut.
Since the density of poles on the imaginary axis diverges in the weak coupling limit, it is reasonable to consider another limit, in which we define the rescaled imaginary frequencyñ = n/v and keep this fixed as v → 0. If we do this, the hypergeometric functions 2 F 1 (a, b, c; z) in (1.9) can be approximated by 1, since their z argument goes to 0, while their a, b, c arguments remains finite. We get the simple approximation: Similarly to (7.5), this expression is useful, as the movement of poles and pole skipping can be easily analyzed.

Assorted comments on the SYK chain
At the level of the fermion two point function, the SYK chain is ultralocal in space, i.e. χ i,x (τ )χ i,x (0) ∝ δ x,x . It is locally critical [8], since in the vacuum the correlator decays as a power law, and at small temperature it takes the form dictated by finite temperature conformal quantum mechanics (2.7). 17 If we move on to the study of the energy density, we saw that our state of matter exhibits diffusive energy transport. 18 Next, we ask, if our system has a continuum limit. We have been working with dimensionless positions and momenta, x and p. If we worked with dimensionful versions, y ≡ xa, k ≡ p/a, we would find that the diffusion constant is (for βJ → 0) .

(7.8)
The only meaningful way to take the continuum limit of the system is to implement the double scaling limit β a → ∞ , βJ → ∞ , D phys β = πγ 6 (βJ ) a β 2 = fixed . (7.9) We find that for all other O(1) values of the dimensionless coupling constant βJ , i.e. all values of v not infinitesimally close to 1, we get a lattice scale diffusion constant. Let us now examine the propagation of chaos in the model. The physical butterfly velocity corresponding to the spatial coordinate y is u B,phys = 2π a β u B , which in the βJ → ∞ (v → 1) limit takes the value u B,phys = π 2 γ 3 (βJ ) a β 2 , (7.10) which is exactly the same combination as appearing in (7.9), and hence remains finite in the continuum limit. This is how it had to be since the remarkable relation D phys = β 2π u 2 B,phys has to be satisfied in the strong coupling limit as discussed in Sec. 7.1.
In the continuum limit the theory is maximally chaotic, with λ(u) = 2π β 1 − |u| u B , see (3.8). The energy density Green function given by (7.6) is extremely simple: it only has a pole with dispersion relation ω = −iD phys k 2 and an infinite radius of convergence in k.
That the continuum limit is so simple, clearly demonstrates that momentum is not an approximately conserved quantity at any scale in the SYK chain. While the collision of the diffusion and a non-hydrodynamic pole on the imaginary ω axis shown on Fig. 8 is reminiscent of the scenario articulated in [30], whereby at large k the collision of poles create two sound modes with ω = ±c s k +. . . , here the poles only depart the imaginary ω axis for a short while, they always have comparable real and imaginary parts, and return to the imaginary ω axis for larger values of the momentum. At fixed momentum, tuning the coupling from weak to strong has a dramatic effect in SYM theory: the closely spaced branch cuts at zero 't Hooft coupling [21] break up into families of poles that form multiple branches of a "Christmas tree" at strong coupling [55], only for the top branch to remain at infinite coupling. The motion of poles in our model is less rich, but fully calculable and should complement the recent studies of the analytic structure of thermal correlators at small but finite coupling [22][23][24][25]: changing the coupling at fixed p leads to occasional collisions of poles that move them out to the complex ω plane in some window of v, and sometimes pole skipping happens, when a line of zeros intersects with a line of poles.
A Four point function for p = 0 As mentioned in the main text, the four point function (4.3) must agree with that of the SYK dot for p = 0 (h = 2). Here we confirm this by doing the n sum. For h = 2 the zero modes (2.17)  The advantage of scaling out the n dependence from the denominator is that now these modes depend polynomially on n apart from the Fourier modes in θ. That is, one can write Φ ± n (θ) = p ± + (n, θ)e i n 2v θ + p ± − (n, θ)e −i n 2v θ , (A. 3) where p ± ± (n, θ) are degree two polynomials in n. Now (4.3) is invariant under rescaling the modes, so we have the expression N q 2 G 0,n (θ 1 , θ 2 ) = v 3 cot πv 2 8 n 2 (n 2 − v 2 ) (−1) n Φ − n (θ 1 , v)Φ + n (θ 2 , v)Θ(−θ 1 + θ 2 ) where we used the Wrosnkian Using the property (A.3), we may pull out the polynomial dependence on n as derivatives and write the position space expression as N q 2 n e inY G 0,n (θ 1 , θ 2 ) = 8v 3 cot πv 2 Θ(π − θ 1 − θ 2 ) × a,b=± and the sum for F (Y ) was done by Sommerfeld-Watson resummation. Since all the p ± a in (A.6) are degree two polynomials in the derivatives, it is easy to now explicitly evaluate (A.6) and confirm that it agrees with the results in [33,34].