Solvable Sachdev-Ye-Kitaev models in higher dimensions: from diffusion to many-body localization

TL;DR

This paper introduces a higher-dimensional generalization of the SYK model that exhibits a transition from a diffusive metal to a many-body localized phase, providing insights into MBL transitions and their critical properties.

Contribution

It proposes a new higher-dimensional SYK model demonstrating a direct transition from diffusion to MBL, with numerical evidence and analysis of critical behavior.

Findings

Diffusive phase with maximal chaos for r<r_c

Vanishing diffusion constant at transition point

Level statistics change from Wigner-Dyson to Poisson

Abstract

Many aspects of many-body localization (MBL) transitions remain elusive so far. Here, we propose a higher-dimensional generalization of the Sachdev-Ye-Kitaev (SYK) model and show that it exhibits a MBL transition. The model on a bipartite lattice has $N$ Majorana fermions with SYK interactions on each site of the $A$ sublattice and $M$ free Majorana fermions on each site the of $B$ sublattice, where $N$ and $M$ are large and finite. For $r$$\equiv$$M/N\!<\!r_c$=1, it describes a diffusive metal exhibiting maximal chaos. Remarkably, its diffusive constant $D$ vanishes [$D$$\propto$$ (r_c-r)^{1/2}$] as $r$$\rightarrow$$r_c$, implying a dynamical transition to a MBL phase. It is further supported by numerical calculations of level statistics which changes from Wigner-Dyson ($r$$<$$r_c$) to Poisson ($r$$>$$r_c$) distributions. Note that no subdiffusive phase intervenes between diffusive and MBL phases. Moreover, the critical exponent $\nu$=0, violating the Harris criterion. Our higher-dimensional SYK model may provide a promising arena to explore exotic MBL transitions.

Figures

Equations

$$H$$

$$\frac{S}{N}$$

$$G_{\psi ,x}^{-1}(i\omega )=-i\omega -{\Sigma }_{\psi ,x}(i\omega ),G_{\eta ,x}^{-1}(i\omega )=-i\omega -{\Sigma }_{\eta ,x}(i\omega ),$$

$${\Sigma }_{\psi ,x}(\tau )=J^2{G}_{\psi ,x}^3(\tau )+\sqrt{r}[t^2{G}_{\eta ,x}(\tau )+t^{\prime 2}{G}_{\eta ,x-1}(\tau )],$$

$${\Sigma }_{\eta ,x}(\tau )=[t^2{G}_{\psi ,x}(\tau )+t^{\prime 2}{G}_{\psi ,x+1}(\tau )]/\sqrt{r},$$

$${\tilde{G}}_{a,x}({\tau }_1,{\tau }_2)$$

$$G_{\psi ,x}^s(\tau )$$

$$G_{\eta ,x}^s(\tau )$$

$$\frac{S_{\text{eff}}}{N}$$

$$D=\frac{\pi }{2}\frac{r\sqrt{1-r}J{t}^{\prime 2}}{(1-r)(t^2+t^{\prime 2})+r^{\frac{3}{2}}{J}^2}.$$

$$F_{\psi \psi ,x{x}^{\prime }}({\tau }_1{\tau }_2{\tau }_3{\tau }_4)=\frac{1}{N^2}⟨T_{\tau }\sum _{ij}{\psi }_{i,x}({\tau }_1){\psi }_{i,x}({\tau }_2){\psi }_{j,x^{\prime }}({\tau }_3){\psi }_{j,x^{\prime }}({\tau }_4)⟩$$

$$=G_{\psi }^s({\tau }_1{\tau }_2)G_{\psi }^s({\tau }_3{\tau }_4)+\frac{1}{N}{\mathcal{F}}_{\psi \psi ,x{x}^{\prime }}({\tau }_1{\tau }_2{\tau }_3{\tau }_4),$$

$$\frac{{\mathcal{F}}_{ab,xy}({\tau }_1{\tau }_2{\tau }_3{\tau }_4)}{G_a^s(\frac{\beta }{2})G_b^s(-\frac{\beta }{2})}$$

$$G_{\psi }^{-1}=-i\omega -{\Sigma }_{\psi },G_{\eta }^{-1}=-i\omega -{\Sigma }_{\eta },$$

$${\Sigma }_{\psi }=\sqrt{r}{\tilde{t}}^2{G}_{\eta },{\Sigma }_{\eta }={\tilde{t}}^2{G}_{\psi }/\sqrt{r},$$

$$G_{\eta }^{-1}$$

$$G_{\psi }^{-1}$$

$$G_{\psi }(\tau )=\frac{\sqrt{r}}{2(1-\gamma )(r-1)}\frac{\text{sgn}(\tau )}{|\tau {|}^2},$$

$$H$$

$$\displaystyle H=\sum_{x=1}^{L}\Big{[}\frac{1}{4!}\sum_{ijkl}J_{ijkl,x}\psi_{i,x}\psi_{j,x}\psi_{k,x}\psi_{l,x}+\sum_{i\alpha}\Big{(}t_{i\alpha,x}i\psi_{i,x}\eta_{\alpha,x}+t_{i\alpha,x}^{\prime}i\eta_{\alpha,x}\psi_{i,x+1}\Big{)}\Big{]},$$

$$S$$

$$\frac{S}{N}$$

$$G_{\psi }(i\omega )=i[(1-r)\pi {]}^{1/4}\frac{\mathrm{sgn}(\omega )}{\sqrt{J|\omega |}},G_{\eta }(i\omega )=\frac{i\sqrt{r}\text{sgn}(\omega )}{(t^2+t^{\prime 2})[(1-r)\pi {]}^{1/4}}\sqrt{J|\omega |}.$$

$$\displaystyle\omega_{1}\approx\frac{J}{2\sqrt{\pi}}(1-r)^{3/2}+\frac{J}{2\sqrt{\pi}}\frac{r^{2}}{(1-r)^{1/2}},~{}~{}~{}~{}~{}~{}\omega_{2}\approx\Big{(}\frac{\sqrt{\pi}}{2}\frac{(t^{2}+t^{\prime 2})^{2}}{J}\frac{\sqrt{1-r}}{r}\Big{)}^{1/3}.$$

$$\frac{I}{NL}$$

$$G_{\psi }\ast {\Sigma }_{\psi }=-1,G_{\eta }\ast {\Sigma }_{\eta }=-1,{\Sigma }_{\psi }=J^2{G}_{\psi }^{q-1}+\sqrt{r}{\tilde{t}}^2{G}_{\eta },{\Sigma }_{\eta }=\frac{1}{\sqrt{r}}{\tilde{t}}^2{G}_{\psi },$$

$$G_{\psi }(\tau )=b\frac{\text{sgn}(\tau )}{|\tau {|}^{2\Delta }},G_{\eta }=\frac{\sqrt{r}{J}^2}{(1-r){\tilde{t}}^2}{G}_{\psi }^{q-1},{\Sigma }_{\psi }=\frac{J^2}{1-r}{G}_{\psi }^{q-1},{\Sigma }_{\eta }=\frac{{\tilde{t}}^2}{\sqrt{r}}{G}_{\psi }$$

$$\displaystyle S=-\frac{\partial F}{\partial T}=-I-T\sum_{\alpha}\Big{(}\frac{\delta I}{\delta G_{\alpha}}\frac{\partial G_{\alpha}}{\partial T}+\frac{\delta I}{\delta\Sigma_{\alpha}}\frac{\partial G_{\alpha}}{\partial T}+\frac{\partial I}{\partial T}\Big{)},$$

$$\displaystyle\frac{1}{2}J^{2}\Big{[}\frac{1}{1-r}-\frac{1}{q^{2}}\Big{]}\int d\tau_{1}d\tau_{2}G^{q}_{\psi}(\tau_{1},\tau_{2})=\frac{1}{2}\beta b^{q}J^{2}\Big{[}\frac{1}{1-r}-\frac{1}{q^{2}}\Big{]}\int_{0}^{\beta}d\tau\frac{1}{\big{(}\frac{\beta}{\pi}\sin\frac{\pi\tau}{\beta}\big{)}^{2}}=0.$$

$$S=\frac{NL}{2}[\text{tr}\log (-{\Sigma }_{\psi })+r\text{tr}\log (-{\Sigma }_{\eta })]=\frac{NL}{2}[\sum _n\log (-{\Sigma }_{\psi }(i{\omega }_n))+r\sum _n\log (-{\Sigma }_{\eta }(i{\omega }_n))].$$

$${\Sigma }_{\psi }(w)=C_{\psi }i\text{sgn}(\omega )|\omega {|}^{1-2\Delta },{\Sigma }_{\eta }(w)=C_{\eta }i\text{sgn}(\omega )|\omega {|}^{2\Delta -1}$$

$$\frac{\partial {\mathcal{S}}_{\text{SYK}}(\Delta )}{\partial \Delta }=-\pi (\frac{1}{2}-\Delta )\tan \pi \Delta .$$

$$\frac{1}{NL}\frac{\partial S_{\psi }(\Delta )}{\partial \Delta }=-\pi (\frac{1}{2}-\Delta )\tan \pi \Delta ,\frac{1}{NL}\frac{\partial S_{\eta }(\Delta )}{\partial \Delta }=-r\pi (\frac{1}{2}-\Delta )\tan \pi (1-\Delta ),\frac{\partial \mathcal{S}(\Delta )}{\partial \Delta }=-\frac{1-r}{1+r}\pi (\frac{1}{2}-\Delta )\tan \pi \Delta$$

Full text

As a foundation of equilibrium statistical mechanics, quantum thermalization and the eigenstate thermalization hypothesis (ETH) for closed quantum systems deutsch1991 ; srednicki1994 ; rigol2008 has received a surge of attention recently. A closed system satisfies ETH and acts as its own heat bath if interactions can thermalize all its own subsystems after long-time dynamic evolution. Nonetheless, it was conceived decades ago and was shown more recently that interacting systems with strong disorders may fail to thermalize but are many-body localized (MBL) anderson1980 ; giamarchi1987 ; altshuler1997 ; polyakov2005 ; altshuler2006 ; huse2007 ; huse2010 ; nayak2013 ; huse2015review ; altman2015review ; moore2016review . Although properties of MBL phases are largely understood now sondhi2013prb ; sondhi2014prb ; demler2014prx ; altman2014prl ; abanin2013 ; huse2014prb ; pollmann2012 ; cenke2016prb ; sid2016arxiv ; dassarma2017 ; xiechen2016 ; huizhai2016 ; swingle2016 ; fradkin2016 ; rqhe2016 ; yuchen2016 , dynamical phase transitions between MBL and thermal phases remain elusive despite the tremendous progress that has been achieved in understanding them lev2015 ; knap2016 ; varma2016 ; altman2015prx ; sid2015prx .

Here, we propose a solvable higher-dimensional generalization of the Sachdev-Ye-Kitaev (SYK) model kitaev2015 ; sachdev1993 and show that it features a dynamical phase transition between a thermal diffusive metal and a MBL phase. The original SYK model consisting of large-$N$ disordered Majorana fermions in zero space dimension was proposed by Kitaev kitaev2015 , which is a generalization of the disordered spin model by Sachdev and Ye sachdev1993 . The SYK model is almost solvable in the large-$N$ limit with an approximate conformal or reparametrization symmetry in low temperature. Moreover, its Lyapunov exponent defined in out-of-time correlations (OTOC) larkin1969 ; kitaev2014 ; shenker2014 saturates the upper bound stanford2016b , implying a holographic dual kitaev2015 to a dilaton gravity theory in nearly AdS2 geometry stanford2016c ; jensen2016 ; verlinde2016a . Various aspects polchinski2016 ; you2016 ; yoon2016 ; stanford2016a ; tezuka2016 ; verbaarschot2017 ; tanasa2017 ; gross2017 ; altland2017 ; franz2017 ; antonio2016 and interesting generalizations fu2016 ; altman2017 ; gu2016 ; gross2016 ; simon2016 ; sachdev2016 ; sachdev2017 ; witten2016 ; klebanov2017 ; volovich2016 ; verlinde2017 ; cenke2017 ; zhou2017 ; gu2017 ; krishnan2016 ; chen2017 ; song2017 of the SYK model have been studied so far.

Our generalized SYK model is defined on bipartite lattices. We focus on the case of one-dimensional lattices, while the model can be generalized to any dimensions. As shown in Fig. 1(a), each unit cell consists of two sites: one site hosting $N$ Majorana fermions with SYK interactions and the other hosting $M$ free Majorana fermions. Two sublattices are coupled via random hopping. The fermion number ratio is denoted as $r$=$M/N$. Here, we consider the case that both $N$ and $M$ are large but finite while the ratio $r$ is fixed. For r$$\ll$$1, SYK physics dominates such that this phase exhibits a finite diffusive constant $D$ and maximal chaos with the Luyapunov exponent satisfying the upper bound $\lambda_{L}$=$2\pi/\beta$, where $\beta$ is the inverse temperature. It is a diffusive metal, similar to the one studied by Gu et al. gu2016 . For r$$\gg$$1, the “free” Majorana fermions on the 1D lattice dominate over the SYK fermions such that weak SYK interactions are irrelevant around the Anderson-localization “fixed” point of free disordered Majorana fermions footnote2 , leading to MBL.

Consequently, we expect that there should be a dynamic phase transition between a thermal (diffusive) phase and a MBL phase as the ratio $r$ varies from small to large. Indeed, for r$$<$$1, our analytical calculations show that the diffusion constant vanishes as D$$\propto$$(1-r)^{1/2} when r$$\to1. This implies that a dynamical phase transition to a MBL phase should occur at $r$=$r_{c}$=$1$. The MBL nature for r$$>$$r_{c} is further supported by our numerical calculations of the many-body level statistics, which qualitatively changes around $r$=$r_{c}$: it follows Poisson distribution for r$$>$$r_{c} but Wigner-Dyson for r$$<$$r_{c}. To the best of our knowledge, it is the first time that a MBL transition is evidenced in a nearly solvable model.

The MBL transition in our generalized SYK models looks qualitatively different from previously studied cases. First, the MBL transition in our generalized SYK model on the 1D lattice occurs between diffusive and MBL phases. This is qualitatively different from previously studied 1D cases, where it was shown by exact diagonalization and real-space renormalization group analysis that a MBL transition can only occur between a subdiffusive phase and a MBL one, both of which have vanishing diffusive constant lev2015 ; knap2016 ; varma2016 ; altman2015prx ; sid2015prx . Second, because of the local criticality in the generalized SYK models, the critical exponent $\nu$=[math] at the MBL transition since the spatial correlation length keeps finite at the transition, which seemly violates the Harris criterion $d\nu\!>\!2$ in systems of $d$ spatial dimensions harris1974 ; chayes1986 ; chandran2015 .

SYK models on 1D lattices: We first introduce the generalized model on 1D lattices, as shown in Fig. 1(a), and consider the cases of more than 1D later. The Hamiltonian of the generalized SYK model in 1D reads:

[TABLE]

where $\psi_{i,x}$ and $\eta_{\alpha,x}$ are SYK Majorana fermions and free Majorana fermions residing on the $A$ site and the $B$ site of the unit cell $x$, respectively, with $i$=1,$\cdots$,$N$ and $\alpha$=1,$\cdots$,$M$. The number of unit cells in the chain is $L$, and the periodic boundary condition is assumed. The SYK fermions on the $A$ sublattice have on-site all-to-all random four-fermion interactions $J_{ijkl,x}$ with mean zero and variance $\langle J_{ijkl,x}^{2}\rangle\!=\!J^{2}3!/N^{3}$. Here, $t$ and $t^{\prime}$ are nearest neighbor random hopping of Majorana fermions within the same unit cell and between neighboring unit cells, respectively, with mean zero and variance $\langle t_{i\alpha,x}^{2}\rangle=t^{2}/\sqrt{MN}$ and $\langle t_{i\alpha,x}^{\prime 2}\rangle=t^{\prime 2}/\sqrt{MN}$. Hereafter, we assume t^{\prime}$$\ll$$t. When we take the large-$N$ limit, we keep the ratio $r\equiv\frac{M}{N}$ fixed. Note that the time-reversal symmetry (\psi$$\rightarrow$$\psi, \eta$$\rightarrow$$-\eta and i$$\rightarrow$$-i) is assumed for the generalized model such that hopping between the same type of fermions is forbidden.

Like in the original SYK model, we use a replica trick to get an effective disorderless model (see the Supplemental Material for details) and introduce bilocal variables: $G^{mm^{\prime}}_{\psi,x}(\tau_{1},\tau_{2})$=$\frac{1}{N}\sum_{i=1}^{N}\psi^{m}_{i,x}(\tau_{1})\psi^{m^{\prime}}_{i,x}(\tau_{2})$ and $G^{mm^{\prime}}_{\eta,x}(\tau_{1},\tau_{2})\!=\!\frac{1}{M}\sum_{\alpha=1}^{M}\eta^{m}_{\alpha,x}(\tau_{1})\eta^{m^{\prime}}_{\alpha,x}(\tau_{2})$, as well as $\Sigma^{mm^{\prime}}_{\psi,x}(\tau_{1},\tau_{2})$, $\Sigma^{mm^{\prime}}_{\eta,x}(\tau_{1},\tau_{2})$ as Legendre multipliers to implement the above identities, where $m,m^{\prime}$ are replica indices. At the large-$N$ limit, different replicas do not interact, so the bilocal fields are diagonal in replica indices, i.e., $G^{mm^{\prime}}=G\delta^{mm^{\prime}}$ and $\Sigma^{mm^{\prime}}=\Sigma\delta^{mm^{\prime}}$. We obtain the following effective action:

[TABLE]

where $G$ and $\Sigma$ are collective bosonic modes and $\iint\equiv\int d\tau_{1}d\tau_{2}$ (integration over two times appears because the replica trick couples fields at different times). The large-$N$ structure is manifest in the effective action above. The saddle-point equations obtained by varying these collective modes are

[TABLE]

where $\tau\!=\!\tau_{1}\!-\!\tau_{2}$. These saddle-point equations are equivalent with Schwinger-Dyson equations obtained from diagrammatic methods stanford2016a ; gu2016 .

Diffusive metals: For r$$\ll1, it is expected that the SYK fermions dominate over the free Majorana fermions in the infrared altman2017 . Similar to features of the original SYK model, the time-derivative terms in Eq. (2) or the $-i\omega$ terms in Eq. (3) are irrelevant in low energy. Remarkably, Eqs.(3-5) in the infrared limit of $\omega\to 0$ are invariant under global (site-independent) reparametrization of time $\tau\rightarrow f(\tau)$,

[TABLE]

where $f^{\prime}(\tau)$=$\frac{df}{d\tau}$, $a$=$\psi$ or $\eta$, and the scaling dimensions $\Delta_{\psi}$=$\frac{1}{4}$, $\Delta_{\eta}$=$\frac{3}{4}$. Like in the SYK model, this is an emergent time reparametrization symmetry at low energy that is explicitly broken by high-energy degrees of freedom in the microscopic model [or the time derivative-terms in the effective action Eq. (2)].

Helped by the emergent reparametrization symmetry, we obtain the following solutions of the Schwinger-Dyson equations in the infrared [Eqs.(3-5) in the limit of $\omega\to 0$]:

[TABLE]

The solutions above are spatially uniform while nontrivial in the time direction, exhibiting local criticality gu2016 ; si2001 ; faulkner2011 . Note that the saddle-point solutions are valid below a cutoff frequency $\omega_{c}$ which scales as \omega_{c}$$\sim$$(1-r)^{1/6} when r$$\rightarrow$$1 (see the SM for details).Using the saddle-point solutions, we find that the zero-temperature entropy per unit cell per Majorana fermion is given by (see the SM for details) $\mathcal{S}=\frac{1-r}{1+r}\mathcal{S}_{\text{SYK}}$, where $\mathcal{S}_{\text{SYK}}=\frac{4\mathcal{C}+\pi\log 2}{8\pi}\approx 0.232$, and $\mathcal{C}\approx 0.916$ is Catalan constant. When $r\rightarrow 1$, the zero temperature entropy vanishes which implies a hint that there is a phase transition at $r_{c}$=$1$.

Note that the saddle-point solutions of Eqs.(7-8) spontaneously break the continuous reparametrization symmetry to SL(2,$R$). Owing to the spontaneous and explicit breaking pattern, site-dependent reparametrization modes $\epsilon_{x}\!=\!f_{x}(\tau)-\tau$ would contribute dominant low-energy fluctuations on top of the saddle-point one, which determine the low-energy physics especially dynamics like transport and the butterfly effect. Note that because $t^{\prime}\ll t$ the relative reparametrization fluctuation within each unit cell (namely, f_{\psi,x}$$-$$f_{\eta,x}) is at high energy and does not affect the physics in the low energy we consider here.

The effective action for the reparametrization modes is given by fluctuations around the saddle-point one, i.e., $S_{\text{eff}}[f]=S[\tilde{G}(f)]-S[G(\tau)]$, where $\tilde{G}_{a,x}(\tau_{1},\tau_{2})=f^{\prime\Delta_{a}}_{x}(\tau_{1})f^{\prime\Delta_{a}}_{x}(\tau_{2})G_{a}(f_{x}(\tau_{1}),f_{x}(\tau_{2}))$ is the Green’s function of $a$=$\psi,\eta$ fermions associated with the spatially dependent time reparametrization $f_{x}(\tau)$=$\tau+\epsilon_{x}$. Note that though the saddle-point solution of $G_{a,x}$ in Eqs.(7-8) is homogenous, its fluctuation associated with the reparametrization modes is generically inhomogeneous. By assuming weak reparametrization $\epsilon_{x}$ as well as performing $\varepsilon$ expansion and series summation (see the Supplemental Material for details), we obtain the effective action up to the quadratic in $\epsilon$,

[TABLE]

where $\omega_{n}$=$2\pi n/\beta$ is the Matsubara frequency, $p$ is momentum, and $\epsilon_{x}(\tau)$=$\frac{1}{\sqrt{L}\beta}\sum_{n,p}\epsilon_{\omega_{n},p}e^{-i\omega_{n}\tau+ipx}$. As shown in the Supplemental Material, $\alpha_{1}$=\frac{1}{64\pi^{2}}\Big{(}\frac{\sqrt{1-r}}{J}+\frac{J}{t^{2}+t^{\prime 2}}\sqrt{\frac{r^{3}}{1-r}}\Big{)} and $\alpha_{2}$=$\frac{1}{128\pi}\frac{rt^{\prime 2}}{t^{2}+t^{\prime 2}}$. Since $J$ and $t$ are both relevant at the UV Gaussian point, they increase as energy scales lower. Thus, $\alpha_{1}$ becomes extremely small due to the emergent reparametrization symmetry, while $\alpha_{2}$ is also small in the homogenous limit, i.e., $t^{\prime}\ll t$. These lead to strong fluctuations of reparametrization modes which dominate the low-energy dynamics.

Having obtained the effective action for the reparametrization modes, we are ready to calculate their contributions to energy transport and OTOC in the limit of N$$\gg$$\beta J$$\gg1. The energy density for small momentum is given by $T_{\omega_{n},p}$=\frac{iN\alpha_{1}}{4\pi}\omega_{n}[\omega_{n}^{2}-\big{(}\frac{2\pi}{\beta}\big{)}^{2}]\epsilon_{\omega_{n},p}. Using the effective action for reparametrization modes, the real-frequency correlator (see the Supplemental Material for details) $\langle T_{-\omega,-p}T_{\omega,p}\rangle=\frac{N\alpha_{1}}{8\pi\beta^{2}}\frac{Dp^{2}}{-i\omega+Dp^{2}}$, where the diffusive constant $D$ is

[TABLE]

Some remarks come with this expression for diffusive transport of energy. First, when $t^{\prime}=0$, different unit cells decouple from each other, and the diffusive constant $D$ vanishes as expected. On the other hand, when r$$\rightarrow[math], the free Majorana fermions vanish, and the system becomes decoupled islands of SYK Majorana fermions and cannot conduct energy. A more interesting observation is that when r$$\rightarrow$$1, the diffusive constant scales as $D\propto(1-r)^{1/2}$, and we expect the system enters a localized phase.

We are now in a position to calculate the OTOC. Consider the following four-point correlation function

[TABLE]

where $T_{\tau}$ denotes imaginary time ordering, $G^{s}_{\psi}$ is given by the saddle-point solutions in Eqs.(7-8), and $\mathcal{F}_{\psi\psi,xx^{\prime}}$ is the connected part coming from the fluctuations around the saddle-point, $\mathcal{F}_{\psi\psi,xx^{\prime}}(\tau_{1}\tau_{2}\tau_{3}\tau_{4})\equiv\langle\delta G_{\psi,x}(\tau_{1}\tau_{2})\delta G_{\psi,x^{\prime}}(\tau_{3}\tau_{4})\rangle$ and is dominated by the reparametrization modes. Similar calculations apply to the OTOC of other operators. In order to evaluate the OTOC, let $\tau_{1}$=$\beta$+$it$, $\tau_{4}$=$\frac{3\beta}{4}$, $\tau_{2}$=$\frac{\beta}{2}$+$it$, $\tau_{3}$=$\frac{\beta}{4}$, and we arrive at (see the Supplemental Material for details)

[TABLE]

with $v_{B}^{2}$=$\frac{2\pi}{\beta}D$ and $a,b$=$\psi,\eta$. We first note that the quantum analog of the Lyapunov exponent defined by OTOC in this phase still saturates the bound $\lambda_{L}$=$\frac{2\pi}{\beta}$ stanford2016b . Second, the butterfly velocity, Lyapunov exponent, and diffusive constant here satisfy a simple and elegant relation: $D$=$\frac{v_{B}^{2}}{\lambda_{L}}$ hartnoll ; blake2016prl . Such a relation was previously obtained in incoherent black holes blake2016 ; blake2017 and higher-dimensional generalizations of the SYK model gu2016 ; sachdev2016 ; footnote . As the butterfly velocity v_{B}$$\propto$$(1-r)^{\frac{1}{4}} is vanishing for r$$\to1, it further indicates that the system shall undergo a localization transition as $r$ crosses the critical value $r_{c}$=1.

MBL phase: For $r\gg 1$, it is expected that the Anderson localization of “free” Majorana fermions for large but finite $N$ dominate in determining low-energy physics and the SYK interaction $J$ is irrelevant. Consequently, the system should fall into a localized phase footnote2 . Similar to the case of $r\ll 1$, we also make a translational invariant ansatz for $r\gg 1$, with which the saddle-point equation can be approximated by

[TABLE]

where ${\tilde{t}}^{2}\equiv t^{2}+t^{\prime 2}$. The exact solutions of the above Schwinger-Dyson equations are obtained in the Supplemental Material. Here, let’s explicitly expand the inverse propagators $G_{a}^{-1}$ around small frequency:

[TABLE]

Although in $G^{-1}_{\eta}$ the bare term $\propto-i\omega$ is renormalized by a factor $\frac{r}{r-1}$, its self-energy is subdominant at low energy, indicating the free Gaussian fixed point of $\eta$ Majorana fermions is stable. (In the limit of r$$\rightarrow$$\infty, G_{\eta}^{-1}$$\rightarrow$$-i\omega, as expected from a free theory.) However, for the $\psi$ fermions, the self-energy actually dominates the behavior of $G_{\psi}$ in low energy, which generates a large anomalous dimension to $\psi$. For simplicity, we keep the leading term in Eq. (16) and make a Fourier transformation,

[TABLE]

where $\gamma\approx 0.577$ is the Euler-Gamma constant. From the propagator of $\psi$ fermions, one deduces its scaling dimension $[\psi]$=$1$, as expected from the r$$\gg1 free fixed point.

Now we explore effects of the SYK interaction $J$. Including $J$ terms leads to a correction to the self-energy, $\delta\Sigma_{\psi}(\tau)$=$J^{2}G_{\psi}(\tau)^{3}$=$\frac{r^{3/2}J^{2}}{8(1-\gamma)^{3}(r-1)^{3}}\frac{\text{sgn}(\tau)}{|\tau|^{6}}$. By a Fourier transformation, $\delta\Sigma_{\psi}(i\omega)\propto\omega^{5}$, which is subdominant in low energy, compared with leading terms in Eq. (16). The same is true for $\delta\Sigma_{\eta}(i\omega)$. Thus, we can conclude the free fixed point with $[\eta]$=[math] and $[\psi]$=$1$ is stable against weak interaction $J$, which self-justifies the assumption we have made. One important consequence is that, as all levels of the free Majorana fermions with random hopping are localized for large and finite $N$ footnote2 ; herbert1971 , MBL emerges in the presence of the weak but irrelevant SYK interaction $J$ altman2013 . It is consistent with vanishing diffusive constant for r$$>$$1. Note that the system has $(M-N)L$ single-particle zero modes localized on the B sites for $r>1$. However, these localized states do not change the MBL phase because they are isolated from the rest of the many-body states and only cause macroscopic degeneracies. One way of removing these extensive localized zero modes but preserving the low-energy physics is to add weak quadratic couplings on B sites, as we show in the Supplemental Material.

Numerical evidences of MBL transitions: We now show numerical evidence of such a phase transition between the thermal and MBL phases. For a MBL phase, its level statistics satisfies the Poisson distribution according to the Berry-Tabor conjecture berry1977 while a thermal phase’s level statistics follows the Wigner-Dyson (WD) distribution. Suppose $\{E_{n}\}$ denotes many-body eigenstate energies in an ascending order and the level spacings between adjacent eigenstates are $\Delta_{n}$=E_{n+1}$$-$$E_{n} with \Delta_{n}$$\geq[math]. The ratio between two consecutive gaps $s_{n}$=$\frac{\Delta_{n+1}}{\Delta_{n}}$ can be employed to characterize the level statistics huse2007 ; atas2013 . The distribution of ratios in MBL phases follows Poisson level statistics $p(s)$=$\frac{1}{(1+s)^{2}}$, while in thermalized phases, it follows WD level statistics $p(s)$=$\frac{81\sqrt{3}}{4\pi}\frac{(s+s^{2})^{2}}{(1+s+s^{2})^{4}}$ (assuming Gaussian unitary ensemble).

Following Ref. you2016 , we plot the distribution of $\log s$, i.e., $P(\log s)=p(s)s$, as shown in Fig. S1. The data are obtained from exactly diagonalizing the model with $J$=$t$=1, $t^{\prime}$=0.5 on a six-site chain with $N$+$M$=10 Majorana fermions per unit cell. The distribution for $(N,M)$=(6,4), (5,5), (4,6) is shown in Fig. 2(a,b,c), respectively. When N$$>$$M (namely, r$$<$$1), the distribution in Fig. S1 follows that of WD; when N$$<$$M (namely, r$$>$$1), the distribution in Fig. S1 follows that of Poisson. When $N$=$M$ (namely, $r$=$r_{c}$=$1$), the distribution in Fig. S1 is in transition between Poisson and WD. Our numerical results imply that a dynamic transition from a thermal to a MBL phase occurs around $r$=$1$. As mentioned before, for $r>1$, each many-body energy level has an extra degeneracy due to the presence of the single-particle zero modes localized on B sites. In the calculation of energy level statistics, we have ignored these trivial degeneracy. The degeneracy can be lifted by adding weak quadratic couplings on B sites, and for such modified case we also calculated the level statistics and obtained the qualitatively same results, as shown in the Supplemental Material.

SYK model on 2D lattices: Our construction of the SYK models in 1D can be straightforwardly generalized to more than 1D. For instance, we consider the generalization to the square lattice as shown in Fig. 3. Each unit cell consists of two sites represented by a square and a disk, where $N$ SYK Majorana fermions and $M$ free Majorana fermions reside, respectively. The model is given by

[TABLE]

where ${\bf x}$ represents unit cells, and ${\bf b}_{a}$ label the vectors connecting neighboring unit cells with ${\bf b}_{1}$=${\bf e}_{1}$=$(1,0)$, ${\bf b}_{2}$=${\bf e}_{2}$=$(0,1)$, and ${\bf b}_{3}$=$(1,1)$. Similarly, $\langle J_{ijkl,{\bf x}}^{2}\rangle$=$3!J^{2}/N^{3}$, $\langle t_{i\alpha,{\bf x}}^{2}\rangle$=$t^{2}/\sqrt{MN}$, and $\langle t_{i\alpha,{\bf x}a}^{\prime 2}\rangle$ =$t^{\prime 2}_{a}/\sqrt{MN}$. (Note that the limit of $t^{\prime}_{3}$=[math] corresponds to the honeycomb lattice). The analysis of the generalized model on 2D and higher-dimensional lattices goes like the 1D chain case. For r$$<$$1, the generalized models on 2D lattices possess similar features, including diffusive energy transport, zero-temperature entropy, and maximum quantum chaos, which are the same as the model on a 1D chain. For instance, the diffusive constant in 2D as a function of $r$ is also given by Eq. (10), which is plotted in Fig. 3(b). For r$$\to[math], $D\propto r$ because in diffusive metal, the SYK Majorana fermions diffuse via free Majorana fermions; while for r$$\rightarrow$$1, $D\propto(1-r)^{1/2}$, which indicates that the system could undergo a dynamical transition into a MBL phase.

Concluding remarks: We have shown that the MBL transition in the generalized SYK models is qualitatively distinct from previously studied ones in other models like the XXZ model in various ways. Intuitively, we think that the qualitative differences are mainly due to the large-$N$ degrees of freedom on each site in the generalized SYK models. In the large-$N$ limit, due to the all-to-all interactions, we can define an effective dimensions d_{\text{SYK}}$$\to$$\infty such that the effective dimensions of the generalized model on the $d$-dimensional lattice is $d^{\ast}$=d_{\text{SYK}}$$+$$d, which approaches infinity. As a consequence, for the SYK model on the $d$=$1$ lattice, there is no subdiffusive phase around the MBL transition because its effective space dimension $d^{\ast}$ is much larger than 1. Moreover, the Harris criterion is not violated by $\nu$=[math] when $d^{\ast}$ is considered as the effective space dimension.

Note that there are questions that remain open. To inspire readers, we provide a few here. First, what is the critical theory governing this MBL transition? Our analysis cannot be applied directly at $r$=$1$, and the critical theory remains unknown. Second, is time-reversal symmetry spontaneously broken in the MBL phase (r$$\gtrsim$$1)? Although we have shown that the $J$ term is irrelevant when r$$\gg$$1, it is possible to be dangerously irrelevant for r$$\gtrsim$$1. Third, how robust is the critical point when other types of interactions are included in the model?

Acknowledgements: We would like to thank Xin Dai, Yingfei Gu, David Huse, Xiaoliang Qi, Cenke Xu, and Shixin Zhang for helpful discussions. This work is supported in part by the MOST of China under Grant No. 2016YFA0301001 (H. Y.) and by the NSFC under Grant No. 11474175 (S.-K. J. and H. Y.).

Note added: After the completion of the present work, we became aware of an upcoming work cenke studying a different generalization of the SYK model with a zero-temperature insulating phase (but not MBL).

I Supplemental Material

I.1 A. Replica action

The Hamiltonian in main text is given by

[TABLE]

where $\psi$ and $\eta$ are SYK Majorana fermions and free Majorana fermions, $J_{ijkl,x}$, $t_{i\alpha,x}$ and $t_{i\alpha,x}^{\prime}$ are independent random couplings with zero mean and variance $\langle J_{ijkl,x}^{2}\rangle=J^{2}3!/N^{3}$, $\langle t_{i\alpha}^{2}\rangle=t^{2}/\sqrt{MN}$, and $\langle t_{i\alpha}^{\prime 2}\rangle=t^{\prime 2}/\sqrt{MN}$. Note that $N$, $M$ are the numbers of SYK Majorana fermions and free Majorana fermions in each sub-lattice respectively, while $L$ is the number of unit cells in the chain. Replica trick utilizes the identity $\log Z=\lim_{n\rightarrow 0}\frac{e^{n\log Z}-1}{n}=\lim_{n\rightarrow 0}\frac{Z^{n}-1}{n}$. Instead of disorder averaging the logarithm of partition function which is difficult to do, one averages over $n$ copies of the system, then take $n\rightarrow 0$ limit. After the disorder average, the replica action is

[TABLE]

where $m,m^{\prime}$ is replica index. As explained in main text, here we consider the diagonal parts. Introduce two collective modes, $G_{\psi,x}(\tau_{1},\tau_{2})=\frac{1}{N}\sum_{i}\psi_{i,x}(\tau_{1})\psi_{i,x}(\tau_{2})$, $G_{\eta,x}(\tau_{1},\tau_{2})=\frac{1}{M}\sum_{\alpha}\eta_{\alpha,x}(\tau_{1})\eta_{\alpha,x}(\tau_{2})$, and corresponding Legendre multipliers $\Sigma_{\psi}$, $\Sigma_{\eta}$, one arrives at

[TABLE]

where Majorana fermions are integrated out.

I.2 B. Two-point functions in diffusive metal

The solutions of saddle point function in frequency domain are given by Fourier transforming the time domain solutions, i.e.,

[TABLE]

While finite temperature solutions are obtained by substitution $\tau\rightarrow\tan\frac{\pi\tau}{\beta}$, i.e., $G_{a}(\tau)=\Lambda_{a}\frac{\text{sgn}(\tau)}{|\frac{\beta}{\pi}\sin\frac{\pi\tau}{\beta}|^{\Delta_{a}}}$. Here \Lambda_{\psi}\equiv\Big{(}\frac{1-r}{4\pi J^{2}}\Big{)}^{1/4}, \Lambda_{\eta}=\frac{1}{2(t^{2}+t^{\prime 2})}\Big{(}\frac{r^{2}J^{2}}{4\pi^{3}(1-r)}\Big{)}^{1/4}, $\Delta_{\psi}=\frac{1}{4}$ and $\Delta_{\eta}=\frac{3}{4}$. Following the same lines in altman2017 , the cutoff frequency for validity of these conformal solutions are given by $\omega_{c}=\text{min}(\omega_{1},\omega_{2})$, where

[TABLE]

When $r$ approach the transition point $r_{c}=1$, $\omega_{c}\sim(1-r)^{1/6}\rightarrow 0$, the saddle point solutions break down.

I.3 C. Zero temperature entropy

In this subsection, we use $I$ to denote action while $S$ denotes entropy to avoid confusion. The homogenous saddle-point action generalized to $q$-body interaction is given by

[TABLE]

where time derivative terms are neglected since they are irrelevant in zero temperature. Equations of motion are

[TABLE]

from which we have $J^{2}G_{\psi}\ast G_{\psi}^{q-1}=-(1-r)$ with the solutions

[TABLE]

where $J^{2}b^{q}\pi=(1-r)(\frac{1}{2}-\Delta)\tan\pi\Delta$, and $\Delta=\frac{1}{q}$. Free energy is given by $F=TI$, from which we can get entropy through

[TABLE]

The second term of Eq. (S9) vanishes when one plugs the solutions; the zero temperature entropy is given by $S_{0}=-\lim_{T\rightarrow 0}I$. Plug the solutions into the second term of the action Eq. (S6), we have

[TABLE]

The vanishing of above integral can be seen by analytical continuing $\tau\rightarrow\tau+it$. Thus only the first term in saddle point action is not vanishing,

[TABLE]

Fourier transforming the self-energy leads to

[TABLE]

where $C_{\psi}$ and $C_{\eta}$ are two constants independent of $\Delta$. These constants are not important since we will take derivative with respect to $\Delta$. According to Ref. stanford2016a , the zero temperature entropy for SYK model satisfies

[TABLE]

where $\mathcal{S}_{\text{SYK}}=\frac{S_{\text{SYK}}}{N}$. Here, analogous to SYK model, we have

[TABLE]

and where $\mathcal{S}=\frac{1}{(N+M)L}(S_{\psi}+S_{\eta})$. The boundary condition can be set by $\mathcal{S}(\frac{1}{2})=0$ because when $q=2$, the system has only quadratic term and unique ground state. Then for our interest case, $\Delta=\frac{1}{4}$, the ground state entropy is given by

[TABLE]

where $\mathcal{S}_{\text{SYK}}=\frac{4\mathcal{C}+\pi\log 2}{8\pi}\approx 0.232$, here $\mathcal{C}\approx 0.916$ is Catalan constant.

I.4 D. Effective action for reparametrization modes in diffusive metal

Inspired by the reparametrization symmetry in infrared, we redefine $\Sigma(\tau_{1},\tau_{2})\rightarrow\Sigma(\tau_{1},\tau_{2})+\delta(\tau_{1}-\tau_{2})\partial_{\tau_{2}}$, and bring the action to

[TABLE]

The redefinitions effectively collect the time derivative term to $S_{UV}$. Moreover, since each unit cell decouples from others in $S_{IR}$, the saddle point equations given by $S_{IR}$ are exactly conformal invariant whose solutions are obtained above. Since $S_{IR}$ has reparametrization symmetry, it vanishes for reparametrization modes. However, this symmetry is explicitly broken by UV part, i.e., $S_{UV}$ will give a small action to reparametrization modes which dominants the four-point correlator in the infrared. The effective action for reparametrization modes is given by fluctuations around saddle point action, i.e., the effective action for reparametrization modes reads ${S[f]}=S_{UV}[\tilde{G}(f)]-S_{UV}[G(t)]$. For the first term in $S_{UV}$, we obtain the effective action using $\varepsilon$-expansion yoon2016 ; verlinde2017 , where $\varepsilon=\frac{1}{2}-\Delta_{\psi}$:

[TABLE]

where \{f,\tau\}=\frac{f^{\prime\prime\prime}}{f^{\prime}}-\frac{3}{2}\Big{(}\frac{f^{\prime\prime}}{f^{\prime}}\Big{)}^{2} denotes Schwartz derivative. For small reparametrization $f_{x}=\tau+\epsilon_{x}(\tau)$,

[TABLE]

where $\epsilon_{x}(\tau)=\frac{1}{2\pi\sqrt{L}}\sum_{n,p}\epsilon_{n,p}e^{-in\tau+ipx}$, and we have set $\beta=2\pi$ for simplicity. For second term in $S_{UV}$, the reparametrization modes are defined as $\tilde{G}_{\psi,x}=G_{\psi}+\delta G_{\psi,x},\tilde{G}_{\eta,x}=G_{\eta}+\delta G_{\eta,x}$. Though the saddle point solutions are uniform in spatial direction, their fluctuations are position-dependent gu2016 . Then the effective action reads

[TABLE]

where $\iint=\int d\tau_{1}d\tau_{2}$. Explicitly, for conformal solutions $G(\tau)=\Lambda\frac{\text{sgn}(\tau)}{|2\sin\frac{\tau}{2}|^{2\Delta}}$, the reparametrization gives rise to

[TABLE]

where $h_{n}(\tau)\equiv\sin\frac{n\tau}{2}\cot\frac{\tau}{2}-n\cos\frac{n\tau}{2}$ and $\tau=\tau_{1}-\tau_{2}$, $\bar{\tau}=\frac{1}{2}(\tau_{1}+\tau_{2})$. Then

[TABLE]

Plug these results in to $S^{(2)}$, we get

[TABLE]

Finally, in terms of the infinitesimal modes, the effective action is given by

[TABLE]

where \alpha_{1}=\frac{1}{64\pi^{2}}\Big{(}\frac{\sqrt{1-r}}{J}+\frac{J}{t^{2}+t^{\prime 2}}\sqrt{\frac{r^{3}}{1-r}}\Big{)}, $\alpha_{2}=\frac{1}{128\pi}\frac{rt^{\prime 2}}{t^{2}+t^{\prime 2}}$. To restore the dimension, note that $\text{dim}[\epsilon_{x}(\tau)]=-1$, and $\text{dim}[\epsilon_{n}]=-2$, where $\text{dim}[...]$ denotes engineer dimension (not to confuse with scaling dimension), we obtain

[TABLE]

where we have expanded for small momentum.

I.5 E. Diffusive constant in diffusive metal

According to Noether’s theorem, the reparametrization modes in time direction couple to energy density, i.e., $\delta S=\int d\tau\epsilon\partial_{\tau}T$, where $\epsilon$ is reparametrization mode and $T$ is energy density. At zero momentum limit, the effective action is given by

[TABLE]

Then one finds that $T=\frac{\alpha_{1}N}{2}(\epsilon^{\prime\prime\prime}+\epsilon^{\prime})$ or in frequency domain, $T_{n}=\frac{i\alpha_{1}N}{4\pi}(n^{3}-n)\epsilon_{n}$. The correlation of energy density is given by

[TABLE]

For small momentum, one can directly generalize this function,

[TABLE]

Plug the propagator from effective action into the correlation, we have

[TABLE]

Note that $\text{dim}[T(\tau)]=1$, we can also restore temperature:

[TABLE]

Analytical continuation from upper half-plane, i.e. $i\omega_{n}\rightarrow\omega+i\delta$, the retarded correlation function is

[TABLE]

where $D\equiv\frac{\alpha_{2}}{\alpha_{1}}$ and we have implicitly extracted a contact term gu2016 ; sachdev2016 , i.e., $\langle T_{-\omega,-p}T_{\omega,p}\rangle-\langle T_{-\omega,0}T_{\omega,0}\rangle$.

I.6 F. Chaos and butterfly velocity in diffusive metal

Four-point functions are defined as

[TABLE]

where $\mathcal{F}_{ab,xy}(\tau_{1},\tau_{2},\tau_{3},\tau_{4})\equiv\langle\delta G_{a,x}(\tau_{1},\tau_{2})\delta G_{b,y}(\tau_{3},\tau_{4})\rangle$ and $a,b=\psi,\eta$. Using the translational symmetry,

[TABLE]

with $\mathcal{F}_{ab,p}(\tau_{1},\tau_{2},\tau_{3},\tau_{4})=L\langle\delta G_{a,-p}(\tau_{1},\tau_{2})\delta G_{b,p}(\tau_{3},\tau_{4})\rangle$. Thus

[TABLE]

For OTOC, let $\tau_{1}=\beta+it,\tau_{4}=\frac{3\beta}{4},\tau_{2}=\frac{\beta}{2}+it,\tau_{3}=\frac{\beta}{4}$, then $x_{12}=\frac{\beta}{2},x_{34}=-\frac{\beta}{2},y_{12}=\frac{3\beta}{4}+it,y_{34}=\frac{\beta}{2}$, and set $\beta=2\pi$ for simplicity, we have

[TABLE]

Plug the propagator obtained from effective action into above equation, we have

[TABLE]

To evaluate above summation, consider the integral

[TABLE]

Making a large semicircle contour to $w\rightarrow+\infty$ of complex plane gu2016 , one finds according to residue theorem

[TABLE]

or

[TABLE]

There is no exponential term in $K$ gu2016 , thus the only exponential growth part is

[TABLE]

Fourier transform to real space, then we have

[TABLE]

The exponent growth of OTOC is now given by (we have restored dimensions)

[TABLE]

Since $x,y$ is dimensionless here (we set the lattice constant to 1), $[v_{B}]=1$ and the butterfly velocity is given by

[TABLE]

When approaching transition point, $r\rightarrow 1$, butterfly velocity vanishes indicating a MBL phase.

I.7 G. Many-body localized phase

Similar to the case of $r\ll 1$, we also make a translational invariant ansatz for $r\gg 1$, with which the saddle point equation can be approximated by

[TABLE]

where ${\tilde{t}}^{2}\equiv t^{2}+t^{\prime 2}$. The exact solutions of the above Schwinger-Dyson equations can be obtained:

[TABLE]

The self-energy part $\Sigma_{a}$ seems to dominate the propagator $G_{a}$ at low frequency or energy due to the terms $\propto\frac{1}{\omega}$ in the denominators of Eqs. (S50-S51). It is true for $\Sigma_{\psi}$; but for $\Sigma_{\eta}$ the $\frac{1}{\omega}$ terms are cancelled in the limit of low frequency.

I.8 G. A modified model

In order to lift the zero modes in MBL phase, we add random quadratic couplings into the Hamiltonian altman2017 ; song2017 . The modified Hamiltonian is

[TABLE]

where $J_{ijkl,x},t_{i\alpha,x},t^{\prime}_{i\alpha,x}$ are random variables same as before, and $V_{\alpha\beta,x}$ refers to onsite random quadratic couplings satisfying Gaussian distribution with $\langle V_{\alpha\beta,x}\rangle=0$, $\langle V_{\alpha\beta,x}^{2}\rangle=\frac{V^{2}}{M}$. We calculate the distribution of level-spacing ratios for the modified Hamiltonian as shown in Fig. S1, from which one can infer that this system also has a transition from diffusive metal to MBL phase.