Quench, thermalization and residual entropy across a non-Fermi liquid to Fermi liquid transition

TL;DR

This paper investigates the thermalization dynamics of a model transitioning from a non-Fermi liquid to a Fermi liquid, revealing rapid thermalization after quenches to NFL and complex slow-quench behavior influenced by residual entropy.

Contribution

It introduces a model coupling SYK fermions to lead fermions, analyzing thermalization and residual entropy effects across the NFL-FL transition with novel non-analytic scaling results.

Findings

Rapid thermalization after sudden quench to NFL.

Multiple prethermal regimes after quench to FL.

Intermediate-$\tau$ non-analytic power-law energy scaling.

Abstract

We study the thermalization, after sudden and slow quenches, of an interacting model having a quantum phase transition from a Sachdev-Ye-Kitaev (SYK) non-Fermi liquid (NFL) to a Fermi liquid (FL). The model has SYK fermions coupled to non-interacting lead fermions and can be realized in a graphene flake connected to external leads. After a sudden quench to the NFL, a thermal state is reached rapidly via collapse-revival oscillations of the quasiparticle residue of the lead fermions. In contrast, the quench to the FL, across the NFL-FL transition, leads to multiple prethermal regimes and much slower thermalization. In the slow quench performed over a time $\tau$, we find that the excitation energy generated has a remarkable intermediate-$\tau$ non-analytic power-law dependence, $\tau^{-\eta}$ with $\eta<1$, which seemingly masks the dynamical manifestation of the initial residual entropy of the SYK fermions. The power-law scaling is expected to eventually break down for $\tau\to\infty$, signaling a violation of adiabaticity, due to the residual entropy present in the SYK fermions.

Equations

$${\mathcal{H}}_c$$

$${\mathcal{H}}_{\psi }$$

$${\mathcal{H}}_{c\psi }(t)$$

$${\mathcal{H}}_c^{(NI)}$$

$$\displaystyle\begin{array}[]{clc}\mathcal{S}=&\int_{\mathcal{C}}\textup{d}z\left[\sum\limits_{i=1}^{{N}}\bar{c}_{i}(z)(\mathfrak{i}\partial_{z}+\mu)c_{i}(z)+\sum\limits_{\alpha=1}^{{M}}\bar{\psi}_{\alpha}(z)(\mathfrak{i}\partial_{z}+\mu)\psi_{\alpha}(z)\right.-\frac{1}{(2{N})^{3/2}}\sum\limits_{\mathclap{ijkl}}J_{ijkl}\bar{c}_{i}\bar{c}_{j}c_{k}c_{l}-{1\over{M}^{1/2}}\sum\limits_{\mathclap{\alpha\beta}}t_{\alpha\beta}^{\psi}\bar{\psi}_{\alpha}\psi_{\beta}\\ &-\left.{f(z)\over({N}{M})^{1/4}}\sum\limits_{\mathclap{i\alpha_{1}}}V_{i\alpha}\bar{c}_{i}\psi_{\alpha}+V_{i\alpha}^{*}\bar{\psi}_{\alpha}c_{i}\right],\end{array}$$

$$⟨{\mathcal{Z}}_{\text{neq}}{⟩}_{dis}$$

$$=\int \mathcal{D}(c,\bar{c},\bar{\psi },\psi )\int d[J,V,t]e^{\mathfrak{i}\mathcal{S}}P[J]P[V]P[t],$$

$$G_c(z_1,z_2)$$

$$G_{\psi }(z_1,z_2)$$

$$\mathcal{S}[\Sigma ,G]=$$

$$\displaystyle\int_{\mathcal{C}}\mathrm{d}z_{1}\mathrm{d}z_{2}\Big{[}+\mathfrak{i}\frac{J^{2}{N}}{4}{G_{c}}(z_{2},z_{1})^{2}{G_{c}}(z_{1},z_{2})^{2}+\mathfrak{i}\frac{t_{\psi}^{2}{M}f(z_{1})f(z_{2})}{2}{G_{\psi}}(z_{2},z_{1}){G_{\psi}}(z_{1},z_{2})$$

$$+\mathfrak{i}{V}^2\sqrt{NM}f(z_1)f(z_2)G_c(z_2,z_1)G_{\psi }(z_1,z_2)+{\Sigma }_c(z_1,z_2)[-\mathfrak{i}N{G}_c(z_2,z_1)]$$

$$\displaystyle+{\Sigma_{\psi}}(z_{1},z_{2})[-\mathfrak{i}{M}{G_{\psi}}(z_{2},z_{1})]\Big{]},$$

$${\Sigma }_c(z_1,z_2)$$

$$+\sqrt{p}{V}^{2}f(z_1)f(z_2)G_{\psi }(z_1,z_2)$$

$${\Sigma }_{\psi }(z_1,z_2)$$

$$(\mathfrak{i}{\partial }_{z_1}+\mu )\delta (z_1-z_2)-{\Sigma }_s(z_1,z_2)=$$

$$(\mathfrak{i}{\partial }_{z_1}+\mu )G_s(z_1,z_2)-$$

$$={\delta }_{\mathcal{C}}(z_1,z_2),$$

$$\left(-\mathfrak{i}{\partial }_{z_2}+\mu \right)G_s(z_1,z_2)-$$

$$={\delta }_{\mathcal{C}}(z_1,z_2),$$

$$G_c^{>}(t_1,t_2)$$

$$G_c^{<}(t_1,t_2)$$

$$G^R(t_1,t_2)$$

$$=\Theta (t_1-t_2)\left[G^{>}(t_1,t_2)-G^{<}(t_1,t_2)\right],$$

$$G^A(t_1,t_2)$$

$$=\Theta (t_2-t_1)\left[G^{<}(t_1,t_2)-G^{>}(t_1,t_2)\right].$$

$$\mathfrak{i}{G}_s^{<}(\omega )/(2\mathrm{Im}{G}_s^R(\omega ))$$

$$G(\omega )={\int }_{-\infty }^{\infty }G(t)e^{\mathfrak{i}\omega t}\text{d}\omega ,$$

$$(\mathfrak{i}{\partial }_{t_1}+\mu )G_s^{>}(t_1,t_2)$$

$$(-\mathfrak{i}{\partial }_{t_2}+\mu )G_s^{>}(t_1,t_2)$$

$$I_{s>}^{(1)}(t_1,t_2)$$

$$+{\int }_{-\infty }^{t_2}{\Sigma }_s^{>}(t_1,t)G_s^A(t,t_2)\text{d}t$$

$$I_{s>}^{(2)}(t_1,t_2)$$

$$+{\int }_{-\infty }^{t_2}{G}_s^{>}(t_1,t){\Sigma }_s^A(t,t_2)\text{d}t,$$

$${\Sigma }^R(t_1,t_2)$$

$${\Sigma }^A(t_1,t_2)$$

$${\Sigma }_c^{>,<}(t_1,t_2)={\Sigma }_c(z_1\to t_{1-},z_2\to t_{2+})$$

$$=J^2{G}_c^{>,<}(t_1,t_2{)}^2{G}_c^{<,>}(t_2,t_1)+\sqrt{p}{V}^{2}f(t_1)f(t_2)G_{\psi }^{>,<}(t_1,t_2)$$

$${\Sigma }_{\psi }^{>,<}(t_1,t_2)={\Sigma }_{\psi }(z_1\to t_{1-},z_2\to t_{2+})=t_{\psi }^2{G}_{\psi }^{>,<}(t_1,t_2)$$

$$+\frac{V^2}{\sqrt{p}}f(t_1)f(t_2)G_c^{>,<}(t_1,t_2).$$

$$(\mathfrak{i}{\partial }_{t_1}+\mu )G_s^{>,<}(t_1,t_2)$$

$$+{\int }_{-\infty }^{t_2}{\Sigma }_s^{>,<}(t_1,t)G_s^A(t,t_2)\text{d}t$$

Full text

I Introduction

One of the major frontiers in condensed matter physics is to describe gapless phases of interacting fermions without any quasiparticles, namely non Fermi liquids (NFL) Sachdev (2011). Recently, new insights about fundamental differences between NFLs and Fermi liquids (FL) have been gained in terms of many-body quantum chaos and thermalization. This new impetus has come from exciting developments in a class of NFLs described by Sachdev-Ye-Kitaev (SYK) model, Sachdev and Ye (1993); Kitaev ; Maldacena and Stanford (2016) and its extensions Gu et al. (2017); Banerjee and Altman (2017); Jian and Yao (2017); Song et al. (2017); Davison et al. (2017); Patel et al. (2018); Chowdhury et al. (2018); Haldar et al. (2018); Haldar and Shenoy (2018), and their connections with black holes in quantum gravity Sachdev (2010); Kitaev ; Sachdev (2015). In particular, the model proposed in ref.Banerjee and Altman (2017) classifies the SYK NFL and a FL as two distinct chaotic fixed points, separated by a quantum phase transition (QPT). In this characterization, the NFL thermalizes much faster than the FL, as quantified by a rate of the onset of chaos or the Lyapunov exponent Maldacena et al. (2016); Kitaev ; Banerjee and Altman (2017).

However, the Lyapunov exponent is computed from an equilibrium dynamical correlation, the so-called out-of-time-ordered correlator Kitaev ; Maldacena and Stanford (2016); Kitaev and Suh (2018). Here, using the model of ref.Banerjee and Altman (2017) as a template, we ask whether such contrast between the NFL and FL persists even for thermalization from a completely out-of-equilibrium situation, e.g. a quantum quench. The model has two species of fermions, interacting SYK fermions coupled to another species of otherwise non-interacting fermions, referred to as lead fermions. A QPT between a strongly interacting NFL and weakly interacting FL phases can be tuned in the model at low energies by varying the ratio, $p$, of numbers of sites on which the two types of fermions reside. Remarkably, the solvable nature of the model allows us to study its full non-equilibrium evolution after a quench exactly. By using non-equilibrium Keldysh field theory in the thermodynamic limit, as well as numerical exact diagonalization (ED) for finite systems, we demonstrate a drastic difference in thermalization rates for the NFL and FL after a sudden quench. In addition, we show that the quasiparticle residue of the lead fermions exhibits a dynamical transition as function of $p$ from collapse-and-revival oscillations to prethermalized plateaus as a function of time. This dynamical transition is similar to that seen in the interaction quench of Hubbard model Eckstein et al. (2009).

Furthermore, the Landau description of a FL is based on the concept of adiabatic time evolution from a non-interacting system under slow switching on of the interaction, without encountering a phase transition. Is it possible to evolve an NFL adiabatically to the FL and vice versa? We argue that such evolution is not possible here due to another intriguing aspect of the SYK NFL, namely the finite zero-temperature residual entropy (density) $S_{0}$ Kitaev ; Sachdev (2015); Maldacena and Stanford (2016). The entropy is related to the Bekenstein-Hawking entropy of the black hole in the dual gravity theory Kitaev ; Sachdev (2015); Maldacena and Stanford (2016); Kitaev and Suh (2018) and has relevance for strange metallic states described by local quantum criticality Gu et al. (2017); Davison et al. (2017); Chowdhury et al. (2018); Haldar and Shenoy (2018). We probe the signature of this entropy in the heat generated during non-equilibrium dynamics and characterize how the putative adiabatic limit is approached in the two phases, and across the QPT, after a slow quench with a finite rate. We show that, remarkably, the heat or energy of excitations $\Delta E$ generated by the quench scales as $\Delta E(\tau)\sim\tau^{-\eta}$ with quench time $\tau$. Moreover, we find a direct manifestation of the equilibrium QPT in the scaling exponent $\eta$. We also contrast all the above results for the interacting model with those for an analogous non-interacting model under identical quench protocols. In particular, we show that the two models show drastically different thermalization behaviors in the thermodynamic limit.

The model studied here, could be realized in a graphene nano flake Can et al. (2019) attached to leads, and under a magnetic field. The QPT also has close parallel in the NFL to FL transition in the multichannel Kondo model Parcollet and Georges (1997). Moreover, the study of dynamics after a quench in our model, where no quasiparticle description exists in one of the phases around the QPT, allows us to probe hitherto unexplored regime of many-body quantum dynamics. This is complementary to the previous studies of dynamics after quench across a QPT in integrable models Essler and Fagotti (2016); D’Alessio et al. (2016) or, weakly interacting systems with well defined quasiparticles Moeckel and Kehrein (2008, 2009); Eckstein et al. (2009). The scaling laws mentioned above can not be explained by the usual Kibble-Zurek scaling Dziarmaga (2010); Polkovnikov et al. (2011), unlike that in integrable or weakly-integrable models Essler and Fagotti (2016); D’Alessio et al. (2016); Moeckel and Kehrein (2008, 2009); Eckstein et al. (2009). Although, there have been a few studies on non-equilibrium dynamics of the original SYK model Eberlein et al. (2017); Sonner and Vielma (2017); Kourkoulou and Maldacena (2017); Bhattacharya et al. (2018), none of them addressed the issue of a quench across a non-trivial QPT.

The remainder of this paper is organized as follows. In sec. II, we introduce the interacting and non-interacting models and describe the quench protocols. Section III discusses the results for the non-equilibrium evolutions after slow and sudden quenches in the large-$N$ limit obtained using non-equilibrium Schwinger-Keldysh method. Some results for finite-$N$ obtained via ED studies in a few limiting cases are also discussed in this section. In sec. IV we conclude with the implications and significance of our results. The details of the non-equilibrium large-$N$ formulations, equilibrium spectral properties of the models and some additional results on the slow quench in the non-interacting model are included in Appendices A, B and C. The analysis of the slow quench using usual adiabatic perturbation theory and the breakdown of adiabaticity in our model are discussed in Appendix D. Additional details of the numerical calculations, ED and the results are given in the Supplementary Material (SM) SM .

II Model

II.1 Interacting model

As described schematically in fig.1, we study a time ($t$)-dependent version of the model in ref.Banerjee and Altman (2017), $\mathcal{H}(t)=\mathcal{H}_{c}+\mathcal{H}_{\psi}+\mathcal{H}_{c\psi}(t)$, where

[TABLE]

The model (fig. 1), has two species of fermions – (1) the SYK fermions ($c$), on sites $i=1,\dots,N$, interacting with random four-fermion coupling $J_{ijkl}$ (eqn. (1a)), drawn from a Gaussian distribution with zero mean and variance $\overline{|J_{ijkl}|^{2}}=J^{2}$; and (2) the lead fermions ($\psi$), on a separate set of sites $\alpha=1,\dots,M$ connected via random all-to-all hopping $t^{\psi}_{{\alpha}{\beta}}$ (eqn. (1b)). The SYK and the lead fermions are quadratically coupled via $V_{i{\alpha}}$; $t^{\psi}_{{\alpha}{\beta}}$ and $V_{i{\alpha}}$ are complex Gaussian random variables with zero mean, and variances $\overline{|t^{\psi}_{{\alpha}{\beta}}|^{2}}=t_{\psi}^{2}$ and $\overline{|V_{i{\alpha}}|^{2}}=V^{2}$, respectively.

The model is exactly solvable for $N,M\to\infty$ with a fixed ratio $p=M/N$, that is varied to go through the QPT between NFL and FL at a critical value $p=p_{c}=1$ Banerjee and Altman (2017). Two crossover scales, $\omega_{NFL}$ and $\omega_{FL}$, approach zero from either sides of the QPT (fig. 1). The residual entropy density $S_{0}(p)$ of the SYK NFL continuously vanishes at the transition Banerjee and Altman (2017). This is one of the unique features of the QPT.

To probe the non-equilibrium dynamics, we make the coupling term in eqn. (1c) time dependent. In particular, we perform geometric quenches (fig.1) by switching on the coupling between the two initially disconnected subsystems (a) suddenly such that $f(t)=\Theta(t)$, the Heaviside step function, and (b) by slowly ramping up the coupling over a time $\tau$, i.e. $f(t)=r(t/\tau)[\Theta(t)-\Theta(t-\tau)]$; $r(x)$ is a ramp function, e.g. $r(x)=x$. Before the quench, the disconnected subsystems, with a preset site-ratio $p$, are at their own thermal equilibria at initial temperatures $T_{i}^{c}$ and $T_{i}^{\psi}$. We take $T_{i}^{c},~{}T_{i}^{\psi}\to 0$ so that SYK and lead fermions belong to the NFL and (non-interacting) FL states, respectively. As shown in fig. 1, for $t\to\infty$, depending on whether $p<1$ or $p>1$, the coupled system eventually is expected to thermalize to either the NFL or the FL state, respectively. In any case, one of the subsystems always undergoes a transition, either from FL to NFL or vice versa, under the quench.

II.2 Non-interacting model

To contrast the behavior of the above interacting model, and to demonstrate the crucial role of interaction in the non-equilibrium dynamics after quench and eventual thermalization, we also consider an analogous non-interacting model. The latter is obtained by replacing the interaction term for the $c$ fermions with a random hopping term similar to the one appearing in the $\psi$-fermion Hamiltonian [eqn. (1)]. To this end, we have $\mathcal{H}^{(NI)}(t)=\mathcal{H}_{c}+\mathcal{H}_{\psi}+\mathcal{H}_{c\psi}(t)$, with

[TABLE]

Here $t^{c}_{ij}$ is a complex Gaussian random variable with zero mean, and variance $\overline{|t^{c}_{ij}|^{2}}=t_{c}^{2}$; $\mathcal{H}_{\psi}$ and $\mathcal{H}_{c\psi}(t)$ are same as in eqn. (1).

Below, we first briefly discuss the method for studying the time evolution of the systems, followed by the results for sudden and slow quenches.

III Results

Non-equilibrium evolution: We use standard Schwinger-Keldysh non-equilibrium Green’s function technique Kamenev (2011); Stefanucci and van Leeuwen (2013) to study the quenches described above. Utilizing the closed-time-contour Schwinger-Keldysh action (see Appendix A) for the model we derive the Kadanoff-Baym (KB) equations for the disorder-averaged non-equilibrium Green’s functions, $G_{s}^{<}(t_{1},t_{2})$, $G_{c}^{>}(t_{1},t_{2})$ ($s=c,\psi$), e.g. $G_{c}^{>}(t_{1},t_{2})=-i\overline{\langle c_{i}(t_{1})c_{i}^{\dagger}(t_{2})\rangle}$; the overline denotes disorder averaging. The KB equations are numerically integrated using a predicator-corrector scheme (see SM , S1) starting from the initial equilibrium Green’s functions (see SM , S2) for the disconnected system. The time-dependence of $\mathcal{H}(t)$ is encoded in KB equations via the local self-energies $\Sigma_{s}$, which could be exactly calculated in the large-$N$ limit for both the interacting and the non-interacting models (see Appendix A).

We obtain the time-dependent expectation value of an observable $\mathcal{O}(t)$, i.e. $\langle\mathcal{O}(t)\rangle\equiv\mathrm{Tr}[\rho(t)\mathcal{O}(t)]$ (see SM , S3), using the Green’s functions. Here $\rho(t)$ is the time-dependent density matrix and $\mathcal{O}(t)$ includes the explicit time dependence, if any, of the observables. To understand thermalization, we track how the contributions of the individual terms in eqn. (1) to the total energy $E(t)=\langle\mathcal{H}(t)\rangle$, e.g. $E_{c\psi}(t)=\langle\mathcal{H}_{c\psi}(t)\rangle$, relax after the quench. Since the whole system is isolated, we estimate the expected temperature $T_{f}$ of the putative thermal state at long times from the total energy $E(t)=E_{f}$, which is conserved after the quench.

III.1 Sudden quench

We first ask whether the contrast between dynamics of the NFL and FL can be seen even when the system is subjected to an abrupt non-equilibrium process. To address this, we study the case when the sub-systems are suddenly connected at $t=0$. We take $J=1$, $t_{\psi}=1$, $V=1$ and low initial temperatures, $T_{i}^{c}=0.05$ and $T_{i}^{\psi}=0$ 111We keep the initial temperature of the SYK subsystem low but finite. Because of the divergent spectral density at $T=0$, the Green’s function for SYK fermions in the initial equilibrium state can not be obtained numerically at strictly zero temperature.. The sudden quench leads to a rather high final temperature $T_{f}\sim 1$ (see SM , S4.1). Before the quench, the lead fermions are non-interacting and the single-fermionic excitations are sharply defined at $T_{i}^{\psi}=0$. To track the quasiparticle evolution, we compute an energy resolved time-dependent occupation, $n_{\psi}(\epsilon,t)=-iG_{\psi}^{<}(\epsilon;t,t)$, for the lead fermions. Here $\epsilon$ are the eigenvalues of the quadratic Hamiltonian in eqn. (1b), i.e. $\mathcal{H}_{\psi}=\sum_{\epsilon}\epsilon\psi_{\epsilon}^{\dagger}\psi_{\epsilon}$ and the Green’s function $G_{\psi}^{<}(\epsilon;t_{1},t_{2})$ is obtained by integrating an appropriate KB equation (see eqn. (A.1)). The quasiparticle residue, $z_{\psi}(t)=n_{\psi}(0^{-},t)-n_{\psi}(0^{+},t)$, is obtained from the occupation discontinuity at $\epsilon=0$. The vanishing of the residue indicates the destruction of the quasiparticles.

Collapse-and-revival oscillations and prethermal plateaus: In our model, $z_{\psi}^{\infty}=z_{\psi}(t\to\infty)$ is expected to vanish for quench to any $p$ since the coupled system either thermalize to NFL or to a finite temperature FL state . As shown in figs. 2(a)-(b), the collapse of the residue happens through two very different routes. First, for $p<p_{c}^{dyn}$, a critical value of $p$ corresponding to a dynamical transition, $z_{\psi}(t)$ undergoes collapse-and-revival oscillations (fig. 2(a)) Second, for $p>p_{c}^{dyn}$, $z_{\psi}(t)$ shows multiple long prethermal plateaus (fig. 2(b)). For $V=1$, we find $p_{c}^{dyn}\approx 1.5$ from $z_{1}(p)\to 0$, where $z_{1}$ is the residue at the first maximum of the oscillations (see curve labeled int. in fig. 2(a) inset). Hence, the critical value $p_{c}^{dyn}$ for this dynamical transition is greater than the ‘equilibrium’ critical ratio $p_{c}=1$. It is encouraging to find that similar oscillations and evidence of a dynamical transition have been observed Eckstein et al. (2009) in the interaction quench across Mott transition in the Hubbard model as well.

However, in contrast to the interaction quench in Hubbard model, where the collapse-and-revival oscillations originate from on-site Hubbard repulsion, in our model, the oscillations are linked to a soft hybridization gap in the lead fermions. The gap appears due to hybridization (eqn. (1c)) between the SYK and lead fermions, and closes at the NFL-FL transition (Appendix B). In fact, we observe a similar dynamical transition in the non-interacting model of eqn. (2) under an analogous sudden quench (see curve labeled n.int. in fig. 2(a) inset). The critical point in this scenario occurs at $p_{c}^{dyn}\approx 0.5$. We emphasize here that although a dynamical transition exists for the non-interacting case, the presence of interactions is crucial for thermalization (discussed in next section), which the non-interacting model fails to achieve for any value of $p$ (see SM , S4.3). The SYK-interactions are also responsible for non-trivially shifting the critical value of $p_{c}^{dyn}$ from $0.5$ to $1.5$ as shown in fig. 2(a) (inset).

Thermalization and long-time steady states: The crucial aspects that distinguish the non-interacting (eqn. (2)) and interacting (eqn. (1)) models, as well as the NFL and FL, are the thermalization process and the long-time steady states. As shown in fig. 2(c) (also see SM , S4.2, fig. S3(a)-(b)), $E_{c\psi}(t)$ reaches the thermal expectation corresponding to the temperature $T_{f}$ very rapidly for $p<p_{c}^{dyn}$, whereas there is a drastically slow, albeit finite, relaxation rate for $E_{c\psi}(t)$ towards the thermal value for $p>p_{c}^{dyn}$. In contrast, $E_{c\psi}(t)$ does not relax to the expected thermal value for any $p$ in the non-interacting model, see fig. 2(g) (also SM , S4.2, fig. S4(a)-(b)). We further analyze the steady state through the Green’s functions $G_{s}(\mathcal{T},\omega)=\int_{-\infty}^{\infty}G_{s}(t_{1}=\mathcal{T}+t/2,t_{2}=\mathcal{T}-t/2)e^{i\omega t}$, where $\mathcal{T}=(t_{1}+t_{2})/2$. In the steady state, $G_{s}(\mathcal{T},\omega)$ becomes independent of $\mathcal{T}$. Moreover, for a thermal steady state at $T_{f}$, the steady-state occupation function $f_{s}^{\infty}(\omega)=\lim_{\mathcal{T}\to\infty}\mathfrak{i}G_{s}^{<}(\mathcal{T},\omega)/(2\mathrm{Im}G_{s}^{R}(\mathcal{T},\omega))$ should be equal to the Fermi function $n_{\mathrm{F}}(\omega,T_{f})=1/(e^{\omega/T_{f}}+1)$, i.e. should satisfy the fluctuation-dissipation theorem (FDT). We find that for the interacting model (eqn. (1)), $f_{s}^{\infty}(\omega)=n_{\mathrm{F}}(\omega,T_{f})$ (see fig. 2(e) top) for the NFL ($p<p_{c}^{dyn}$), whereas FDT gets violated (fig. 2(e) bottom) in the FL regime ($p>p_{c}^{dyn}$) for the largest $\mathcal{T}$ ($\sim 100J^{-1}$, where $J$ is the interaction strength in eqn. (1a)) accessed.

We find that the FDT is never satisfied for the non-interacting model at any $p$, as seen in fig. 2(f) (also see SM , S4.2, fig. S4(c)-(d)). This is expected for the non-interacting case, where the long-time steady state is described by a generalized Gibbs ensemble (GGE) instead of the usual thermal Gibbs ensemble Polkovnikov et al. (2011); Essler and Fagotti (2016); D’Alessio et al. (2016); SM . Nevertheless, even in the interacting model, the FL phase shows an approximate GGE behavior by attaining a prethermal steady state within the time accessible in our numerical calculations. The prethermal GGE will presumably relax to a thermal state over a much longer time scale D’Alessio et al. (2016); Stark and Kollar (2013). Similar behaviors have been seen for quenches to FL phases in other interacting models Manmana et al. (2007). In contrast, the strong interaction leads to rapid thermalization for the NFL phase. Hence the sudden quench in the interacting model demonstrates drastically different thermalization dynamics between the NFL and the FL phases.

It is worthwhile to ask whether the contrast in the thermalization behaviors of NFL and FL persists even at finite $N$. In fig. 2(h), we show (see also SM , S4.2.1) the results for $E_{c\psi}(t)$ obtained from ED studies of the model of eqn. (1) for $N=16$. The ED gives results similar to that at large-$N$, shown in fig. 2(g). Another pertinent question is whether the thermalization times of $E_{c\psi}(t)$ in NFL and FL phase can be directly related to their respective Lyapunov time scales ($\tau_{\mathrm{L}}$) Banerjee and Altman (2017). Here it is important to note that the sudden quench in our model leads to (see SM , S4.1, fig. S2) substantially high temperature $T_{f}\gtrsim J$. As a result, the relaxation of various high-energy modes also influence $E_{c\psi}(t)$, making it hard to isolate $\tau_{\mathrm{L}}$ from the relaxation of high-energy modes.

We would also like to note that even though some specific low-energy features of SYK NFL and FL phases, like the temperature dependence of the Lyapunov time $\tau_{\mathrm{L}}$, cannot be ascertained from the dynamics after sudden quench, our results show that the low-energy fixed points for the initial states still drastically influence the thermalization process. This is despite the fact that the sudden quench leads to substantially high final temperatures.

III.2 Slow quench

We next address the question whether the initial decoupled NFL and FL subsystems can be adibatically evolved to the final states of the coupled system. To this end, we consider the slow quench where the coupling is changed slowly through a ramp, $r(t/\tau)$. Here we keep both the subsystems at some low initial temperature $T_{i}^{c}=T_{i}^{\psi}=T_{i}$ for $t<0$ and define the heat or excitation energy Polchinski and Rosenhaus (2016); Eckstein and Kollar (2010), $\Delta E(\tau)=E(\tau)-\langle\mathcal{H}(\tau)\rangle_{T_{i}}$, produced during the quench; $\langle\mathcal{H}(\tau)\rangle_{T_{i}}$ is the thermal expectation of the final Hamiltonian $\mathcal{H}(\tau)$ at the initial temperature $T_{i}$. As shown in fig. 3(a), remarkably, we find that $\Delta E(\tau)\sim\tau^{-\eta(p)}$ with $\eta(p)<1$, i.e. a non-analytic power-law scaling. The exponent $\eta$ has a strong non-monotonic dependence on $p$ with a minimum around the QPT (fig. 3(b)), revealing signatures of equilibrium QPT in the non-equilibrium evolution. We also find non-analytic scaling for the quench in the non-interacting model (see Appendix C). However, the exponent has a very weak dependence on $p$ (fig. 3(b) inset). The particular non-analytic power laws cannot be explained through a standard adiabatic perturbation theory Polkovnikov (2005); Dziarmaga (2010); De Grandi et al. (2010); Eckstein and Kollar (2010), as we show in Appendix D. Also, a Kibble-Zurek-type argument Dziarmaga (2010); Polkovnikov et al. (2011) cannot be given for such a zero-dimensional system. We find the exponent to depend on ramp shape as well (see SM , S5.1). This is also not expected from adiabatic perturbation theory for an exponent $\eta<1$ Eckstein and Kollar (2010); SM .

One possible promising route to understand the non-standard exponent $\eta(p)$ in the intermediate time window after the quench could be to construct a low-energy theory for the model of eqn.(1) along the line of Schwarzian theory for the pure SYK model Maldacena et al. (2016); Kitaev and Suh (2018). A recent study Almheiri et al. (2019) analyzes the quench dynamics using the Schwarzian theory for a SYK model suddenly coupled to a large thermal bath made out of another SYK model. However, the situation is somewhat more involved in our model due the strong back action of SYK fermions on the lead fermions in the NFL phase ($p<1$) and that of the lead fermions on the SYK fermions in the FL phase ($p>1$). Such back action is absent in the model of ref.Almheiri et al., 2019 since the bath is infinitely larger than the system, and since both the bath and the system are described by SYK model. In our case, one needs to start with two uncoupled low-energy theories, one corresponding to the Schwarzian action for the SYK fermions with scaling dimension $\Delta=1/4$, and the other for the non-interacting fermions with scaling dimension $\Delta=1/2$. Hence the resulting theory after the quench, will not be that of the standard Schwarzian mode, but a different and more complicated one that includes the strong back actions among the subsystems. We would discuss this effective theory elsewhere Banerjee et al. .

As alluded earlier, for the quench to any finite $p$, e.g. from the NFL to FL, the residual entropy $S_{0}$ of the NFL, implies a violation of adiabaticity even for an arbitrary slow quench. The excitation energy $\Delta E(\tau)$ characterizes how $S_{0}$ metamorphoses into thermal excitations in the FL. The latter has $S_{0}=0$, and hence even an arbitrary slow quench must lead to $\Delta E(\tau\to\infty)\neq 0$ and a $T\neq 0$ state, having a thermal entropy that at least accounts for the $T=0$ entropy of the initial NFL state. Hence, the observed powerlaw, implying asymptotic approach to the adiabatic limit $\Delta E(\tau\to\infty)=0$, is surprising. It would suggest that $S_{0}$ is not manifested as thermal excitations in the final state for $\tau\to\infty$. Hence, we do not expect the power law to eventually persist for any finite $p$ for $\tau\to\infty$. We see an indication of only a weak deviation from the intermediate-$\tau$ power law for $p=10$ around $\tau\sim 30-40$ (fig. 3(a)(inset)). From the intermediate-$\tau$ power law, we can estimate a time scale, much longer than presently accessible in our calculations, where the scaling is expected to be violated due to $S_{0}$ (see Appendix D). This mechanism of the violation of adiabaticity due to $S_{0}$ in the large-$N$ limit is different from the hitherto known routes Polkovnikov et al. (2011) of adiabaticity breaking. As we show in Appendix D, physics beyond the large-$N$ limit Bagrets et al. (2016) suggests that the limits $\tau\to\infty$ and $N\to\infty$ do not commute, also indicating the absence of the adiabatic limit Polkovnikov et al. (2011).

As shown in fig. 3(c), a steady-state value of $E_{c\psi}(t)$, consistent with the thermal value is attained very rapidly within the NFL for any $\tau$. In contrast, a ‘glassy’ behavior is seen within the FL, where $E_{c\psi}(t)$ relaxes very slowly for small $\tau$, but relaxes almost instantaneously for larger $\tau$ values (fig. 3(d)).

IV Conclusions

In conclusion, our study of sudden and slow quenches in a large-$N$ model of NFL-FL transition reveal a remarkably rich non-equilibrium phase diagram and sharp contrasts between non-interacting, FL and NFL phases. The sudden quench allows us to track the distinct evolutions of initially prepared well defined quasipartcile state in the NFL and FL phases and establish the existence of a dynamical phase transition which is different from the equilibrium NFL-FL quantum phase transition. In the context of slow quenches, unique features of the NFL-FL QPT and the low-temperature state of SYK model, such as strongly interacting fermionic excitations and residual zero-temperature entropy, allows us to probe completely unexplored regime of out-of-equilibrium quantum many-body dynamics compared to previous studies of integrable and weakly-integrable systems. These unusual features lead to remarkable intermediate non-analytic scaling of excitation-energy production with the quench-duration and the eventual breakdown of quantum adiabaticity. A natural future extension would be to go beyond large-$N$ to study evolution for longer times $\sim N$.

Acknowledgements.

We thank Ehud Altman, Subroto Mukerjee, Diptiman Sen, Joel Moore, Vijay B. Shenoy, Emil Uzbashyan, Soumen Bag, Renato Dantas and Paul A. McClarty for useful discussions. SB acknowledges support from The Infosys Foundation, India.

Appendix A Non-equilibrium Green’s functions and Kadanoff-Baym (KB) equations

A.1 Interacting model

We find the non-equilibrium Green’s functions and the corresponding Kadanoff-Baym equations using the Schwinger-Keldysh closed contour formalism. To do this we write the Schwinger-Keldysh actionKamenev (2011) for the time-dependent Hamiltonian ${\cal{H}}$ in eqn. (1), i.e.

[TABLE]

The contour variable $z$ lies on the usual Keldysh contour Kamenev (2011) (fig. 4) with the forward ($+$) or backward ($-$) branches.

The explicit time dependence in the action is introduced via the function $f(z)=f(t)$. The non-equilibrium generating functional for obtaining time-dependent expectation values is defined as $\mathcal{Z}_{\textup{neq}}=\int\mathcal{D}(c,\bar{c},\bar{\psi},\psi)\ e^{\mathfrak{i}\mathcal{S}[\bar{c},c,\bar{\psi},\psi]}$, under the two usual assumptions. First, the initial density-matrix $\hat{\rho}$ is independent of any disorder, and all the disorder dependence has been pushed into the time evolution operators. Second, the disorder is switched on sometime in the infinitely long past so that the system has enough time to equilibrate to the conditions created by the disorder dependent Hamiltonian. These assumptions allow us to implement the averaging of $\mathcal{Z}_{\textup{neq}}$ over all disorder realizations as follows

[TABLE]

where $P[.]$ s are the Gaussian probability distributions for the couplings $J_{ijkl},~{}t_{\alpha\beta}^{\psi}$ and $V_{i\alpha}$ appearing in eqn. (1). We perform the Gaussian integrals over the disorder distributions and define the large-$N$ fields,

[TABLE]

that live on the contour, and the corresponding Lagrange multipliers $\Sigma_{c,\psi}(z_{1},z_{2})$. Finally, after integrating out the fermions we end up with the action

[TABLE]

where the matrix $(\mathfrak{i}\partial_{1}+\mu)\mathbf{1}$ has elements of the form $(\mathfrak{i}\partial_{z_{1}}+\mu)\delta(z_{1}-z_{2})$. The elements for matrices $\mathbf{{\Sigma_{c}}}$, $\mathbf{{\Sigma_{\psi}}}$ are ${\Sigma_{c}}(z_{1},z_{2})$, ${\Sigma_{\psi}}(z_{1},z_{2})$ respectively. The action $\mathcal{S}$ is extremized with respect to $G$ and $\Sigma$ to produce the large-$N$ saddle point equations

[TABLE]

and

[TABLE]

where $G_{s}^{-1}$ is the inverse of the matrix $G_{s}$ with its elements given by $G_{s}(z_{1},z_{2})$ and $s=c,\psi$. We rewrite eqn. (10) by multiplying with $G_{s}$ from the right and the left, respectively,

[TABLE]

which are integro-differential equations satisfied by $G_{s}(z_{1},z_{2})$, and where $\delta_{c}$ is the Dirac-delta function defined on the contour.

From the contour-ordered Green’s function $G_{s}(z_{1},z_{2})$, we obtain the disorder averaged real-time non-equilibrium Green’s functions, greater ($>$), lesser ($<$), and retarded $R$, from $G_{s}(z_{1},z_{2})$ at the saddle point, e.g.

[TABLE]

and similarly for $G_{\psi}(z_{1},z_{2})$. The first (second) sign in the suffix of $G_{s}(t_{1\pm},t_{2\mp})$ indicates whether the $z_{1}$($z_{2}$) coordinate lies in the forward (+ve) or backward (-ve) branch of the contour.

Steady state: In a steady-state the Green’s functions are invariant under time translational, e.g., $G^{R}(t_{1},t_{2})=G^{R}(t_{1}-t_{2})=G^{R}(t)$.

Thermal equilibrium: For thermal equilibrium at a temperature $T$, in addition to the above steady-state condition, the Green’s functions satisfy the fluctuation dissipation theorem (FDT) Kamenev (2011), e.g.

[TABLE]

where $n_{\mathrm{F}}(\omega,T)=1/(e^{\omega/T}+1)$ is the Fermi function and $G_{<,>,R,A}(\omega)$ are Fourier transforms of $G_{<,>,R,A}(t)$ defined by

[TABLE]

The above conditions allow us to test whether a system, after undergoing a non-equilibrium process, has reached a steady state or not, and whether the steady state is consistent with thermal equilibrium.

The Kadanoff-Baym Equations From eqn. (11) and eqn. (12), we obtain the time-evolution equations for $G_{s}^{>,<}$, e.g.

[TABLE]

where we have used the fact $\delta_{\mathcal{C}}(z_{1}\to t_{1-},z_{2}\to t_{2+})=0$. Finally, using Langreth rulesLangreth and Wilkins (1972) we get

[TABLE]

where $\Sigma^{R}$ ($\Sigma^{A}$) is the retarded (advanced) self energies, given by

[TABLE]

where, using eqn. (A.1),

[TABLE]

Similarly, one can repeat the above analysis for the $G_{s}^{<}$ case. The integro-differential equations, for $G_{s}^{>}$ and $G_{s}^{<}$ can be stated together as

[TABLE]

This set of equations are called the Kadanoff-Baym (KB) equations, which along with the relations given in eqn. (A.1) and eqn. (A.1), set up a closed system of equations that can be time evolved in the $t_{1}-t_{2}$ plane (see fig. S1), starting from an initial condition for $G$s and $\Sigma$s.

A.2 Non-interacting model

Following procedure similar to that discussed above for the interacting model, we obtain the disorder averaged Schwinger-Keldysh action for the non-interacting model of (2) and the saddle point equations.

The large-$N$ self-energies are given by

[TABLE]

Consequently, a set of contour Kadanoff-Baym equations similar to the ones given in eqn. (11) and eqn. (12) is obtained. Using the Langreth rules the above equations involving contour indices $z_{1}$,$z_{2}$ can be changed to real time variables $t_{1}$, $t_{2}$ to give the final Kadanoff-Baym equations (see eqn. (A.1)) for the system with the following expressions for the self-energies

[TABLE]

Appendix B Gap closing transitions in the interacting model and non-interacting model

The spectral functions obtained from the equilibrium Green’s functions (see SM , S2) for the connected system, i.e. $V=1$, for various $p$-values are discussed below. Fig. 5(a) and (b) show the results for the non-interacting model, whereas fig. 6(a) and (b) show the spectral functions for the interacting model. In the non-interacting case, we find that the spectral function for the $\psi$ fermions have a soft-gap for smaller values of $p$, which closes completely around $p=0.5-0.7$. This is consistent with the dynamical transition that we observe for the sudden quench of the non-interacting model, which we discuss in SM , S4.3. The $c$-fermion spectral function does not have a soft-gap for small values of $p$, but a soft-gap begins to form near $p\sim 1.9$ as we increase $p$. This is expected since, the non-interacting model has an additional symmetry under $p\to 1/p$ and $c\leftrightarrow\psi$. Moving on to the interacting model, we find a similar soft-gap closing scenario taking place for the $\psi$-fermion spectral functions, as shown in fig. 6(b). However, this time the gap closes completely around $p=1.6-2.6$, which is far away from the equilibrium NFL to FL transition point of $p=1.0$. This is again consistent with the dynamical transition critical point (see fig. 2(a) inset in the main text) that we find in the sudden quench of the interacting model. The $c$-fermion spectral functions have a highly peaked form around $\omega=0$, for smaller $p$ values, due to the presence of a divergent $T=0$ spectral function coming from the NFL fixed point. At higher values of $p$, the peak subsides and a gap begins to form in the spectral function.

Appendix C Excitation energy as a function of quench time for the non-interacting model

We also study the dependence of excitation energy $\Delta E$ on the quench duration $\tau$ for the slow quench of the non-interacting model of eqn. (2). We find that there exists a powerlaw relationship between $\Delta E$ and $\tau$ here as well, as shown in fig. 7, which we report in the main text. However, the dependence of the powerlaw exponent $\eta$ on $p$ is qualitatively different from the interacting case as shown in fig. 3(b) of the main text.

Appendix D Breakdown of adiabatic perturbation theory and the absence of adibatic limit

Here we elaborate on the connection between the zero-temperature residual entropy of the SYK NFL and the absence of the adiabatic limit for the slow quenches described in the main paper. In the large-$N$ limit, $S_{0}$ reflects the exponentially dense many-body energy spectrum near the ground state in the NFL phase Maldacena and Stanford (2016), and the QPT in our model marks a transition from exponentially small many-body level spacing, $\Delta\sim e^{-S_{0}(p)N}$, in the NFL to $\Delta\sim 1/N$ in the FL. Hence, the residual entropy $S_{0}$ cannot be thought of as a thermodynamic entropy strictly at $T=0$, i.e., when the $T\to 0$ limit is taken first, keeping $N$ finite and then $N\to\infty$ limit is taken, $S_{0}=0$. In the large-$N$ description, the limit is taken the other way around, and it captures the exponentially dense many-body level spectrum near the ground state in the NFL phase. However, at any non-zero temperature (\ {\raise-2.15277pt\hbox{\buildrel>\over{\sim}}}\ e^{-S_{0}N}), which could be infinitesimally small in the large-$N$ limit, $S_{0}$ is the true thermodynamic entropy. Hence, we expect this entropy to be manifested in the large-$N$ non-equilibrium dynamics during a slow quench, implying the absence of the adiabatic limit. However, as mentioned in the main paper, surprisingly we find the intermediate-$\tau$ non-analytic power law scaling, that seems to mask the effect of residual entropy. In the following, we first try to estimate the power law from the so-called adiabatic perturbation theory Polkovnikov (2005); Dziarmaga (2010); De Grandi et al. (2010); Eckstein and Kollar (2010), that has been previously used for non-interacting and weakly-interacting systems. We show that the adibatic perturbation theory cannot explain the $p$-dependent exponent directly obtained from the direct non-equilibrium evolution, discussed in the main text. Furthermore, we discuss the possible modes of violation of the adiabaticity in the large-$N$ limit, as well as, beyond it.

D.0.1 Adiabatic perturbation theory in the large-$N$ limit

In the adiabatic perturbation theory Polkovnikov (2005); Dziarmaga (2010); De Grandi et al. (2010); Eckstein and Kollar (2010), the time-dependent term (eqn.(1c)) is treated as a perturbation assuming a weak strength of the ramp in eqn.(1), i.e. $f(t)=\Delta fr(t/\tau)$ with $\Delta f\ll 1$. To this end, we obtain the energy generated during quench Eckstein and Kollar (2010) as,

[TABLE]

Here $R(x)=|\int_{0}^{1}dsr^{\prime}(s)e^{ixs}|^{2}$, $r^{\prime}(x)$ is derivative of the ramp function. $\mathcal{A}_{\mathcal{V}}(\omega)$ is a $T=0$ spectral function corresponding to the disordered-averaged imaginary-time correlation function of the operator $\mathcal{V}=(NM)^{-1/4}\sum_{i\alpha}(V_{i\alpha}c_{i}^{\dagger}\psi_{\alpha}+\mathrm{h.c.})$, i.e. $-\overline{\langle\mathcal{T}_{\tau}\mathcal{V}(\tau)\mathcal{V}(0)\rangle}$. This can be computed in the large-$N$ limit, and we obtain,

[TABLE]

where $\rho_{c}(\omega)$ and $\rho_{\psi}(\omega)$ are the spectral functions of the SYK and lead fermions, respectively, for the uncoupled systems before the quench. For long quench time $\tau\gg\Omega^{-1}$, where $\Omega\approx J,t_{\psi}$, we can use the low-energy forms, $\rho_{c}(\omega)\sim|\omega|^{-1/2}$ and $\rho_{\psi}(\omega)\sim$ constant, for $\omega\ll\Omega$. These give $\mathcal{A}_{\mathcal{V}}(\omega)\sim\omega^{1/2}$. It can be shown Eckstein and Kollar (2010), $R(x\to\infty)\sim x^{-2n}$, where $n\geq 1$ is the order of the derivative discontinuity in the ramp (see SM , S5.1). As a result,

[TABLE]

Moreover, the integral above is convergent for $\tau\to\infty$ since $2n>1/2$. Hence, the adiabatic perturbation theory predicts a non-analytic power law with exponent $\eta=1/2$ for any $p$. This, of course, does not agree with the strongly non-monotonic $\eta(p)$ obtained from the direct non-equilibrium calculations (fig.3(b)). Hence, the adiabatic perturbation theory does not work in our case. The theory assumes non-degenerate ground state Eckstein and Kollar (2010), and it is an interesting question whether such theory could at all be applied to a phase with exponentially dense many-body spectrum near ground state, even though some of the effect of the dense spectrum is incorporated in the large-$N$ single-particle density of states.

D.0.2 Breakdown of adiabaticity in the large-$N$ limit

As discussed in the main-text, the presence of a residual entropy in the SYK fermions, prior to the quench, implies that a finite amount of excitation energy $\Delta E>0$ must be produced even in the $\tau\to\infty$ limit. Therefore, the powerlaw relationship $\Delta E\sim\tau^{-\eta}$, in principle, should break down at some large $\tau$. We now provide a route to find an estimate for the time $\tau_{\textup{break}}$ at which we can expect the powerlaw behavior to break down. This can be done by assuming an iso-entropic (the initial and final entropies are taken to be equal) limit of the quench process. We first calculate the temperature $T_{f}$ necessary for the final Hamiltonian $H_{f}$ to hold the initial entropy $S_{i}$ by solving the equilibrium problem and demanding

[TABLE]

Using the definition of $\Delta E$ given in the main text we can estimate $\tau_{\textup{break}}$ from the condition, $\Delta E(\tau_{\mathrm{break}})\approx\Delta E_{T_{f}}=\langle\mathcal{H}(\tau_{\mathrm{break}})\rangle_{T_{f}}-\langle\mathcal{H}(\tau_{\mathrm{break}})\rangle_{T_{i}}$. The latter is the excitation energy produced by converting the initial entropy to thermal excitations, and the excitations generated by the quench cannot be lower than $\Delta E_{T_{f}}$. This leads to,

[TABLE]

where $\alpha$ can be extracted from powerlaw fits to $\Delta E$ using

[TABLE]

Performing such an estimate for $p=1.5$ case, we find $\tau_{\textup{break}}\sim 200$, a time which is not easily accessible using the numerical algorithm of SM , S1.

D.0.3 The adiabatic perturbation theory beyond large-$N$ and the absence of adiabatic limit

Within the adiabatic perturbation theory discussed above, we can go beyond the large-$N$ theory by incorporating some finite-$N$ corrections, at least due to the single-particle level spacing, following ref.Bagrets et al., 2016. It is still not known how to incorporate the effects of many-body level spacing. Nevertheless, it has been shown in ref.Bagrets et al., 2016 that the SYK spectral function $\rho_{c}(\omega)$ changes from the divergent $|\omega|^{-1/2}$ behavior to $\rho_{c}(\omega)\sim\Delta_{s}^{-1}|\omega|^{1/2}$ for $\omega\ll\Delta_{s}\sim J/(N\ln{N})$. The large prefactor $N\ln{N}$ in the $\sqrt{|\omega|}$ dependence presumably arises from the dense spectrum, even though the density of states is suppressed at low energies. We obtain $\mathcal{A}_{\mathcal{V}}(\omega)\sim\omega^{3/2}$ for $\omega\ll\Delta_{s}$, giving

[TABLE]

Hence, keeping $N$ fixed, we obtain $\Delta E(\tau)\sim\tau^{-3/2}N\ln{N}$ for $\tau\to\infty$. Clearly, the limits $\tau\to\infty$ and $N\to\infty$ do not commute. This indicates that the adiabatic limit can not be reached.

S1 Numerical integration of KB equations: Predictor corrector algorithm

The time propagation of KB equations (see eqn. (A.1)) is done numerically on the $t_{1}$-$t_{2}$ (see fig. S1) time plane using a predictor-corrector algorithm. The algorithm works by guessing a value for the Green’s function ($G^{<}(t_{1},t_{2})$, $G^{<}(t_{2},t_{1})$) at a new time step by using a predictor scheme and then uses the said value along with past values to correct the guess by using a corrector scheme. The initial conditions for the Green’s functions are encoded in quadrant $C$ as shown in fig. S1 using the equilibrium solutions for the functions obtained by using the initial Hamiltonian from before the quench. We have discussed the details of constructing the equilibrium Green’s functions in sec. S2. Before we discuss the rest of the algorithm, we define a new center of time variable $\mathcal{T}$ and a relative time variable $\mathbf{t}$ as

[TABLE]

Since, $\frac{\partial}{\partial\mathcal{T}}=\frac{\partial}{\partial t_{1}}\frac{\partial t_{1}}{\partial\mathcal{T}}+\frac{\partial}{\partial t_{2}}\frac{\partial t_{2}}{\partial\mathcal{T}}=\frac{\partial}{\partial t_{1}}+\frac{\partial}{\partial t_{2}},$ we derive yet another differential equation in the variable $\mathcal{T}$ by subtracting the two equations in eqn. (A.1), and get the following

[TABLE]

which we shall use later on to time step the Green’s functions when $t_{1}=t_{2}$. We now discuss the two parts of the algorithm– part a) is the predictor scheme, which we implement using Euler discretization, and part b) a corrector scheme which will be implemented implicitly.

S1.1 Predictor scheme

We define $F_{1,2}^{>,<}$, such that

[TABLE]

We discretize the time domain and represent the discretized times with indexed symbols $t_{n},t_{m},t_{i}$ etc., see fig. S1, also $\Delta t$ is the step size of our time domain and $-\infty$ ($\infty$) is a large negative (positive) value set by $-T_{max}$ ($T_{max}$). Using the said definitions the $F_{1,2}$ terms can be approximated as

[TABLE]

The Green’s functions at $(t_{n+1}=t_{n}+\Delta t,t_{m})$ and $(t_{n},t_{m+1})$ can then be predicated using

[TABLE]

in the directions indicated by the blue and red arrows in fig. S1. Suppose we have all the information, i.e. $G^{>,<},\Sigma^{>,<}$, in the grid $[-\infty,t_{k}]\times[-\infty,t_{k}]$ and want to extend it to $t_{k+1}=t_{k}+\Delta t$, then $F_{1,2}^{>,<}$, for $n,m\leq k$, can be easily calculated since the information needed to evaluate $F_{1,2}^{>,<}$ is readily available. Although the equations in eqn. (S1.1) are perfectly valid for obtaining a prediction for both $G^{>}$, $G^{<}$ along the horizontal and vertical lines in fig. S1, they are computationally expensive. We can cut down on the cost by using the following property of the nonequilibrium Green’s functions

[TABLE]

To this end, we time step one of the Green’s functions, for e.g. $G^{>}$, along the vertical line in the direction of the blue arrows shown in fig. S1, i.e. $(t_{n},t_{m})\to(t_{n+1},t_{m})$, and the lesser Green’s function $G^{<}$ along the horizontal line in the direction of the red arrows, i.e. $(t_{n},t_{m})\to(t_{n},t_{m+1})$, and then use eqn. (S1.1) to find $G^{>}$ ($G^{<}$) on the horizontal (vertical) line as follows

[TABLE]

Note that we have not yet obtained a prediction for the diagonal point $(t_{k+1},t_{k+1})$, i.e, $G^{>,<}(t_{k+1},t_{k+1})$. This can be done using eqn. (LABEL:eq:KB_eq_diag_1_app), which gives us the prediction formula

[TABLE]

S1.2 Corrector scheme

In order to correct the values for the Green’s function at $(t_{n=k+1},t_{-\infty\leq m\leq k})$, $(t_{-\infty\leq n\leq k},t_{m=k+1})$ and $(t_{k+1},t_{k+1})$ we need to find $F_{1}^{>,<}(t_{n=k+1},t_{m})$, $F_{2}^{>,<}(t_{n},t_{m=k+1})$ and $F_{2}^{>,<}(t_{n=k+1},t_{m=k+1})$ which can be done by substituting the predicted values for $G^{>,<}$, derived earlier, into eqn. (S1.1) and obtain

[TABLE]

A prediction for the self-energies appearing above can be obtained by using the closure relation given in eqn. (A.1). Having calculated the new values for $F_{1,2}$, we can now get a better estimate for the Green’s functions as follows

[TABLE]

Using these corrected values we can find the corrected value for $G^{>,<}(t_{k+1},t_{k+1})$ as well, which is

[TABLE]

The above correction process can be carried out arbitrary number (we have found that one corrective iteration is sufficient to get good results) of times till one finds that the new values for $[-\infty,t_{k+1}]\cup[-\infty,t_{k+1}]$ have converged to a desired accuracy.

S1.3 Evaluation schemes for the integrands

In the above discussion we chose a very simple integration scheme for the integrals appearing in the KB equations for the sake of clarity. We now slightly modify the scheme to get a better error performance. We approximate the integrands in eqn. (S1.1) using the trapezoidal rule, such that

[TABLE]

and use these new expressions at every place where the evaluation of $F_{1,2}^{>,<}$ were performed in the previous discussion.

S2 Equilibrium Green’s functions Calculations

To obtain the equilibrium Green’s functions we follow a similar approach to Banerjee and Altman (2017), and incorporate the conditions of equilibrium into the saddle point equations eqn. (A.1). First, we assume time translation symmetry, which makes every two-point correlator a function of $t=t_{1}-t_{2}$. Therefore, we can write eqn. (A.1) as

[TABLE]

where we have made the $V$-term time independent, and introduced the symbols – $J_{s}=J^{2}\ (t_{\psi}^{2})$, $q_{s}=2\ (1)$, $\bar{s}=\psi\ (c)$ when $s=c\ (\psi)$. The exponent of the factor $\sqrt{p}$ takes the value $+1$ when $s=c$ and $-1$ when $s=\psi$. To calculate the Green’s functions for the (dis)connected system we set $V=(1)\ 0$. Using the relations between $G^{>}$, $G^{<}$ and $G^{K},\ G^{R},\ G^{A}$ (see sec. S2) and the condition eqn. (17), we have

[TABLE]

where $\rho_{s}(\omega)$ are the spectral functions for the respective fermion flavors. Writing $G^{>,<}(t)$ in terms of $G^{>,<}(\omega)$ in eqn. (S2) and using the definition for retarded function $\Sigma^{R}(t)=\Theta(t)[\Sigma^{>}(t)-\Sigma^{<}(t)]$, we find $\Sigma^{R}(\omega)$ to be

[TABLE]

where $n_{(1-4)s}(t)$ are defined as

[TABLE]

The integro-differential equations in eqn. (A.1), under the assumptions of equilibrium, reduces to the Dyson equations for both the fermion flavors and takes the form

[TABLE]

where $G_{s}^{R}$ is the retarded Green’s function. The spectral function appearing in eqn. (S2) and eqn. (S2.47) is related to the retarded Green’s function as

[TABLE]

Using eqn. (S2.47), eqn. (S2.51) and eqn. (S2.52) we iteratively solve for the spectral function $\rho_{s}(\omega)$. The iterative process is terminated when we have converged to a solution for $\rho_{s}(\omega)$ with sufficient accuracy. Using the converged value of $\rho_{s}(\omega)$ the equilibrium versions of $G^{>}(t)$ and $G^{<}(t)$ can be obtained by using the left most equations in eqn. (S2) and then taking a Fourier transform.

S3 Time dependent expectation values

The time dependent expectation values for the energy components $E_{c}$,$E_{\psi}$, $E_{c\psi}$ etc. discussed in the main text are derived using the identityKamenev (2011)

[TABLE]

where $Z_{neq}$ is the generating functional, introduced in sec. A, now containing an additional source field $\eta(z)\hat{O}(t)$ with $\eta(z)$ defined as

[TABLE]

and $\hat{O}$ is an operator whose expectation value we want to evaluate. The following source fields are used to derive the time dependent energy components

$(1+\eta(z)){1\over(2{N})^{3/2}}\sum\limits_{\mathclap{ijkl}}J_{ijkl}\bar{c}_{i}\bar{c}_{j}c_{k}c_{l}$, $(1+\eta(z)){1\over{M}^{1/2}}\sum\limits_{\mathclap{\alpha\beta}}t_{\alpha\beta}^{\psi}\bar{\psi}_{\alpha}\psi_{\beta}$ for evaluating $\langle{\cal{H}}_{c}(t)\rangle$ $\equiv E_{c}(t)$ and $\langle{\cal{H}}_{\psi}(t)\rangle$ $\equiv E_{\psi}(t)$ respectively. 2. 2.

$f(z)(1+\eta(z))\sum\limits_{\mathclap{i\alpha}}V_{i\alpha}\bar{c}_{i}\psi_{\alpha}+\mbox{h.~{}c.}$, for evaluating $\langle{\cal{H}}_{c\psi}(t)\rangle$ $\equiv E_{c\psi}(t)$.

Using the identity in eqn. (S3.53) we calculate the disorder averaged time dependent expectation values for the individual parts of the Hamiltonian defined in eqn. (1), and find

[TABLE]

Using the above quantities, an excitation energy produced during the quench can be defined as

[TABLE]

where $\langle\mathcal{H}(\tau)\rangle_{T_{i}}$ is the thermal expectation of energy for the final Hamiltonian $\mathcal{H}(\tau)$ at the initial temperature $T_{i}$ from which the quench began.

An additional observable can be defined corresponding to the occupation density $n_{\psi}(\epsilon,t)$ of a $\psi$-fermion having energy $\epsilon$. This is possible since, the $\psi$ fermions in the non-interacting model are connected via random hoppings, (see eqn. (2)) and hence can be diagonalized to yield a set of eigenstate fermions which we label by $\epsilon$. The Kadanoff-Baym equations satisfied by the Green’s functions for the $\psi$ fermions are

[TABLE]

where

[TABLE]

The number density for the $\psi$ fermions can be obtained from the Keldysh Green’s function using

[TABLE]

where $G^{K}_{\psi}(\epsilon;t,t)=G^{>}_{\psi}(\epsilon;t,t)+G^{<}_{\psi}(\epsilon;t,t)$.

S4 Sudden Quench

We now briefly discuss some of the results obtained from the sudden quench of the interacting as well as the non-interacting models, and also provide a comparison between the two cases.

S4.1 Energy time-series

We evaluate the total energy for the system, i.e.

[TABLE]

as a function of time for various values of site-fraction $p$. The time-series data is shown in fig. S2(a). Interestingly, the sudden change of $f(t)$ in eqn. (1), from 0 to 1, is an iso-energetic process and occurs without any work input to the system. The tiny fluctuations near $t=0$ are an artifact of discretization and reduces as the resolution of the time grid (see fig. S1) is increased. We verify the iso-energetic nature of the quench by comparing the final temperature obtained from the FDT relation in eqn. (17), with the temperature obtained from equilibrium calculations. The temperature in the latter case was found by solving for the value that produced the same energy for the connected $c-\psi$ system equal to the initial disconnected one. As shown in fig. S2(b), the temperature vs. $p$ curves obtained for the above two cases are qualitatively similar proving that the sudden quench of the interacting model is indeed iso-energetic. The sudden quench of the non-interacting model is an iso-energetic process as well. However, unlike the interacting case, the non-interacting energy time-series is independent of $p$. This is consistent with the model, since in the initial disconnected system the $c$ and $\psi$ fermions are both described by a semi-circular density of states scaled appropriately by factors of $1/(1+p)$ and $p/(1+p)$ respectively, which cancel out when the energies of the two fermion flavors are added.

S4.2 Thermalization of the occupation function and $E_{c\psi}$

The thermalization behavior for the interacting model is shown in fig. S3. The energy associated with the bonds between the $c$ and $\psi$ sites are shown as a function of time for various values of site-fraction $p$. We find for small values of $p$, i.e. $p<1$, the energy reaches the equilibrium value (shown with arrow heads, see fig. S3(a)) very rapidly. This is expected, since the $p<1$ phase after the quench is a SYK-type NFL, which is known to be a fast thermalizer. On the other hand, for values of $p>3.0$ the approach to the equilibrium state slows down rapidly, as demonstrated by the negative slope of $E_{c\psi}-t$ curves, for $t>0$, in fig. S3(b). Again, this behavior is consistent with the equilibrium model since at larger values of $p$ a FL-state is expected to dominate the thermalization time scales of the system. The same phenomenon is observed in the distribution functions of the $\psi$ fermions as well and which is reported in the main text. For small values of $p$, like $p=0.1$, (see fig. 2(e) top) the distribution function, evaluated sometime after the quench, matches well with that obtained from a equilibrium calculation using the temperature determined via energy conservation. However, at a large value of $p=8.0$, see fig. 2(e) bottom, the match is poor due to the slowed down approach towards equilibrium.

Moving on to the sudden quench of the non-interacting model, we find a drastically different behavior, in that the system completely fails to thermalize. Fig. S4(a) demonstrates this feature very clearly, the energy $E_{c\psi}$ reaches a steady state value, right after the quench, which is very different from the expected equilibrium value. Unlike the interacting model for large $p$ values, where approach to equilibrium was simply slowed down, in this situation we have a complete halt on equilibration. The above fact is further supported by the form of the $\psi$-fermion distribution function as shown in fig. S4(b) and (c). Irrespective of the value of site-fraction $p$, the distribution function fails to converge towards its equilibrium counterpart. This is in sharp contrast to the interacting model, where we had a regime of $p$ values for which the distribution function converged really well with the equilibrium results.

In summary, we find that the interacting model can be distinguished from the non-interacting model by studying their thermalization behavior. In the case of interacting model, the system either thermalizes rapidly or slowly depending on the regime of $p$ values we are looking at. On the other hand, the non-interacting model ceases to thermalize at all and reaches a steady state, far away from equilibrium, irrespective of the value of site fraction $p$.

S4.2.1 Exact diagonalization study of thermalization of $E_{c\psi}$ at finite $N$

As discussed in the main text, one can ask whether the crucial features of the large-$N$ non-equilibrium dynamics, e.g. the fast and slow thermalization in the NFL and FL, respectively, persist at finite $N$. To answer this, we perform numerical exact diagonalization (ED) studies of the sudden quench in our model. We take the $T=0$ direct product ground state $|\Psi_{0}\rangle=|NFL\rangle\otimes|FL\rangle$ as the initial state of the initially decoupled system, where $|NFL\rangle$ is the ground state of the SYK model (eqn. (1a)) with $N$ sites and $|FL\rangle$ the ground state of the lead fermions described by eqn. (1b) for a particular disorder realization. As in the large-$N$ calculations, we switch on the coupling term $\mathcal{H}_{c\psi}$ (eqn. (1c)) at $t=0$. After the quench, $|\Psi_{0}\rangle$ is no longer an eigenstate of the Hamiltonian $\mathcal{H}$ (eqn. (1)) and its time evolution is given by $|\Psi(t)\rangle=e^{-\mathfrak{i}\mathcal{H}t}|\Psi_{0}\rangle$. The disorder-averaged $E_{c\psi}(t)=\overline{\langle\Psi_{0}(t)|\mathcal{H}_{c\psi}|\Psi_{0}(t)\rangle}$ is shown in fig. S3(c) and in fig. 2(h) (main text). The disorder average is taken over $300$ disorder realizations. To obtain the time-evolution, the initial state $|\Psi_{0}\rangle=\sum_{n}c_{n}|\phi_{n}\rangle$ is written in terms of eigenstate $|\phi_{n}\rangle$ of the post quench Hamiltonian $\mathcal{H}$, so that $|\Psi(t)\rangle=\sum_{n}c_{n}e^{-\mathfrak{i}E_{n}t}|\phi_{n}\rangle$, where $c_{n}=\langle\phi_{n}|\Psi_{0}\rangle$ and $E_{n}$’s are eigen energies of $\mathcal{H}$.

The ED results are obtained at half filling for total system size $N+M=16$. The plots (fig. S3(c) and in fig. 2(h)) are shown for $p=0.33$ (i.e. $N=12,M=4$) and for $p=3.0$ (i.e. $N=4,M=12$). According to eigenstate thermalization hypothesis (ETH) D’Alessio et al. (2016), if the coupled system thermalizes, then, in the long-time limit, $E_{c\psi}$ is expected to reach the value, $E^{d}_{c\psi}=\overline{\mathrm{Tr}(\rho_{d}H_{c\psi})}$, given by the diagonal ensemble corresponding to the density matrix $\rho_{d}=\sum_{n}|c_{n}|^{2}|\phi_{n}\rangle\langle\phi_{n}|$. The expectation value in the diagonal ensemble matches with the thermal expectation, that is obtained from microcanonical ensemble, within a sub-extensive correction. We have checked that there is indeed a sub-extensive difference ($\sim 1\%$ ) between the values of $E_{c\psi}$ in the diagonal and microcanonical ensembles . Hence, to avoid finite size effects we compare long time $E_{c\psi}(t)$ with the diagonal expectation value. As shown in fig. S3(c), for $p=0.33$, $E_{c\psi}(t)$ reaches diagonal ensemble expectation value much faster than that for $p=3.0$. Unlike the large-$N$ results in Fig.2(e), we could not study the thermalization deeper in the FL phase due to the limitation of system sizes ($N$ and $M$) in ED. Nevertheless, the difference of thermalization rates in the NFL and FL is quite evident even from the comparison of $E_{c\psi}(t)$ for $p=0.33$ and $p=3$.

S4.3 Collapse-revival to prethermal transition in the non-interacting model

In this subsection we discuss the sudden quench of the non-interacting defined in sec. A.2. The results of this exercise are given in fig. S5. Indeed, we find the same qualitative behavior as the interacting case, however with rescaled values for $p_{c}^{dyn}$. In fact, the way the oscillations disappear as a function of $p$ (see fig. S5(c)) are very similar to the interacting case (see inset in fig. 2(a)) with $p_{c}^{dyn}\approx 0.5$ in this case. The oscillations in $z_{\psi}$ that appear in fig. 2(a) are also reproduced (see fig. S5(a)) and occur when $p\leq{p_{c}^{dyn}}$. The prethermal plateaus and long thermalization times shown in fig. 2(b) are obtained at higher but different values of $p$, see fig. S5(b). Hence overall, the relaxation features of the fermion distribution function associated with the $\psi$ fermions, for the interacting model, are also observed in the sudden quench of the non-interacting model. Furthermore, the value of $p_{c}^{dyn}$ for this case can be explained by closing of a soft-gap in the spectral function of the fermions, see sec. B for details.

S5 Slow quench

In this section, we discuss the details of the results that we obtained from the slow quench studies of both the non-interacting and the interacting models.

S5.1 Effect of ramp shapes on the interacting model quench

We join the SYK $c$-fermions with $\psi$-fermions over a time-period $\tau$ using ramps having various shapes and heights.

The ramp shapes that we use, and shown in fig. S6, are chosen following ref.Eckstein and Kollar, 2010, and are classified in the order of their smoothness. In particular, $r_{n}(t)$ will have its $n$-th derivative to be discontinuous at the start and end of the ramp (fig. S6(b)-(d)). The ramp functions for $n=1,2,3$ are given below

[TABLE]

The $\eta-p$ curves, calculated at an initial temperature $T_{i}=0.05$, for the ramp shapes in eqn. (S5.61) are shown in fig. S7(a). We find that there exists a dependence on ramp-shapes. However, this dependence diminishes drastically around the transition point $p=1$, indicating that the intrinsic properties of the critical point are becoming more prevalent. Deep within the phases, the dependence on ramp-shape is the strongest, with the the stronger effect manifesting in the FL limit($p\to\infty$) (see fig. S7(b),(c)).

S5.1.1 Effect of temperature on the interacting model quench

In the case of the interacting model the initial temperature $T_{i}=0$ is not readily accessible, therefore we ask how does the exponent $\eta$ change as a function $T_{i}$ for any given $p$ value? The behavior of $\eta(p,T_{i})$ as a function of $1/T_{i}$, for $p=0.1$, is shown in fig. S8(a), and suggests that $\eta(p,T_{i})$ converges to a finite value at low temperatures. Further, we obtain the $\eta-p$ curves (see fig. S8(b)) at various temperatures, and find that they indeed converge to a finite value as $T_{i}\to 0$. This suggests that the powerlaw behavior that we observe is a genuine property of the quantum ground states involved in the quench.