The chiral SYK model

TL;DR

This paper introduces a 1+1D chiral SYK model with all-to-all interactions, showing exact solvability at large N, integrability at small N, and a transition to chaos as N increases, revealing rich quantum dynamics and chaos signatures.

Contribution

It generalizes the SYK model to a chiral 1+1D setting, providing exact solutions, integrability results, and insights into the transition from integrability to chaos with increasing N.

Findings

Exact Schwinger-Dyson solutions at large N

Integrability for N ≤ 6 via bosonization

Chaos signatures in OTOCs and Lyapunov exponents

Abstract

We study the generalization of the Sachdev-Ye-Kitaev (SYK) model to a $1+1$ dimensional chiral SYK model of $N$ flavors of right-moving chiral Majorana fermions with all-to-all random 4-fermion interactions. The interactions in this model are exactly marginal, leading to an exact scaling symmetry. We show the Schwinger-Dyson equation of this model in the large $N$ limit is exactly solvable. In addition, we show this model is integrable for small $N\le6$ by bosonization. Surprisingly, the two point function in the large $N$ limit has exactly the same form as that for $N=4$, although the four point functions of the two cases are quite different. The ground state entropy in the large $N$ limit is the same as that of $N$ free chiral Majorana fermions, leading to a zero ground state entropy density. The OTOC of the model in the large $N$ limit exhibits a non-trivial spacetime structure reminscent of that found by Gu and Kitaev for generic SYK-like models. Specifically we find a Lyapunov regime inside an asymmetric butterfly cone, which are signatures of quantum chaos, and that the maximal velocity dependent Lyapunov exponent approaches the chaos bound $2\pi/\beta$ as the interaction strength approaches its physical upper bound. Finally, the model is integrable for (at least) $N\le6$ but chaotic in the large $N$ limit, leading us to conjecture that there is a transition from integrability to chaos as $N$ increases past a critical value.

Equations

$$S_{\text{SYK}}=\int dt\left(\frac{i}{2}\sum _i{\psi }_i{\partial }_t{\psi }_i-i^{\frac{q}{2}}\sum _{i_1\cdots i_q}{J}_{i_1\cdots i_q}{\psi }_{i_1}\cdots {\psi }_{i_q}\right),$$

$$S=\int dtdx\mathcal{L}=\int dtdx\left(\frac{i}{2}\sum _{i=1}^N{\psi }_i({\partial }_t+{\partial }_x){\psi }_i+\sum _{1\le i<j<k<l\le N}{J}_{ijkl}{\psi }_i{\psi }_j{\psi }_k{\psi }_l\right),$$

$$\begin{array}{rl}& ⟨J_{ijkl}⟩=0,\\ & ⟨J_{ijkl}{J}_{i^{\prime }{j}^{\prime }{k}^{\prime }{l}^{\prime }}⟩=\frac{3!J^2}{(N-1)(N-2)(N-3)}{\delta }_{i{i}^{\prime }}{\delta }_{j{j}^{\prime }}{\delta }_{k{k}^{\prime }}{\delta }_{l{l}^{\prime }}(i<j<k<l,i^{\prime }<j^{\prime }<k^{\prime }<l^{\prime }),\end{array}$$

$$i({\partial }_t+{\partial }_x){\psi }_i+\sum _{1\le j<k<l\le N}{J}_{ijkl}{\psi }_j{\psi }_k{\psi }_l=0.$$

$$\begin{array}{rl}& [({\psi }_i{\psi }_j)(t,x),({\psi }_k{\psi }_l)(t,x^{\prime })]\\ =& \frac{i}{2\pi }({\delta }_{ik}{\delta }_{jl}-{\delta }_{il}{\delta }_{jk}){\partial }_x\delta (x-x^{\prime })-\delta (x-x^{\prime })({\delta }_{ik}{\psi }_j{\psi }_l-{\delta }_{il}{\psi }_j{\psi }_k+{\delta }_{jl}{\psi }_i{\psi }_k-{\delta }_{jk}{\psi }_i{\psi }_l),\end{array}$$

$$H=\int dx\mathcal{H},\mathcal{H}=-\frac{i}{2}\sum _{i=1}^N{\psi }_i{\partial }_x{\psi }_i-\sum _{1\le i<j<k<l\le N}{J}_{ijkl}{\psi }_i{\psi }_j{\psi }_k{\psi }_l,$$

$$T_0^x=\frac{i}{2}\sum _{i=1}^N{\psi }_i\left({\partial }_t-\frac{J^2}{4{\pi }^2}{\partial }_x\right){\psi }_i.$$

$$T_x^0=-\frac{i}{2}\sum _{i=1}^N{\psi }_i{\partial }_x{\psi }_i,T_x^x=\frac{i}{4}\sum _{i=1}^N{\psi }_i({\partial }_t-{\partial }_x){\psi }_i.$$

$$S=\int dtdx\left(i\sum _{\sigma =\uparrow ,\downarrow }{c}_{\sigma }^{†}({\partial }_t+{\partial }_x)c_{\sigma }-J{n}_{\uparrow }{n}_{\downarrow }\right),$$

$$c_{\sigma }(t,x)={\eta }_{\sigma }{e}^{i{\varphi }_{\sigma }(t,x)},c_{\sigma }^{†}(t,x)={\eta }_{\sigma }^{†}{e}^{-i{\varphi }_{\sigma }(t,x)}(\sigma =\uparrow ,\downarrow ),$$

$$[{\varphi }_{\sigma }(t,x),{\varphi }_{{\sigma }^{\prime }}(t,x^{\prime })]=i\pi {\delta }_{\sigma {\sigma }^{\prime }}\text{sgn}(x-x^{\prime }),$$

$$n_{\sigma }=c_{\sigma }^{†}{c}_{\sigma }=\frac{{\partial }_x{\varphi }_{\sigma }}{2\pi },-i{c}_{\sigma }^{†}{\partial }_x{c}_{\sigma }=\frac{({\partial }_x{\varphi }_{\sigma }{)}^2}{4\pi }.$$

$$S=-\int dtdx\left(\frac{1}{4\pi }\sum _{\sigma =\uparrow ,\downarrow }{\partial }_x{\varphi }_{\sigma }({\partial }_t{\varphi }_{\sigma }+{\partial }_x{\varphi }_{\sigma })+\frac{J}{4{\pi }^2}{\partial }_x{\varphi }_{\uparrow }{\partial }_x{\varphi }_{\downarrow }\right).$$

$$S=-\int dtdx\left(\sum _{\alpha =\pm }\frac{1}{4\pi }{\partial }_x{\varphi }_{\alpha }({\partial }_t{\varphi }_{\alpha }+u_{\alpha }{\partial }_x{\varphi }_{\alpha })\right),$$

$$u_{\pm }=1\pm \frac{J}{2\pi },$$

$$⟨{\varphi }_{\alpha }(t,x){\varphi }_{{\alpha }^{\prime }}(0,0)⟩=-{\delta }_{\alpha {\alpha }^{\prime }}\log \left[2\pi i(u_{\alpha }t-x-i{0}^{+})\right],$$

$$\begin{array}{rl}⟨{\varphi }_{\alpha }(t,x){\varphi }_{{\alpha }^{\prime }}(0,0){⟩}_{\beta }& =\sum _{m=-\infty }^{\infty }⟨{\varphi }_{\alpha }(t+im\beta ,x){\varphi }_{{\alpha }^{\prime }}(0,0){⟩}_{\beta =\infty }\\ & =-{\delta }_{\alpha {\alpha }^{\prime }}\log \left\{2i\beta u_{\alpha }\sinh \left[\frac{\pi }{\beta }\left(t-u_{\alpha }^{-1}x-i{0}^{+}\right)\right]\right\}.\end{array}$$

$$T_0^0=\frac{1}{4\pi }\sum _{\alpha =\pm }{u}_{\alpha }{\partial }_x{\varphi }_{\alpha }{\partial }_x{\varphi }_{\alpha },T_0^x=-\frac{1}{4\pi }\sum _{\alpha =\pm }{u}_{\alpha }{\partial }_t{\varphi }_{\alpha }{\partial }_x{\varphi }_{\alpha },$$

$$\mathcal{E}=\langle T^{0}_{\ 0}\rangle_{\beta}=-\frac{1}{4\pi}\sum_{\alpha=\pm}u_{\alpha}\partial_{x}^{2}\langle\phi_{\alpha}(t,x)\phi_{\alpha}(0,0)\rangle_{\beta}\Big{|}_{t\rightarrow 0,x\rightarrow\epsilon}=\frac{-1}{2\pi\epsilon^{2}}+\frac{\pi}{12\beta^{2}}(u_{+}^{-1}+u_{-}^{-1})+\mathcal{O}(\epsilon^{2})\ ,$$

$$j_{\mathcal{E}}=\langle T^{x}_{\ 0}\rangle_{\beta}=\frac{1}{4\pi}\sum_{\alpha=\pm}u_{\alpha}\partial_{t}\partial_{x}\langle\phi_{\alpha}(t,x)\phi_{\alpha}(0,0)\rangle_{\beta}\Big{|}_{t\rightarrow 0,x\rightarrow\epsilon}=-\frac{u_{+}^{2}+u_{-}^{2}}{4\pi\epsilon^{2}}+\frac{\pi}{6\beta^{2}}+\mathcal{O}(\epsilon^{2})\ .$$

$$\mathcal{E}=⟨T_0^0{⟩}_{\beta }=\frac{\pi }{12{\beta }^2}(u_{+}^{-1}+u_{-}^{-1}),j_{\mathcal{E}}=⟨T_0^x{⟩}_{\beta }=\frac{\pi }{6{\beta }^2}.$$

$$G_{\beta }(t,x)=\frac{1}{4}\sum _{i=1}^4⟨{\psi }_i(t,x){\psi }_i(0,0){⟩}_{\beta }=\frac{1}{4}\sum _{\sigma =\uparrow ,\downarrow }⟨c_{\sigma }(t,x)c_{\sigma }^{†}(0,0)+c_{\sigma }^{†}(t,x)c_{\sigma }(0,0){⟩}_{\beta }.$$

$$\begin{array}{rl}{G}_{\beta }(t,x)& =\frac{1}{2}\sum _{\sigma =\uparrow ,\downarrow }{e}^{⟨{\varphi }_{\sigma }(t,x){\varphi }_{\sigma }(0,0){⟩}_{\beta }}=\prod _{\alpha =\pm }{e}^{\frac{1}{2}⟨{\varphi }_{\alpha }(t,x){\varphi }_{\alpha }(0,0){⟩}_{\beta }}\\ =& \frac{1}{2i\beta \sqrt{u_{+}{u}_{-}}}\frac{1}{\sqrt{\sinh \left[\frac{\pi }{\beta }(t-u_{+}^{-1}x-i{0}^{+})\right]\sinh \left[\frac{\pi }{\beta }(t-u_{-}^{-1}x-i{0}^{+})\right]}},\end{array}$$

$$⟨T_0^0{⟩}_{\beta }=⟨-\frac{i}{2}\sum _{i=1}^4{\psi }_i{\partial }_x{\psi }_i-J{\psi }_1{\psi }_2{\psi }_3{\psi }_4{⟩}_{\beta }=\frac{i}{4}\sum _{i=1}^4⟨{\psi }_i({\partial }_t-{\partial }_x){\psi }_i{⟩}_{\beta }.$$

$$\mathcal{E}=\langle T^{0}_{\ 0}\rangle_{\beta}=-i(\partial_{t}-\partial_{x})G_{\beta}(t,x)\Big{|}_{t\rightarrow 0,x\rightarrow\epsilon}=\frac{1}{\pi\epsilon^{2}}+\frac{\pi}{12\beta^{2}}(u_{+}^{-1}+u_{-}^{-1})+\mathcal{O}(\epsilon^{2})\ .$$

$$j_{\mathcal{E}}=\langle T^{x}_{\ 0}\rangle_{\beta}=-2i\left(\partial_{t}-\frac{J^{2}}{4\pi^{2}}\partial_{x}\right)G_{\beta}(t,x)\Big{|}_{t\rightarrow 0,x\rightarrow\epsilon}=\frac{u_{+}^{-2}+u_{-}^{-2}}{2\pi\epsilon^{2}}+\frac{\pi}{6\beta^{2}}+\mathcal{O}(\epsilon^{2})\ ,$$

$${\mathcal{F}}_4(t_1,x_1,\cdots ,t_4,x_4)\equiv \frac{1}{4^2}\sum _{i,j=1}^4⟨{\psi }_j(t_1,x_1){\psi }_i(t_2,x_2){\psi }_j(t_3,x_3){\psi }_i(t_4,x_4){⟩}_{\beta }$$

$${\mathcal{F}}_4(t_1,x_1,\cdots ,t_4,x_4)=-\frac{1}{4}{G}_{\beta }^{13}{G}_{\beta }^{24}\left(\frac{G_{\beta }^{12}{G}_{\beta }^{34}}{G_{\beta }^{14}{G}_{\beta }^{23}}+\frac{G_{\beta }^{14}{G}_{\beta }^{23}}{G_{\beta }^{12}{G}_{\beta }^{34}}+\frac{{\tilde{G}}_{\beta }^{12}{\tilde{G}}_{\beta }^{34}}{{\tilde{G}}_{\beta }^{14}{\tilde{G}}_{\beta }^{23}}+\frac{{\tilde{G}}_{\beta }^{14}{\tilde{G}}_{\beta }^{23}}{{\tilde{G}}_{\beta }^{12}{\tilde{G}}_{\beta }^{34}}\right),$$

$${\tilde{G}}_{\beta }(t,x)=\frac{e^{\frac{1}{2}⟨{\varphi }_{+}(t,x){\varphi }_{+}(0,0){⟩}_{\beta }}}{e^{\frac{1}{2}⟨{\varphi }_{-}(t,x){\varphi }_{-}(0,0){⟩}_{\beta }}}=\sqrt{\frac{u_{-}\sinh \left[\frac{\pi }{\beta }(t-u_{-}^{-1}x-i{0}^{+})\right]}{u_{+}\sinh \left[\frac{\pi }{\beta }(t-u_{+}^{-1}x-i{0}^{+})\right]}}.$$

$${\mathcal{F}}_4(t,x)=\frac{1}{2{\beta }^2{u}_{+}{u}_{-}}\tanh \left[\frac{2\pi }{\beta }(t-u_{+}^{-1}x)\right]\tanh \left[\frac{2\pi }{\beta }(t-u_{-}^{-1}x)\right].$$

Full text

1. 1 Introduction

2. 2 The 1+1 dimensional chiral SYK model

3. 2.1 The Lagrangian

4. 2.2 Point splitting regularization and the energy momentum tensor

5. 3 Exact solution for $N=4$ via bosonization

6. 3.1 Bosonization of the action

7. 3.2 Fermion two point function for $N=4$

8. 3.3 Fermion four point function for $N=4$

9. 4 The chiral SYK model in the large $N$ limit

10. 4.1 Schwinger-Dyson equation and two point function

11. 4.1.1 The zero temperature solution

12. 4.1.2 Discussion on the real space UV regularization

13. 4.1.3 The finite temperature solution

14. 4.2 Thermal quantities and the ground state entropy density

15. 4.3 The OTOC, chaos exponent and butterfly cone

16. 5 The model at other small finite $N$

17. 5.1 Exact solution for $N=5$

18. 5.2 Exact solution for $N=6$

19. 5.3 $N\geq 7$ and thermalization of chiral Majorana fermions

20. 6 In closing

21. 7 Acknowledgments

22. A Bosonization and point splitting

23. A.1 Bosonization and normal ordering

24. A.2 Operator product expansion and commutation relation

25. A.3 Energy momentum tensor

26. B Saddle Point Approximation for the OTOC

1 Introduction

The last few years have seen an enormous amount of activity around the Sachdev-Ye-Kitaev (SYK) model [2, 3, 4, 5]. This activity is a consequence of the confluence of several themes. First, the model is partially solvable in a tractable large $N$ limit which it shares with the class of so-called tensor models [6, 7, 8] and the literature on dynamical mean field theory. This large $N$ limit is distinct from the large $N$ limits of vector and matrix models. Second, the SYK model exhibits many body quantum chaos or ergodicity as can be seen by computing the four point out of time order correlator (OTOC) in the same large $N$ limit. This has fed into the current intense interest in such foundational questions in quantum statistical mechanics. Finally, features of the model suggest that it is the holographic dual to some theory of quantum gravity in $1+1$ dimensions, and the prospect of a solvable model of holography has naturally excited enormous interest.

The SYK Model lives in $0+1$ dimensions—it is essentially a model on a quantum dot—and various authors have considered extensions to higher dimensions. One family of extensions builds higher dimensional models out of lattices of SYK dots, and this has led to interesting results on transport and on non-Fermi liquid behavior [9, 10, 11, 12, 13, 14, 15, 16]. A second family of extensions involves nonchiral homogeneous models with interactions in the new large $N$ limit which can be viewed either as arising from random couplings or from tensorial fields [17, 18, 19, 20, 21, 22]. In particular, the random extension to $1+1$ dimensional nonchiral fermions shows marginally irrelevant SYK interactions [19]. In this paper we present a third extension, this time to a homegenous chiral system in $1+1$ dimensions, where we study the resulting chiral SYK model at large $N$.

Before turning to our model, we briefly recapitulate some aspects of the original Sachdev-Ye-Kitaev (SYK) model [2, 3, 4, 5]. This model contains $N$ Majorana fermions $\psi_{i}$, which have all-to-all $q$-fermion random interactions, where $q\geq 2$ is an even number. The $0+1\equiv 1$d action of the SYK model is

[TABLE]

where $J_{i_{1}\cdots i_{q}}$ are a set of random Gaussian couplings satisfying $\langle J_{i_{1}\cdots i_{q}}^{2}\rangle=J^{2}(q-1)!/N^{q-1}$. In the large $N$ limit, the SYK model has a large ground state entropy due to the existence of highly degenerate spin glass like states. At low temperatures, an asymptotic conformal symmetry emerges from the collection of those states. The conformal symmetry is weakly broken and the corresponding soft mode physics is captured by a Schwarzian action. The low temperature dynamics governed by this action has been argued to be maximally chaotic or scrambled. In this context this means that the out-of-time-order correlators (OTOCs) grow as $1-{1\over N}e^{\lambda(t-t_{0})}$ with time $t$ and that the Lyapunov exponent $\lambda$ approaches the so-called “chaos bound”111In systems with locally bounded Hilbert spaces there does not appear to be such a Lyapunov regime, hence our qualifier. $2\pi/\beta$, where $\beta$ is the inverse of temperature [23, 24, 25]. It turns out that the Schwarzian action can also be viewed as emergent from the description of Near Extremal Black Holes when we consider gravitational dynamics near the AdS2 throat [26, 27, 28, 29, 30]. In the gravitational system, the black hole horizon provides a large ground state degeneracy and the conformal symmetry is the asymptotic symmetry of AdS2. With the Schwarzian action, the SYK model provides a concrete example of holographic dual system for near extremal black holes, although it is not clear whether the full model has an exact gravitational dual. Many striking conclusions about near extremal black holes have been drawn from this duality, such as the construction of a traversable wormhole [31, 32, 33] and black hole level statistics [34, 35, 36].

Turning now to our paper, we present a generalization of the SYK model to $N$ flavors of $1+1$ dimensional chiral Majorana fermions with all-to-all interactions (see action (2.1)), which we call the chiral SYK model. The interactions are uniform in the spacetime and random with respect to fermion flavors. Since fermions in $1+1$ dimensions have a scaling dimension $1/2$, any $q>4$ fermion interactions will be irrelevant, so we only focus on the $q=4$ fermion interactions. From the perspective of Hibert space dimension, the model we study here is the minimal generalization of the SYK model into higher spacetime dimensions. In condensed matter physics, chiral fermions can exist on the edge of $2+1$ dimensional gapped chiral topological states, and can interact with each other on the edge. For instance, chiral complex fermions occur on the edge of integer quantum Hall states [37, 38] or Chern insulators [39, 40], each of which is equivalent to two chiral Majorana fermions. Odd number of chiral Majorana fermions can be realized on the edge of non-Abelian fractional quantum Hall states [41], $p+ip$ chiral topological superconductors [42, 43] or chiral spin systems [44].

Unlike nonchiral systems, the 4-fermion interactions of $1+1$ dimensional chiral fermions are generically exactly marginal. Besides, such 4-fermion interactions break the Lorentz symmetry explicitly, since the product of 4 chiral fermion fields has conformal spin $2$ instead of zero. For example, in the celebrated chiral Luttinger model [45, 46, 47] of two complex chiral fermions, the 4-fermion interaction leads to a spin-charge velocity separation independent of energy scale, which is therefore exactly marginal; meanwhile, the Lorentz symmetry is absent since spin and charge excitations have different velocities. In our $1+1$ dimensional chiral SYK model in the large $N$ limit, the 4-fermion interactions are also exactly marginal and break the Lorentz symmetry. We show this explicitly by solving the two point function of the model exactly. Therefore, all the features of our model are energy scale independent.

The exact marginality of 4-fermion interactions in our chiral SYK model leads to a scaling symmetry, which enables us to solve the model exactly. In polar coordinates $re^{i\theta}=\tau+ix$ of imaginary time $\tau$ and spatial position $x$, the scaling symmetry restricts the zero temperature two point function to be proportional to $r^{-1}$, therefore the $1+1$ dimensional Schwinger-Dyson equations in the large $N$ limit reduce to $1$ dimensional equations in the boost angle $\theta$, which are exactly solvable. The exact two point function at finite temperature $\beta^{-1}$ can also be obtained by observations based on the zero temperature solution. We find that the two point function in real space takes a form of the product of two $0+1$ dimensional SYK propagators in the conformal limit moving at two different velocities $u_{\pm}$ determined by the interaction strength ($u_{+}>u_{-}>0$). With the exact solution, we calculate the ground state entropy by directly evaluating the large $N$ effective action (4.42), and find it is the same as the ground state entropy of $N$ flavors of free chiral Majorana fermions. Depending on the spatial boundary condition, the ground state entropy could be either [math] (anti-periodic boundary condition) or $N\log 2$ (periodic boundary condition, which has $N$ zero modes). In any case, the ground state entropy density vanishes, so there are no spin-glass like states in the model. We also investigate the behavior of the OTOC using the retarded kernel method in [4, 5], and identify a tilted butterfly cone where the OTOC grows exponentially along a space-time trajectory with fixed velocity. The butterfly cone exactly coincides with the causality cone, which has edges at velocities $u_{+}$ and $u_{-}$ (Figure 10(e)). The butterfly cone can be further divided into three regions which behave differently: for the region with small velocity $x/t$ near the left butterfly cone edge, the $1/N$ piece of the OTOC behaves as $\exp\left[{2\pi\over\beta}(u_{-}^{-1}x-t)\right]$. For the region in the middle of butterfly cone, the $1/N$ piece of the OTOC grows as $\exp\left[\overline{\lambda}t-{\#(t-v_{c}x)^{2}\over t}\right]$, where $0<\overline{\lambda}<{2\pi\over\beta}$, and $u_{-}<v_{c}<u_{+}$. When the interaction strength approaches its allowed upper bound (where $u_{-}\rightarrow 0$), the Lyapunov exponent $\overline{\lambda}$ approaches the chaos bound $2\pi/\beta$. For the large velocity region near the right butterfly cone edge, the $1/N$ piece of OTOC grows as $\exp\left[{2\pi\over\beta}(t-u_{+}^{-1}x)\right]$, where the exponent in the $t$ direction saturates the chaos bound $2\pi/\beta$. A similar saturation of the chaos bound near the butterfly edge was observed in [9], and recently explained on generic grounds by [1] using the ladder identity. However, we note that this is not the case near the left butterfly edge of our model, where the exponent in the $t$ direction is $-2\pi/\beta$, i.e., opposite to the chaos bound value. It is interesting that our model has both vanishing ground state entropy density and maximal Lyapunov exponent when the interaction strength approaches its upper bound, which is not a feature expected from black holes.

We also study the model at small finite $N$. In particular, in the cases $N=4,5$ and $6$, the model is seen to be integrable by bosonization. The $N=4$ model is nothing but the chiral Luttinger model with two complex fermions, which has a spin-charge separation. Surprisingly, we find that the $N=4$ fermion two point function is exactly of the same form as the large $N$ two point function. We conjecture that this coincidence may be due to the fact that the model has an SO($N$) symmetry for both $N=4$ and in the large $N$ limit (upon averaging over the random couplings). Together with the scaling symmetry, the SO($N$) symmetry may already pin down the form of two point function. Nevertheless, the four point function for $N=4$ is completely different from that in the large $N$ limit. This is because the $N=4$ case is integrable while the large $N$ case is chaotic. As a result, the four point function for $N=4$ only has power-law or exponential decay (up to a background constant from conserved charges), and is insensitive to time orderings. For $N=5$, upon a flavor basis rotation, the model decouples into an $N=4$ model plus a free chiral Majorana fermion. For $N=6$, by a flavor basis rotation and bosonization, we can rewrite the model as three free chiral bosons with three distinct velocities. For finite $N\geq 7$, the model becomes much more complicated and probably unsolvable. We conjecture that the model has a transition from integrablity to chaos when $N$ increases, and this transition probably happens between $N=6$ and $N=7$. This would also be the transition from non-thermalization to thermalization among the $N$ chiral Majorana flavors. In chiral topological condensed matter systems, thermalization of chiral states on the edge could affect the thermal measurements in experiments [48, 49, 50, 51].

While the bosonization of the $0+1$ dimensional SYK model has been attempted [52], we note that our chiral SYK model for any $N$ can always be written straightforwardly in the bosonized representation. For $N\leq 6$, it is advantageous to solve the model in the bosonized representation; while in the large $N$ limit, the fermion representation is more convenient.

The paper is organized as follows: In Sec. 2, we introduce the $1+1$ dimensional chiral SYK model for general number of chiral Majorana fermion flavors $N$, and find the energy-momentum tensor of the model. In Sec. 3, we solve the model for $N=4$ via bosonization. We then calculate the thermal quantities and the four point function for $N=4$, and show the OTOC in this case only has exponential decay. In Sec. 4, we investigate the chiral SYK model in the large $N$ limit by solving the Schwinger-Dyson equations exactly. Based on the exact two point function, we discuss its thermal properties and ground state entropy, and examine the large $N$ OTOC for signatures of chaos and operator spreading. In section 5, we give the exact solutions for $N=5$ and $N=6$, and discuss what might happen at $N\geq 7$ and its consequence for thermalization. Finally, we briefly summarize our results in Sec. 6.

2 The 1+1 dimensional chiral SYK model

In this section, we describe the $1+1$ dimensional chiral SYK model of chiral Majorana fermions which we will study in this paper. We will introduce the Lagrangian and the energy momentum tensor of the model. In particular, the model requires a proper choice of point splitting regularization for the operator product expansion (OPE) of the chiral Majorana fermion fields, which we will specify.

2.1 The Lagrangian

In $1+1$ dimensions, a free chiral Majorana fermion field $\psi$ has a scaling dimension $1/2$, so a $q$-fermion interaction term has a scaling dimension $q/2$, which is marginal if $q=4$, and irrelevant if $q>4$. This was considered as a major obstacle in the generalization of the SYK model to $1+1$ dimensions, since it is important to have the interaction relevant in the original $0+1$ dimensional SYK model. However, it is still worthwhile to examine the effect of the marginal random $4$-fermion interaction in $1+1$ dimensions, and whether some features of the SYK model are preserved.

In this paper, we study the direct generalization of the $q=4$ SYK model to $N$ flavors of $1+1$ dimensional right-moving chiral Majorana fermions $\psi_{i}$ ($1\leq i\leq N$), which has the action

[TABLE]

where $\psi_{i}$ satisfies the anticommutation relation $\{\psi_{i}(t,x),\psi_{j}(t,x^{\prime})\}=\delta_{ij}\delta(x-x^{\prime})$, and $\mathcal{L}$ denotes the Lagrangian density. The couplings $J_{ijkl}$ are real, and antisymmetric with respect to any two indices. In the large $N$ case, we assume $J_{ijkl}$ obey the random Gaussian distribution with

[TABLE]

where the coupling strength $J\geq 0$ is kept constant when we take the large $N$ limit. Note that here $J_{ijkl}$ is only random with respect to Majorana fermion flavor indices, while is uniform in the entire spacetime. This is fundamentally different from the chiral interactions studied in [10] which are spatially disordered. As we shall show, the spacetime uniform 4-fermion interaction $J_{ijkl}$ we study here is exactly marginal, thus is important at all energy scales.

The chiral 4-fermion interactions $J_{ijkl}$ explicitly break the Lorentz symmetry, and thus also break the conformal symmetry. This is because the chiral Majorana fermion $\psi_{i}$ carries a conformal spin $1/2$, thus the 4-fermion interaction term has a nonzero conformal spin $2$ and is not Lorentz invariant. The absence of Lorentz symmetry is, however, not a problem in condensed matter systems. For example, the $N$ chiral Majorana fermions can live on the edge of $N$ copies of the $2+1$d $p+ip$ chiral topological superconductor, which has no intrinsic Lorentz symmetry.

The Euler-Lagrange equation of motion of the action (2.1) can be easily found to be

[TABLE]

Throughout this paper, we will work in both real time $t$ and imaginary time $\tau=it$, depending on which one is more convenient.

2.2 Point splitting regularization and the energy momentum tensor

In quantum field theories (QFTs), the product of quantum fields at the same spacetime coordinate needs to be regularized carefully. The regularization is usually done by finding the operator product expansion (OPE) of two fields at split spacetime points $(t,x)$ and $(t^{\prime},x^{\prime})$, taking the limit $t^{\prime}\rightarrow t$ and $x^{\prime}\rightarrow x$, and keeping the normal ordered part of the limit. In QFTs without Lorentz invariance, the regularization may depend on the direction of point splitting, so one needs to specify the point splitting direction.

In this paper, we specify the point splitting to be in the $x$ direction whenever we regularize the product of field operators at the same spacetime point. Namely, for instance, we define the product of two chiral Majorana fermion fields as $\psi_{i}(t,x)\psi_{j}(t,x)=\lim_{x^{\prime}\rightarrow x}:\mathrel{\psi_{i}(t,x)\psi_{j}(t,x^{\prime})}:$, where $:\mathrel{\mathcal{O}}:$ stands for the normal ordering of operator $\mathcal{O}$. This point splitting choice is conventional in condensed matter systems, where there is generically an ultraviolet (UV) cutoff in the $x$ direction (lattice constant, etc) playing the role of the point splitting. Such a point splitting in the $x$ direction leads to the commutation relation on a constant time slice

[TABLE]

which can be derived from the OPE of free chiral Majorana fields $\psi_{i}$ (Appendix A). Note that $i\psi_{i}\psi_{j}$ ($1\leq i<j\leq N$) form the representation of the current operators of SO$(N)_{1}$, thus Eq. (2.4) is nothing but the commutation relation of the SO$(N)_{1}$ Kac-Moody algebra. For the model (2.1) to be self-consistently defined, we require this commutation relation on a constant time slice to persist when the chiral Majorana fermions $\psi_{i}$ are no longer free ($J>0$).

The Hamiltonian of the model can be obtained via the Legendre transformation as

[TABLE]

where $\mathcal{H}$ is the energy density. Note that the kinetic term can also be written into the Sugawara form, $-\frac{i}{2}\sum_{i}\psi_{i}\partial_{x}\psi_{i}=-\frac{\pi}{N-1}\sum_{i<j}:\mathrel{\psi_{i}\psi_{j}}::\mathrel{\psi_{i}\psi_{j}}:$ (Appendix A). Therefore, the kinetic term and the 4-fermion interaction term $\psi_{i}\psi_{j}\psi_{k}\psi_{l}$ are of the same 4-fermion form (with scaling dimension $2$ and conformal spin $2$). This suggests that the interaction couplings $J_{ijkl}$ are exactly marginal, which we will verify by calculations in the rest of our paper. This is in contrast to the nonchiral interaction in the Thirring model, which is either marginally irrelevant (for large $N$) [10] or marginally relevant [53].

We can also derive the energy momentum tensor $T^{\mu}_{\ \nu}$ of the model. By definition, $T^{0}_{\ 0}$ has the physical meaning of the energy density, so we have $T^{0}_{\ 0}=\mathcal{H}$ as given by Eq. (2.5). The energy current $T^{x}_{\ 0}$ can be derived from the energy conservation law $\partial_{t}T^{0}_{\ 0}+\partial_{x}T^{x}_{\ 0}=i[H,\mathcal{H}]+\partial_{x}T^{x}_{\ 0}=0$. Making use of the equation of motion (2.3), the commutation relation (2.4) and the correlation (2.2), we can obtain the energy current averaged over all random couplings $J_{ijkl}$ (Appendix A)

[TABLE]

In a similar way, we can derive the momentum density and the momentum current components of $T^{\mu}_{\ \nu}$ as (Appendix A)

[TABLE]

Due to the absence of Lorentz symmetry at $J>0$, the energy momentum tensor $T^{\mu}_{\ \nu}$ is not a symmetric tensor. We note that the expression of $T^{\mu}_{\ \nu}$ depends on the choice of the point splitting direction; here our results are for point splitting along the $x$ direction.

In a gapped $2+1$ dimensional condensed matter system with a chiral edge theory on its $1+1$ dimensional edge, the energy current $j_{\mathcal{E}}=\langle T^{x}_{\ 0}\rangle_{\beta}$ as a function of temperature $\beta^{-1}$ gives a thermal Hall conductance $\kappa_{xy}=\partial j_{\mathcal{E}}/\partial(\beta^{-1})$. For scaling invariant edge theories, the energy current is of the form $j_{\mathcal{E}}=\pi c/12\beta^{2}$, and yields a thermal Hall conductance $\kappa_{xy}=\pi c/6\beta$, where $c$ is a constant. For conformal field theories (CFTs), $c=c_{R}-c_{L}$ is the right-moving central charge $c_{R}$ minus the left-moving central charge $c_{L}$ [54]. In our case, the chiral model (2.1) is a CFT of $N$ free chiral Majorana fermions when $J=0$, which has $c_{R}=N/2$ and $c_{L}=0$. Therefore, we conclude that our model has $\kappa_{xy}=N\pi/12\beta$ when $J=0$. When $J>0$, the model is no longer a CFT but is still scaling invariant, and we need to calculate the coefficient $c$ using Eq. (2.6).

3 Exact solution for $N=4$ via bosonization

The $1+1$ dimensional chiral SYK model (2.1) with $N=4$ flavors of chiral Majorana fermions is integrable via bosonization, which is equivalent to the celebrated chiral Luttinger model with spin-charge separation in condensed matter physics. $N=4$ is also the minimal number of flavors which allows a 4-fermion interaction. In this section, we review this exact solution for $N=4$, and examine the two point function and OTOC of $\psi_{i}$, as a representative of non-chaotic integrable models. Furthermore, in Sec. 5 we shall show the model is also integrable for $N=5$ and $N=6$, but is probably no longer integrable for all $N\geq 7$.

Since there is only one coupling parameter $J_{1234}$ for $N=4$, it is meaningless to talk about the distribution of $J_{ijkl}$ as defined by Eq. (2.2). Instead, we simply set $J_{1234}=J\geq 0$ as a constant. Note that this definition of $J$ is consistent with that defined in Eq. (2.2) for $N=4$.

It is convenient to define two chiral complex fermion fields $c_{\uparrow}=(\psi_{1}+i\psi_{2})/\sqrt{2}$ and $c_{\downarrow}=(\psi_{3}+i\psi_{4})\sqrt{2}$, distinguished by a spin index $\uparrow$ or $\downarrow$. The action (2.1) for $N=4$ can then be rewritten as

[TABLE]

where $n_{\sigma}=c^{\dagger}_{\sigma}c_{\sigma}$ is the spin $\sigma$ fermion density. In this form, Eq. (3.1) is simply the chiral Luttinger model, which can be exactly solved by bosonization and exhibits a spin-charge separation physics. Physically, it could describe the edge of a Chern number $\nu=2$ quantum Hall state with a short-range repulsive interaction.

We comment that the $N=4$ action (3.1) has an SO(4) symmetry upon rotations among the four chiral Majorana fields $\psi_{i}$ [55, 56]. In contrast, there is generically no SO($N$) symmetry for $N\geq 5$ due to the presence of couplings $J_{ijkl}$ which are non-symmetric under SO($N$) rotations. Therefore, the $N=4$ model is more symmetric.

3.1 Bosonization of the action

Here we briefly review the bosonization procedure of the model action (3.1). The bosonization is done by redefining the fermion operators in terms of vertex operators

[TABLE]

where $\phi_{\sigma}(t,x)$ are two scalar boson fields satisfying the commutation relation

[TABLE]

with $\text{sgn}(x)$ denoting the sign of $x$. The coefficients $\eta_{\sigma}$ and $\eta_{\sigma}^{\dagger}$ are the Klein factors [57, 58, 59], which recover the anti-commutation relation between different flavors of fermion fields, and do not appear in fermion bilinears. By imposing a point splitting regularization in the $x$ direction, one can derive the following normal ordered identities in an OPE expansion (Appendix A):

[TABLE]

The action (3.1) can then be bosonized into

[TABLE]

Note that the action $S$ only contains bilinears of $\phi_{\sigma}$, so the boson fields $\phi_{\sigma}$ are free. By redefining two new boson fields $\phi_{\pm}=(\phi_{\uparrow}\pm\phi_{\downarrow})/\sqrt{2}$, we can rewrite the action as

[TABLE]

where $\phi_{+}$ and $\phi_{-}$ are two independent chiral boson fields decoupled from each other. The velocities $u_{+}$ and $u_{-}$ of the two chiral boson fields are given by

[TABLE]

respectively. From this solution, we see that the Lorentz symmetry is indeed explicitly broken when $J>0$, since there are two different “speeds of light” $u_{+}$ and $u_{-}$. Besides, we also see that the interaction $J$ is exactly marginal as we expected, which simply modifies the boson velocities $u_{\pm}$ in an energy scale independent way.

Eq. (3.7) also indicates a physical bound $0\leq J<2\pi$ for the interaction strength $J$. This ensures that both velocities $u_{\pm}$ are positive, and thus the total chirality of the model is kept invariant, as shown in Fig. 1(a). If $J>2\pi$, one would have the chirality of the boson field $\phi_{-}$ flipped as shown in Fig. 1(b), which is unphysical. This is because the action (3.1) is derived based on the normal ordering assumption that all the $k<0$ ($k>0$) fermion states are occupied (empty); while dispersions as shown in Fig. 1(b) would have both $k<0$ and $k>0$ occupied states and invalidate the normal ordering assumption.

We may still try to understand what happens when $J>2\pi$ if the $1+1$ dimensional chiral SYK model describes the edge theory of a bulk gapped $2+1$ dimensional condensed matter system. Physically, as long as the bulk gap does not close, the chirality of the edge should be preserved, as protected by the (Bogoliubov-de Gennes) Chern number of the ground state in the bulk. In fact, to understand the dispersion of $J>2\pi$, one has to take into account the irrelevant nonlinearities at high energies in a physical condensed matter system. For instance, a condensed matter system always has a spatial UV cutoff $\ell_{0}$ (lattice constant, etc), so one would expect the delta function interaction $Jn_{\uparrow}n_{\downarrow}$ to be smeared into a finite range $\ell_{0}$ interaction $\bar{J}(x-x^{\prime})n_{\uparrow}(x)n_{\uparrow}(x^{\prime})$, and $J=\int dx\bar{J}(x)$ only describes the interaction at the low energy limit. The model is then no longer scaling invariant. This will modify the velocities of the two boson fields $\phi_{\pm}$ into $u_{\pm}(k)=1\pm\tilde{J}(k)/2\pi$, where $k$ is the momentum, and $\tilde{J}(k)$ is the Fourier transform of $J(x)$. Due to the cutoff $\ell_{0}$, one has $\tilde{J}(k)\rightarrow J$ when $|k|\ll\ell_{0}^{-1}$, and $\tilde{J}(k)\rightarrow 0$ when $|k|\gtrsim\ell_{0}^{-1}$. Therefore, when $J>2\pi$, the dispersions of the edge chiral bosons up to very high energies will be as shown in Fig. 1(c), where the shaded area stands for the bulk states. In particular, the velocity of $\phi_{-}$ is negative near $k=0$, but goes back to positive when $|k|\gtrsim\ell_{0}^{-1}$. Therefore, the ground state of $J>2\pi$ is severely changed from that of $J<2\pi$, and the low energy theory becomes three right-moving chiral bosons and one left-moving chiral boson, which keeps the total chirality preserved but is no longer purely chiral.

In this paper, we shall only discuss interaction strengths $J$ inside the physical bound $0\leq J<2\pi$, so that we do not need to worry about any UV nonlinearities, and the model is purely chiral and exactly scaling invariant. At zero temperature, the bosonic two point function (without time ordering) of $\phi_{\alpha}$ from the action (3.6) can be derived as

[TABLE]

where $\alpha,\alpha^{\prime}=\pm$, and $0^{+}$ stands for a real positive infinitesimal number. Since the bosons are free, their two point function at finite temperature $\beta^{-1}$ can be derived by summing over all the zero temperature two point functions at time $t+im\beta$ ($m\in\mathbb{Z}$):

[TABLE]

All the $n$ point functions of boson fields $\phi_{\pm}$ can then be obtained from the above two point function and the Wick’s theorem.

The bosonic action (3.6) allows us to calculate the energy momentum tensor $T^{\mu}_{\ \nu}$ and thermal quantities in terms of boson fields $\phi_{\alpha}$. For instance, the energy density and energy current are given by

[TABLE]

which can be derived from Noether’s theorem. With the bosonic two point function in Eq. (3.9), we can calculate the values of $T^{0}_{\ 0}$ and $T^{x}_{\ 0}$ at temperature $\beta^{-1}$ by performing a point splitting $\epsilon$ in the $x$ direction, and then take the limit $\epsilon\rightarrow 0$. This leads to an energy density

[TABLE]

and an energy current

[TABLE]

The $\epsilon^{-2}$ terms in both expressions come from the vacuum expectation values, thus should be subtracted. Therefore, we find the physical energy density and energy current are

[TABLE]

In the bosonized representation here, different choice of point splitting direction in the above calculations will only affect the $\epsilon^{-2}$ terms, and does not affect the physical result of $\mathcal{E}$ and $j_{\mathcal{E}}$.

Eq. (3.13) leads to a thermal Hall conductance $\kappa_{xy}=\partial j_{\mathcal{E}}/\partial(\beta^{-1})=\pi/3\beta$. This is in agreement with the expectation that a gapped condensed matter system with $N$ flavors of chiral Majorana fermions on the edge has a quantized thermal Hall conductance $\kappa_{xy}=N\pi/12\beta$ [54], where $N=4$ here. We can also derive the thermal entropy density $\mathcal{S}$ from $\beta^{-1}[\partial\mathcal{S}/\partial(\beta^{-1})]=\partial\mathcal{E}/\partial(\beta^{-1})$, which gives $\mathcal{S}=(\pi/6\beta)(u_{+}^{-1}+u_{-}^{-1})$. The zero temperature entropy density of the free bosons should be zero, which is used here.

3.2 Fermion two point function for $N=4$

Since we are mainly interested in fermion correlation functions in this paper, it is useful to refermionize the above bosonized picture for $N=4$. First, we derive the average fermion two point function (without time ordering) at temperature $\beta^{-1}$, which is defined as

[TABLE]

By Eq. (3.2), we can replace $c_{\sigma}$ and $c_{\sigma}^{\dagger}$ by vertex operators $e^{i\phi_{\sigma}}$ and $e^{-i\phi_{\sigma}}$, respectively. Then by Wick’s theorem, one has $\langle e^{\pm i\phi_{\sigma}(t,x)}e^{\mp i\phi_{\sigma}(0,0)}\rangle_{\beta}=e^{\langle\phi_{\sigma}(t,x)\phi_{\sigma}(0,0)\rangle_{\beta}}$, so the fermion two point function can be derived as

[TABLE]

where we have used the two point function of boson fields $\phi_{\pm}=(\phi_{\uparrow}\pm\phi_{\downarrow})/\sqrt{2}$ in Eq. (3.9).

To verify Eq. (3.15) is the correct fermion two point function, we can calculate the finite temperature energy density $\mathcal{E}=\langle T^{0}_{\ 0}\rangle_{\beta}$ and energy current $j_{\mathcal{E}}=\langle T^{x}_{\ 0}\rangle_{\beta}$ in the fermion basis using $G_{\beta}(t,x)$ in Eq. (3.15). By the expression of $\langle T^{0}_{\ 0}\rangle_{\beta}=\mathcal{H}$ in Eq. (2.5) and the equation of motion (2.3), we can rewrite the energy density as

[TABLE]

This expression involves products of fermion fields at the same spacetime point, which therefore requires a point-splitting regularization, which is chosen to be the $x$ direction in our paper. By performing a point splitting $\epsilon$ in the $x$ direction, the expression becomes a derivative of the two point function $G_{\beta}$, and we find the energy density given by

[TABLE]

Similarly, the energy current $j_{\mathcal{E}}=\langle T^{x}_{\ 0}\rangle$ defined by Eq. (2.6) can also be rewritten as a derivative of $G_{\beta}$ via a point splitting $\epsilon$ in the $x$ direction, and can then be computed as

[TABLE]

where the definition of $u_{\pm}$ in Eq. (3.7) is used. As one can see, the physical pieces of $\mathcal{E}$ and $j_{\mathcal{E}}$ proportional to $\beta^{-2}$ are consistent with the results in Eq. (3.13) derived in the boson picture, although their vacuum $\epsilon^{-2}$ terms do not agree with those in Eqs. (3.11) and (3.12). This is because fermions and bosons have different vacuum energies. We note that in the fermionic picture here, it is important to have the point splitting direction in the $x$ direction.

3.3 Fermion four point function for $N=4$

We can further calculate any fermionic $n$ point functions for $N=4$ by rewriting them in terms of vertex operators of $\phi_{\pm}$. Here we calculate the averaged four point function

[TABLE]

for $N=4$. One can rewrite the fermion fields $\psi_{i}$ as vertex operators $(e^{i\phi_{\sigma}}\pm e^{-i\phi_{\sigma}})/\sqrt{2}$ of the boson fields $\phi_{\uparrow,\downarrow}$. By Wick’s theorem, the $n$ point correlation of vertex operators $e^{iq_{\alpha}\phi_{\alpha}}$ satisfies $\langle\prod_{j=1}^{n}e^{iq_{j}\phi_{\alpha_{j}}}\rangle=\delta_{\sum_{j}q_{j},0}\exp\left(-\sum_{\langle ij\rangle}q_{i}q_{j}\langle\phi_{\alpha_{i}}\phi_{\alpha_{j}}\rangle\right)$. With this formula and some rearrangements, we arrive at

[TABLE]

where we have defined the abbreviated notation $G_{\beta}^{ij}=G_{\beta}(t_{i}-t_{j},x_{i}-x_{j})$ for the two point function given by Eq. (3.15), while the function $\widetilde{G}_{\beta}^{ij}=\widetilde{G}_{\beta}(t_{i}-t_{j},x_{i}-x_{j})$ is defined by

[TABLE]

One can easily check that Eq. (3.20) reduces to the free fermion four point function when $J=0$, in which case $\widetilde{G}_{\beta}(t,x)\equiv 1$. If we set $(t_{2}-i\frac{\beta}{2},x_{2})=(t_{4},x_{4})=(0,0)$ and $(t_{1}-i\frac{\beta}{2},x_{1})=(t_{3},x_{3})=(t+i\frac{\beta}{4},x)$, Eq. (3.20) becomes the definition of the regularized equal thermal circle separation OTOC (see Eq. (4.46)), which can be simplified as

[TABLE]

Fig. 2(a) shows normalized $\mathcal{F}_{4}(t,x)$ along the $x=0$ direction, the $t=0$ direction and the $t=x$ line when $J=1.2\pi$. Fig. 2(b) shows the snapshot of normalized $\mathcal{F}_{4}(t,x)$ as a function of $x$ at fixed time $t=2\beta/\pi$ (green dashed-dotted line) and $t=10\beta/\pi$ (blue solid line), respectively. At late time $t$, one can see that $\mathcal{F}_{4}(t,x)$ decays to $-1$ in the region $u_{-}t<x<u_{+}t$ as denoted by the shaded area in Fig. 2(c), while decays to $1$ in the regions $x<u_{-}t$ and $x>u_{+}t$. Along a direction $x=vt$ with velocity away from $u_{-}$ and $u_{+}$, the OTOC $\mathcal{F}_{4}(t,x)$ always decays exponentially in $t$ towards $1$ or $-1$ for $t\gtrsim\beta$. This can be seen more explicitly by noting that $\mathcal{F}_{4}(t,x)$ can be expanded in nonnegative powers of the exponentially decaying function $e^{-2\pi t/\beta}$. For integrable models, the OTOCs have no fundamental difference from the TOCs.

We note that the shaded region $u_{-}t<x<u_{+}t$ in Fig. 2(c) is exactly the causal cone of the model, which can be seen from the retarded fermion Green’s function

[TABLE]

obtained by analytical continuation of the two point function (3.15). Therefore, the behavior of the OTOC function $\mathcal{F}_{4}(t,x)$ is intuitively understandable. Outside the causal cone in the large $t$ or large $x$ limit, the two $\psi_{i}$ fields do not communicate with the two $\psi_{j}$ fields, so $\mathcal{F}_{4}(t,x)$ tends to two decoupled two point functions $-G^{13}_{\beta}G^{24}_{\beta}$, which is a constant. Inside the causal cone, one would expect the correlations between the two $\psi_{i}$ fields and the two $\psi_{j}$ fields to be large, and the OTOC $\mathcal{F}_{4}(t,x)$ should strongly deviate from the decoupled value outside the causal cone.

4 The chiral SYK model in the large $N$ limit

Now we proceed to study the $1+1$ dimensional chiral SYK model in the large $N$ limit, and examine its analogy and difference compared with the $0+1$ dimensional SYK model. We show that the two point function of chiral Majorana fermions is exactly solvable in the large $N$ limit, and the interaction strength $J$ is exactly marginal. We will also show the model with a nonzero interaction strength within the physical parameter range $0<J<2\pi$ is chaotic in the large $N$ limit, and its maximal velocity dependent Lyapunov exponent in the OTOC approaches the maximal chaos bound $2\pi/\beta$ when $J\rightarrow 2\pi$. The butterfly cone of the OTOC is asymmetric. For any $0<J<2\pi$, the time direction Lyapunov exponent at the right butterfly edge always saturates the chaos bound $2\pi/\beta$, while at the left butterfly edge the time direction Lyapunov exponent is always $-2\pi/\beta$.

4.1 Schwinger-Dyson equation and two point function

In this subsection, we compute the two point function of the $1+1$ dimensional chiral SYK model in the large $N$ limit by exactly solving the Schwinger-Dyson equation. For convenience, we will work in Euclidean spacetime where the time $\tau=it$ is imaginary. Since the model is translationally invariant in both space and time, its two point function $G(-i\tau_{1},x_{1};-i\tau_{2},x_{2})=G[-i(\tau_{1}-\tau_{2}),(x_{1}-x_{2})]$ will only depend on the spacetime coordinate difference of the two points. Here we define $G(-i\tau,x)$ to be the imaginary time ordered two point Green’s function averaged over all Majorana fermion flavors. For zero temperature, it is defined by

[TABLE]

where $T$ stands for time ordering, $\tau=it$ is the imaginary time, $\psi_{i}(\tau,x)=e^{H\tau}\psi_{i}(x)e^{-H\tau}$ is the Majorana field at imaginary time $\tau$, and $\Theta(\tau)$ is the Heaviside step function. For finite temperature $\beta^{-1}$, the definition becomes

[TABLE]

where $Z=\text{Tr}e^{-\beta H}$ is the partition function, and $T$ still represents time ordering. Due to the anticommutation of fermions, the two point function satisfies the Kubo-Martin-Schwinger condition $G(-i\tau,x)=-G(-i\tau-i\beta,x)$.

The analytical continuation of $G(-i\tau,x)$ allows us to derive different kinds of real time Green’s functions. The real time two point function without time ordering $G(t,x)$ can be obtained by substituting $\tau\rightarrow it+0^{+}$, namely, $G(t,x)\equiv G(t-i0^{+},x)$. In Sec. 4.3, we will need to use the retarded Green’s function $G_{R}(t,x)=N^{-1}\sum_{i}\langle\{\psi_{i}(t,x),\psi_{i}(0,0)\}\rangle=\Theta(t)[G(t-i0^{+},x)-G(t+i0^{+},x)]$, and the Wightman function with half thermal circle separation defined by $G_{lr}(t,x)=G(t-i\frac{\beta}{2},x)$. Besides, the real-time ordered Green’s function (i.e., the Feynman propagator) is given by $G_{T}(t,x)=G(t-i0^{+}\text{sgn}(t),x)$.

4.1.1 The zero temperature solution

We first solve the Schwinger-Dyson equation for the two point function $G(-i\tau,x)$ at zero temperature. As we will see, the scaling invariance allows us a short cut to obtain the zero temperature solution exactly.

The imaginary-time ordered Green’s function of free chiral Majorana fermions at zero temperature can be easily derived to be

[TABLE]

which has a Fourier transform

[TABLE]

Here $\omega_{\tau}$ is the frequency along the imaginary time $\tau$, and $k$ is the momentum along the spatial direction.

When the interaction strength $J>0$, similar to the $0+1$ dimensional SYK model, the two point function $G$ to the leading order of $1/N$ is contributed by the melon diagrams as shown in Fig. 3. The solid lines and dashed lines in Fig. 3 denote the free two point function $G^{f}$ and the average over $J_{ijkl}$, respectively. By defining a translationally invariant self energy $\Sigma(-i\tau,x)$, we can summarize the Feynman diagams into the Schwinger-Dyson equation for the imaginary time ordered two point function $G$ and the self energy $\Sigma$:

[TABLE]

The scaling invariance of the action (2.1) tells us that $G(-i\omega_{\tau},k)$ and $\Sigma(-i\omega_{\tau},k)$ at zero temperature have scaling dimensions $-1$ and $+1$, respectively. Therefore, one cannot ignore the $i\omega_{\tau}+k$ term in the first equation of (4.5) in any cases, which has the same scaling dimension as $\Sigma(i\omega_{\tau},k)$. As a consequence, the Schwinger-Dyson equation of the $1+1$ dimensional chiral SYK model here does not have a reparametrization symmetry, which is present at low energies in the $0+1$ dimensional SYK model Schwinger-Dyson equation because of the irrelevance of the $i\omega_{\tau}$ term therein.

The exact scaling invariance, however, helps us constrain the form of $G$ and $\Sigma$. By introducing the polar coordinates $(\kappa,\theta_{k})$ in momentum space

[TABLE]

we can write down the following ansatz for $G$ and $\Sigma$ based on their scaling dimensions:

[TABLE]

where $f(z)=-f(-z)$ is an odd function (corresponding to fermionic statistics) defined on the unit complex circle $|z|=1$. The unknown part of the two point function $G$ and self energy $\Sigma$ is only its angular dependence on $\theta_{k}$, which comes form the Lorentz symmetry breaking by the interaction $J$. Note that the form of this ansatz already satisfies the first equation of Eq. (4.5).

The remaining task is to solve for the function $f(z)$ from the second equation of Eq. (4.5). This is an equation in real space, so we need to Fourier transform our ansatz (4.7). For convenience, we also use polar coordinates $(r,\theta)$ for real space defined by

[TABLE]

By integrating along the radial coordinate $\kappa$ first, we can write the Green’s function $G$ in real space as

[TABLE]

By defining $w=e^{-i\theta}$ and denoting $z=we^{i\theta_{k}}$, we can view the above integral as a contour integral along the circle $|z|=1$ in the complex $z$ plane. By assuming $f(z)$ can be analytically continued to the complex $z$ plane as a holomorphic function, and has no branch points or intrinsic singularities for $|z|\leq 1$, we can perform the contour integral to obtain

[TABLE]

where $f^{\prime}(z)$ denotes the first derivative of $f(z)$. The integral picks up the residues at the pole $z=i-i0^{+}$ as shown in Fig. 4 and other possible poles $\xi$ satisfying $\xi f(\xi)=i$ inside the unit circle. Whether or not such additional poles $\xi$ exist depends on the unknown function $f(z)$. The pole at $z=-i-i0^{+}$ is outside the contour and thus has no contribution.

The Fourier transform of the self energy $\Sigma$ to the real space can be calculated in a similar way. First, we can integrate over the radial momentum $\kappa$ to rewrite the Fourier transform of $\Sigma$ as a contour integral:

[TABLE]

where we again defined $w=e^{-i\theta}$ and $z=we^{i\theta_{k}}$, and the integral contour is the unit circle $|z|=1$ as shown in Fig. 4. In this case, the integrand has two poles at $z=\pm i-i0^{+}$, and other possible poles from the unknown function $f(z)$. By adding up the residues of poles inside the contour, we find the contour integration leads to a real space self energy function

[TABLE]

where $f^{\prime\prime}$ and $f^{\prime}$ are the second and first derivatives of $f$, respectively, and $\xi$ runs over all poles of the function $f(z)$ in the unit disk $|z|<1$.

For now we assume the function $f(z)$ in the unit disk $|z|<1$ satisfies $|zf(z)|<1$ and has no intrinsic singularities or branch points. We will come back to verify this after we obtain the solution of $f(z)$. This assumption ensures that $f(z)$ has at most one pole at $z=0$, and the last residue terms contributed by $|\xi|<1$ in Eq. (4.10) and in Eq. (4.12) all vanish. Therefore, the expressions of $G(-i\tau,x)$ in Eq. (4.10) and $\Sigma(-i\tau,x)$ in Eq. (4.12) are greatly simplified, and the second Schwinger-Dyson equation in Eq. (4.5) becomes a second order ordinary differential equation (ODE) for $f(z)$:

[TABLE]

By changing the variable $z$ to $s=z^{2}$, and defining a new function $g(s)=-izf(z)$, we can further simplify the above differential equation into

[TABLE]

Such a second order ODE can be solved by two steps of integration. First, we note that the second derivative of $g$ can be rewritten as \frac{d^{2}g}{ds^{2}}=\frac{dg}{ds}\frac{d}{dg}\Big{(}\frac{dg}{ds}\Big{)}=\frac{1}{2}\frac{d}{dg}\Big{[}\Big{(}\frac{dg}{ds}\Big{)}^{2}\Big{]}. Therefore, we can integrate Eq. (4.14) over $g$ to obtain

[TABLE]

where $c_{0}$ is the constant of integration. The square root of the above equation is a simple first order ODE. The second step is then to solve this first order ODE by direct integration, after which we arrive at a general solution

[TABLE]

where $s_{0}$ is the constant of integration in the second step. The function $f(z)$ is then given by

[TABLE]

Next, we need to determine the constants $c_{0}$ and $s_{0}$. When $J=0$, the system is free, so $f(z)$ should vanish identically, which is proportional to the self energy $\Sigma(-i\omega_{\tau},k)$. This fixes the constant $c_{0}=1$. The other constant $s_{0}$ is determined by the choice of point splitting. As we mentioned in Sec. 2.2, we have chosen the point splitting in the spatial $x$ direction when regularizing the OPE of $\psi_{i}$, which leads to the commutation relation (2.4). The commutation relation then requires the Green’s function $G(-i\tau,x)$ on the constant time slice $\tau=0$ to be the same as the free Green’s function, namely, $G(0,x)=G^{f}(0,x)=\frac{i}{2\pi x}$. From Eq. (4.10), this requires $f(z)=0$ when $z=\pm 1$. Therefore, we find the constant of integral $s_{0}=1$. This fixes the form of the solution to

[TABLE]

In particular, our assumption in solving the Schwinger-Dyson equation, that $|zf(z)|<1$ and $f(z)$ has no intrinsic singularities or branch points when $|z|<1$, is satisfied for $0\leq J<2\pi$. Therefore, $2\pi$ is the upper bound of $J$ for the above solution (4.18) to be self-consistent. We will discuss more on the bound of $J$ at the end of this subsection.

From Eqs. (4.10) and (4.12), we can then obtain the real space two point Green’s function in the large $N$ limit as

[TABLE]

and the real space self energy

[TABLE]

where the two velocities are $u_{\pm}=1\pm J/2\pi$. The above expressions of $G$ and $\Sigma$ have branch cuts so the function is not uniquely defined unless specified. Here we specify the branch of their definitions by requiring $G(-i\tau,x)\rightarrow{i\over 2\pi x}$ and $\Sigma(-i\tau,x)\rightarrow-{iJ^{2}\over 8\pi^{3}x^{3}}$ when $|x|\gg|\tau|$, and arranging the branch cut to be the straight line segment from $-iu_{-}\tau$ to $-iu_{+}\tau$. This procedure uniquely determines the functions $G$ and $\Sigma$ (except for the singular point $x=\tau=0$).

Remarkably, the two point function in the large $N$ limit in Eq. (4.19) has exactly the same form as that of $N=4$ given in Eq. (3.15) at zero temperature, except that the meanings of $J$ are different: here $J$ is the Gaussian average amplitude of the random couplings $J_{ijkl}$, while for $N=4$ $J$ is defined as the only coupling $J_{1234}$. The form of this solution also indicates that the coupling strength $J$ is exactly marginal, which determines the two velocities $u_{\pm}$ irrespective of the energy scale.

Besides, similar to the $N=4$ case, the large $N$ solution here also requires a physical bound of the interaction strength $0\leq J<2\pi$, which ensures the two velocities $u_{\pm}$ in Eq. (4.19) are positive. The physical reason for this bound is analogous to that we explained for $N=4$ in Fig. 1, namely, the total chirality of the system would no longer be preserved for $J>2\pi$ unless certain UV nonlinearities are taken into account. More explicitly, this can be seen from the zero temperature spectral weight obtained by

[TABLE]

where $G(\omega+i0^{+},k)$ is the momentum space retarded Green’s function obtained by doing an analytical continuation $i\omega_{\tau}\rightarrow\omega+i0^{+}$ of the imaginary time ordered Green’s function $G(i\omega_{\tau},k)$. Fig. 5 shows the spectral weight $A(\omega,k)$ at a fixed momentum $k>0$, which indicates the eigenstate energies of the model at momentum $k$ are distributed within the energy range $[u_{-}k,u_{+}k]$. Therefore, if $J>2\pi$, one would have $u_{-}<0$, and part of the eigenstates will have a negative velocity and reversed chirality. This would lead to a severe change of ground state similar to that illustrated by Fig. 1(c), which cannot be resolved unless certain UV nonlinearities are considered. Here we shall only restrict ourselves in the parameter range $0\leq J<2\pi$, so that we do not need to consider any UV nonlinearities.

4.1.2 Discussion on the real space UV regularization

in Sec. 4.1.1, we derived the zero temperature solution (4.19,4.20) by performing Fourier transforms of $G$ and $\Sigma$ from the momentum space to the real space. In particular, we determined the constant $s_{0}$ in Eq. (4.17) by considering the $x$ direction point splitting regularization we chose. This indicates that the real space UV regularization of this model has important physical consequences on correlation functions, and thus has to be handled correctly. An immediate example is to transform the real space $G$ and $\Sigma$ back to the momentum space: since the self energy $\Sigma$ diverges as $1/r^{3}$ in the real space, its Fourier transform requires a UV regularization, and different UV regularization schemes will give different results as we will see below. In this subsection, we will show that the correct UV regularization scheme for doing real space integrations in consistency with point splitting in $x$ direction is to take a UV cutoff of time $|\tau|\geq\epsilon$. The readers who are solely interested in the main results of the paper can skip this subsection and go directly to Sec. 4.1.3.

One may also have noticed that the solution (4.19) factorizes into the product of two conformal propagators of dimension ${1\over 4}$ and velocities $u_{\pm}$. Naively, this looks like two decoupled $0+1$ dimensional SYK systems moving at different velocities $u_{\pm}$. However, this cannot be true, since otherwise they will satisfy the low-energy approximate Schwinger-Dyson equation of the $0+1$ dimensional SYK model, rather than the exact $1+1$ dimensional Schwinger-Dyson equation (4.5) here. In fact, it is exactly the UV regularization we will discuss below that invalidates such a factorization into two $0+1$ dimensional SYK systems.

Let us start by considering the Fourier transformations of the real space $G$ and $\Sigma$ in Eqs. (4.19,4.20) back into the momentum space. As we have said, the Fourier transformation of $\Sigma$ requires a UV regularization. If the UV regularization is correct, we should be able to find that the Fourier transformed $G$ and $\Sigma$ satisfy the first Schwinger-Dyson equation in Eq. (4.5).

We first examine the regularization scheme of taking a UV cutoff $|\tau|\geq\epsilon$ in the time direction, where $\epsilon>0$ is a small number. In this case, the Fourier transformations of both $G$ and $\Sigma$ are proportional to the generic form

[TABLE]

where $\sigma$ is a constant. For $\sigma=-1/2$ and $\sigma=-3/2$, one obtains the momentum space two point function $G(-i\omega_{\tau},k)=\frac{1}{2\pi}W_{-\frac{1}{2}}(-i\omega_{\tau},k)$ and self energy $\Sigma(-i\omega_{\tau},k)=\frac{J^{2}}{8\pi^{3}}W_{-\frac{3}{2}}(-i\omega_{\tau},k)$, respectively. Without loss of generality, here we assume $k>0$. The integration of $x$ is then along the contour of the real $x$ axis closing in the lower half plane as shown in Fig. 6. When $\tau<0$, the contour integral is zero. When $\tau>0$, the contour encloses a branch cut between $-iu_{+}\tau$ and $-iu_{-}\tau$, thus can be deformed into a contour closely surrounding to the branch cut in Fig. 6. This enables us to rewrite the Fourier transformation as

[TABLE]

where we have defined $y=ix$. The integration can be further simplified by changing the variables of integration to $y_{+}=\frac{u_{+}\tau-y}{u_{+}-u_{-}}$ and $y_{-}=\frac{y-u_{-}\tau}{u_{+}-u_{-}}$, after which the region of integration is bound by $y_{+}\geq 0$, $y_{-}\geq 0$ and $y_{+}+y_{-}\geq\epsilon$. The integral then takes the form

[TABLE]

where we have defined $k_{\pm}=u_{\pm}k+i\omega_{\tau}$. The real parts of $k_{\pm}$ are always positive, so the integral is well convergent in the IR. We have two possible UV divergences, one is from the $\epsilon$ cutoff, and the other is from the end point [math] of the integral. For the latter, we can remind ourselves that the integral is actually a contour integral, and therefore we should take the principle value (or equivalently by analytic continuing $\sigma$). While for the former, if we neglect the $\epsilon$ piece, the two integrals of $y_{+}$ and $y_{-}$ would decouple, and we would get an answer factorized into the product of two decoupled $0+1$ dimensional SYK Green’s functions or self energies in the conformal limit along $y_{+}$ and $y_{-}$ direction. However, as we shall see, the $\epsilon$ piece is important to produce an additional contribution to the self energy $\Sigma$, which enables $\Sigma$ and $G$ to satisfy the $1+1$ dimensional Schwinger-Dyson equation exactly.

When $y_{+}<\epsilon$, the integral of $y_{-}$ yields a result $k_{+}^{-1-\sigma}\Gamma(1+\sigma,k_{+}(\epsilon-y_{+}))$, where $\Gamma(\sigma,\epsilon)=\int_{\epsilon}^{\infty}y^{\sigma-1}e^{-y}dy$ is the incomplete Gamma function. While when $y_{+}>\epsilon$, the integral of $y_{-}$ simply gives the usual Gamma function $k_{+}^{-1-\sigma}\Gamma(1+\sigma)$. By expanding in $\epsilon-y_{+}$, we have

[TABLE]

where we have used the definition of beta function $B(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}=\int_{0}^{1}y^{a-1}(1-y)^{b-1}dy$. For $\sigma=-1/2$, the integral (4.25) is regular as $\epsilon\rightarrow 0$, so the $\epsilon$ dependent terms can be ignored, and we find the Fourier transform of the two point function to be

[TABLE]

as expected from Eq. (4.7). This also indicates that the Fourier transform of $G$ does not require a UV regularization. For $\sigma=-3/2$, the leading $1/\epsilon$ divergent term of Eq. (4.25) vanishes, and we arrive at the Fourier transform of the self energy

[TABLE]

again in agreement with Eq. (4.7). We see that the UV cutoff $\epsilon$ gives rise to a regular piece $i\omega_{\tau}+k$, which is important for the answers to exactly satisfy the first Schwinger-Dyson equation in Eq. (4.5). This shows that UV cutoff $|\tau|\geq\epsilon$ is the correct regularization scheme, which we will give an understanding later in this subsection.

We could also try to use a different regularization scheme, for example, take a UV cutoff $|x|\geq\epsilon$ in the spatial direction. Following a similar calculation, one will arrive at an expression analogous to Eq. (4.24), except that the region of integration is bounded by $u_{-}y_{+}+u_{+}y_{-}\geq\epsilon$ and $y_{\pm}\geq 0$. This does not affect the two point function $G$ and the $(k_{+}k_{-})^{1/2}$ piece of the self energy $\Sigma$, but does change the UV contributed bare piece. More explicitly, it is easy to check that the resulting Fourier transforms in this case are

[TABLE]

[TABLE]

which obviously do not satisfy the Schwinger-Dyson equation (4.5). In fact, if we instead set $s_{0}=-1$ in the solution (4.17), which corresponds to a time $\tau$ direction point splitting, one will find that the spatial UV cutoff $|x|\geq\epsilon$ becomes the correct UV regularization.

The above discussion shows that the real space UV regularization scheme depends on the choice of point splitting direction, and incorrect UV regularization will lead to answers inconsistent with the Schwinger-Dyson equation. Here we give a physical understanding why UV cutoff $|\tau|\geq\epsilon$ is the correct regularization scheme for point splitting in the $x$ direction. A heuristic understanding is the following: taking the cutoff $|\tau|\geq\epsilon$ for the self energy $\Sigma(-i\tau,x)$ is equivalent to setting $\Sigma(0,x)=0$ on the constant time slice $\tau=0$. Recall that point splitting in the $x$ direction leads to the commutation relation (2.4), and requires the two point function $G(-i\tau,x)$ to reduce to the free two point function $G(0,x)=G_{f}(0,x)=\frac{i}{2\pi x}$ on the constant time slice $\tau=0$ (which is how $s_{0}=1$ in Eq. (4.17) is determined). This free nature of $G(0,x)$ therefore intuitively agrees with the vanishing of the self energy $\Sigma(0,x)=0$ on the time slice $\tau=0$.

A more intuitive explanation is illustrated by Fig. 7. In real space, the self energy $\Sigma(-i\tau,x)\sim J^{2}\langle(\psi_{j}\psi_{k}\psi_{l})(0,0)(\psi_{j}\psi_{k}\psi_{l})(\tau,x)\rangle$ is the correlation of the operator product $\psi_{j}\psi_{k}\psi_{l}$, the definition of which needs a point splitting regularization. For point splitting along the $x$ direction we have chosen, the three Majorana fields in $\psi_{j}\psi_{k}\psi_{l}$ should be splitted by a spatial spacing $\epsilon_{x}$ as shown in Fig. 7(a). Since each $\psi_{i}$ is a point-like operator, one could view the regularized $\psi_{j}\psi_{k}\psi_{l}$ as a “flat-shape” operator with size $\epsilon_{x}$ in the $x$ direction and $\epsilon_{\tau}$ in the $\tau$ direction, where $\epsilon_{\tau}/\epsilon_{x}\rightarrow 0$. With such a point splitting, the self energy correlation $\Sigma(-i\tau,x)$ is only meaningful when the two operators $(\psi_{j}\psi_{k}\psi_{l})(0,0)$ and $(\psi_{j}\psi_{k}\psi_{l})(\tau,x)$ do not overlap with each other, namely, when $(x/\epsilon_{x})^{2}+(\tau/\epsilon_{\tau})^{2}\geq 1$. Therefore, when calculating the Fourier transformation of $\Sigma(-i\tau,x)$, one should take a UV cutoff $(x/\epsilon_{x})^{2}+(\tau/\epsilon_{\tau})^{2}\geq 1$ with $\epsilon_{\tau}/\epsilon_{x}\rightarrow 0$, i.e., perform the Fourier integration outside the shaded area as shown in Fig. 7(b). This UV cutoff can be further simplified by noting that the zero temperature solution in Eqs. (4.19) and (4.20) is scaling invariant. Via a scaling transformation, we can scale $\epsilon_{x}\rightarrow\infty$ while keeping $\epsilon_{t}=\epsilon$ small, which still satisfies $\epsilon_{t}/\epsilon_{x}\rightarrow 0$. The UV cutoff is then simplified to $|\tau|\geq\epsilon$, i.e., outside the shaded area of Fig. 6(c). Therefore, this is the correct UV cutoff corresponding to point splitting in the $x$ direction.

It is interesting that in our case, even though the magnitude of the UV cutoff $\epsilon$ does not affect the final result, the spacetime direction along which we put the cutoff $\epsilon$ does have a physical consequence. This is not very common in high energy theory, since there we usually study interactions preserving the Lorentz symmetry, which guarantees that the direction of UV cutoff does not matter. In our particular example, the interaction breaks Lorentz invariance, and such a UV cutoff direction dependence arises.

4.1.3 The finite temperature solution

We now proceed to solve the Schwinger-Dyson equation of our $1+1$ dimensional chiral SYK model at finite temperature $\beta^{-1}$, which is given by Eq. (4.5) with $i\omega_{\tau}$ replaced by the Matsubara frequency $i\omega_{n}=\frac{(2n+1)\pi}{\beta}$ ($n\in\mathbb{Z}$). Unlike the zero temperature case, the temperature $\beta^{-1}$ sets an energy scale for the problem, and the Schwinger-Dyson equation is no longer scaling invariant. Therefore, the scaling invariant ansatz method in Sec. 4.1.1 no longer applies. Besides, the $1+1$ dimensional Schwinger-Dyson equation (4.5) has no reparameterization symmetry, so we cannot do conformal transformation from zero temperature to finite temperature. This makes the finite temperature problem much more difficult.

However, there are two hints from the zero temperature solution (4.19,4.20) which help us guess the finite temperature solution. The first hint is that the large $N$ zero temperature two point function $G$ has exactly the same form as that of $N=4$ (see Eq. (3.15) with $\beta\rightarrow\infty$). The second hint is that the large $N$ zero temperature $G(-i\tau,x)$ and $\Sigma(-i\tau,x)$ factorize into the product of two zero temperature two point functions and self energies of the $0+1$ dimensional SYK model in the conformal limit along $\tau-iu_{+}^{-1}x$ and $\tau-iu_{-}^{-1}x$ directions, respectively. Although in Sec. 4.1.2 we showed such a factorization is invalidated in the momentum space by an indispensable UV cutoff regularization contribution, we could expect it is valid in the IR (since the model is scaling invariant, any finite length scale is IR if no UV cutoff is put in by hand). This motivates us to guess that the above two facts still hold at finite temperature $\beta^{-1}$, so that the finite temperature two point function in the large $N$ limit is given by

[TABLE]

and the finite temperature self energy is $\Sigma_{\beta}(-i\tau,x)=J^{2}G_{\beta}(-i\tau,x)^{3}$ by the second Schwinger Dyson equation in Eq. (4.5). Namely, we guess the large $N$ finite temperature two point function $G_{\beta}$ is equal to the $N=4$ finite temperature two point function as well, and also factorizes into the product of two finite temperature $0+1$ dimensional SYK two point functions in the conformal limit along $\tau-iu_{+}^{-1}x$ and $\tau-iu_{-}^{-1}x$ directions. In the remainder of this subsection, we shall verify that Eq. (4.30) is indeed the large $N$ finite temperature solution.

We shall verify the validity of the two point function $G_{\beta}$ in Eq. (4.30) and the corresponding self energy $\Sigma_{\beta}$ by Fourier transforming them into the momentum space, and check whether they satisfy the first Schwinger-Dyson equation in Eq. (4.5). For convenience, we set $\beta=2\pi$ in the Fourier calculations below, while the answers for generic temperature $\beta^{-1}$ can be restored by a rescaling $(i\omega_{n},k)\rightarrow\frac{\beta}{2\pi}(i\omega_{n},k)$ in the end. As we have shown in Sec. 4.1.2, the correct UV regularization for the Fourier transform is to take a UV cutoff $|\tau|\geq\epsilon$. Therefore, we calculate the following general function

[TABLE]

with a cutoff $\epsilon$ for small $\tau$, where $\sigma$ is a constant. The Fourier transformations of the two point function and the self energy are by definition given by $G(i\omega_{n},k)={1\over 2\sqrt{u_{+}u_{-}}}W_{-1/2}(i\omega_{n},k)$ and $\Sigma(i\omega_{n},k)={J^{2}\over 32\pi^{2}(u_{+}u_{-})^{3/2}}W_{-3/2}(i\omega_{n},k)$, respectively.

We first change $\tau$ to the complex variable $z=e^{i\tau}$, so that the integration over $\tau$ in Eq. (4.31) becomes an integration of $z$ along the unit circle contour $|z|=1$ counterclockwise from $1+i\epsilon$ to $1-i\epsilon$ as shown in Fig. 8. At a given $x$, the integrand has three branch points [math], $e^{-u_{-}^{-1}x}$, and $e^{-u_{+}^{-1}x}$ on the nonnegative real axis of $z$. For half-integral $\sigma$ which we are interested in, $z=0$ is not a branch point. Therefore, we can deform the contour along $|z|=1$ into a contour closely bypassing the nonnegative real axis of $z$ from $1+i\epsilon$ to [math] and then to $1-i\epsilon$, as shown in Fig. 8. Eq. (4.31) can then be rewritten as

[TABLE]

where the integration over $z$ is from $1+i\epsilon$ to [math] above the real axis, and then from [math] to $1-i\epsilon$ below the real axis. For convenience, we denote $z=e^{-t}$, where $t=-i\tau$ is the Wick rotation of the imaginary time $\tau$. The deformed integral contour of $z$ is then equivalent to a contour of $t$ going from $-i\epsilon$ to $+\infty$ and then back to $+i\epsilon$ in straight lines.

For a given $z=e^{-t}$ (and $\sigma$ being a half integer), the integrand has branch cuts for $x$ along the line segments from $u_{-}(t+i{2\pi m\over\beta})$ to $u_{+}(t+i{2\pi m\over\beta})$ ($m\in\mathbb{Z}$). Without loss of generality, we assume the momentum $k>0$, so that the integration over $x$ is along a contour of the real $x$ axis and closed in the lower half plane at infinity. Therefore, the integration of $x$ picks up the contributions of branch cuts from $u_{-}(t+i{2\pi m\over\beta})$ to $u_{+}(t+i{2\pi m\over\beta})$ with $m\leq 0$ when the imaginary part of $t$ is infinitesimally negative (i.e., when $t$ is integrated from $-i\epsilon$ to $+\infty$), and picks up the contributions of branch cuts from $u_{-}(t+i{2\pi m\over\beta})$ to $u_{+}(t+i{2\pi m\over\beta})$ with $m<0$ when the imaginary part of $t$ is infinitesimally positive (i.e., when $t$ is integrated from $+\infty$ to $i\epsilon$). After $t$ is integrated, the contributions of all $m<0$ branch cuts will cancel, and only the contribution of the $m=0$ branch cut from $u_{-}t$ to $u_{+}t$ remains. Considering these facts, we can deform the $x$ contour into a contour closely surrounding the branch cut, and write Eq. (4.32) as

[TABLE]

Since the physical result should not depend on whether the cutoff $\epsilon$ is real or not, by adiabatical continuation we can replace the lower bound $-i\epsilon$ of the $t$ integral by $\epsilon>0$. To further simplify the expression, we define two new variables $y_{+}=t-u_{+}^{-1}x$ and $y_{-}=u_{-}^{-1}x-t$. The integral then becomes

[TABLE]

where we have defined $k_{\pm}={u_{\pm}\over u_{+}-u_{-}}\omega_{n}-i{u_{+}u_{-}\over u_{+}-u_{-}}k$, and $\Delta(\epsilon,y_{+})=\max[(\frac{u_{+}}{u_{-}}-1)\epsilon-\frac{u_{+}}{u_{-}}y_{+},0]$. Here the integration of $y_{\pm}$ are all done along the real axis. If the $\epsilon$ piece is ignored, by formula $\int_{0}^{\infty}dye^{ky}(1-e^{-y})=B(1+\sigma,-k)=\frac{\Gamma(1+\sigma)\Gamma(-k)}{\Gamma(1+\sigma-k)}$, we see that $W_{\sigma}$ factorizes into the product of two Beta functions. With a finite $\epsilon>0$, we can expand in $\Delta(\epsilon,y_{+})$, which follows exactly what we did from Eq. (4.24) to Eq. (4.25) in Sec. 4.1.2. This gives a result

[TABLE]

where we have defined $\tilde{\epsilon}=(\sqrt{\frac{u_{+}}{u_{-}}}-\sqrt{\frac{u_{-}}{u_{+}}})\epsilon$. For $\sigma=-1/2$, the $\epsilon$ pieces vanishes as $\epsilon\rightarrow 0$. For $\sigma=-3/2$, the leading $1/\epsilon$ term vanishes, while there is a constant piece contributed by the $\epsilon$ expansion. Therefore, we find the Fourier transformation $G(i\omega_{n},k)={1\over 2\sqrt{u_{+}u_{-}}}W_{-1/2}$ and $\Sigma(i\omega_{n},k)={J^{2}\over 32\pi^{2}(u_{+}u_{-})^{3/2}}W_{-3/2}$ to be given by

[TABLE]

and

[TABLE]

where we have recovered the answer to arbitrary temperature $\beta^{-1}$. One can then easily see that they satisfy the Schwinger-Dyson equation (4.5) at finite temperature (with $i\omega_{\tau}$ replaced by $i\omega_{n}$ therein). Therefore, we have proved that Eq. (4.30) is the correct finite temperature two point function in the large $N$ limit. One can also check that in the limit $\beta\rightarrow\infty$, Eqs. (4.36, 4.37) reduce to the zero temperature solution (4.26, 4.27) in momentum space.

4.2 Thermal quantities and the ground state entropy density

The thermal quantities of the large $N$ chiral SYK model can be derived from the exact finite temperature solution we obtained in Sec. 4.1.3. In fact, since the large $N$ two point function $G_{\beta}$ in Eq. (4.30) is of exactly the same form as that of $N=4$ in Eq. (3.15), the energy density $\mathcal{E}$ and the energy current $j_{\mathcal{E}}$ in the large $N$ limit can be calculated in exactly the same way as we did for $N=4$ in Eqs. (3.17) and (3.18). This gives an energy density

[TABLE]

and an energy current

[TABLE]

where we have eliminated the unphysical $1/\epsilon^{2}$ terms. This gives a thermal Hall conductance $\kappa_{xy}=\partial j_{\mathcal{E}}/\partial_{\beta}=N\pi/12\beta$ if the system is an edge of a $2+1$ dimensional bulk gapped topological state. In particular, $\kappa_{xy}$ does not depend on $J$ and is quantized at order $\mathcal{O}(N)$, which agrees with the general belief that $\kappa_{xy}$ is a topologically invariant quantity. The entropy density $S$ at temperature $\beta$ can be calculated from $\beta^{-1}[\partial\mathcal{S}/\partial(\beta^{-1})]=\partial\mathcal{E}/\partial(\beta^{-1})$, which gives $\mathcal{S}=\mathcal{S}_{0}+(N\pi/24\beta)(u_{+}^{-1}+u_{-}^{-1})$, where $\mathcal{S}_{0}$ denotes the ground state entropy density. In the below, we discuss the ground state entropy density $\mathcal{S}_{0}$ of our $1+1$ dimensional chiral SYK model.

A remarkable feature of the $0+1$ dimensional SYK model is that it has a large ground state entropy of order $N$, which indicates an unusual ground state degeneracy of order $e^{N}$. This implies a duality between a $0+1$ dimensional SYK model on the boundary of a $1+1$ dimensional AdS spacetime and a black hole in the bulk, which has a large entropy proportional to the event horizon area. As a comparison, in the $1+1$ dimensional chiral SYK model here we can examine the ground state entropy density $\mathcal{S}_{0}$. However, unlike the $0+1$ dimensional SYK model where the Hilbert space dimension $2^{N/2}$ is finite, the model here as a field theory has an infinite dimensional Hilbert space per unit length, which therefore makes the ground state entropy density possibly divergent and ill-defined. Therefore, we need to impose a physical spatial UV cutoff $\epsilon_{x}$ to make the Hilbert space dimension per unit length finite.

The entropy density at temperature $\beta^{-1}$ is defined by $\mathcal{S}=L^{-1}(1-\beta\partial_{\beta})\log Z$, where $L\rightarrow\infty$ is the spatial size of the system, and $Z$ is the partition function. By taking the limit $\beta\rightarrow\infty$, we obtain the ground state entropy density $\mathcal{S}_{0}$. In the large $N$ limit, we can take a Gaussian average over the random interactions $J_{ijkl}$, and write the partition as 222Here we write the disorder average over $J_{ijkl}$ as the annealed disorder $Z=\langle Z\rangle_{J}$, however, for large $N$ it is equal to the quenched disorder $Z=e^{\langle\log Z\rangle_{J}}$ to the leading order [5]. This is because in the replica treatment $\langle Z^{n}\rangle_{J}$ is equal to $\langle Z\rangle_{J}^{n}$ up to order $1/N^{3}$, which leads to $\langle\log Z\rangle_{J}=\lim_{n\rightarrow 0}\frac{\langle Z^{n}\rangle_{J}-1}{n}=\lim_{n\rightarrow 0}\frac{\langle Z\rangle^{n}_{J}-1}{n}=\log\langle Z\rangle_{J}$.

[TABLE]

One can first integrate out $J_{ijkl}$. Then, by inserting $1=\int\mathcal{D}\bar{G}\mathcal{D}\bar{\Sigma}\exp[-\frac{1}{2}\int d\tau dx\int d\tau^{\prime}dx^{\prime}\bar{\Sigma}(N\bar{G}-\sum_{i}\psi_{i}\psi_{i})]$, and integrating out the Majorana fermion fields $\psi_{i}$, one can rewrite the action as a functional of two functions $\bar{G}(-i\tau,x;-i\tau^{\prime},x^{\prime})$ and $\bar{\Sigma}(-i\tau,x;-i\tau^{\prime},x^{\prime})$ in the below:

[TABLE]

One can easily check that the Euler-Lagrange equation for $\bar{G}$ and $\bar{\Sigma}$ is exactly the Schwinger-Dyson equation satisfied by the two point function $G$ and the self energy $\Sigma$. Therefore, the on-shell values of the functions $\bar{G}$ and $\bar{\Sigma}$ are simply equal to the two point function $G$ and the self energy $\Sigma$ we solved in Sec. 4.1.3. To the leading order of $1/N$, the partition function $Z$ is just given by fixing $\bar{G}$ and $\bar{\Sigma}$ to their the on-shell values $G$ and $\Sigma$, namely,

[TABLE]

Instead of calculating the above complicated expression directly, we can compare its value with the free case $J=0$. This can be done by evaluating the derivative $\partial_{J}(\log Z/L)$ with respect to $J$. Making use of the Schwinger-Dyson equations $\partial_{\tau}-i\partial_{x}-\Sigma=G^{-1}$ (which should be understood as matrices in the Hilbert space) and $\Sigma(-i\tau,x)=J^{2}G(-i\tau,x)^{3}$, one would find that all the terms proportional to $\partial_{J}G$ and $\partial_{J}\Sigma$ vanish, and one arrives at

[TABLE]

In order for the above integral to be convergent, we need to impose a UV cutoff $(x/\epsilon_{x})^{2}+(\tau/\epsilon_{\tau})^{2}\geq 1$ with $\epsilon_{\tau}/\epsilon_{x}\rightarrow 0$ as shown in Fig. 7(b), where $\epsilon_{x}>0$ is now a physical spatial cutoff for the Hilbert space dimension per unit length not to diverge. Importantly, note that $\sin^{2}[{\pi\over\beta}(\tau-iu_{+}^{-1}x)]\sin^{2}[{\pi\over\beta}(\tau-iu_{-}^{-1}x)]$ is invariant under $\beta\rightarrow-\beta$, therefore the above expression is an odd function of $\beta$. As a result, we would find the following to be also an odd function of $\beta$:

[TABLE]

where $a_{j}$ ($j\geq 0$) are constants depending on $\epsilon_{x},\epsilon_{\tau}$ and $J$. The coefficient $a_{j}$ is of order $\mathcal{O}(\epsilon_{\tau,x}^{2j-2})$, and in particular, $a_{1}(J)$ does not depend on $\epsilon_{\tau}$ and $\epsilon_{x}$. In the free case $J=0$, the entropy density is easily known to be given by $\mathcal{S}_{f}=\mathcal{S}_{f0}+N\pi/(12\beta)$, where $\mathcal{S}_{f0}$ denotes the ground state entropy density of $N$ flavors of free chiral Majorana fermions. Therefore, we find the entropy density at interaction $J$ to be given by

[TABLE]

It can be directly calculated from Eq. (4.43) that $\partial_{J}a_{1}(J)=\frac{N}{96}(u_{-}^{-2}-u_{+}^{-2})$ 333since $\partial_{J}a_{1}(J)$ is $\epsilon_{x}$ and $\epsilon_{\tau}$ independent, it can be calculated using UV cutoff $|\tau|\geq\epsilon$ as illustrated by Fig. 7(b)-(c), instead of the more complicated cutoff $(x/\epsilon_{x})^{2}+(\tau/\epsilon_{\tau})^{2}\geq 1$ with $\epsilon_{\tau}/\epsilon_{x}\rightarrow 0$.. Therefore, the coefficient of the thermal part of entropy density is $\frac{N\pi}{12}+2a_{1}(J)=\frac{N\pi}{24}(u_{+}^{-1}+u_{-}^{-1})$, which is in agreement with our calculations below Eq. (4.39). The ground state entropy density $\mathcal{S}_{0}=\lim_{\beta\rightarrow\infty}\mathcal{S}$ for any $0\leq J<2\pi$ is therefore equal to that in the free fermion case, namely, $\mathcal{S}_{0}=\mathcal{S}_{f0}$. It is known that $N$ flavors of free chiral Majorana fermions has a unique ground state for anti-periodic spatial boundary condition, and has $2^{N}$ degenerate ground states for periodic spatial boundary condition due to $N$ zero modes. This gives a total ground state entropy either $\mathcal{S}_{f0}L=0$ or $\mathcal{S}_{f0}L=N\log 2$. Therefore, in the infinite spatial size limit $L\rightarrow\infty$, the ground state entropy density of the free fermions is $\mathcal{S}_{f0}\rightarrow 0$. We then conclude that the ground state entropy density for any interaction strength $0\leq J<2\pi$ is $\mathcal{S}_{0}=\mathcal{S}_{f0}=0$. Basically, in order to have $\mathcal{S}_{0}$ different from $\mathcal{S}_{f0}$, one has to have a constant $1/\epsilon_{x}$ piece independent of $\beta$ in $L^{-1}\log Z$, which is absent in Eq. (4.44) here.

The fact that the ground state entropy density $\mathcal{S}_{0}=\mathcal{S}_{f0}=0$ for any $0\leq J<2\pi$ can also be physically understood from the spectral weight $A(\omega,k)$ in the large $N$ limit plotted in Fig. 5. Since the spectral weight $A(\omega,k)$ is nonzero only when $u_{-}\leq\omega/k\leq u_{+}$, this suggests that any many-body eigenstate $|n,k\rangle$ with momentum $k$ (we use $n$ to denote different states) has an energy bounded by $E_{n,k}\in[u_{-}k,u_{+}k]$, and all the eigenstates have $k\geq 0$ (with both the energy and momentum of the ground state energy defined as zero). This means all the many-body state energy eigenvalues $E_{n,k}$ at momentum $k$ of the chiral SYK system are no smaller than that of $N$ flavors of free chiral Majorana fermions with velocity $u_{-}$, and no larger than that of $N$ flavors of free chiral Majorana fermions with velocity $u_{+}$. Accordingly, the partition function $Z$ of the chiral SYK model at any temperature $\beta^{-1}$ will also be bounded by $Z_{+}\leq Z\leq Z_{-}$, where $Z_{\pm}$ is the partition function of $N$ flavors of free chiral Majorana fermions with velocity $u_{\pm}$, respectively. Since the ground state entropy density of free chiral Majorana fermions $\mathcal{S}_{f0}=0$ is independent of their velocity, the above bound of partition function $Z$ indicates that the ground state entropy density of the chiral SYK model has to be $\mathcal{S}_{0}=\mathcal{S}_{f0}=0$ for any $0\leq J<2\pi$. From this argument, we can see that the chiral nature of the model (and translational invariance) protects its ground state entropy density to be at zero.

4.3 The OTOC, chaos exponent and butterfly cone

In this section, we derive the chaos exponent and butterfly velocity in the OTOC of the model. Following [4], we define the regularized OTOC four point function in real time as

[TABLE]

where $y=e^{-\beta H/4}=\rho(\beta)^{1/4}$ separates evenly the four fermion fields by a quarter of the thermal circle. The leading contribution in the early time OTOC comes from the contractions between two $\psi_{i}$ and between two $\psi_{j}$, which gives an order $1$ piece $-G(-i\frac{\beta}{2},0)G(t_{1}-t_{2}-i\frac{\beta}{2},x_{1}-x_{2})$. The next order contribution comes from contraction of $\psi_{i}$ with $\psi_{j}$, and is of order $1/N$ at the early time. Therefore, we can separate the order $1$ and order $1/N$ pieces of the OTOC as

[TABLE]

where the function $\delta\mathcal{F}$ is then of order $1$ at the early time by definition.

To the leading order of $1/N$, the Feynman diagrams contributing to $\delta\mathcal{F}$ are the ladder diagrams as shown in Fig. 9 [4, 5]. For chaotic systems where there is a Lyapunov regime (which usually requires unbounded local Hilbert spaces), when $t_{1}=t_{2}=t$, one expects $\delta\mathcal{F}$ to grow exponentially with respect to time $t$ in the Lyapunov regime. For large $N$ models like ours, the Lyapunov regime is in the interval $\beta\lesssim t\lesssim\beta\log N$. In this case, at late time $t$ when $\delta\mathcal{F}$ becomes large, it will approximately satisfy the self-consistent equation

[TABLE]

where $K_{R}$ is the retarded kernel defined as

[TABLE]

as illustrated by the second line of Fig. 9. Here we have denoted for short $t_{ij}=t_{i}-t_{j}$, and $x_{ij}=x_{i}-x_{j}$. The functions $G_{R}$ and $G_{lr}$ are the retarded Green’s function and the Wightman correlator with half thermal circle separation, respectively, which can be explicitly derived via the analytic continuation of Eq. (4.30) as

[TABLE]

and

[TABLE]

Therefore, the function $\delta\mathcal{F}$ is an eigenfunction of kernel $K_{R}$ with eigenvalue $1$. To derive the exponential growth of $\delta\mathcal{F}$, we first need to compute all the eigenfunctions of kernel $K_{R}$ which have eigenvalue $1$.

To simplify the computation of the eigenfunctions, we make a spacetime coordinate transformation from $(t,x)$ to a new basis

[TABLE]

where we have used the relation $u_{\pm}=1\pm J/2\pi$. The transformation is nonsingular as long as the interaction strength $J>0$. After this change of coordinates, Eq. (4.48) can be rewritten as

[TABLE]

One can see the equation does not have an explicit dependence on the interaction strength $J$ any more. In addition, the kernel $K_{R}$ factorizes into the product of a function $k_{R}^{+}(t_{i}^{+})$ of coordinates $t_{i}^{+}$ multiplied by a function $k_{R}^{-}(-t_{i}^{-})$ of $-t_{i}^{-}$ ($1\leq i\leq 4$), each of which up to a factor is nothing but the kernel of a $1$d $q=4$ SYK model along the $t^{+}$ or the $-t^{-}$ time direction in the conformal limit. Therefore, an eigenfunction of the $1+1$ dimensional kernel $K_{R}$ here can also be factorized into the product of the eigenfunctions of two $0+1$ dimensional SYK model along $t^{+}$ and $-t^{-}$ directions, respectively. More explicitly, the eigenfunction of the kernel $K_{R}$ can be written into the form

[TABLE]

where $h_{+}$ and $h_{-}$ are constants to be determined. The corresponding eigenvalues of $K_{R}$ are

[TABLE]

The above eigenfunction (4.54) and eigenvalue (4.55) can be directly verified by substituting them into Eq. (4.53). Since there is no UV divergence in the above calculations, we do not have the point-splitting issue that arises in the calculation of the self energy $\Sigma$. This is consistent with our expectation that the chaos is controlled by the IR physics.

When the spacetime coordinates of the two $\psi_{j}$ fields in Eq. (4.46) coincide, namely, $t_{1}=t_{2}=t$ and $x_{1}=x_{2}=x$, the eigenfunction can be rewritten as

[TABLE]

Since the four point function $\delta\mathcal{F}$ at any fixed time $t$ cannot diverge when $x\rightarrow\pm\infty$, we need to constrain $p$ to be real, and $f_{h_{+},h_{-}}$ is then simply a plane wave solution with momentum $2\pi p/\beta$. By solving Eq. (4.55), we can derive $\varkappa$ as a function of $p$ to be

[TABLE]

The function is defined in the branch where $\chi(p)\rightarrow-i(1-\mathcal{J})p-\frac{\mathcal{J}}{1+\mathcal{J}}$ as $p\rightarrow\infty$. Note that $\varkappa(p)$ in general takes complex values, so the Lyapunov exponent of the eigenfunction $f_{h_{+}(p),h_{-}(p)}$ is given by $\lambda(p)=\text{Re}\left[{2\pi\over\beta}\varkappa(p)\right]$ (where $p$ is real), which is an even function of $p$ since $\varkappa(p)=\varkappa(-p)^{*}$. Its maximum value $\overline{\lambda}$ is always achieved at $p=0$, which has a magnitude

[TABLE]

As $J$ grows from [math] to $2\pi$ (i.e., $\mathcal{J}$ grows from [math] to $1$), the maximal Lyapunov exponent $\overline{\lambda}$ grows from [math] to the maximal chaos bound $2\pi/\beta$, as shown in Fig. 10(a).

We now discuss the chaotic behavior and butterfly cone of the function $\delta\mathcal{F}(t,x)=\delta\mathcal{F}(t,x;t,x)$ with $t_{1}=t_{2}=t>0$ and $x_{1}=x_{2}=x$ set in Eq. (4.46), namely, the $1/N$ piece of the OTOC of our chiral SYK model. For this purpose, we need to know the weight coefficient $\rho(p)$ of the eigenfunction $f_{h_{+}(p),h_{-}(p)}$ in $\delta\mathcal{F}$. Due to the translation symmetry, one can check that each total momentum $2\pi p/\beta$ component of the four point function $\delta\mathcal{F}(p,x_{12},t_{1},t_{2})=\int dxe^{i\frac{2\pi}{\beta}px}\delta\mathcal{F}(t_{1},x_{1}+x;t_{2},x_{2}+x)$ satisfies the ladder diagram of Fig. 9 by itself. This allows us to employ the remarkable identity between magnitude and the chaos exponent derived in [1] to determine the coefficient $\rho(p)$. Up to multiplication of some mild function, the weight coefficient is dominated by the factor

[TABLE]

Therefore, to a good approximation the OTOC function $\delta\mathcal{F}$ is given by the integral

[TABLE]

As discussed in details in [1], the integral at large $x$ and $t$ can be done by using saddle point approximation for $p$. The saddle point $p_{v}$ of the integrand of Eq. (4.60) in the complex plane of $p$ is given by

[TABLE]

which depends on the velocity $v=x/t$ we look at. Since $\varkappa(p)$ is real on the imaginary $p$ axis, $\varkappa^{\prime}(p)$ is purely imaginary on the imaginary $p$ axis. Besides, $\varkappa(p)$ has two Riemann surface branches connected by a branch cut. For velocity $u_{-}\leq v\leq u_{+}$, this yields two saddle points $p_{v}$ and $\tilde{p}_{v}$ on the imaginary axes of the two Riemann surface branches, as shown in Fig. 10(c). The original integration contour along the real axis in Eq. (4.60) can be deformed into two steepest descent contours passing the two saddle points $p_{v}$ and $\widetilde{p}_{v}$, respectively, as shown in Fig. 10(c). If the contour deformation does not cross any poles, the integral of Eq. (4.60) will be dominated by the contribution of the first saddle point $p_{v}$ as (Appendix B)

[TABLE]

where the velocity $v$ dependent Lyapunov exponent $\lambda_{v}=2\pi[\varkappa(p_{v})+ip_{v}v]/\beta$ is real. The maximal value of $\lambda_{v}$ in Eq. (4.62) is equal to $\overline{\lambda}$ given by Eq. (4.58), which is reached at velocity

[TABLE]

namely, when the saddle point $p_{v}=0$. For velocity $v$ in the vicinity of $v_{c}$, the $v$ dependent Lyapunov exponent is approximately given by $\lambda_{v}\approx\overline{\lambda}-\frac{2\pi\xi}{\beta}(v-v_{c})^{2}$, and thus

[TABLE]

where $\xi=\frac{(3-2\mathcal{J}^{2})^{3/2}}{6\mathcal{J}(1-\mathcal{J}^{2})^{2}}$. Besides, $\lambda_{v}$ drops to zero at a velocity $v_{s}$ satisfying $v_{s}=i\varkappa^{\prime}(p_{s})=i{\varkappa(p_{s})/p_{s}}$, where $p_{s}$ is the saddle point position for velocity $v_{s}$. Solving this equation gives us two velocities $v_{s}^{\pm}$ satisfying $\lambda_{v_{s}^{\pm}}=0$, which have the expressions

[TABLE]

When the velocity $v<v_{s}^{-}$ or $v>v_{s}^{+}$, we have $\lambda_{v}<0$, so the saddle point contribution (4.62) no longer gives an exponential growth in $t$, and is not valid for large $t$. In particular, we note that for any $0<J<2\pi$ we have $v_{s}^{-}>u_{-}$ and $v_{s}^{+}<u_{+}$ (Fig. 10(e)), so $v_{s}^{\pm}$ are within the causality cone of the system between $u_{-}$ and $u_{+}$ (where the retarded Green’s function $G_{R}$ is nonzero).

If the contour deformation for Eq. (4.60) crosses a pole of the integrand, the integral will pick up both the pole contribution and the saddle point contribution. The weight coefficient $\rho(p)$ has two poles $p_{1}$ and $p_{2}$ satisfying $\varkappa(p_{1})=\varkappa(p_{2})=1$ on the positive imaginary axes of the two Riemann surface branches:

[TABLE]

and two poles $\widetilde{p}_{1}$ and $\widetilde{p}_{2}$ satisfying $\varkappa(\widetilde{p}_{1})=\varkappa(\widetilde{p}_{2})=-1$ on the imaginary axes of the two branches

[TABLE]

The poles $p_{2}$ and $\widetilde{p}_{2}$ can be shown to never give a dominant contribution (Appendix B).

As the velocity $v$ increases to $v>v_{*}=i\varkappa^{\prime}(p_{1})={2-2\mathcal{J}^{2}\over 2-\mathcal{J}}$, the saddle point $p_{v}$ will move upwards and pass the pole $p_{1}$, and (i.e., $v>v_{*}$) the integral (4.60) will pick up and be dominated by the residue of the pole $p_{1}$. In particular, one always has $v_{c}<v_{*}<v_{s}^{+}$ for $0<J<2\pi$ (Fig. 10(f)), so this transition from saddle to pole contribution at $v=v^{*}$ happens before $\lambda_{v}$ given by the saddle point contribution decreases to zero. Therefore, for $v=x/t>v_{*}$, the OTOC function has the form

[TABLE]

and the velocity dependent Lyapunov exponent $\lambda_{v}$ will be given by the pole contribution as $\lambda_{v}=2\pi(1-vu_{+}^{-1})/\beta$. In particular, one sees that $\lambda_{v}$ in this spacetime region decreases to zero exactly at velocity $v=u_{+}$, i.e., at the right causality edge. For velocity $v>u_{+}$, the growth behavior of $\delta\mathcal{F}$ disappears.

As the velocity $v$ decreases to $v<v_{*}^{\prime}=i\varkappa^{\prime}(\widetilde{p}_{1})={2-2\mathcal{J}^{2}\over 2+\mathcal{J}}$, the saddle point $p_{v}$ will move downwards and pass the pole $\widetilde{p}_{1}$. Accordingly, the integral (4.60) will be dominated by the residue of the pole $\widetilde{p}_{1}$. Similarly, one has $v_{s}^{-}<v_{*}^{\prime}<v_{c}$ for $0<J<2\pi$ (Fig. 10(f)), so this transition from saddle to pole contribution at $v=v_{*}^{\prime}$ also happens before $\lambda_{v}$ decreases to zero. So the OTOC for $v=x/t<v_{*}^{\prime}$ is given by

[TABLE]

and the velocity dependent Lyapunov exponent is $\lambda_{v}=2\pi(vu_{-}^{-1}-1)/\beta$. Therefore, one sees $\lambda_{v}$ also decreases to zero at velocity $v=u_{-}$, i.e., the left causality edge. For velocity $v<u_{-}$, the growth of $\delta\mathcal{F}$ disappears.

The numerical accuracy of the above saddle point approximation is demonstrated in Appendix B. Fig. 10(f)-(g) shows the velocity dependent Lyapunov exponent $\lambda_{v}$ as a function of velocity $v=x/t$ for two cases $J=1.2\pi$ and $J=1.9\pi$, respectively. We can therefore define the butterfly cone of our model as $u_{-}<x/t<u_{+}$ as shown by the shaded regions in Fig. 10(d), inside (outside) which $\lambda_{v}>0$ ($\lambda_{v}<0$). Roughly speaking, the butterfly cone characterizes the range of spreading of local chiral operators as a function of time [60]. It coincides exactly with the causality cone. We can summarize the behavior of the OTOC function $\delta\mathcal{F}$ as follows:

1. $v_{-}<{x/t}<v_{*}^{\prime}$, in which case the pole $\widetilde{p}_{1}$ contribution dominates, and $\delta\mathcal{F}$ grows for $\beta\lesssim t\lesssim\beta\log N$ according to Eq. (4.69). This is illustrated by the left dashed shaded region in Fig. 10(d). The velocity dependent Lyapunov exponent in this region is given by $\lambda_{v}=2\pi(vu_{-}^{-1}-1)/\beta$. Note that the $t$ direction Lyapunov exponent in this region is $-2\pi/\beta<0$, although the velocity dependent Lyapunov exponent is positive.

2. $v_{*}^{\prime}<{x/t}<v_{*}$, where the saddle point $p_{v}$ contribution dominates, and $\delta\mathcal{F}$ grows exponentially in $t$ as given by Eq. (4.62), which is valid for time $\beta\lesssim t\lesssim\beta\log N$. This butterfly region is illustrated by the solid shaded region in Fig. 10(d). The velocity dependent Lyapunov exponent $\lambda_{v}$ reaches its maximal value $\overline{\lambda}$ at velocity $v_{c}$ in this range (Eq. (4.63)). Fig. 10(b) illustrates the OTOC $\mathcal{F}(t,v_{c}t)$ along the $x=v_{c}t$ direction, which is proportional to $1-e^{\overline{\lambda}(t-t_{sc})}$, where $t_{sc}\sim\mathcal{O}(\beta\log N)$ is the scrambling time.

3. $v_{*}<{x/t}<u_{+}$, where the pole $p_{1}$ contribution dominates, and the exponential growth of $\delta\mathcal{F}$ for $\beta\lesssim t\lesssim\beta\log N$ is given by Eq. (4.68). This butterfly region is illustrated by the right dashed shaded region in Fig. 10(d). The velocity dependent Lyapunov exponent in this region is given by $\lambda_{v}=2\pi(1-vu_{+}^{-1})/\beta$. In addition, it is also worthwhile to note that for fixed $x>0$ in this spacetime region, $\delta\mathcal{F}$ grows in $t$ with an exponent $2\pi/\beta$, saturating the maximal chaos bound.

4. $x/t>u_{+}$ or $x/t<v_{s}^{-}$. In this region, $\delta\mathcal{F}$ no longer shows a growing behavior in $t$. Accordingly, the normalized OTOC $\mathcal{F}\sim 1-\frac{1}{N}\delta\mathcal{F}$ will tend to $1$ at large time $t$.

As one can see from the constant time $t$ slices of the (normalized) OTOC $\mathcal{F}(t,x)$ shown in Fig. 10(h) and (i) for $J=1.2\pi$ and $J=1.9\pi$, respectively, the exponential growth is peaked at velocity $v_{c}$ inside butterfly cone represented by the green shaded area. The peak becomes sharper and tends to $u_{-}$ as $J$ increases. In some sense, the behavior of the chiral SYK model looks like a nonchiral chaotic model with its edges of butterfly cone tilted to velocities $u_{\pm}$, and its time direction tilted to velocity $v_{c}$.

5 The model at other small finite $N$

We have seen in Sec. 3 that the $1+1$ dimensional chiral SYK model (2.1) is integrable for $N=4$, and in Sec. 4 that the model is chaotic in the large $N$ limit. Therefore, we conjecture there exists a critical number of Majorana flavors $N_{c}$, below which ($N<N_{c}$) the model is integrable, while equal to or above which ($N\geq N_{c}$) the model is non-integrable and chaotic. We note that here this is a well-defined statement, since the Hilbert space dimension of the $1+1$ dimensional model (2.1) is infinite for any $N$, and the energy spectrum cannot be exactly solved (i.e., integrable) unless the model has infinite number of conserved quantities defined in terms of local operators. In contrast, the $0+1$ dimensional SYK model with $N$ Majorana fermions has a finite Hilbert space at any finite $N$, therefore it can always be exactly solved by diagonalizing a finite dimensional Hamiltonian matrix. In this sense, the $0+1$ dimensional SYK model is “integrable” at any finite $N$, and there is no clear boundary between integrable and non-integrable during the increase of $N$.

In this section, we show that the $1+1$ dimensional chiral SYK model is also integrable by bosonization for $N=5$ and $N=6$. However, the interactions of the model becomes increasingly complicated for $N\geq 7$. This leads us to conjecture that the boundary between integrable and chaotic is $N_{c}=7$.

The difficulty in solving the chiral SYK model (2.1) for a finite $N>4$ comes from the fact that there are large number of interaction terms $J_{ijkl}\psi_{i}\psi_{j}\psi_{k}\psi_{l}$ which do not mutually commute. However, we can always rotate the chiral Majorana basis $\psi_{i}$ to a new basis $\psi_{i}^{\prime}=O_{ij}\psi_{j}$, where $O$ is an SO(N) matrix, under which $J_{ijkl}$ transforms as a rank-4 tensor. In some cases, one can find a proper basis to reduce the number of interaction terms. This is the basic idea how we can solve the model for $N=5$ and $N=6$.

5.1 Exact solution for $N=5$

For the chiral model with $N=5$, there are $5$ independent interactions $J_{ijkl}$, which forms a total antisymmetric rank-4 tensor. We can redefine them using their Hodge dual as $I_{i}=\frac{1}{4!}\epsilon_{ijklm}J_{jklm}$ (repeated indices are summed), where $\epsilon_{ijklm}$ is the Levi-Civita symbol. $I_{i}$ is nothing but a rewriting of the $5$ interactions $J_{ijkl}$. The Lagrangian density of the model then becomes

[TABLE]

Under an SO(5) rotation $\psi_{i}^{\prime}=O_{ij}\psi_{j}$, the interactions $I_{i}$ transform as an SO(5) vector. Therefore, we can find a proper SO(5) matrix $O_{ij}$ to rotate the vector $I_{i}$ to $I_{i}^{\prime}=O_{ij}I_{j}=(0,0,0,0,I)$, where $I=\sqrt{\sum_{i=1}^{5}I_{i}^{2}}$. Under this new basis, the Lagrangian becomes

[TABLE]

One easily sees that $\psi_{5}^{\prime}$ is a free chiral Majorana fermion decoupled from the other four fermions $\psi_{1\leq i\leq 4}^{\prime}$. The Lagrangian for $\psi_{1\leq i\leq 4}^{\prime}$ simply reduces to the $N=4$ model (with interaction strength $I$) we discussed in Sec. 3, which is integrable by bosonization. Therefore, we see the $N=5$ model is integrable, which is equivalent to two free bosons with velocities $u_{\pm}$ and a free chiral Majorana fermion with velocity $1$. We can then easily obtain the averaged two point function for $N=5$ at temperature $\beta^{-1}$ as

[TABLE]

where $u_{\pm}=1\pm\frac{I}{2\pi}$. Any higher $n$ point function can be calculated using the Wick’s theorem. From the two point function $G(-i\tau,x)$ or from the bosonic picture, we can derive the finite temperature energy density to be $\mathcal{E}=\frac{\pi}{12\beta^{2}}(u_{+}^{-1}+u_{-}^{-1}+\frac{1}{2})$, and the energy current $j_{\mathcal{E}}=5\pi/24\beta^{2}$. This yields a quantized thermal Hall conductance $\kappa_{xy}=5\pi/12\beta$, and a thermal entropy density $\mathcal{S}=\frac{\pi}{6\beta}(u_{+}^{-1}+u_{-}^{-1}+\frac{1}{2})$.

5.2 Exact solution for $N=6$

Now we turn to the $N=6$ case. Following a similar idea, we can redefine the interaction rank-4 tensor $J_{ijkl}$ using its Hodge dual $I_{ij}=\frac{1}{4!}\epsilon_{ijklmn}J_{klmn}$, which is antisymmetric in $i$ and $j$. Under an SO(6) rotation $\psi_{i}^{\prime}=O_{ij}\psi_{j}$, the interactions $I_{ij}$ transform as an antisymmetric rank-2 tensor. Therefore, there exists an SO(6) transformation $O_{ij}$ which brings $I_{ij}$ into the standard $2\times 2$ block diagonal form

[TABLE]

After changing into this new basis of $\psi_{i}^{\prime}$, the Lagrangian of the model at $N=6$ becomes

[TABLE]

We can then bosonize the chiral Majorana fermions as $\psi_{1}^{\prime}+i\psi_{2}^{\prime}=e^{i\phi_{1}}$, $\psi_{3}^{\prime}+i\psi_{4}^{\prime}=e^{i\phi_{2}}$ and $\psi_{5}^{\prime}+i\psi_{6}^{\prime}=e^{i\phi_{3}}$, respectively, after which the Lagrangian becomes

[TABLE]

where the velocity matrix is a real symmetric matrix given by

[TABLE]

Note that we have only quadratic terms of $\phi_{i}$ in the above Lagrangian. The velocity matrix can then be diagonalized into $V_{ij}^{\prime}=\text{diag}(u_{1},u_{2},u_{3})$ by an SO(3) transformation $\phi_{i}^{\prime}=Q_{ij}\phi_{j}$, where $Q_{ij}$ is a $3\times 3$ special orthogonal matrix. After the transformation, one gets $3$ decoupled free boson fields $\phi_{i}^{\prime}$ with velocities $u_{i}$ ($i=1,2,3$), respectively, and the Lagrangian becomes

[TABLE]

Therefore, the model is also integrable via bosonization for $N=6$. One can then show the averaged fermion two point function for $N=6$ at temperature $\beta^{-1}$ to be given by

[TABLE]

Any higher $n$ point function can be calculated using the Wick’s theorem. It is also straightforward to derive the finite temperature energy density $\mathcal{E}=\frac{\pi}{12\beta^{2}}(u_{1}^{-1}+u_{2}^{-1}+u_{3}^{-1})$, and the energy current $j_{\mathcal{E}}=\pi/4\beta^{2}$. This yields a quantized thermal Hall conductance $\kappa_{xy}=\pi/2\beta$, and a thermal entropy density $\mathcal{S}=\frac{\pi}{6\beta}(u_{1}^{-1}+u_{2}^{-1}+u_{3}^{-1})$.

5.3 $N\geq 7$ and thermalization of chiral Majorana fermions

The interaction terms $J_{ijkl}$ of the chiral SYK model become increasingly complicated when $N\geq 7$. For instance, for $N=7$, the rank-4 total antisymmetric tensor $J_{ijkl}$ has $35$ parameters in total. One can try to make some of them zero by an SO($7$) rotation $\psi_{i}^{\prime}=O_{ij}\psi_{j}$. However, since the SO(7) group has only $21$ generators, the SO(7) rotation matrix $O_{ij}$ can at most reduce $21$ or less of the interaction parameters $J_{ijkl}$ to zero, while there are still $14$ or more nonzero parameters $J_{ijkl}$ remaining. After bosonization, these remaining $14$ or more interaction terms cannot all be bilinears of boson fields $\phi_{i}$ of the form $\partial_{x}\phi_{i}\partial_{x}\phi_{j}$ but necessarily contain many nonlinear terms of the form $\cos(\phi_{i}\pm\phi_{j}\pm\phi_{k}\pm\phi_{l})$. Therefore, the model is no longer free bosons for $N=7$, which makes the model extremely complicated. This is also the case for all $N>7$. Therefore, we conjecture the model is no longer integrable for $N\geq N_{c}$ and starts to show chaotic behaviors in some aspects, where the critical number $N_{c}\geq 7$, and is probably equal to $7$.

One of the physical differences between $N<N_{c}$ and $N\geq N_{c}$ will be the thermalization of the model, which is closely related to integrability. For instance, consider the spatial configuration of a set of chiral Majorana edge states as shown in Fig. 11, where the $N$-th mode $\psi_{N}$ is initially spatially separated from the other $N-1$ modes $\psi_{j\leq N-1}$, and merges with them after propagating a distance. Such a configuration can be realized by using the edges of $N$ copies of the $p+ip$ chiral topological superconductor. Assume $\psi_{j\leq N-1}$ and $\psi_{N}$ are connected with thermal baths of temperature $T_{1}$ and $T_{2}$, respectively. For $N<N_{c}$, the model is integrable, so one expects the $N$ chiral Majorana fermion modes do not thermalize with each other during propagation. Instead, the model with $N\geq N_{c}$ is non-integrable, and will have all the chiral Majorana fermion modes thermalize among themselves after they merge together.

The thermalization of the chiral edge states of a topological condensed matter system may affect the experimental measurement of its thermal Hall conductance. For instance, the recent experiment observed a thermal Hall conductance of the filling $\nu=5/2$ fractional quantum Hall state disagreeing with theoretical predictions [48, 61, 62, 63], and one of the possible explanations is suggested to be the lack of thermal equilibrium among the chiral edge modes [49, 50, 51]. Our model here might be too idealized to describe a realistic edge of condensed matter system, which always have spatial disorders and coupling to phonons, etc. However, our results still suggest that thermalization is more difficult for fewer number of chiral edge modes.

6 In closing

We have introduced a new chiral model in $1+1$ dimensions which is solvable in the large $N$ limit like its $0+1$ dimensional namesake and tensor models. Indeed, the exact scaling symmetry of our $1+1$ dimensional chiral SYK model makes it exactly solvable at all energy scales, thus makes it a “simpler” chaotic model than the the $0+1$ dimensional SYK model which is only exactly solvable at low energies. We have calculated the two point function of the model, and verified that the thermal Hall conductance of the chiral model at nonzero interaction strength $J>0$ in the large $N$ limit is quantized at $\kappa_{xy}=N\pi/12\beta$, which is generically believed to be a topological invariant. We also find the ground state entropy at nonzero interaction strength $J$ is no different from that of $N$ free chiral Majorana fermions, which leads to a vanishing ground state entropy density.

The OTOC in the large $N$ limit generically has a Lyapunov regime and an asymmetric butterfly cone. The velocity dependent Lyapunov exponent $\lambda_{v}$ along $x=vt$ has its maximum reaching the chaos bound $2\pi/\beta$ when the interaction strength $J$ approaches the upper bound $2\pi$. Therefore, the model can approach the maximally chaotic limit without having a nonzero ground state entropy density, which is not expected from a black hole picture. Besides, for any interaction strength $J>0$, the $t$ direction Lyapunov exponent $\lambda$ for a fixed $x$ near the right butterfly edge always saturates the chaos bound $2\pi/\beta$, while near the left butterfly edge the $t$ direction Lyapunov exponent $\lambda$ is always $-2\pi/\beta$ (although the velocity dependent Lyapunov exponent is positive). This feature near the left butterfly edge is different from that in nonchiral models. Furthermore, the model is integrable for small $N\leq 6$ as can be seen by bosonization, but becomes increasingly complicated for $N\geq 7$. Therefore, we conjecture there is a transition from integrablity to chaos as $N$ increases to some number $N_{c}$, and the transition is probably at $N_{c}=7$.

Somewhat surprisingly we find that the two point function of the chiral SYK model in the large $N$ limit is exactly the same as that for $N=4$, although the four point functions of the two cases are completely different. We conjecture such a coincidence of two-point functions is because of the presence of the SO($N$) symmetry in both the large $N$ case (effectively) and the $N=4$ case.

It will also be interesting to investigate whether this model has Schwarzian like soft modes and an effective Schwarzian action as in [17]. If so, we could ask whether there is a holographic interpretation. Besides, in this paper we only calculate the OTOC in the late time chaos regime. To get more information about the operators in the model, one can study the OPE expansion of the four point function. Furthermore, we have been restricting ourselves to interaction strengths within the range $0\leq J<2\pi$ in this paper. One may ask what happens for large $N$ if $J>2\pi$, with proper UV nonlinearities included to preserve the total chirality (as shown in Fig. 1(c)). In this case, the model is no longer a purely chiral model. We shall leave these questions to future studies.

We briefly discuss how spatial disorder may affect the model, which always exists in practical condensed matter systems. The most relevant disordered term one can add to the Hamiltonian (2.5) is a quadratic fermion term

[TABLE]

where $A_{ij}(x)$ is an $x$ dependent real antisymmetric matrix. However, such a term can be eliminated by a spatially dependent SO(N) rotation $\psi_{i}\rightarrow(e^{-\int^{x}_{-\infty}A(x^{\prime})dx^{\prime}})_{ij}\psi_{j}$. In addition, the 4-fermion coupling $J_{ijkl}$ may also has an $x$ dependence due to disorder or the above spatially dependent SO(N) rotation. However, as long as the couplings $J_{ijkl}$ still have long range spatial correlations, namely, $\langle J_{ijkl}(x)J_{ijkl}(x^{\prime})\rangle\rightarrow 3!J^{2}/N^{3}$ as $|x-x^{\prime}|\rightarrow\infty$, we expect the physics in this paper remains unchanged in the IR.

Lastly, we comment that interacting chiral models with a large number of degrees of freedom (either Majorana fermions or complex fermions) arise naturally in two settings. The first is, naturally enough, that of “infinite layer” quantum Hall systems [64, 65, 66, 67] or $p+ip$ chiral superconductors. The second is the low energy theory of a piece of Fermi surface (of complex fermions with a conserved U(1) charge) in spatial dimension $d_{s}\geq 2$ with forward scatterings [68, 69]. In this case, the fermion dispersion $\omega=v_{F}k_{\perp}$ is linear in the momentum $k_{\perp}$ normal to the Fermi surface, where $v_{F}$ is the Fermi velocity. Therefore, one can treat the fermion at each fixed transverse momentum $\mathbf{k}_{\parallel}$ tangent to the Fermi surface as a chiral fermion in the $k_{\perp}$ direction, and regard $\mathbf{k}_{\parallel}$ as a flavor index which runs over large $N$ number of values. We expect our chiral SYK model to further shed lights on the general understanding of such chiral systems.

7 Acknowledgments

We thank Y. Gu, F.D.M. Haldane, I. Klebanov, J. Maldacena, X.L. Qi, N. Seiberg, D. Stanford, H. Verlinde, J.Y. Chen and X.Q. Sun for helpful discussions. We also thank G. Sárosi and M. Mezei for valuable comments on our OTOC calculations [70]. BL is supported by Princeton Center for Theoretical Sciences at Princeton University. SLS acknowledges support from the United States Department of Energy via grant No. DE-SC0016244. ZY is supported by Charlotte Elizabeth Procter Fellowship from Princeton University.

Appendix A Bosonization and point splitting

In this appendix, we show how to deal with the $x$ direction point splitting regularization of chiral Majorana fermion fields $\psi_{i}$ and $\psi_{j}$ via bosonization, and derive the commutation relation (2.4) and the energy momentum tensor $T^{\mu}_{\ \nu}$.

A.1 Bosonization and normal ordering

We shall follow the procedure of bosonization to regularize the point splitting. Given two chiral Majorana fermion fields $\psi_{i}$ and $\psi_{j}$ ($i\neq j$), we can define a complex fermion $c$ and its bosonization

[TABLE]

where $\phi$ is a bosonic field, and $:\mathrel{\mathcal{O}}:$ stands for the normal ordering of operator $\mathcal{O}$. Since all the point splitting we are considering is along the $x$ direction, we assume all the fields in the below are on the same time slice, i.e., at the same $t$. For simplicity, we shall therefore omit the time variable $t$ of all fields.

On the constant time slice, the chiral Majorana fermion fields satisfy the anticommutation relation $\{c(x),c^{\dagger}(x^{\prime})\}=\{c(x),c(x^{\prime})\}=\{c^{\dagger}(x),c^{\dagger}(x^{\prime})\}=0$ for $x\neq x^{\prime}$, which requires the boson commutation relation

[TABLE]

and thus

[TABLE]

To correctly deal with normal ordering, we can separate the boson field $\phi$ into

[TABLE]

where $\varphi(x)$ and $\varphi^{\dagger}(x)$ are the collection of annihilation and creation operators in the mode expansion, respectively. The normal ordering rule is then that $\varphi^{\dagger}(x)$ should always be on the left of $\varphi(x)$. They satisfy $[\varphi(x),\varphi(x^{\prime})]=[\varphi^{\dagger}(x),\varphi^{\dagger}(x^{\prime})]=0$. Using Eq. (A.3) and the fact that

[TABLE]

where $0^{+}$ represents positive infinitesimal, one finds that

[TABLE]

Integrating the two equations yields

[TABLE]

where the constant of integration is chosen for normalization reason. In particular, from Eq. (A.7) one can verify that $[\phi(x),\phi(x^{\prime})]=-\log\left(\frac{x-x^{\prime}+i0^{+}}{x^{\prime}-x+i0^{+}}\right)=i\pi\text{sgn}(x-x^{\prime})$, in agreement with Eq. (A.2).

A.2 Operator product expansion and commutation relation

We can now calculate the operator product expansion of chiral Majorana fermion fields, with a point splitting in the $x$ direction. First, we consider the bilinear $c^{\dagger}(x)c(x^{\prime})$ with small $x$ direction separation $x-x^{\prime}$:

[TABLE]

The first constant term $\frac{i}{2\pi(x-x^{\prime}+i0^{+})}$ is the vacuum term. Therefore, by taking the limit $x-x^{\prime}\rightarrow 0$, we find the normal ordered density operator is

[TABLE]

In particular, the density operator is proportional to $\partial_{x}\phi$ because the point splitting is in the $x$ direction. Then using Eq. (A.3), we can also derive the following expected commutation relation:

[TABLE]

Next, we consider the kinetic term by point splitting:

[TABLE]

where we have split the two fermion fields at points $x$ and $x^{\prime}$, respectively, and used the OPE in Eq. (A.8). By taking the limit $x-x^{\prime}\rightarrow 0$, we find the normal ordered kinetic term is

[TABLE]

When there are more than $2$ chiral Majorana fermion fields $\psi_{i}$, we can define multiple boson fields, for instance

[TABLE]

where $1\leq i\leq[N/2]$, and $\eta_{i}$ are the Klein factors satisfying $\{\eta_{i},\eta_{j}\}=\{\eta_{i}^{\dagger},\eta_{j}\}=\{\eta_{i}^{\dagger},\eta_{j}^{\dagger}\}=0$ for $i\neq j$, and $\eta_{i}^{\dagger}\eta_{i}=\eta_{i}\eta_{i}^{\dagger}=1$ [59]. Different flavors of boson $\phi_{i}$ commute with each other. The use of the Klein factors is only to ensure the anticommutation relation between different fermion fields, and can usually be neglected in calculations where only bosonic terms are involved.

We can now derive the commutation relation (2.4) of fermion bilinears, $[(\psi_{i}\psi_{j})(x),(\psi_{k}\psi_{l})(x^{\prime})]$. First, the commutation is trivially zero when $i=j$ or $k=l$. Secondly, when $i,j,k,l$ are all distinct, the commutation is also obviously zero. Thirdly, when $i=k\neq j=l$, we can define $(\psi_{i}+i\psi_{j})/\sqrt{2}=:\mathrel{e^{i\phi}}:$, and use Eqs. (A.3) and (A.9) to find

[TABLE]

Lastly, when $i=k\neq j$ and $j\neq l\neq i$, we can define two boson fields $\phi$ and $\phi^{\prime}$ by $(\psi_{i}+i\psi_{j})/\sqrt{2}=\eta:\mathrel{e^{i\phi}}:$ and $(\psi_{l}+i\psi_{m})/\sqrt{2}=\eta^{\prime}:\mathrel{e^{i\phi^{\prime}}}:$, where $\psi_{m}$ is another arbitrary fermion field, while $\eta$ and $\eta^{\prime}$ are the Klein factors. The commutator in this case is then

[TABLE]

We can then summarize these cases into the commutation relation (2.4). Note that since all $i\psi_{i}\psi_{j}$ ($1\leq i<j\leq N$) form the generators of SO$(N)$ group, so Eq. (2.4) is nothing but the SO$(N)_{1}$ Kac-Moody algebra.

Besides, for each pair of fermion fields $\psi_{i}$ and $\psi_{j}$ ($i\neq j$), from Eqs. (A.9) and (A.11) we can formally rewrite their kinetic term as

[TABLE]

Therefore, we can rewrite the total kinetic term of $N$ flavors of $\psi_{i}$ as

[TABLE]

which is simply the Sugawara construction of the kinetic term.

A.3 Energy momentum tensor

Here we derive the energy momentum tensor of the $1+1$ dimensional chiral SYK model. The energy density can be simply derived from Legendre transformation of the Lagrangian, as shown in Eq. (2.5):

[TABLE]

In general, the full energy momentum tensor can be derived from Noether’s theorem. However, to avoid complications from the point-splitting regularization, we will not use Noether’s theorem, instead we will directly employ the energy momentum conservation law.

To derive the energy current $T^{x}_{\ 0}$, we consider the energy conservation law $\partial_{t}T^{0}_{\ 0}+\partial_{x}T^{x}_{\ 0}=i[H,\mathcal{H}]+\partial_{x}T^{x}_{\ 0}=0$, namely, $\partial_{x}T^{x}_{\ 0}(x)=-i\int dx^{\prime}[\mathcal{H}(x^{\prime}),\mathcal{H}(x)]$. To calculate the right hand side, we first note that from the bosonized expressions above, one can verify that

[TABLE]

[TABLE]

and

[TABLE]

We also need to know the commutator $[(\psi_{i^{\prime}}\psi_{j^{\prime}}\psi_{k^{\prime}}\psi_{l^{\prime}})(x^{\prime}),(\psi_{i}\psi_{j}\psi_{k}\psi_{l})(x)]$. The general expression is rather complicated, which we shall not derive here. Instead, we note that the coefficient of the commutator $[(\psi_{i^{\prime}}\psi_{j^{\prime}}\psi_{k^{\prime}}\psi_{l^{\prime}})(x^{\prime}),(\psi_{i}\psi_{j}\psi_{k}\psi_{l})(x)]$ in $[\mathcal{H}(x^{\prime}),\mathcal{H}(x)]$ is $J_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}J_{ijkl}$. Here we shall only derive the energy current $T^{x}_{\ 0}$ averaged over $J_{ijkl}$ in the large $N$ limit, so one expects this commutator has a nonzero contribution only if its coefficient $J_{i^{\prime}j^{\prime}k^{\prime}l^{\prime}}J_{ijkl}$ does not average to zero, namely, only if $(i^{\prime}j^{\prime}k^{\prime}l^{\prime})=(ijkl)$ up to permutations. The commutator in this case can be easily derived to be

[TABLE]

The other terms will take the form such as $J_{ijk^{\prime}l^{\prime}}J_{ijkl}\partial_{x}(\psi_{k^{\prime}}\psi_{l^{\prime}}\psi_{k}\psi_{l})(x)$, which we expect will have a zero average value, since $J_{ijk^{\prime}l^{\prime}}J_{ijkl}$ is random and has no correlation with $\psi_{k^{\prime}}\psi_{l^{\prime}}\psi_{k}\psi_{l}$.

Therefore, one finds the energy current $T^{x}_{\ 0}$ averaged over $J_{ijkl}$ satisfying

[TABLE]

where in the last line we have used the equation of motion (2.3). Therefore, we arrive at the average energy current

[TABLE]

We note that this expression is exact instead of an average when $N=4$, in which case there is only one interaction $J_{1234}=J$.

We now turn to the momentum density $T^{0}_{\ x}$ and momentum current $T^{x}_{\ x}$. The momentum density can be directly written down by its definition as

[TABLE]

The energy current can then be derived from the momentum conservation law:

[TABLE]

where again we have used the equation of motion (2.3). Therefore, we find

[TABLE]

Note that the momentum density $T^{x}_{\ x}$ is in fact the same as the energy density $T^{0}_{\ 0}$.

Appendix B Saddle Point Approximation for the OTOC

In this section, we give more details on the saddle point approximation for the OTOC integral (4.60) in the large $N$ limit. In particular, the contour deformation requires a careful consideration of both branches of the Riemann surface of the function $\varkappa(p)$ defined in Eq. (4.57).

The function $\varkappa(p)$ in Eq. (4.57) has two branch points in the upper half complex plane at

[TABLE]

where $\mathcal{J}=J/2\pi\in[0,1)$ as we have defined previously. As a result, $\varkappa(p)$ lives on a Riemann surface with two branches. To be specific, we define a branch cut along the straight line segment between the two branch points $p_{b1}$ and $p_{b2}$ as shown by the wavy line in Fig. 12(a), which separate the two branches of Riemann surface. We then denote the branch where the original contour of the OTOC integral (4.60) lies in as branch $1$ (Fig. 12(a)), and the other branch as branch $2$. More specifically, on branch $1$ we have $\varkappa(0)=\frac{\mathcal{J}(\sqrt{3-2\mathcal{J}^{2}}-1)}{1-\mathcal{J}^{2}}=\frac{\beta}{2\pi}\overline{\lambda}$ (see Eq. (4.58)), while on branch $2$ we have $\varkappa(0)=\frac{\mathcal{J}(-\sqrt{3-2\mathcal{J}^{2}}-1)}{1-\mathcal{J}^{2}}$.

On branch $1$, the function $\chi(p)\rightarrow-i(1-\mathcal{J})p-\frac{\mathcal{J}}{1+\mathcal{J}}$ as $p\rightarrow\infty$. Therefore, the numerator $\exp\left[\frac{2\pi t}{\beta}(\varkappa+ipv)\right]$ of the integrand in Eq. (4.60) at velocity $x/t=v$ tends to $\exp\left[i\frac{2\pi t}{\beta}(v-u_{-})p\right]$ at large $p$. Since we are interested in velocities $v$ within the causality cone $v\in[u_{-},u_{+}]$, the original integral contour of Eq. (4.60) along the real axis should always be closed in the upper half complex plane of branch $1$ (counterclockwise), as shown in Fig. 12(a) and (e).

The weight coefficient $\rho(p)=1/\cos[\pi\varkappa(p)/2]$ in the integrand of Eq. (4.60) has two poles

[TABLE]

on the imaginary axis satisfying $\varkappa(p_{1})=\varkappa(p_{2})=1$, and another two poles

[TABLE]

satisfying $\varkappa(\widetilde{p}_{1})=\varkappa(\widetilde{p}_{2})=-1$. The two poles $p_{1}$ and $\widetilde{p}_{2}$ are symmetric about the branch cut, and so do the two poles $p_{2}$ and $\widetilde{p}_{1}$. For interaction strenth $0<J<\pi$, poles $p_{2}$ and $\widetilde{p}_{1}$ are located in branch $1$ above and below the branch cut, respectively, while poles $p_{1}$ and $\widetilde{p}_{2}$ are located in branch $2$ above and below the branch cut, respectively. Such an example at $J=0.6\pi$ is given in Fig. 10(a), where we use solid dots for poles lying in branch 1, and hollow dots for poles lying in branch $2$. For interaction strength $\pi<J<2\pi$, all the poles $p_{1}$, $p_{2}$, $\widetilde{p}_{1}$ and $\widetilde{p}_{2}$ are located in branch $1$. Fig. 10(e) shows such an example at $J=1.2\pi$.

Meanwhile, for a given velocity $x/t=v$, the numerator $\exp\left[\frac{2\pi t}{\beta}(\varkappa+ipv)\right]$ has two saddle points $p_{v}$ and $\widetilde{p}_{v}$ on the imaginary axis satisfying $\varkappa^{\prime}(p_{v})+iv=\varkappa^{\prime}(\widetilde{p}_{v})+iv=0$, which are symmetric about the branch cut:

[TABLE]

The first (second) saddle point $p_{v}$ ($\widetilde{p}_{v}$) moves upwards (downwards) on the imaginary axis as $v$ increases from $v_{-}$ to $v_{+}$. When the velocity $u_{-}<v<1$, both of the two saddle points $p_{v}$ and $\widetilde{p}_{v}$ are located in branch $1$ (Fig. 12(b)&(c)). While when the velocity $1<v<u_{+}$, the both saddle points $p_{v}$ and $\widetilde{p}_{v}$ will move to branch $2$ (Fig. 12(d)).

The saddle point approximation generically works in the following way: if the original integration contour can be deformed into some deepest decent (stationary phase) contour passing a saddle point of the integrand without crossing any poles, the integral can be evaluated by the value of the integrand at the saddle point. If the deformed contour contains multiple deepest decent contours passing several saddle points, the contributions of these saddle points will add, and usually one of the saddle points will dominate the integral. If the contour deformation crosses some pole of the integrand, the integral will be contributed by both the saddle point values and the pole residues.

We first discuss the contour deformation when the velocity $u_{-}<v<1$, in which case both saddle points $p_{v}$ and $\widetilde{p}_{v}$ are in Riemann branch $1$. As shown in Fig. 12(b), the steepest decent contour of saddle point $p_{v}$ hits the other saddle point $\widetilde{p}_{v}$, forming an ellipse surrounding the branch cut in branch $1$; then from $\widetilde{p}_{v}$ the steepest decent contour extends along the imaginary axis to $+i\infty$ of branch $1$ and $-i\infty$ of branch $2$ (across the branch cut). In the figure, we use solid lines to represent contours lying in branch $1$, and dashed lines for contours lying in branch $2$. Hitting the other saddle point $\widetilde{p}_{v}$ makes the deepest decent contour of saddle point $p_{v}$ ambiguous, which is accidental. To resolve this ambiguity, we can add a infinitesimal imaginary part to the velocity $v$, so that $p_{v}$ and $\widetilde{p}_{v}$ slightly deviates from the imaginary axis. After such a procedure, the steepest decent contours passing $p_{v}$ and $\widetilde{p}_{v}$ no longer coincide, which are illustrated by Fig. 12(c): the contour passing $p_{v}$ extends from $-i\infty$ of branch $2$ to the counterclockwise ellipse surrounding the branch cut in branch $1$, and then extends to $+i\infty$ of branch $1$. While the contour passing $\widetilde{p}_{v}$ goes directly from $i\infty$ of branch $1$ to $-i\infty$ of branch $2$. One can easily check that the original contour in Fig. 12(a) can be continuously deformed into the sum of the two steepest decent contours passing $p_{v}$ and $\widetilde{p}_{v}$. Therefore, both saddle points $p_{v}$ and $\widetilde{p}_{v}$ contributes to the integral. If $-i\widetilde{p}_{v}<-ip_{2}$ or $-i\widetilde{p}_{v}<-i\widetilde{p}_{2}$, the deformation will cross the pole $p_{2}$ or $\widetilde{p}_{2}$; while if $-ip_{v}>-ip_{1}$ or $-ip_{v}<-i\widetilde{p}_{1}$, the contour deformation will cross the pole $p_{1}$ or $\widetilde{p}_{1}$. In these cases, the residues of the poles crossed during the deformation will then also contribute.

Second, we discuss the case when the velocity $1<v<u_{+}$, where both saddle points $p_{v}$ and $\widetilde{p}_{v}$ move to Riemann branch $2$. Fig. 12(d) shows such an example, we we use hollow dots for $p_{v}$ and $\widetilde{p}_{v}$ to remind the reader that they are in branch $2$. Now the steepest decent contour passing $p_{v}$ extends from $-i\infty$ of branch $2$ to a clockwise ellipse surrounding the branch cut in branch $2$, and then extends to $i\infty$ of branch $1$ across the branch cut. The other contour passing $\widetilde{p}_{v}$ still extends from $i\infty$ of branch $1$ to $-i\infty$ of branch $2$. In this case, one may find that deforming the original contour of Fig. 12(a) into the sum of the two steepest decent contours passing $p_{v}$ and $\widetilde{p}_{v}$ requires crossing the two branch points, and may wonder whether this is legal. In fact, it is indeed legal to cross the branch points in this configuration. This is because here the Riemann surface with two branches connected by a branch cut is topologically equivalent to a two dimensional sphere, so the branch points are topologically nonsingular and allow a contour to cross through. By drawing the Riemann surface as a sphere, one can show the original contour in Fig. 12(a) can be smoothly deformed into the sum of the two steepest contours of $p_{v}$ and $\widetilde{p}_{v}$ in Fig. 12(d). Therefore, the integral is still contributed by the two saddle points $p_{v}$ and $\widetilde{p}_{v}$, and residues of any of the poles $p_{1}$, $p_{2}$, $\widetilde{p}_{1}$ and $\widetilde{p}_{2}$ crossed during the contour deformation.

All together, in the entire velocity range $u_{-}<v<u_{+}$, to calculate the integral (4.60) we only need to add the contributions from saddle points $p_{v}$ and $\widetilde{p}_{v}$, and the residues of poles $p_{1}$, $p_{2}$, $\widetilde{p}_{1}$ and $\widetilde{p}_{2}$ if crossed during the contour deformation. The evaluation at large $t$ can be further simplified by noting the following two facts. First, for velocity $u_{-}<v<u_{+}$, the residues of the two poles $p_{2}$ and $\widetilde{p}_{2}$ give velocity dependent Lyapunov exponents

[TABLE]

respectively. Therefore, the residues of poles $p_{2}$ and $\widetilde{p}_{2}$ never give an exponential growth in $t$ along velocity $v$, and can always be ignored. Secondly, we can compare the exponent in $t$ contributed by the two saddle points $p_{v}$ and $\widetilde{p}_{v}$. Since $p_{v}$ and $\widetilde{p}_{v}$ are symmetric about the branch cut, we can denote them as $p_{v}=\frac{i\mathcal{J}}{1-\mathcal{J}^{2}}-i\eta$ and $\widetilde{p}_{v}=\frac{i\mathcal{J}}{1-\mathcal{J}^{2}}+i\eta$, respectively, where $\eta$ is real (Eq. (B.4)). Accordingly, the exponents they contribute at velocity $v$ are given by

[TABLE]

[TABLE]

respectively. For velocity $u_{-}<v<u_{+}$ which we are interested in, one can easily see that $\varkappa(p_{v})+ivp_{v}>\varkappa(\widetilde{p}_{v})+iv\widetilde{p}_{v}$, namely, the exponent contributed by saddle point $p_{v}$ is always larger. Therefore, we conclude that the first saddle point $p_{v}$ always dominates over the second saddle point $\widetilde{p}_{v}$.

Given the above two facts, we find that we only need to consider the contributions of the saddle point $p_{v}$ and the two poles $p_{1}$ and $\widetilde{p}_{1}$. When $-i\widetilde{p}_{1}<-ip_{v}<-ip_{1}$, the contour deformation does not pass either pole $p_{1}$ or pole $\widetilde{p}_{1}$, so the integral is given by the saddle point $p_{v}$ contribution (Eq. (4.62)). This is the case in the examples of Fig. 12 (c) and (e), where the former (latter) has $0<J<\pi$ ($\pi<J<2\pi$) and thus pole $p_{1}$ located in branch $2$ (branch $1$). When $-ip_{v}>-ip_{1}$, the OTOC integral is contributed by both the saddle point $p_{v}$ and the residue of pole $p_{1}$, and the exponent contributed by pole $p_{1}$ is always larger in this case. Therefore, the integral is dominated by pole $p_{1}$ (Eq. (4.68)). An example of this case is shown in Fig. 12(d). Lastly, when $-ip_{v}<-i\widetilde{p}_{1}$, the OTOC integral is contributed by both the saddle point $p_{v}$ and the residue of pole $\widetilde{p}_{1}$, and is always dominated by pole $\widetilde{p}_{1}$.

Fig. 13 shows a comparison between the velocity dependent Lyapunov exponent from saddle point approximation and that from the direct numerical integration of Eq. (4.60) for $J=1.2\pi$. This example demonstrates the accuracy of the above saddle point approximation. 444We thank G. Sárosi and M. Mezei for pointing out the critical velocity $v_{*}^{\prime}$ and the verification via numerical integration [70].