Local criticality, diffusion and chaos in generalized Sachdev-Ye-Kitaev models

The Sachdev-Ye-Kitaev model is a $(0+1)$-dimensional model describing Majorana fermions or complex fermions with random interactions. This model has various interesting properties such as approximate local criticality (power law correlation in time), zero temperature entropy, and quantum chaos. In this article, we propose a higher dimensional generalization of the Sachdev-Ye-Kitaev model, which is a lattice model with $N$ Majorana fermions at each site and random interactions between them. Our model can be defined on arbitrary lattices in arbitrary spatial dimensions. In the large $N$ limit, the higher dimensional model preserves many properties of the Sachdev-Ye-Kitaev model such as local criticality in two-point functions, zero temperature entropy and chaos measured by the out-of-time-ordered correlation functions. In addition, we obtain new properties unique to higher dimensions such as diffusive energy transport and a"butterfly velocity"describing the propagation of chaos in space. We mainly present results for a $(1+1)$-dimensional example, and discuss the general case near the end.


Introduction
Holographic duality, also known as the anti-de-Sitter space/conformal field theory (AdS/CFT) correspondence, refers to the duality between quantum many-body systems in d spatial dimensions and gravitational theories in (d + 1)-dimensional asymptoticaly anti-de Sitter geometries [1][2][3]. Various pieces of evidence suggest that holographic duality is a generic phenomenon that applies beyond the super Yang-Mills theories where the original conjecture was proposed. However, it is difficult to find concrete models for which the duality can be verified directly by comparing bulk and boundary calculations. A very interesting development in this direction is the Sachdev-Ye-Kitaev model [4], which is a (0+1)-dimensional model describing random interactions between N Majorana fermions. When N is large and the temperature is low, the model has an emergent approximate time reparameterization symmetry which is weakly broken. In this limit, two-point functions and four-point functions can be computed. The results suggest that the model has a weakly coupled holographic dual, which includes dilaton gravity in an approximate AdS 2 geometry [5][6][7][8][9][10], weakly coupled to an infinite number of matter fields [11].
The SYK model is a modification of a quantum spin model proposed by Sachdev and Ye more than 20 years ago [12] (which was also related to holographic duality in Ref. [13]). In the large N and low temperature limit N βJ 1 (with β the inverse temperature and J the average coupling strength), the behavior of the model is controlled by a large N saddle point, with the fermion two-point function G(τ, τ ) = 1 N j χ j (τ )χ j (τ ) playing the role of a semi-classical "order parameter" with small fluctuations suppressed by 1 N . At low temperature, G(τ, τ ) has a power law dependence on τ − τ in the infared, suggesting that the low energy dynamics of this model might be conformally invariant [14]. In fact, the (0 + 1)-d conformal symmetry is only approximate, and in fact the low-temperature dynamics are dominated by the specific way in which the symmetry is broken [4,11].
A particularly interesting type of four-point function is the out-of-time-order correlation function (OTOC) F (t) = χ j (t)χ k (0)χ j (t)χ k (0) , which has been proposed as a measure of chaos in quantum systems [15][16][17][18][19][20]. Physically, the decrease of this four-point function measures the increase of the size 1 of the operator anticommutator {χ j (t), χ k (0)}, which indicates how sensitive the system is to an initial perturbation created by acting with the fermion operator χ j (0). The exponential time dependence of the connected part of the OTOC F (t) conn. ∝ e λ L t defines an inverse time scale λ L which can be considered a quantum analog of Lyapunov exponent. Ref. [21] proved a general upper bound λ L 2π β (for a regularized form of OTOC with imaginary time evolution), which is saturated for theories with an Einstein gravity dual. Interestingly, the SYK model also saturates this upper bound [4,22,11]. Many other aspects of the SYK model (and a similar model for complex fermions) have been investigated recently [23, 22, 24-27, 8, 10, 28].
Given the interesting properties of the (0 + 1)-dimensional SYK model, it is natural to look for higher dimensional cousins. In this paper, we propose a family of higher dimensional variants of the SYK model, which remain solvable in the large N limit. Our model can be defined on an arbitrary discrete lattice in arbitrary spatial dimensions. There are N fermions on each site with an SYK Hamiltonian. Different sites have independent SYK couplings, and neighboring sites are coupled by random four-fermion terms. Using the same techniques as those applied to the SYK model, we can study two-point functions and four-point functions in space-time. The spatial locality of the model allows us to study transport properties and propagation of quantum chaos in space-time. We find that the disorder-averaged two-point functions vanish between different lattice sites, and have the same local critical behavior as in the SYK model within each site. Our model also has the same zero temperature entropy at each site as the SYK model. Correlation between different sites begins at the level of four-point functions. Similar to the SYK model case, one can consider the fermion four-point function as being mediated by a series of collective fields, with the leading contribution coming from energy fluctuations. In our generalized models, the four-point function allows us to study the dynamics of collective fields in spacetime. In particular, we find a diffusive dynamics of energy density, which means this model describes a diffusive strongly correlated metal phase. The OTOC can also be studied, which now has both spatial and temporal dependence. We show that at low temperature the Lyapunov exponent still saturates the chaos bound 2π β . The propagation of quantum chaos in space-time can be characterized by a "butterfly effect velocity" x = v B t, which means that the anticommutator {χ jx (t), χ jy (0)} becomes significant at t = |x − y|/v B + (const). Interestingly, the diffusion constant D and butterfly velocity v B in our generalized SYK model satisfy a simple relation D = v 2 B 2πT , which realises a bound conjectured on diffusion in incoherent metal and agrees with the holographic calculation on incoherent black hole [29][30][31].
The remainder of the paper is organized as follows. In section 2 we will briefly review the properties of the original SYK model. In section 3, we define the generalized SYK model in higher dimensions. For concreteness, we work on a (1 + 1)-dimensional example and study its correlation functions and thermodynamic properties in detail. In section 4 we study the dynamics of collective fields in this system based on an operator-product expansion of fermion four-point functions. In particular, we show that one of the collective fields describes the diffusion of energy in this system. In section 5, we study the OTOC of the (1 + 1)-dimensional model and obtain the Lyapunov exponent and the butterfly velocity. In section 6 we discuss the general form of the model in generic dimensions and graphs. In section 7 we end the paper with a summary and discussion of further topics.

Review of the SYK model
In this section, we briefly review some basic facts about the SYK model. The SYK model consists of N Majorana fermions χ j , j = 1, 2, . . . , N with a random four-fermion interaction [4] H = 1 j<k<l<m N J jklm χ j χ k χ l χ m , {χ j , χ k } = δ jk (1) where {J jklm } are independent random couplings with zero mean J jklm = 0, and variance of individual elements is given by Here, the symbol J without indices is a dimensionful constant that sets the scale of the Hamiltonian, and the factors 1 3! N 3 are for later convenience. The interaction is all-to-all, so that there is no spatial locality, and this model should be considered as a (0 + 1)-d quantum mechanical system. The model is solvable at large N , and exhibits holographic behavior at strong coupling N βJ 1. In this limit, the (finite temperature) two-point function [12,14] where we denote τ 1 − τ 2 by τ 12 for simplicity here and below. One can also compute the four-point function in the holographic limit [4,11]. The interesting piece of the four-point function is the connected part, which begins at order 1 At leading order of 1/βJ and 1/N , when we consider a configuration of imaginary (Euclidean) times with the ordering τ 1 > τ 2 > τ 3 > τ 4 , then the correlator F factorizes: where α K ≈ 2.852. The factorization indicates that the OPE of two fermion operators is dominated by the conserved quantity in this model -the Hamiltonian itself, and the fluctuation of energy determine the four point function in this configuration.
Quantum chaos can be diagnosed by the out-of-time-ordered correlator (OTOC). The OTOC is calculated by starting with the the Euclidean correlator with a "j-k-j-k" configuration of times, e.g. τ 1 > τ 3 > τ 2 > τ 4 , and then giving the times τ 1 and τ 2 large real-time parts. Making the particular choice that all four points are evenly spaced on the imaginary time circle, we have [11] which exhibits an exponential growth for real time t greater than the dissipation time t d ∼ β, i.e., F ∼ −βJ exp 2π β t for t β. This growth was recognized as a signature of chaos [15,16,18], moreover, the corresponding growth exponent λ L = 2π β saturates the chaos bound proposed in Ref. [21].
The two-point function and four-point function shown above can be calculated via standard Feynman diagrams with large N simplification. Another approach is to use the disorder-averaged effective action. After introducing a pair of auxiliary bi-local fields "Green's function" G(τ 1 , τ 2 ) and "self-energy" Σ(τ 1 , τ 2 ) and carrying out the disorder average of the partition function, the fermions can be integrated out, leaving [32,4]: The large N prefactor implies that the problem is essentially classical. The saddle point reproduces the Schwinger-Dyson equations that can also be derived via Feynman diagrams: where we assumed time translation symmetry to simplify the equations. One can also obtain the connected part of the four point function by considering quantum fluctuations around the saddle point. Among all the quantum fluctuations, there is a special class induced by the reparametrization of the time circle f ∈ Diff(S 1 ), which contributes the leading piece at strong coupling βJ 1. The effective description of this part turns out to be well-captured by a local action proportional to the Schwarzian derivative [4,11]. Remarkably, the same form of effective action also appears in the AdS 2 Einstein-dilaton theory [5][6][7][8][9][10].
One can also derive the thermodynamic properties from the large-N saddle point free energy: In the second line we write the free energy in a low temperature expansion, 3 where U ≈ −0.0406J is the ground state energy, S 0 ≈ 0.232 is the zero temperature entropy [32,4], βJ is the specific heat [11]. The entropy term can be derived by inserting the conformal saddle point solution (2) in the effective action. The specific heat can be derived from knowledge of the leading (in 1/βJ) correction to the conformal saddle, but the energy requires the exact (numerical) finite βJ solution.

The generalized SYK model
In this section, we will present a simple way to generalize the SYK model to higher dimensions while keeping the solvable properties of the model in the large-N limit. For concreteness of the presentation, in this section we focus on a (1 + 1)-dimensional example, which describes a one-dimensional array of SYK models with coupling between neighboring sites. It should be clear how to generalize, and we will discuss more details of the generalization to arbitrary dimensions and generic graphs in section 6.  In (1 + 1)-d, our model describes a coupled array of SYK model sites, as shown in figure 1. Each site containins N 1 Majorana fermions with SYK interactions drawn independently for each site. Each pair of neighboring sites are then further coupled via random four fermion interactions with two of the fermions from each site. The Hamiltonian has the following form:

Definition of the chain model
where x labels the lattice sites, and {χ j,x } j=1,2,...,N ; x=1,2,...,M are the Majorana fermion operators satisfying anti-commutation relations and periodic boundary condition: {χ j,x , χ k,y } = δ xy δ jk , χ j,0 ≡ χ j,M . N is the number of Majorana fermions on each site and M is the number of the sites, or equivalently, the length of the chain. In this expression, we restrict the range of indices in the sum such that each term only appears once. The first term describes the on-site SYK interaction, while the second term is the nearest neighbor random four fermion coupling. The random couplings {J jklm,x } and {J jklm,x } are drawn independently for each value of x, from a distribution with zero mean and variance defined in the following way: The normalization, especially the factors of N , are chosen to make the large N limit uniform and ensures the dimension 1 coupling constants J 0 and J 1 represent the average strength of the thermal bath seen by each fermion field. Comparing to the original SYK model, our model clearly has spatial locality. One can view our model as either M coupled SYK sites, or equivalently as a big SYK site with N M Majorana fermions but with inhomogeneous coupling strength-the non-local couplings (between sites |x − y| > 1) are suppressed. As will be discussed in section 6, it is not essential to have the 2-2 coupling between neighboring sites. Introducing more generic couplings such as χ j,x χ k,x+1 χ l,x+2 χ m,x+3 is straightforward as long as the coefficients of these terms are all independent variables. However, in this section we will focus on the 2-2 coupling case for simplicity.
A first observation of the model is that the Hamiltonian doesn't contain any quadratic term, so that the free propagator is diagonal not only in flavor indices but also in spatial coordinate T τ χ j,x (τ 1 )χ k,y (τ 2 ) free = 1 2 δ jk δ xy sgn(τ 12 ). Furthermore, at leading order the j, x j, x interaction vertices under random average only contribute to the diagonal part, see figure 2. Therefore, the dressed Green's function T τ χ j,x (τ 1 )χ k,y (τ 2 ) is also diagonal in both flavor and spatial coordinates 4 . In addition, the diagonal Green's functions T τ χ j,x (τ 1 )χ j,x (τ 2 ) are independent of j due to an SO(N ) symmetry of the model after averaging over disorder. Therefore we have T τ χ j,x (τ 1 )χ k,y (τ 2 ) = δ jk δ xy In this section, we employ the large N effective action approach to analyze the model. In principle, to study the quenched problem, one should introduce n replicas and analyze the disorder averaged partition function Z n . However, as in the original SYK model, our chain model self-averages at large N . One can verify this by checking the replicon off-diagonal contribution to the averaged partition function Z n . Figure (3) shows the leading replicon off-diagonal diagram, which is suppressed by 1/N 3 relative to the connected diagonal diagrams. We will be interested in correlators at order 1/N , so we can assume Z n = Z n and directly work with the disorder averaged partition function and correlators. In other words, quenched equals annealed at the order we work, so we can study annealed quantities.

The effective action and the saddle point
To derive the effective action, we begin by integrating the partition function Z[{J jklm , J jklm }] over disorder {J jklm } and {J jklm } with Gaussian distributions. This will produce eightfermion interaction terms that are non-local in time: The expression of the averaged partition function can be further simplified by introducing bilocal auxiliary fields G x (τ 1 , τ 2 ) and Σ x (τ 1 , τ 2 ). Σ x (τ 1 , τ 2 ) are Lagrange multipliers which impose the constraints Integrating out the fermion fields χ jx (τ ) one obtains the effective action of the bilocal fields: All terms except the last term ∝ J 2 1 describe decoupled SYK models at each site, and the J 2 1 term couples the bilocal fields G x (τ 1 , τ 2 ) on neighboring sites. In the large N limit, the saddle point of this action determines the leading order behavior of the two-point function G x (τ 1 , τ 2 ). The saddle point equation δS eff δG = 0 and δS eff δΣ = 0 produce the Schwinger-Dyson equations (assuming time translation symmetry): with G x (τ ) = G x (τ + τ 1 , τ 1 ) and G(iω) its fourier transformation, and similarly for Σ x (τ ) and Σ x (iω). This set of equations can be equivalently derived via Feynman diagrams: where the thick lines represent dressed Green's functions and the gray disk represents the self-energy.
Notice that the chain model under disorder average is translation invariant and also has a translation invariant free propagator G free x (τ ) = 1 2 sgn(τ ). Therefore we consider translation invariant solutions of the Schwinger-Dyson equations: Comparing to equation (7), one sees the Schwinger-Dyson equations reduce to exactly the same form as those of a (0 + 1)-d SYK model with the coupling constant J = J 2 0 + J 2 1 . Therefore, we can directly apply the SYK model results in Ref. [4,11]. In particular, we immediately know that the solution in the conformal limit N βJ 1 has the following familiar form: Here and below, we will denote the effective coupling J 2 0 + J 2 1 in our model as J. In summary, the two-point function in our chain model is local in space and powerlaw decaying (at low temperature) in time, a behavior known as local quantum criticality [33,34]. This saddle point solution (and the finite βJ corrections to it) is the starting point for studying other properties of the model.
For example, inserting the solution into the action, we derive the saddle point approximation to the partition function. This gives the order N term in the free energy: Here we have used space translation symmetry to simplify the notation. (M is the number of lattice sites, i.e. the spatial volume of the system.) Using the saddle point solution one can see that the free energy density F N M agrees exactly with that in the SYK model Eq. (8) with J = J 2 0 + J 2 1 . Therefore, in the large N limit our model has exactly the same zero temperature entropy per fermion S 0 ≈ 0.232, and specific heat c v ≈ 0.396 βJ per fermion. It should be noted that the zero temperature entropy is extensive.
We should remark here that the exact agreement on thermodynamic properties with the SYK model only holds at leading order in the 1/N expansion. This is similar and related to the discussion of the two-point function, where the agreement with the SYK model only holds at leading order. We will see below that 1/N effects, which are the quantum fluctuations around the large N saddle point, are different in our chain model than in the SYK model.

Fluctuations of the collective fields and four-point functions
The fermion four-point function can be determined from the fluctuations about the saddle point just discussed. At large N the saddle is sharp, and the connected four point function is small, of order 1 N . This is determined by the Gaussian fluctuations of bilocal fields G x (τ 1 , τ 2 ) and Σ x (τ 1 , τ 2 ) about the saddle.
It is convenient to expand about the saddle using variables g x , σ x defined by where we have rescaled the fluctuation fields g x , σ x by prefactors |G s | −1 and |G s | for convenience. 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. Using Eq. (13), the connected, averaged four-point function of the fermions (defined anal-ogously to Eq. (3)) can be written as the connected two-point function of g x (τ 1 , τ 2 ): More precisely, F N is the connected part of the fermion four point function. We will see that F is of order one at large N , so the connected correlator is of order 1 N . To compute the g x g y correlator, 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 where we continue to use the notation that J = J 2 0 + J 2 1 . It is straightforward to 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 [4,11]: The effective action of g x is which determines the fermion four-point function: Here we have written the correlator in a compact matrix form.
Comparing Eq. (29) with the four-point function of the SYK model, the only difference is with the spatial kernel S. Replacing S xy by δ xy (i.e. taking J 1 = 0) reproduces the SYK result, as expected.
The behavior of the four-point function (29) can be analyzed by diagonalizing the kernel ( K −1 − S). First of all, due to translation symmetry it is straightforward to do a spatial Fourier transformation to (x − y) and define the Fourier component Then one can diagonalize the temperal kernel K in the same way as in the SYK model. We write the antisymmetric eigenfunctions Ψ h,n (τ 1 , τ 2 ) where n labels the fourier mode for the sum of the two times, and h specifies the dependence on the difference of the times. Writing the corresponding eigenvalues of K as k(h, n), the four-point function can be expressed as (where we have used the fact that the symmetrized kernel is Hermitian). The only difference from the original SYK model is the factor of s(p) 1 in the denominator. 5 The details of the eigenvectors and eigenvalues of the temporal kernel, k(h, n) and Ψ h,n have been worked out in Ref. [11], and we will use them in the following sections. It should be noted that Eq. (31) returns to the SYK result at zero momentum p = 0, as one expects from the form of the effective action.

Symmetry breaking and pseudo-Goldstone mode
Although the general discussion above is sufficient for calculating four-point functions, it is helpful to gain more physical understanding by analyzing symmetry properties of the saddle point solution. As in the SYK model, our effective action (14) admits an approximate reparametrization symmetry of time in the IR limit ω → 0, where one can ignore the −iω in the first term. The reparameterization transformation is defined by  (20) breaks the emergent reparametrization and lifts the Goldstone modes to pseudo-Goldstone modes. (c) The situation in the SYK chain model is similar to a ferromagnetic spin chain with a small pinning field, where the SU(2) symmetry is "almost spontaneously" broken to U(1), leading to a pseudo-Goldstone mode.
The symmetry Diff(S 1 ) is broken spontaneously to PSL 2 (R) by the the solution in equation (21). If this symmetry were exact, the spontaneous symmetry breaking would have produced infinite number of Goldstone modes corresponding to the spatially dependent reparametrizations f x ∈ Diff(S 1 )/ PSL 2 (R). The −iω term in the effective action (20) plays the role of a small explicit symmetry breaking field, which selects the solution in Eq. 21 and turns the Goldstone bosons to pseudo-Goldstone bosons. This is similar to the situation in a ferromagnetic spin chain with small external magnetic field (see figure 4). The effect of this symmetry breaking term is small at large βJ, which means that the pseudo-Goldstone bosons have large fluctuations and make the most important contributions to the long-wavelength dynamics. In particular, they are responsible for the diffusion of energy that we will analyze in section 4 and the chaos characteristics we will study in section 5.
From the view point of the effective action, the pseudo-Goldstone modes are those fluctuations along the "nearly-flat" direction around the saddle point in the potential of S eff [G]. As we have discussed, these fluctuations correspond to the spontaneously and ex-plicitly broken reparametrization symmetry, and therefore can be parametrized by residue target space Diff(S 1 )/ PSL 2 (R) at each point in space. More explicitly, f x acts on the saddle point in the following way For small deformations of time f x (τ ) = τ + x (τ ), the quadratic action for x (τ ) can be determined by diagonalizing the kernel K and using (28). Building on the results of [11], we will find below that to quadratic order in the spatial momentum p, this leads to the action The first term is familiar from the SYK model. It is local in time, and can be interpreted as a quadratic approximation to the Schwarzian derivative action at each site. The coefficient 1 βJ tells us that at large βJ limit the reparameterization fields are soft modes, due to the approximate reparameterization symmetry. The second term describes a simple coupling of the reparameterization modes at different sites, but with a nonlocal form as a function of time. 6 We will see that together, these two terms determine both the energy diffusion dynamics and chaos behavior. 7 In the following two sections, we will analyze two different regions of the four-point function with different ordering of the time variables τ 1 , τ 2 , τ 3 , τ 4 as shown in figure (5). The four-point functions with ordering τ 1 > τ 2 > τ 3 > τ 4 , which we will sometimes refer to as j-j-k-k order, determines the dynamics of collective fields in our model, while the order τ 1 > τ 3 > τ 2 > τ 4 which we will refer to as j-k-j-k determines OTOC after analytic continuation, and characterizes chaos in our system.

The OPE region and the Energy Transport
As we have seen from the two-point functions and four-point functions, single fermion fields χ jx (τ ) do not propagate spatially in our model. The only fields that have nontrivial dy-space imaginary time namics are the collective bilocal fields g x (τ 1 , τ 2 ), which are singlet in the flavor SO(N ) symmetry. In the four-point function , and take the limit |τ 2 −τ 1 |, |τ 4 −τ 3 | |τ |, we are effectively taking an operator product expansion (OPE) of the fermion fields χ j,x (τ 1 )χ j,x (τ 2 ) and similarly for χ k,y (τ 3 )χ k,y (τ 4 ). The four-point function becomes a sum over the propagators of different collective fields that appear in the OPE. Roughly speaking, different collective fields correspond to different families of eigenvalues in Eq. (31).
The eigenvalues and eigenfunctions of the temporal kernel have been studied in Ref. [4,11]. In the strong coupling limit βJ 1, the two-point functions take the conformal form (21). This correlation function is covariant under Möbius transformations (PSL 2 (R)) of the time coordinate, which is used in diagonalizing the temporal kernel. It turns out [4,22,11] that the eigenvalues of the kernel are given by the simple formula The label set I consists of a discrete series and a continuum [22,11]. For small momentum , we will see that 1 βJ corrections to this formula are important. First, however, we study the case of finite p where the conformal large βJ limit is straightforward.

The OPE at finite p
The fact that the eigenvalues (36) don't depend on n is a consequence of the conformal symmetry of the eigenvalue problem at large βJ. It makes it possible to do the sum over The sum and integral over the remaining eigenvector index h ⊂ I can then be done following Ref. [11] by deforming the contour of integration over the continuum portion of I. As the contour is pushed to infinity, one encounters two sets of poles: one set conveniently cancels the contribution from the discrete h ⊂ I and the other gives the answer (in the region τ 1 >τ 2 >τ 3 >τ 4 where η < 1) The only difference in (38) from the regular SYK case [11] is that we have the factor of s(p).
It is also useful to consider the fourier transform of F p back to position space. 8 One can show that c 2 m (p) is analytic in a strip surrounding the real p axis, so the fourier transform will decay exponentially in x. More precisely, c 2 m (p) has both a pole and a branch cut for complex p. The pole is at the location where h m (p) becomes an even integer, and the tan(πh m /2) in the denominator diverges. The branch cut is due to the multivaluedness of the solution to k(h m ) = s(p) −1 and occurs at the location where k (h m (p)) = 0. For large m, one can show that both singularities are at a distance of order log m from the real p axis, so the contribution of the m mode will decay in space as e − log(m)|x| for large x and large m. Here x is measured in lattice units, so the modes only propagate a few lattice sites. However, note that the decay factor in the exponential, log(m)|x| ∼ log(h)|x|, grows less rapidly with h than in a conformally invariant (1 + 1) dimensional theory, where we would have h|x|.

The h = 2 contribution for small p
For small momentum p, the first solution to k(h m )s(p) = 1 approaches h 0 (0) = 2. Because of the factor of tan πhm 2 in the denominator of Eq. (40), this leads to a large OPE coefficient c 2 0 (p) that diverges like p −2 as p → 0. It will be important to understand how this is cut off at finite βJ.
The would-be divergence can be traced to the h = 2 ⊂ I eigenfunctions of K. In this section we study the h = 2 piece in more detail, which requires taking into account the non-conformal corrections to k(h, n). The correction can be derived by including the 1/βJ correction to the conformal form of the two-point function, as has been discussed in Ref. [11]. In the following we will apply the SYK results to our generalized model.
To start, we note that due to time translation symmetry, the eigenfunctions Ψ 2 (τ 1 , τ 2 ) for h = 2 have the form Ψ 2 (τ 1 , τ 2 ) = Ψ 2,n (τ 12 )e −i 2π β n τ 1 +τ 2 2 . In the conformal limit βJ 1, Ψ 2,n (τ 12 ) has the following explicit form: In the above expression, γ 2 n = 3 π 2 |n|(n 2 −1) is the prefactor to normalize Ψ 2,n . In the first order perturbation theory, one can use these zeroth order (conformally covariant) eigenfunctions to obtain the first order corrections to the eigenvalues, which are given in Ref. [11] by where α K ≈ 2.852. Summing over the eigenfunctions with this improved formula for the eigenvalue, we find the h = 2 contribution In the second step we have expanded s(p) in the long wavelength limit s(p) 1 − J 2 1 3J 2 p 2 , and defined a constant D which, as will be shown below, is the energy diffusion constant of the system: Notice that Eq. (44) has a smooth p → 0 limit, which reduces to the corresponding formula of the SYK model in Ref. [11].
OPE from the correction of h = 2 parts Figure 6: (a) The conformal limit of the four-point function corresponds to collective fields that are locally critical and short-range correlated in space (see equation (41)). (b) The non-conformal corrections of the four-point function corresponds to the reparameterization field, which has a diffusive dynamics in space-time.

Energy transport
In subsections 4.1 and 4.2 we discussed contributions to the fermion four-point functions that are naturally organized in terms of the OPE. The operator product of two fermion fields χ j,x (τ + τ 12 2 ), χ j,x (τ − τ 12 2 ) generates an infinite family of collective fields, including the reparameterization field (h = 2 contribution) and the local conformal fields φ m . In the strong coupling limit βJ 1, the dominant contribution for small p is the h = 2 piece, which has diffusive behavior as can be seen from Eq. (44) and (34). It is natural to relate the h = 2 contribution to energy transport properties, since the time reparameterization field is related to diffeomorphisms.
To understand the relation to energy transport more explicitly, we start from p = 0, where the h = 2 four-point function (44) reduces to the SYK result. As shown in Ref. [11], F p=0,h=2 (τ, τ 12 , τ 34 ) has a factorized form in the j-j-k-k order region (τ 1 >τ 2 >τ 3 >τ 4 ), which satisfies the following equality [9]: Here c v is the specific heat per fermion, such that is the thermal fluctuation of the total energy. Here H is the Hamiltonian of the chain model. Thus Eq. (46) can be rewritten as Eq. (48) implies that the h = 2 piece of the fermion four-point function at p = 0 describes the thermal fluctuation of the two-point function G s (τ 12 ) induced by the fluctuation of total energy. If we define the energy density of the chain model as T 00 (x), and define its Fourier component as T 00 (p) = M −1/2 x T 00 (x)e −ipx , we have T 00 (p = 0) = H/ √ M . It is natural to generalize Eq. (48) to finite (small) momentum and express the energy density correlation function as T τ T 00 (−p, τ )T 00 (p, 0) conn. F p,h=2 (τ, τ 12 , τ 34 ) Eq. (49) is expected to hold in the limit τ 12 , τ 34 → 0 (which guarantees the j-j-k-k order for any finite τ ). In this limit, by a Taylor expansion of F p,h=2 (τ, τ 12 , τ 34 ) in τ 12 and τ 34 one obtains Evaluating ∂ β G(τ ) for small τ and using c v = πα K 16 √ 2βJ and Eq. (49), we find with ω n = 2πn β . This equation directly gives the Matsubara correlator C τ 00 (p, iω n ). By analytically continuing iω n → ω + iδ in C τ 00 (p, iω n ) and subtracting a contact term, one obtains the retarded energy-density correlation function in the (p, ω) space 9 : Here we have taken the small frequency limit and omitted the higher order terms in ω.
Eq. (51) tells us that energy density perturbations in our system satisfy a diffusion equation with diffusion constant D given by Eq. (45). From energy density correlations one can also derive the thermal conductance: For p = 0 this formula reduces to κ = N c v D = C v D, with C v = N c v the specific heat per site.
In summary, in this section we analyzed the fermion four-point function in the limit τ 12 , τ 34 → 0, which determines the behavior of SO(N ) singlet collective fields in the chain model. The only collective field with nontrivial spatial dynamics is the time reparameterization field, the dynamics of which describes diffusion of energy in this disordered system, with a temperature independent diffusion constant. Although there are infinite number of other collective fields in this system, all of them are locally critical and decay exponentially in space with order 1 correlation length.

Chaos and the butterfly velocity
Another interesting aspect of the model is the OTOC, which have been studied as a diagnostic of chaos. In addition to the exponent λ L that determines the exponential growth of (anti)commutators in the large N dynamics within a single site, the spatial locality in our generalized model allows us to study the dynamics of chaos in space. An important parameter is the speed v B [16,35] defined as the speed of growth of the "filled light cone" that marks the region where operators have large (anti)commutators with an initial operator. It has the interpretation of the speed at which the butterfly effect spreads in space, and has also been related to the Lieb-Robinson velocity [36] recently [37,38]. In this section we will evaluate v B . Figure 7: Double Keldysh-Schwinger contour with operators equally spaced in imaginary time. (Note that here we use the convention that the real/imaginary part of t is the real/imaginary time, which is different from our convention in previous section.) The OTOC can be computed by analytic continuation of the imaginary time ordered correlation function. A convenient special case to study is the "regularized OTOC," where the operators are equally spaced in imaginary time as shown in figure 7: Tr (rχ j,x (t)rχ k,0 (0)rχ j,x (t )rχ k,0 (0)) , with r = ρ(β) 1/4 = e − β 4 H /Z 1/4 . The regularization does not change the qualitative behavior of the OTOC. For fixed x we expect the following behavior for F (x, t): at early times it should be approximately equal to the disconnected product −G( β 2 ) 2 (the minus sign is because we have fermions), and at late times it should be small, indicating a large anticommutator between χ j,x (t) and χ k,0 (0). The butterfly velocity v B is defined as the rate at which the region where F is small expands outwards as we increase t. Of course, to actually see that F becomes small, we would have to sum all orders in the 1 N expansion. We can only compute the first 1 N term, but we assume that the exact F becomes small at around the time where this correction becomes comparable to the order one disconnected contribution −G( β 2 ) 2 . In the imaginary time configuration, the 1 N part of the four point function is given by − F N , where F is the function studied in the previous section. More precisely, we studied the spatial fourier transform F p . To compute the 1 N term in F (x, t), we will continue F p to the configuration in Eq. (53), and then finally fourier transform back to position space. We will see that there are some subtleties involved in getting an expression for F p that is accurate for the small momenta that dominate this fourier transform.
To begin, we will warm up by studying the case where p 2 J βJ 2 1 so that we can use the conformal limit of the kernel. From Eq. (37), we have that the cross ratio of the times in the configuration (53) is For large t, this is small and purely imaginary, but we begin the continuation from t = 0, where η is greater than one. The four point function and this continuation are discussed in detail for the SYK model in Ref. [11]. The only difference in our case is that we have to insert a factor of s(p). After the contour manipulation and expansion of the relevant hypergeometric functions [11], one finds that the growing term at small η is Using Eq. (54), we see that this implies an exponential growth F p ∼ e 2π β (h * (p)−1)t , so we have a momentum-dependent chaos exponent λ L (p) = 2π This leads to a p −2 divergence in the prefactor coming from the tan(πh/2) in the denominator of (55). As before, this divergence can be traced to the h = 2 eigenfunctions, which we can treat more accurately by directly continuing Eq. (44) to the configuration (53). In appendix B, we show that the growing term after this continuation is This expression has a smooth p → 0 limit. However, notice that the exponential growth is e 2π β t independent of p, whereas we have argued that at finite p the exact answer should be modified to e 2π β (h * (p)−1)t . Also, at finite βJ, we expect a modification of the growth exponent proportional to 1 βJ . Both of these modifications must come from the sum over h = 2 eigenfunctions. We can perturbatively compute both of them as follows. We have two small parameters, p 2 and 1 βJ , and we consider both to be first order quantities, of order . If we expand in , we expect to find At order −1 , we have 1 b 1 e 2π β t . Comparison with (57) then determines b 1 . At order 0 we expect a term 2πλ 1 βb 1 te 2π β t coming from the small shift in λ L . Since this is at order 0 , we can compute it in the theory with p = 0 and βJ = ∞. Indeed, Ref. [11] did find such a term in this limit, which for our case evaluates to 6π β te 2π β t . Matching to the expected term, we find λ 1 = 3b 1 , giving the formula The exponent in Eq. (59) is correct to order p 2 and order 1 βJ and the inverse of the prefactor is correct to the same accuracy. One can check that when βJ = ∞ it agrees with the exact result Eq. (55) to the claimed order in p 2 .
We now have an expression that is sufficiently accurate at small p, so we can do the final step and transform Eq. (59) to position space to compute F (x, t). Approximating the discrete fourier transform as an integral, we have The integrand has a pole at momentum p = i 2π βD , which dominates the behavior for large x, leading to 10 This is the main result of this section, giving the butterfly velocity v B . The order 1 N term competes with the order one term, indicating that the anticommutator has become large, when t = β 2π log N a 1 + x v B . The formula for v B given in Eq. (61) agrees with the relation identified in holographic theories, see Ref. [31].
It is remarkable to find a simple relation between the diffusion constant and butterfly velocity of this type. In this model, it is a consequence of the fact that the same reparameterization degrees of freedom are responsible both for energy diffusion dynamics and the OTOC chaos behavior. This is a property that the model shares with conventional 10 Some constants that appear in the following equations are holographic theories, where the the gravitational field in the bulk determines both of these phenomena. It would be interesting to work out whether v 2 B = 2πD β continues to hold at higher orders in 1 βJ . At least naively, because other modes besides reparameterizations will become important, one would not expect this equality to persist beyond infinite βJ.
There is an interesting subtlety in the fourier transform that we glossed over above. In fact, the pole dominates only for x v B t βJ . This means that for x v B J log N , the pole analysis will not be correct at the time when the anticommutator becomes large. For such x, we can approximate (60) another way by replacing b(p) in the denominator by b(0) and doing the Gaussian integral. This leads to which is accurate for x v B t βJ . In this region, we find that the "butterfly cone" is rounded out, see figure 8. It is rather striking that these two different regions of behavior characterized by (62) and (61) were also identified in the analysis of stringy corrections to the holographic computation of F (x, t) in [39].
We will close this section with one further comment. A surprising feature of the large x behavior (61) is that the growth as a function of time is given by e 2π β t , despite the finite p and 1 βJ corrections to λ L in Eq. (59). These corrections, which are both negative for real momenta, cancel against each other when we evaluate at the pole at imaginary p. It would be interesting to know if this persists at higher orders in 1 βJ . Note that in the small x region (62), the growth as a function of time is decreased by a 1 βJ correction, with the same coefficient as in the original SYK model.

General construction and higher dimensions
In previous sections, we have focused on a (1 + 1)-dimensional chain example of the generalized SYK model, but our construction actually applies to generic dimensions and more general forms of interactions. In this section, we will discuss the general form of our model. We will start from the simple case of higher dimensional regular lattices and then discuss the even more general cases beyond that.

Generalized SYK model on higher dimensional lattices
In the chain model example, the different sites are only coupled by a 2-2 coupling with two fermions from each site. This restriction is chosen for simplicity, which is not necessary. In general, our model can be defined on arbitrary graphs, including higher dimensional translation invariant lattices and non-translation invariant graphs (see figure (9) for an illustration on square lattice). We denote the set of sites in the graph as Γ, and label the sites by x, y, z....
here in the sum over x, y, z, w, one can always define an order of the sites, and avoid duplication due to different order of x, y, z, w. For fixed sites x, y, z, w, we can further J uuwz w z u x y J xxyy Figure 9: An example of two-dimensional square lattice model. The Hamiltonian could contain all kinds of random four-fermion terms, for example, J jklm,uuwz χ j,u χ k,u χ l,w χ m,z and J jklm,xxyy χ j,x χ k,x χ l,y χ m,y as shown in the figure. where J xyzw is fixed in the large N limit. This model includes SYK model and the chain model discussed above as two special cases. The SYK model corresponds to taking J xyzw = J uniform for all x, y, z, w, which gives a completely non-local Hamiltonian (or equivalently it can be treated as the case Γ only contains a single site). The chain model is obtained by setting J xyzw to zero except {J xxxx } x∈Γ and {J xxyy } x∈Γ for y = x+1.
Remarkably, this general model can still be solved in the same way as the chain model. From the Feynman diagrams shown in figure 10 one can directly see that the disorder averaged two-point function is still diagonal between different sites, even if there is no symmetry reason for it to vanish. In path integral approach, we can write down the 11 For example: if x = y = z = w, we restrict 1 j < k < l < m N ; and if x = y = z = w, we restrict 1 j < k N and 1 l < m N , as we did in the chain model. effective action of the same collective fields G x (τ 1 , τ 2 ) and Σ x (τ 1 , τ 2 ): where C xyzw are combinatorial factors which depend on how many sites in xyzw coincide. For example C xxxx = 1 2·4! , C xxyy = 1 2·2!·2! for x = y. 12 In the large N limit, the corresponding saddle point equation always admits a translation invariant solution G x (τ 1 , τ 2 ) = G s (τ 1 , τ 2 ), which is identical to the two-point function of SYK model with effective coupling Here we denote |Γ| as the total number of sites. Similar to the chain model, we can expand G x (τ 1 , τ 2 ) around the saddle point by defining G x (τ 1 , τ 2 ) = G s (τ 12 ) + |G s (τ 12 )| −1 g x (τ 1 , τ 2 ). Expanding the effective action to quadratic order of g x (τ 1 , τ 2 ) leads to δS eff [g] = 3J 2 4 x,y d 4 τ g x (τ 1 , τ 2 ) K −1 (τ 1 , τ 2 ; τ 3 , τ 4 )δ xy − S xy δ(τ 13 )δ(τ 24 ) g y (τ 3 , τ 4 ) (67) The effective action has the same form as equation 28 except for a different spatial kernel S xy : The discussion so far does not rely on translation symmetry, and applies to general graphs. If the graph is a d-dimensional lattice and the coupling has translation symmetry, the spatial kernel will also have the same symmetry, so that S xy = S(x − y) where the labels x, y should be considered as d-dimensional vectors. In this case the spatial kernel can be diagonalized by Fourier transformation. If S xy is short-ranged, the Fourier 12 In general, a group of sites xyzw defines a partition (n 1 , n 2 , ..., n k ) of 4: k i=1 n i = 4 with n i the multiplicity of site i. For example, xxxx corresponds to partition (4) and xxyy corresponds to the partition (2,2). The combinatorial factor is given by C xyzw = 1 2 1 i ni! transformation s( p) is a smooth function of p, which can be expanded at small p as s( p) 1 − i=1,2,...,d a i p 2 i . In the same way as in the chain model, we obtain an energy diffusion, and a i determines the diffusion constant of energy along the i-th direction.
Following the same approach as in the chain model case, we can also study the OTOC in the general model. Similar to the (1 + 1)-d case, one finds a Lyapunov exponent 2π β saturating the chaos bound, and a butterfly velocity v B (for translation invariant systems). When the spatial kernel S xy is short-ranged, the universal relation D = v 2 B 2πT still holds. More details of the higher-dimensional calculation is given in Appendix C.

Models with global symmetry
As we learned from Ref. [12,23,25], a complex fermion version of the SYK model can be defined, which has similar properties such as local critical two-point functions. The Hamiltonian of the complex fermion at chemical potential µ = 0 is H = j,k,l,m J jklm c † j c † k c l c m with 1 j < k N , 1 l < m N . J jklm are also independent random variables. Since one can always write complex fermion operators in Majorana operators, the complex SYK model can be viewed as a Majorana model with 2N Majorana fermions c j = 1 2 (χ j1 + iχ j2 ). The coefficients in this Majorana model are different from that of an SYK model because of the U(1) symmetry requirement.
It is natural to generalize the single-site complex fermion SYK model to models defined on a general graph with a generic global symmetry. The most general form of the Hamiltonian is given by H = x,y,z,w∈Γ N j,l,k,m=1 a,b,c,d;P J P jklm,xyzw η abcd P χ j,a,x χ k,b,y χ l,c,z χ m,d,w with j = 1, 2, ..., N , a = 1, 2, ..., L labels N L Majorana fermions {χ j,a,x } at each site x ∈ Γ. The index a carries a representation of a global symmetry group. (For Majorana operators the representation has to be real. In other words, χ i,a,x is transformed under a subgroup of SO(L).) η abcd P for each P is an invariant rank-4 tensor in the symmetry group. For example, for symmetry group U(1) SO(2), there are two possible invariant tensors η abcd 0 = δ ab δ cd , η abcd 1 = ab cd . 13 We assume the couplings J P ijkl are independent variables 13 A natural way to define η abcd P is to write η abcd P = f P m ab f P m cd , with f P m ab the Clebsch-Gordan coefficients that maps the representation carried by ab to the representation P . However, it should be noted that the choice may be redundant. In other words, due to the anti-commutation of Majorana fermion operators, one may not need all representations P to expand an invariant tensor. with J P ijkl,xyzw = 0, J P ijkl,xyzw Once the Hamiltonian is defined one can try to study the effective action in the same way as before. Now we have to introduce the matrix bilocal field and the corresponding Lagrange multiplier Σ ab x (τ 1 , τ 2 ). The effective action after integrating out fermions is Although the effective action is more complicated, the large N saddle point approximation still apply, at least if we keep L finite in the large N limit. 14 The saddle point condition of this effective action gives the Schwinger-Dyson equations of G and Σ. In general there may be saddle points where G ab and Σ ab have off-diagonal matrix elements, which we don't know how to solve analytically. However, there always exists a diagonal saddle point solution G ab x (τ 1 , τ 2 ) = δ ab G s (τ 12 ), for which the effective action reduces to the SYK model with a coupling Therefore one can always take an expansion around this saddle point and study its stability using our knowledge about SYK model solution. If this saddle point is stable, these models have the same locally critical two-point functions as the SYK model, but have different 1 N fluctuations since the fluctuation g ab x (τ 1 , τ 2 ) is a matrix field. This will lead to different four-point functions. We will leave more systematic investigation of these models to future works.

General q
Another type of generalization one can consider is to include interactions that involve q fermions at a time, rather than just four. For q 4, the detailed analysis of the chain model will be very similar to q = 4, with a few minor changes in numerical coefficients. The relation v 2 B = 2πD/β will remain correct at large coupling. One can also consider including terms with different values of q in the same model. In general, the terms with higher values of q are more RG-irrelevant, so the terms with the lowest values of q will dominate in the infrared. Nevertheless, the subleading terms can have interesting effects. Suppose we have terms with q 1 , q 2 , with q 1 < q 2 , and we use the same dimensionful coupling J for both interactions. Then a perturbative analysis of the Schwinger Dyson equations indicates that for τ β, τ J 1 we will have (omitting coefficients) The first term is the naive conformal limit in the pure q 1 theory. The second term is the correction to the conformal limit, again in the pure q 1 theory; this term leads to the α K correction to the kernel that eventually gives the action for the reparameterization modes. The third term is the new feature of the theory with both q 1 and q 2 . It reflects the contribution of the irrelevant q 2 -fermion operator, which affects the correlator at quadratic order. The dimension of this operator near the IR is ∆ = q 2 q 1 , and when ∆ < 3 2 we have p < 2 q 1 + 1, so this third term will dominate over the second in (74). In such a case, the analysis of the reparameterization action would have to be redone, following e.g. appendix D of [9]. It would be interesting to compute the energy and chaos dynamics in the resulting model.

Conclusion and discussion
In this paper, we generalize the (0 + 1)-dimensional SYK model of Majorana fermions to higher spatial dimensions. The generalized model retains many interesting properties of the SYK model, such as local criticality, extensive zero temperature entropy and maximal chaos. On top of that, the spatial locality of our generalized model leads to many new physical properties. We find that single Majorana fermions in our model do not propagate between different sites, but collective modes made by pairs of fermions have nontrivial spatial dynamics. In particular, the most important collective mode in the low-energylong-wavelength limit is the time reparameterization field, the dynamics of which describes the diffusion of energy in this system, with a temperature independent diffusion constant D. This result tells us that our model describes a strongly correlated diffusive metal. The dynamics of the same reparameterization field also determines OTOC of fermion operators, from which we can also study the butterfly effect in this model. Our result shows that chaos spreads in space with a "butterfly velocity" v B . Remarkably, at strong coupling, the diffusion constant and the butterfly velocity satisfies a simple relation D = v 2 B /2πT , in consistency with the proposal in the literature about incoherent metals [29][30][31].
Our model pointed out a new class of solvable interacting lattice models in condensed matter physics. Usually, solvable models are mapped to weakly interacting theories such as mean field theories, so that they are not "chaotic", while the interesting phenomena in chaotic systems cannot be studied in solvable models. The generalizations of SYK model is a rare example of solvable but still chaotic systems. Therefore this model provides an interesting platform for studying various properties of strongly correlated systems, such as thermalization, entanglement propagation, dissipative transport, etc. It is also natural to ask whether further generalizations of these models allow us to investigate the possibility of many-body localization and phase transition between localized and delocalized phases.
From the perspective of holographic duality, the generalized SYK models might be considered as models that are dual to some kind of incoherent black hole (see [31] and references therein), but the details of this duality require further work. At strong coupling, the models do share a key property with conventional holographic systems, which is that a single set of degrees of freedom dominate and describe both energy diffusion and the chaos behavior. In a holographic theory, the relevant degrees of freedom are the bulk gravitational field. Here, it is a reparameterization of time that can vary from place to place. It would be interesting to derive the action (34) for these degrees of freedom from a subset of the metric degrees of freedom on some black hole background, in a similar way as the derivation of (0+1)-dimensional Schwarzian action from the Einstein-dilaton theory in approximately AdS 2 background [7,9]. It would also be interesting to understand the relationship of this action to recent work on hydrodynamic actions [40,41]. Another natural question is whether there are higher-dimensional translation invariant generalizations of SYK models which are dual to weakly coupled gravity theories in the bulk.
cussions. We especially acknowledge Subir Sachdev for helpful discussions and comments on the draft. This work is supported by the National Science Foundation through the grant No. DMR-1151786 (YG and XLQ). D.S. is supported by the Simons Foundation grant 385600.

A Diagrammatic derivation for four-point functions
In this section, we present a diagrammatic derivation of the fermion four-point functions. The four-point functions are the leading correlation functions which couple different sites. In the language of collective field G x (τ 1 , τ 2 ), the four-point function comes from its quantum fluctuations. The four-point functions are essential for understanding transport properties and also the OTOC measure of chaos, as discussed in section 4 and 5.
The connected four-point function is defined as We follow the main text to use G s to denote the saddle point of bilocal field, which is the Green's function here. Note that the four-point function is non-vanishing after disorder average only if the spatial and flavor indices appear in pairs, i.e. SO(N ) singlet. Similar to the SYK model, in the large N limit, the leading contributions to the connected piece of four-point functions are of order 1/N and consist of the ladder diagrams: where the thick lines are dressed Green's functions solved in the previous sections. The interaction vertices are paired by random disorder fields J x and/or J x . Comparing to those in the SYK model [4,22,11], the ladder diagrams in the chain model have extra labels indicating the spatial coordinates ( figure 11). Each ladder can either couple fermions on the same site, or bring two fermions at site x to the neighboring site x ± 1, as is shown in figure 12. The on-site terms are contributed by both J x and J x terms in the Hamiltonian, while the nearest neighbor terms are only from the J x term. In general, we can complete the summation of ladder diagrams (equation 76) using the Schwinger-Dyson equation: x, τ 1 y, τ 3 x, τ 2 y, τ 4 . . . Figure 11: Ladder diagrams connecting two sites far apart. To connect two sites with distance n = |x − y|, one need at least n ladders. And the four-point function includes all possible such diagrams and the partner terms with (τ 3 ↔ τ 4 ). x x y x y y x J (c) J -J contraction, y = x ± 1 Figure 12: Three types of the "rungs": type (a) is the same as the one appears in SYK, induced by the interactions between fermions at same sites; type (b) comes from the interactions between fermions at site x and nearest neighbor y = x ± 1, but the "rails" carry the same site indices, therefore, the effect of interaction doesn't propagate to next site; type (c) also comes from the interactions between nearest neighbor sites, and the rails get shifted by ±1.
where the gray box represents the dressed interaction vertex. In terms of algebraic formula, this equation is written as In this equation, the ladder kernel K = K xy (τ 1 , τ 2 ; τ 3 , τ 4 ) is treated as an operator acting on functions of one spatial coordinate x and two time variables. There are two kinds of nonzero matrix elements of K, given by the diagrams in figure 12 (a) (b) and those in (c) respectively: For the translation invariant saddle point solution, G s x (τ ) is independent from x, so that K 1 and K 2 are only different by the coefficient in front. Therefore K xy (τ 1 , τ 2 ; τ 3 , τ 4 ) has the separable form The spatial kernel S xy is a simple tight-binding hopping matrix ( i.e. an identity matrix plus a lattice Laplacian), and the temperal kernel K(τ 1 , τ 2 ; τ 3 , τ 4 ) is identical to that of the (0 + 1)-d SYK model with coupling constant J. (One should be reminded that we denote J = J 2 0 + J 2 1 .) The separable form of the kernel K xy (τ 1 , τ 2 ; τ 3 , τ 4 ) allows us to directly apply the results in the SYK model [4,11] to diagonalize the kernel and solve the four-point functions in our model. Assuming a formal diagonalization K xy = h,n,p k(h, n, p)|Ψ h,n,p Ψ h,n,p |, we can express the four-point function by the inner product: F = h,n,p |Ψ h,n,p 1 1−k(h,n,p) Ψ h,n,p |F 0 where |Ψ h,n,p = Ψ h,n,p (τ 1 , τ 2 , x) is some antisymmetric eigenfunctions in time which n labels the fourier mode for the sum of the two times, and h specifies the dependence on the difference of the times, and p labels the fourier mode for space. Technically, one can further simplify the calculation using the symmetrized kernel 15 [11] : where the symmetrized temporal kernel is defined as K(τ 1 , τ 2 , τ 3 , τ 4 ) = 3J 2 G s (τ 13 )·|G s (τ 34 )|· G s (τ 42 )·|G s (τ 21 )|. The simplification works by the following steps: (1) add two rungs (with absolute value for convenience) to the ladders in F (one on the left 16 , one on the right), see figure 13; (2) we get a multiple of "curved boxes" , then express it as a power of the 15 Roughly speaking, this trick is used to avoid the computation of the inner product Ψ|F 0 16 the extra factor 3J 2 S xx is needed to construct the kernel · . . . · = . . . Figure 13: A trick [11] to simplify the expressions by using the symmetrized kernel. symmetrized kernel; (3) in the end, sum over all ladders, which is now a geometric series: x 3J 2 S xx |G s (τ 12 )|F x y (τ 1 , τ 2 ; τ 3 , τ 4 )|G s (τ 34 | = 2 . . . where the factor of 2 comes from the counting for extra term with 3 ↔ 4. This is the general formal expression for the connected Euclidean four-point function. We can proceed by diagonalizing the spatial kernel via plane waves e ipx , i.e., k(h, n, p) = s(p)k(h, n) and write the general expression above in momentum space: where k(h, n) is the eigenvalue of temporal kernel. This agrees with equation 31 derived from effective action.

B Summation trick and the prefactor
In this appendix we determine the large t behavior of the contribution to the OTOC from the h = 2 modes. We choose a special configuration: τ 12 = τ 34 = β/2, | Re(τ )| < π/2 (87) where τ = τ 1 +τ 2 −τ 3 −τ 4 2 is the center of mass time separation, and the requirement | Re(τ )| < β/4 ensures the out-of-time order. The sum over h = 2 modes in equation (44) is We would like to continue τ = it and determine the large t behavior of this sum. To do this we consider the following integral Convergence is guaranteed by | Re(τ )| < β/4. The integrand has poles at ω = 1 and 2, 4, 6. . . . on R + . When we deform the contour to the right (see figure 14), we find Re ω Im ω The key point is that when we continue to large real time τ = it, the integral I remains convergent and non-growing in time, so we must have This directly gives (57).

C Diffusion and the butterfly velocity in general dimensions
In this section, we sketch the computation relevant to the diffusion and the butterfly velocity in general dimensional models with translation symmetry. We won't derive the exact formula for most general case, but will instead present key steps that determines the diffusion constant and butterfly velocity.
In the mode with transnational symmetry, the four-point function has simple expression in terms of momentum eigenvalue: F p( τ 1 , τ 2 ; τ 3 , τ 4 ) = 2 3J 2 G s (τ 12 )G s (τ 34 ) h,n k(h, n) 1 − k(h, n)s( p) Ψ h,n (τ 1 , τ 2 )Ψ * h,n (τ 3 , τ 4 ) Notice here the momentum p represents a general high dimensional vector. Further restrict the model to be local, we have a small p expansion for eigen-value s( p) 1 − j a j p 2 j . As we have noticed in the SYK and SYK chain model, the contribution relevant to the energy fluctuation arises from h = 2 modes, which corresponds to the reparametrization fields (pseudo-Goldstone modes). At strongly coupling limit, the eigenvalue for h = 2 receives a correction depends on n. k(h = 2, n) 1 − √ 2α K βJ |n| + . . .. Therefore, the pole that determines the diffusion constant has the simple form: this is the formula in Matsubara frequency ω n = 2π β n, after rotating to real time, we have diffusion pole: Next is to compute the butterfly velocity, which can be extract from the OTO correlation function. Again, we choose the time configuration which placed four time equally spacing around imaginary time circle, i.e., we are computing: χ j, x (t + i 3β 4 )χ k,0 (i β 2 )χ j, x (t + i β 4 )χ k,0 (0) Tr (rχ j, x (t)rχ k,0 (0)rχ j, x (t )rχ k,0 (0)) (95) Here the spatial coordinate x represents a high dimensional vector. Analogous to the computation in the chain model, we goes to momentum space, and plug in the h = 2 eigen-functions, which is the leading contribution: F ( p, t) h=2 ∝ n 2 even (−1) n/2 cos(n 2π β τ ) n 2 − 1 n 2πn β + j D j p 2 j (96) Using the trick in appendix B and analytic to real time, we have a formula in momentum space: Here we get a Lyapunov exponent λ L from h = 2 contribution. In general, the exponent will receive a 1/βJ correction as well as small p 2 correction from h = 2 part. Following the argument in section 5, we have In general dimension, we need to compute the following fourier transformation to get the butterfly velocity in x j direction: Notice t ∼ β log N is a large parameter. Therefore, we can integral over all direction except p j by saddle point approximation, where the saddle point is p k = 0 k = j: Then this fourier transformation essentially goes back to the 1-dimensional case, for which we know that at large x j , the integral is dominated by the pole 2π β + D j p 2 j = 0. The pole determines the exponential decaying profile for F (x j , t) at real space, with butterfly velocity v 2 B,j = 2πT D j (102)