Energy diffusion and the butterfly effect in inhomogeneous Sachdev-Ye-Kitaev chains

We compute the energy diffusion constant $D$, Lyapunov time $\tau_{\text{L}}$ and butterfly velocity $v_{\text{B}}$ in an inhomogeneous chain of coupled Majorana Sachdev-Ye-Kitaev (SYK) models in the large $N$ and strong coupling limit. We find $D\le v_{\text{B}}^2 \tau_{\text{L}}$ from a combination of analytical and numerical approaches. Our example necessitates the sharpening of postulated transport bounds based on quantum chaos.

Introduction 1 A few years ago it was noted that many experimentally realized "strange metals" seem to be characterized by Drude "transport time" of order τ * ∼ /k B T [1]. As this time scale was proposed to be the "fastest possible" time scale governing quantum dynamics [2,3], it was conjectured by Hartnoll [4] that these strongly interacting strange metals remained metallic due to a fundamental bound on diffusion: D v 2 τ * . In strongly interacting systems without quasiparticles, many-body quantum chaos provides a natural velocity v and time scale τ * for such a diffusion bound [5,6]. In some large-N models, it is natural to define a Lyapunov time τ L and butterfly velocity v b as [7,8] V (x, t)W (0, 0)V (x, t)W (0, 0) for rather general Hermitian operators V and W . τ l is an analogue of the Lyapunov time from classical chaos: it describes the rate at which quantum coherence is lost. Similarly, v b is called the butterfly velocity: it governs the speed of chaos propagation, and is somewhat analogous to a state-dependent Lieb-Robinson velocity [9]. Since we now have a velocity and time scale which can be defined in any strange metal, [5,6] noted that D was naturally related to v 2 b τ l in some simple holographic settings. A natural question to ask is whether Hartnoll's conjecture becomes Originally (2) was observed to hold for charge diffusion constant [5], with theory-dependent O(1) prefactors, but there are now multiple known counterexamples [10,11,12]. It is more compelling that (2) hold for the energy diffusion constant: as argued in [11,13,14], v b characterizes the loss of quantum coherence, a process related to quantum "phase relaxation" which should also characterize energy fluctuations and diffusion. Furthermore, additional evidence for an energy diffusion bound of the form (2) has arisen in holographic models in the low temperature limit [6,15], and models of Fermi surfaces coupled to gauge fields [13]: at weak coupling, [16,17] have also proposed a relation between diffusion and chaos.
Much of the recent literature focuses on "homogeneous" models of disorder: these can crudely be thought of as models where momentum is not a conserved quantity (hence, by Noether's Theorem, microscopic translation symmetry has been broken), yet the effective equations governing transport remain spatially homogeneous. However, transport coefficients can be sensitive to how translation symmetry has been broken [14], so it is important to test the robustness of any transport bound in inhomogeneous models. Such a test was performed for charge diffusion in a family of holographic models [10], and the inequality in (2) was found to be reversed. In this paper, we test (2) for the energy diffusion constant in an inhomogeneous system: an inhomogeneous analogue of a Sachdev-Ye-Kitaev (SYK) chain of Majorana fermions. We will describe this solvable model of a "strange metal" without quasiparticles in more detail in Section 2. Our main result is that in this model the energy diffusion constant is upper bounded by chaos: We will prove this in the limit where the inhomogeneity is parametrically slowly varying, and provide examples where the ratio D/v 2 b τ l is arbitrarily small, in Section 3. Thus, we do not expect (2) to be generically true in disordered strange metals: holding as a strict inequality up to a finite O(1) prefactor which may be theory-dependent 1 but robust to disorder. Furthermore, in a generic, nearly translation invariant 1 + 1 dimensional field theory without a global U(1) symmetry, the natural diffusion constant is the diffusion of energy. As noted in [6], this diffusion constant will be parametrically large due to the fact that translations are only weakly broken [14]. Hence, we now have examples of strange metals where D is either much larger or much smaller than v 2 b τ l . It is still an interesting open question whether transport properties of strange metals are related to quantum chaos. Any such relation may be restricted to particular diffusion constants and/or models. For example, it may be the case that diffusion constants of translation invariant field theories are related to chaos [13,15].
We hope that variational methods, related to those developed in [18,19,20] for hydrodynamic and holographic models, may be useful in providing rigorous bounds on transport and chaos in disordered strange metals.
For the remainder of the paper, we set = k B = 1.

The Inhomogeneous SYK Chain 2
The SYK model is a strongly interacting large-N model in 0 + 1 spacetime dimensions. It was introduced a long time ago as a model of disordered quantum magnets [21,22]; it was revived more recently [23] due to its possible connection to AdS 2 holography [24,25], which is a toy model for quantum gravity [26,27,28,29,30,31]. Although it has since been shown that the SYK model does not admit a simple holographic dual [32], it does share many fascinating properties of holographic theories, including being "maximally chaotic" [33].
The model that we introduce is a generalization of the SYK chain model developed in [34]. 2 Consider a one-dimensional lattice of L sites, where on each lattice site x there exist N Majorana fermions χ i,x (i = 1, . . . , N ), obeying the standard commutation rules {χ i,x , χ j,y } = δ ij δ xy . The Hamiltonian of these fermions is where the couplings {J ijkl,x } and {J ijkl,x } are all assumed to be independent Gaussian random variables drawn from a distribution with zero mean and following variance: The important (and only) difference comparing to [34] is that we do not assume the variances J 2 1,x and J 2 0,x take the same value for each x. We are interested in the thermodynamic limit L → ∞. One can show that the replica-diagonal 3 partition function, at inverse temperature β = 1/T , can be written as a path integral over two bilocal fields: with the Euclidean time action The "Green's function" {G x (τ 1 , τ 2 )} and "self energy" {Σ x (τ 1 , τ 2 )} are functions of two time variables. The product Σ x G x in above formula is abbreviation for Σ 0,x and J 2 1,x show up to exactly quadratic order in (7) because the random couplings J ijkl,x and J ijkl,x were Gaussian random variables. As the manipulations in this section are essentially identical to [32,34], we only present the few steps where they differ in an important way (due to the absence of translational invariance on average).
It is convenient to rewrite the interaction term (G 4 -term) into the following form: If one chooses J 0,x and J 1,x such that for each x, is a constant independent of x, then the effective on-site coupling is easily seen to be x-independent. The saddle point equations of S eff become and they admit an x-independent approximate solution: which becomes exact at βJ → ∞ (conformal) limit. The system also has a uniform specific heat per site c ≈ 0.396 βJ . Thus, as in [34], this saddle point is identical to the 0 + 1-dimensional SYK model of [23,32] at coupling constant J. If the choice (9) is not made, then the saddle point equations do not admit a homogeneous solution, and it is unclear what the effective theory is.
In the strong coupling limit, N βJ 1, and long wavelength limit, the physics of interest to us is governed by the fluctuations induced by reparametrization modes f x ∈ Diff(S 1 ), which act as ). To quadratic order of the infinitesimal fluctuations, and leading order in 1/βJ expansion, the effective action for the fluctuations has a simple form in Fourier where all x-dependence is contained in the tridiagonal matrix and α = √ 2α k ≈ 12.7 is a constant determined by numerics [32]. A few more steps of this derivation are contained in Appendix A. The long wavelength limit alluded to earlier is the regime when the eigenvalues of C xy are not larger than 1/βJ, which (as we will see) do exist even for the disordered matrix. The derivation of this effective action is identical to [34]: in this previous work, C xy was translation invariant and so (12) was written in momentum space, where the matrix C xy becomes diagonal.
By writing C in the form we immediately recognize that it is positive definite and can be interpreted as the first-order finitedifference discretized version of the differential operator The interpretation of C xy as an approximate differential operator becomes exact when J 2 1,x varies slowly. Letting E denote spatial averages over x, suppose that with M a (large) integer, and f (x) a non-zero function for O(1) argument. To leading order in 1/M , the low-lying spectrum of the discrete operator C xy will be identical to the continuum differential operator (15). In our model, we will take J 1,x to be an arbitrary function of x, simply constrained to 0 J 2 0,x as defined in (9) would be negative). The properties of the matrix C xy will then depend on the inhomogeneity that we encode through x-dependent J 1,x .

Diffusion and Chaos 3
From the effective action (12), we are able to extract the thermal response functions. The procedure is identical to [34] and a diffusion pole is found in the energy density (T tt ) two-point function: T tt x,n T tt y,−n T ∼ |ω n |δ xy + Upon proper analytic continuation to real time, we interpret (17) as having diffusive poles (on the negative imaginary axis) whenever iω is an eigenvalue of C xy . C xy is analogous to a tight-binding-model hopping matrix. If C xy commutes with a discrete translation operator, then we expect plane wave eigenstates, the lowest-lying of which will have an eigenvalue ∼ L −2 in a chain of length L. If C xy is random, then strictly speaking all eigenstates of C xy at fixed ω are localized in the continuum. However, because C xy is analogous to a discretized differential operator (15) which has an exact delocalized zero mode, the low-lying spectrum of C xy will look diffusive on length scales L. In other words, the localization length grows faster than ω −1/2 [37], and the smallest nontrivial eigenvalue scales as L −2 in this case as well.
Hence, the diffusion constant D is finite.
In fact, the lowest-lying non-trivial eigenvectors u x of C xy are well approximated by plane waves: which can be verified numerically. See Appendix B for more comments on this equation. The eigenvalue of such a u x will be Dq 2 , with D an effective diffusion constant. In the large M limit, we may compute D by solving the following differential equation as q → 0: The constant D is computed in Appendix B: This equation has a straightforward physical interpretation. Because the specific heat in our model is x-independent to leading order [34], D is proportional to the thermal conductivity. One can approximate our inhomogeneous chain by joining together homogeneous SYK chains of length L M . Within each of these chains, the thermal conductivity is proportional to the diffusion constant D(x). Joining together these segments of length L , we find a resistor network: hence, the thermal resistivity spatially averages. This leads to (20). In Appendix B, we argue that (20) is applicable even beyond the large M limit, under some assumptions which work relatively well in practice numerically, so long as finite size effects are small. Now we study the butterfly velocity, defined by out-of-time-ordered correlation functions of spatially separated operators. In order to extract v b , we study the (properly regularized) connected out-of-timeordered correlation function. One finds [34], in the region β t β log N that Comparing to (1), we observe that in this model, as in the usual SYK model, This matrix inverse is the discrete analogue of the Green's function In the long range disorder limit, when at each point x the solution of this equation is exponentially decaying [10]: This equation can be derived by noting that, away from the points x = y, we can change coordinates from x to s, defined by D(x)∂ x ≡ ∂ s . One then finds  (20) and (26)) and found numerically, in an inhomogeneous chain with J 2 1 = (1 + J 2 ) −1 , with J = a 0 + a 1 cos(2πx/M ) for a 0 = 0.5 and a 1 = 0.6. The trend of v 2 b τ l to decrease at smaller M towards the "limit" set by diffusion is evident. Right: disordered profiles with J = a 0 + a 1 X, with X = c n cos(φ n + k n x) and c n and φ n random variables chosen so that X ∼ O(1). We take a 0 = 0.5 and vary a 1 , and study the ratio D/v 2 b τ l as a function of M for different realizations of disorder. Discrepancies between the two are enhanced as M becomes large, as expected. We have set L = 6000 and βJ = 25. and when D(s) varies slowly, one can straightforwardly write Hence, we find that D and v b are not equal. Using the Cauchy-Schwarz inequality it is straightforward to conclude The physics at play here is essentially the same as in holography, where charge diffusion was shown to obey a similar inequality, for the same reasons [10]. We have numerically computed D and v b in SYK chains of finite length L with periodic boundary conditions. D is found by averaging the two smallest non-vanishing eigenvales of C xy . v −1 b is found by computing the typical value of − log G(x; y)/|x − y| for |x − y| ∼ L/2. 4 So long as D(x) > D min > 0, we find that D agrees with the "resistor chain" prediction (20) for any M , so long as L/M 20, to within about 0.1% residual error (which is possibly a numerical finite size effect). Indeed, the derivation of (20) in Appendix B does not rely on the assumption that J 1 is slowly varying, so this is not surprising. As shown in Figure 1, we see that while D agrees very well with the "hydrodynamic" theory, v 2 b agrees with the hydrodynamic theory at large M , while partly approaching D/τ l from above as M → 1. This behavior is not surprising: as M becomes shorter, the four-point function (21) begins to "self-average" over the  1,x < X) = X 1+a Θ(X). The circular data points with error bars are numerical data, and dashed lines are the mean of the hydrodynamic predictions for each chain. Qualitative agreement is observed. Right: the numerically computed ratio D/v 2 b τ l as a function of a. We see that this ratio rapidly drops for a 0, as finite size effects become appreciable. This data is taken at βJ = 25.
inhomogeneity in a manner analogous to diffusion. Nevertheless, D < v 2 b τ l holds as a strict inequality in the inhomogeneous systems that we have studied numerically, even when M = 1 (no correlations among D(x). Our violation of the relation D = v 2 b τ l is not limited to the regime of long range disorder. So far, the discrepancies between D and v 2 b τ l are only on the order of a few percent in our numerical data. Yet (29) implies that there is no possible upper bound on diffusion due to quantum chaos. Defining a natural probability measure on D(x) as we estimate that if then v b > 0 but D = 0. 5 We have looked for this parametric breakdown of the relationship between D and v 2 b in smaller chains with M = 1. As shown in Figure 2, we see qualitative agreement (but quantitative disagreement) with our hydrodynamic predictions for D and v 2 b τ l (accounting for finite size effects). As a 0, we observe that the ratio D/v 2 b τ l becomes dependent on the length of the chain, and decreases for the longer chain. This provides evidence that even when M = 1, this inhomogeneous SYK chain may have D = 0 but v b > 0 in the thermodynamic limit.
Finally, let us comment on the low temperature limit β → ∞ (while, of course, taking N → ∞ as well such that βJ N ). At small enough temperature, keeping the inhomogeneity fixed, (24) will break down. Figure 2 suggests that the large M limit is not required to obtain substantial deviations from A more interesting subtlety that arises at very low temperature is the difference between periodic inhomogeneity and random inhomogeneity. For periodic inhomogeneity, one can diagonalize the matrix C xy as a periodic tight-binding hopping matrix in a larger unit cell, and at low enough temperatures, the solution to (23) can be well-approximated by considering only the physics of the lowest band. In this regime, one will recover D = v 2 b τ l . As the period of the periodic inhomogeneity grows longer, the temperature above which D = v 2 b τ l decreases. We expect that for random inhomogeneity (where the eigenstates do not form bands, but are in fact localized) one finds D < v 2 b τ l at all finite temperatures, and provide some numerical evidence for this in Figure 3.

Outlook 4
We have presented a modification of the SYK chain model of [34], in which there is an upper bound on the diffusion constant: D v 2 b τ l . As we pointed out in the introduction, this suggests that there is no (simple) generic bound relating transport and quantum chaos in all strange metals.
One might ask whether our violation of (2) could be found in a "homogeneous" model that does not rely on explicit translation symmetry breaking in the low energy effective description. In the SYK chains that we have studied, this is easy to accomplish at leading order in 1/N , because there is no difference betwen averages over annealed disorder vs. quenched disorder. Hence, consider the partition function with Z SYK the partition function defined in (6), and P[J 2 1,x ] a translation invariant function where J 2 1,x have support on some finite domain. We may think of P[J 2 1,x ] as either a partition function for the "slow" dynamical variables J 2 1,x , or as accounting for certain correlated non-Gaussian fluctuations in the random variables J ijkl,x and J ijkl,x of the microscopic Hamiltonian (4). 6 So long as D v 2 b τ l for each choice of J 2 1,x , by linearity, this inequality will remain true even in the homogeneous model (32) after averaging over the ensemble P[J 2 1,x ] of random couplings.

Derivation of Effective Action A
It is convenient to expand about the saddle using renormalized variables g x , σ x defined by where we have rescaled the fluctuation fields g x , σ x by prefactors |G s | −1 and |G s |. It should be noticed that although the saddle point is uniform in space and translation invariant in time, the fluctuation fields have generic space-time dependence. Now we expand the effective action to second order in the fluctuation fields g, σ, which leads to The spatial kernel S xy is a tight-binding hopping matrix with C xy defined in (13). Next we integrate out σ x and obtain a quadratic action for g x alone. We define K as the (symmetrized) four-point function kernel of the SYK model: We have defined τ ij ≡ τ i − τ j . The effective action of g x is therefore The leading contribution comes from the soft modes which can be identified as induced by reparametrization f x ∈ Diff(S 1 ). They can be interpreted as energy fluctuations. We can linearized these modes and write them in the fourier space: f x (τ ) = τ + x (τ ), n = β 0 dτ e 2πinτ β (τ ). For such modes, we know from [32] that the kernel obeys Following [32], we now find the following effective action for the linearized modes given in (12). The additional numerical prefactors in (12), including |n|(n 2 − 1), come from properly normalizing the eigenvector n,x .

Eigenvalues of Inhomogeneous Diffusion Matrix B
The following argument is reminiscent of [10]. Let us define the matrix and consider the limit q → 0. We postulate that the lowest lying eigenvalues of C xy are proportional to D eff q 2 : Multiplying on both sides by Q we obtain We now look for a series solution to this equation of the form Qu = u 0 + iqu 1 − q 2 u 2 + · · · .
Such a series expansion is reasonable -we expect in any finite chain that the lowest few eigenvectors are delocalized, and have confirmed this numerically. At O(q 0 ), the Q matrix is simply the identity, and so we must take (u 0 ) x = 1.
At O(q 1 ), we find the equation Q is invertible, and hence we conclude that if u 1 is non-trivial, the left-most D must act on a non-trivial vector. Because D has a one-dimensional kernel we conclude We may fix the constant c as follows. First left-multiply by S −1 , and keep only first order terms in q, to obtain cS −1 zw (u 0 ) w = (DQ −1 ) zw (u 0 ) w + iqD zw (u 1 ) w ≈ −iqQ −1 zw (u 0 ) w + iqD zw (u 1 ) w .
In the second equality, we have exploited the fact that Q −1 u 0 is an eigenvector of D of eigenvalue 1 − e iq . Now, we perform a non-rigorous sleight-of-hand: at leading order in q, we may treat qQ −1 zw (u 0 ) w = q(u 0 ) w . We may then left-multiply by u 0 , to remove only the second term of (46), and obtain c = −iq u 0 · u 0 u 0 · S −1 u 0 .
On physical grounds, this is the statement that the eigenvector looks a lot like a plane wave, and this seems to be true numerically. We now have the O(q) corrections to the eigenvector u y , and to leading order q use the variational principle to 'exactly' compute the effective diffusion constant. We find (u 0 + iqu 1 ) · D T SD · (u 0 + iqu 1 ) = q 2 c 2 u 0 · S −1 u 0 = D eff q 2 u 0 · u 0 (48) which gives us, for a periodic chain of L sites: One way to rigorously obtain (46) is to extend a chain of length L into an infinite chain by tiling the same S x repeatedly: S x+L = S x . Using the discrete translation symmetry and choosing the appropriate definition of q, one can demand u x = u x+L , and then sum over only L sites in (46), killing the right-most term. In practice, we have often found this tiling to be unnecessary numerically.