Quantum quench of the Sachdev-Ye-Kitaev Model

TL;DR

This paper investigates the non-equilibrium dynamics of the Sachdev-Ye-Kitaev models, revealing thermalization behavior and providing exact solutions in the large q limit, advancing understanding of quantum systems without quasiparticles.

Contribution

It offers a detailed analysis of quantum quench dynamics in SYK models, including numerical and exact solutions, highlighting thermalization processes in these complex systems.

Findings

Final state is thermal according to Kadanoff-Baym equations

Thermalization rate is proportional to the absolute temperature

In the large q limit, thermalization of Green's function is instantaneous

Abstract

We describe the non-equilibrium dynamics of the Sachdev-Ye-Kitaev models of fermions with all-to-all interactions. These provide tractable models of the dynamics of quantum systems without quasiparticle excitations. The Kadanoff-Baym equations show that the final state is thermal, and their numerical analysis appears consistent with a thermalization rate proportional to the absolute temperature of the final state. We also obtain an exact analytic solution of the non-equilibrium dynamics in the large $q$ limit of a model with $q$ fermion interactions: in this limit, the thermalization of the fermion Green's function is instantaneous.

Figures

Equations

$$H=(i{)}^{\frac{q}{2}}\sum _{1\le i_1<i_2<...<i_q\le N}{j}_{i_1{i}_2...i_q}{\psi }_{i_1}{\psi }_{i_2}...{\psi }_{i_q}$$

$$\{{\psi }_i,{\psi }_j\}={\delta }_{ij},$$

$$⟨j_{i_1...i_q}^2⟩=\frac{J^2(q-1)!}{N^{q-1}}.$$

$$H(t)=i\sum _{i<j}{j}_{2,ij}f(t){\psi }_i{\psi }_j-\sum _{i<j<k<l}{j}_{4,ijkl}g(t){\psi }_i{\psi }_j{\psi }_k{\psi }_l,$$

$$\displaystyle\mbox{for $t<0$, $f(t)=1$ and $g(t)=1$}\,;$$

$$\displaystyle\mbox{for $t>0$, $f(t)=0$ and $g(t)=1$}\,.$$

$$\mathcal{T}=(t_1+t_2)/2,t=t_1-t_2.$$

$$\underset{\mathcal{T}\to \infty }{\lim }G(t_1,t_2)=G_{{\beta }_f}(t)$$

$${\beta }_{\text{eff}}(\mathcal{T})={\beta }_f+\alpha \exp \left(-\Gamma \mathcal{T}\right).$$

$$\Gamma =C/{\beta }_f,{\beta }_f{J}_4\gg 1,$$

$$G^{>}(t_1,t_2)=-i⟨\psi (t_1)\psi (t_2)⟩=-\frac{i}{2}\left[1+\frac{1}{q}g(t_1,t_2)+\dots \right]$$

$$g(t,t)=0.$$

$$\frac{{\partial }^2}{\partial t_1\partial t_2}g(t_1,t_2)=2{\mathcal{J}}^2{e}^{g(t_1,t_2)},$$

$$g(t_1,t_2)=\ln \left[\frac{-h_1^{\prime }(t_1)h_2^{\prime }(t_2)}{{\mathcal{J}}^2(h_1(t_1)-h_2(t_2){)}^2}\right],$$

$$h_{\alpha }(t)\to \frac{a{h}_{\alpha }(t)+b}{c{h}_{\alpha }(t)+d},$$

$$h_1(t)=\frac{a{e}^{i\theta }{e}^{\sigma t}+b}{c{e}^{i\theta }{e}^{\sigma t}+d},h_2(t)=\frac{a{e}^{-i\theta }{e}^{\sigma t}+b}{c{e}^{-i\theta }{e}^{\sigma t}+d}.$$

$$\sigma =2\mathcal{J}\sin (\theta ).$$

$$g(t_1,t_2)=\ln \left[\frac{-{\sigma }^2}{4{\mathcal{J}}^2{\sinh }^2(\sigma (t_1-t_2)/2+i\theta )}\right].$$

$${\beta }_f=\frac{2(\pi -2\theta )}{\sigma }.$$

$$\sigma =\frac{2\pi }{{\beta }_f}.$$

$$Z=\int \mathcal{D}\psi e^{iS[\psi ]}$$

$$\begin{split}S[\psi]&=\int_{\mathcal{C}}dt\Bigl{\{}\frac{i}{2}\sum_{i}\psi_{i}\partial_{t}\psi_{i}-i\sum_{i<j}j_{2,ij}f(t)\psi_{i}\psi_{j}+\sum_{i<j<k<l}j_{4,ijkl}g(t)\psi_{i}\psi_{j}\psi_{k}\psi_{l}\Bigr{\}}\,.\end{split}$$

$$Z=\int D\psi \int D{j}_{2,ij}\int D{j}_{4,ijkl}P(j_{4,ijkl})P(j_{2,ij})e^{iS[\psi ]},$$

$$\displaystyle\mathcal{P}_{1}(j_{2,ij})=\sqrt{\frac{N}{2J_{2}^{2}\pi}}\operatorname{exp}\Big{(}-\frac{N}{2J_{2}^{2}}\sum_{i<j}j_{2,ij}^{2}\Big{)},$$

$$\displaystyle\mathcal{P}_{2}(j_{4,ijkl})=\sqrt{\frac{N^{3}}{12J_{4}^{2}\pi}}\operatorname{exp}\Big{(}-\frac{N^{3}}{12J_{4}^{2}}\sum_{i<j<k<l}j^{2}_{4,ijkl}\Big{)}.$$

$$\begin{split}Z=\int D\psi\operatorname{exp}\Bigl{\{}&-\frac{1}{2}\int_{\mathcal{C}}dt_{1}\sum_{i}\psi_{i}\partial_{t_{1}}\psi_{i}-\frac{J_{2}^{2}}{4N}\sum_{i,j}\int_{\mathcal{C}}dt_{1}\int_{\mathcal{C}}dt_{2}f(t_{1})f(t_{2})\psi_{i}(t_{1})\psi_{i}(t_{2})\psi_{j}(t_{1})\psi_{j}(t_{2})\\ &-\frac{3J_{4}^{2}}{4!N^{3}}\sum_{i,j,k,l}\int_{\mathcal{C}}dt_{1}\int_{\mathcal{C}}dt_{2}g(t_{1})g(t_{2})\psi_{i}(t_{1})\psi_{i}(t_{2})\psi_{j}(t_{1})\psi_{j}(t_{2})\psi_{k}(t_{1})\psi_{k}(t_{2})\psi_{l}(t_{1})\psi_{l}(t_{2})\Bigr{\}}.\end{split}$$

$$G(t_1,t_2)=-\frac{i}{N}\sum _i{\psi }_i(t_1){\psi }_i(t_2).$$

$$\begin{split}Z=&\int D\psi\int DG\int D\Sigma\operatorname{exp}\Bigl{\{}-\frac{1}{2}\int_{\mathcal{C}}dt_{1}\sum_{i}\psi_{i}\partial_{t_{1}}\psi_{i}+\frac{J_{2}^{2}N}{4}\int_{\mathcal{C}}dt_{1}\int_{\mathcal{C}}dt_{2}f(t_{1})f(t_{2})G(t_{1},t_{2})^{2}\\ &-\frac{3J_{4}^{2}N}{4!}\int_{\mathcal{C}}dt_{1}\int_{\mathcal{C}}dt_{2}g(t_{1})g(t_{2})G(t_{1},t_{2})^{4}+\frac{i}{2}\int_{\mathcal{C}}dt_{1}\int_{\mathcal{C}}dt_{1}\Sigma(t_{1},t_{2})\bigl{[}G(t_{1},t_{2})+\frac{i}{N}\sum_{i}\psi_{i}(t_{1})\psi_{i}(t_{2})\bigr{]}\Bigr{\}}.\end{split}$$

$$\begin{split}Z=&\int DG\int D\Sigma\operatorname{exp}\Bigl{\{}iS\bigl{[}G,\Sigma\bigr{]}\Bigr{\}}\end{split}$$

$$\begin{split}S\bigl{[}G,\Sigma\bigr{]}=&-\frac{iN}{2}\operatorname{tr}\operatorname{log}\bigl{[}-i(G_{0}^{-1}-\Sigma)\bigr{]}-\frac{iJ_{2}^{2}N}{4}\int_{\mathcal{C}}dt_{1}\int_{\mathcal{C}}dt_{2}f(t_{1})f(t_{2})G(t_{1},t_{2})^{2}\\ &+\frac{3iJ_{4}^{2}N}{4!}\int_{\mathcal{C}}dt_{1}\int_{\mathcal{C}}dt_{2}g(t_{1})g(t_{2})G(t_{1},t_{2})^{4}+\frac{iN}{2}\int_{\mathcal{C}}dt_{1}\int_{\mathcal{C}}dt_{2}\Sigma(t_{1},t_{2})G(t_{1},t_{2})\end{split}$$

$$G_0^{-1}(t_1,t_2)-G^{-1}(t_1,t_2)$$

Full text

I Introduction

The non-equilibrium dynamics of strongly interacting quantum many-particle systems have been the focus of much theoretical work Kamenev (2013). These are usually studied within the Schwinger-Keldysh formalism, which can describe evolution from a generic initial state to a final state which reaches thermal equilibrium at long times. The thermodynamic parameters of the final state (e.g. temperature) are determined by the values of the conjugate conserved quantities (e.g. energy). The Kadanoff-Baym equations obtained from this formalism describe the manner and rate by which this final thermal state is reached.

The Kadanoff-Baym equations are usually too difficult to solve in their full generality. Frequently, a quasiparticle structure has been imposed on the spectral functions, so that the Kadanoff-Baym equations reduce to a quantum Boltzmann equation for the quasiparticle distribution functions. Clearly such an approach cannot be employed for final states of Hamiltonians which describe critical quantum matter without quasiparticle excitations. The most common approach is then to employ an expansion away from a regime where quasiparticles exist, using a small parameter such as the deviation of dimensionality from the critical dimension, or the inverse of the number of field components: the analysis is then still carried out using quasiparticle distribution functions Damle and Sachdev (1997); Sachdev (1998).

In this paper, we will examine the non-equilibrium dynamics towards a final state without quasiparticles. We will not employ a quasiparticle decomposition, and instead provide solutions of the full Kadanoff-Baym equations. We can obtain non-equilibrium solutions for the Sachdev-Ye-Kitaev (SYK) models Sachdev and Ye (1993); Kitaev (2015); Maldacena and Stanford (2016) with all-to-all and random interactions between $q$ Majorana fermions on $N$ sites. These models are solvable realizations of quantum matter without quasiparticles in equilibrium, and here we shall extend their study to non-equilibrium dynamics. We shall present numerical solutions of the Kadanoff-Baym equations for the fermion Green’s function at $q=4$, and in Section IV, an exact analytic solutions in the limit of large $q$ (no quasiparticles are present in this limit). The large $q$ solution relies on a remarkable exact SL(2,C) invariance, and has connections to quantum gravity on AdS2 with the Schwarzian effective action Sachdev (2010a, b); Almheiri and Polchinski (2015); Kitaev (2015); Maldacena and Stanford (2016); Jensen (2016); Engelsöy et al. (2016) for the equilibrium dynamics, as we will note in Section V.

The equilibrium SYK Hamiltonian we shall study is (using conventions from Ref. Maldacena and Stanford, 2016)

[TABLE]

where $\psi_{i}$ are Majorana fermions on sites $i=1\ldots N$ obeying

[TABLE]

and $j_{i_{1}i_{2}...i_{q}}$ are Gaussian random variables with zero mean and variance

[TABLE]

Our numerical analysis of the SYK model will be carried out on the case with time-dependent $q=2$ and $q=4$ terms:

[TABLE]

where $j_{2,ij}$ and $j_{4,ijkl}$ are random variables as described above, and $f(t)$ and $g(t)$ can be arbitrary functions of time. We will focus on a particular quench protocol in which

[TABLE]

For $t<0$, we have a thermal initial state with both $q=2$ and $q=4$ interactions present. In this case the $q=2$ free fermion terms dominate at low energies, and hence we have Fermi liquid behavior at low temperatures. For $t>0$, we have only the $q=4$ Hamiltonian as in Eq. (1), which describes a non-Fermi liquid even at the lowest energies Sachdev and Ye (1993). Just after the quench at $t=0$, the system is in a non-equilibrium state. We show from the Kadanoff-Baym equations, and verify by our numerical analysis, that the $t\rightarrow\infty$ state is indeed a thermal non-Fermi liquid state at $q=4$, which is at an inverse temperature $\beta_{f}=1/T_{f}$. The value of $\beta_{f}$ is such that the total energy of the system remains the same after the quench at $t=0^{+}$.

We note in passing that an alternative quench protocol in which $f(t)=0$ for all $t$, while $g(t)$ is time dependent (or the complementary case in which $g(t)=0$ for all $t$, while $f(t)$ is time-dependent), does not yield any non-trivial dynamics (this was pointed out to us by J. Maldacena). This case corresponds to an overall rescaling of the Hamiltonian, which does not change any of the wavefunctions of the eigenstates. So an initial state in thermal equilibrium will evolve simply via a time reparameterization which is given in Eq. (54).

We will describe the time evolution of the system for $t\geq 0$ by computing a number of two-point fermion Green’s functions $G(t_{1},t_{2})$. As the system approaches thermal equilibrium, it is useful to characterize this correlator by

[TABLE]

Then we have

[TABLE]

where $G_{\beta_{f}}$ is the correlator in equilibrium at an inverse temperature $\beta_{f}$. For large $\mathcal{T}$ we will characterize $G(t_{1},t_{2})$ by an effective temperature $\beta_{\rm eff}(\mathcal{T})$ (see below). From our numerical analysis we characterize the late time approach to equilibrium by

[TABLE]

This defines a thermalization rate, $\Gamma$. Our complete numerical results for $\Gamma$ appear in Fig. 7, and they are reasonably fit to the behavior

[TABLE]

where $C$ is a numerical constant independent of the initial state. Thus, at low final temperatures, the thermalization rate of the non-Fermi liquid state appears proportional to temperature, as is expected for systems without quasiparticle excitations Sachdev (1999).

I.1 Large $q$ limit

Further, we will also consider a model in which the $t<0$ Hamiltonian has a $q$-fermion interaction $\mathcal{J}$, and a $pq$ fermion interaction $\mathcal{J}_{p}$. At $t>0$ we then quench to only a $q$ fermion interaction $\mathcal{J}$. We will show that the non-equilibrium dynamics of this model are exactly solvable in the limit of large $q$ taken at fixed $p$.

The large $q$ Green’s function can be written as Maldacena and Stanford (2016)

[TABLE]

where

[TABLE]

For $t_{1}>0$ or $t_{2}>0$, we find that $g(t_{1},t_{2})$ obeys the two-dimensional Liouville equation

[TABLE]

where $\mathcal{J}$ is proportional to the final interaction strength (see Eq. (56)). The most general solution of this equation can be written as Tsutsumi (1980)

[TABLE]

where $h_{1}(t_{1})$ and $h_{2}(t_{2})$ are arbitrary functions of their arguments. A remarkable and significant feature of this expression for $g(t_{1},t_{2})$ is that it is exactly invariant under SL(2,C) transformations

[TABLE]

where $\alpha=1,2$ and $a,b,c,d$ are arbitrary complex numbers. We will use the Schwinger-Keldysh analysis to derive ordinary differential equations that are obeyed by $h_{1,2}(t)$ in Section IV. Naturally, these equations will also be invariant under SL(2,C) transformations. We will show that for generic initial conditions in the regime $t_{1}<0$ and $t_{2}<0$, the solutions for $h_{1}(t_{1})$ and $h_{2}(t_{2})$ at $t_{1}>0$ and $t_{2}>0$ can be written as

[TABLE]

The complex constants $a,b,c,d$, and the real constants $\sigma$, $\theta$ are determined by the initial conditions in the $t_{1}<0$ and $t_{2}<0$ quadrant of the $t_{1}$-$t_{2}$ plane. Note that we pick a particular SL(2,C) orientation in the $t_{1}\leq 0$, $t_{2}\leq 0$ quadrant, and then there is no further SL(2,C) arbitrariness in the $t_{1}>0$, $t_{2}>0$ quadrant. Inserting Eq. (15) into Eqs. (11,13) at $t_{1}=t_{2}$, we obtain

[TABLE]

For general $t_{1}>0$ and $t_{2}>0$, inserting Eq. (15) into Eq. (13) we obtain

[TABLE]

The surprising feature of this result is that it depends only upon $t=t_{1}-t_{2}$, and is independent of $\mathcal{T}=(t_{1}+t_{2})/2$. Indeed Eq. (17) describes a state in thermal equilibrium Parcollet and Georges (1999); Maldacena and Stanford (2016) at an inverse temperature

[TABLE]

So the large $q$ limit yields a solution in which the fermion Green’s function thermalizes instantaneously at $t=0^{+}$ at the final temperature given by Eq. (18). This could indicate that the thermalization rate $\Gamma$ diverges as at $q\rightarrow\infty$. Alternatively, as pointed out to us by Aavishkar Patel, the present large $q$ solution could describe a pre-thermal state, and the $1/q^{2}$ corrections in Eq. (10) have a finite thermalization rate; in such a scenario, thermalization is a two-step process, with the first step occuring much faster than the second. But we will not examine the $1/q^{2}$ corrections here to settle this issue.

Note also that in the limit $\beta_{f}\mathcal{J}\gg 1$, we have $\theta\ll 1$ and then

[TABLE]

Then Eqs. (10) and (17) describe the $1/q$ expansion of the low temperature conformal solution Parcollet and Georges (1999) describing the equilibrium non-Fermi liquid state (see Appendix B). And the value of $\sigma$ in Eq. (19) is the maximal Lyapunov exponent for quantum chaos Maldacena et al. (2016). This chaos exponent appears in the time evolution of $h_{1,2}(t)$. However, in the large $q$ limit, it does not directly control the rapid thermalization rate. We note that recent studies of Fermi surfaces coupled to gauge fields, and of disordered metals, found a relaxation/dephasing rate which was larger than the Lyapunov rate Patel and Sachdev (2017); Patel et al. (2017).

We will begin by setting up the Schwinger-Keldysh formalism for Majorana fermions in the SYK model in Sections II. The numerical solution of the Kadanoff-Baym equations appears in Section III. Finally, the large $q$ limit is described in Section IV.

II Kadanoff-Baym equations from the path integral

We construct the path integral for Majorana fermions from the familiar complex fermion path integral Sedrakyan and Galitski (2011) by expressing the complex fermions $\Psi_{i}$ in terms of two real fermions $\psi_{i}$ and $\chi_{i}$ i.e. $\Psi_{i}=\psi_{i}+i\chi_{i}$. Since $\chi_{i}$ is just a spectator which does not appear in the Hamiltonian, we can disregard its contribution and write down the path integral representation of the partition function as

[TABLE]

with the action

[TABLE]

Here, the Majorana fields $\psi_{i}$ live on the closed time contour $\mathcal{C}$, also known as Schwinger-Keldysh contour. In order to include the disorder, we average the partition function over Gaussian distributed couplings as

[TABLE]

where the probability distributions are given by

[TABLE]

This realization of disorder allows us to perform the integrals with respect to $j_{2,ij}$ and $j_{4,ijkl}$. The Gaussian integrals lead to

[TABLE]

We introduce the bilinear $G$, which must fulfill the constraint

[TABLE]

To implement this constraint, we introduce the Lagrange multiplier $\Sigma$ leading to

[TABLE]

Rescaling integration variables as $\Sigma\rightarrow iN\Sigma$ and integrating over $\psi$ yields

[TABLE]

with

[TABLE]

and the free Majorana Green’s function $G_{0}^{-1}(t_{1},t_{2})=i\partial_{t}\delta_{\mathcal{C}}(t_{1},t_{2})$. ​Varying the action with respect to $\Sigma$ and $G$ yields the Dyson equation and the the self-energy, respectively:

[TABLE]

Note that the time arguments are with respect to the full time contour, and the matrix structure of the Keldysh formalism is implicit. We define the following two Green’s functions

[TABLE]

where $t_{i}^{-}$ lives on the lower contour and $t_{i}^{+}$ lives on the upper contour In the case for Majorana fermions there is only one independent component of the Green’s function, even in non-equilibrium, due to the relation Babadi et al. (2015)

[TABLE]

The bare greater and lesser Green’s functions are given by

[TABLE]

We now use $G^{>}$ and $G^{<}$ to obtain the retarded, advanced and Keldysh Green’s functions:

[TABLE]

Similarly to the Green’s function we introduce retarded and advanced self-energies

[TABLE]

and for more details on the Schwinger-Keldysh formalism and the saddle point approximation we refer to Refs. Babadi et al., 2015; Kamenev, 2013; Maldacena and Stanford, 2016. In order to obtain the Kadanoff-Baym equations, we rewrite the Schwinger-Dyson equations (28) as

[TABLE]

Using the Langreth rules Stefanucci and van Leeuwen (2013)

[TABLE]

we obtain the equations of motion for $G^{>}$ from the last equation

[TABLE]

where

[TABLE]

and $\Sigma^{R}$ and $\Sigma^{A}$ can be obtained from (34) and  (35). In the following sections we will solve this system of equations numerically and also solve the generalization of this model in the large q-limit.

III Kadanoff-Baym equations: Numerical Study

In the following we study the non-equilibrium dynamics described by the time-dependent Hamiltonian in Eq. (4). The functions $f(t)$ and $g(t)$, which are arbitrary so far, specify the quench protocol. We will use them so switch on(off) or rescale couplings. Hence, all considered quench protocols are of the form

[TABLE]

and will be denoted by $(J_{2,i},J_{4,i})\rightarrow(J_{2,f},J_{4,f})$. Specifically, we choose the following four protocols

• A)

$(J_{2,i},0)\rightarrow(J_{2,f},0)$ quench: The energy scale of the random hopping model gets suddenly rescaled.

• B)

$(0,0)\rightarrow(0,J_{4,f})$ quench: The system is quenched from an uncorrelated state by switching on the interaction of the SYK model.

• C)

$(J_{2,i},0)\rightarrow(0,J_{4,f})$ quench: The system is prepared in the ground state of the random hopping model and then quenched to the SYK Hamiltonian.

• D)

$(J_{2,i},J_{4,f})\rightarrow(0,J_{4})$ quench: Sarting from a ground state with the quadratic and the quartic terms present the quadratic term is switched off.

The Kadanoff-Baym equations for the SYK model in Eqs. (40) and (41) are be solved numerically. Due to the absence of momentum dependence we are able to explore the long time regime after the quench. The quench time appears at $t_{1}=t_{2}=0$. For $t_{1}<0$ and $t_{2}<0$, we are in thermal equilibrium and the time dependence of the Green’s function is determined by the Dyson equation. When quenching the ground state of the random hopping model, we can use the exact solution for $G^{>}$ as an initial condition. In all other cases, we solve the Dyson equation self-consistently. For further details we refer to Appendix C. For $t_{1}>0$ or $t_{2}>0$ we solve Kadanoff-Baym equation. Integrals in the Kadanoff-Baym equation are computed using the trapezoidal rule. After discretizing the integrals the remaining ordinary differential equations are solved using a predictor-corrector scheme, where the corrector is determined self-consistently by iteration. For long times after the quench, the numerical effort is equivalent to a second-order Runge-Kutta scheme because the self-consistency of the predictor-corrector scheme converges very fast. Right after the quench, our approach reduces numerical errors significantly and is thus advantageous.

In order to interpret the long time behavior Green’s functions, we briefly comment on properties of thermal Green’s functions. In thermal equilibrium all Green’s function only depend on $\tau$. Further, the imaginary part of the Fourier transformed retarded Green’s function is given by the spectral function

[TABLE]

In thermal equilibrium the Kubo-Martin-Schwinger (KMS) conditions Kamenev (2013) establishes a connection between $G^{>}(\omega)$ and $G^{<}(\omega)$ which is given by

[TABLE]

The last relation can be proved by using the Lehmann representation. Employing the definition of $G^{R}$ we can extract the spectral function

[TABLE]

from the Green’s functions $G^{>}$. Similarly the Keldysh component of the fermionic Green’s function is related to the spectral function,

[TABLE]

Out of equilibrium the Green’s functions depend on $\tau$ and $\mathcal{T}$. However, we can still consider the Fourier transform with respect to $\tau$ given by

[TABLE]

which is also known as the Wigner transform. From this object we obtain a spectral function out of equilibrium at time $\mathcal{T}$ as

[TABLE]

Numerically, we determine $G^{>}(t_{1},t_{2})$ and compute the spectral function abd the Keldysh component. In order to investigate thermalization behaviour, we can use a generalization of Eq. (47),

[TABLE]

where the quantities on the right hand side are obtained from the Green’s functions via a Wigner transformation. In a thermal state this equality holds with a time-independent effective inverse temperature $\beta$. In a non-thermal state, this relation allows one to quantify the deviation from a thermal state and to determine the timescale of thermalization if the final state is indeed thermal.

III.1 $(J_{2,i},0)\rightarrow(J_{2,f},0)$ quench: Rescaling of the Random-Hopping model

Figure 1 shows the spectral function of the random hopping model long before and long after a parameter quench. These results were obtained from a numerical Fourier transformation of the retarded Green’s function as described by Eq. (48) with a broadening of $\delta=0.025$. In Fig. 2, we show the frequency dependence of the Keldysh component of the fermionic Green’s function for the same quench protocol as in Fig. 1.

In Fig. 3, we show the ratio between $iG^{K}(\mathcal{T},\omega)$ and $A(\mathcal{T},\omega)$ long before and after the parameter quench. The data is only reliable for frequencies $|\omega|<J_{2,f}$. For these frequencies, $iG^{K}(\mathcal{T},\omega)/A(\mathcal{T},\omega)$ is almost flat – up to numerical artifacts. We can determine the inverse temperature $\beta$ using Eq. (50), and find that $\beta J_{2}$ is roughly constant during the quench. All three results of this $(J_{2,i},0)\rightarrow(J_{2,f},0)$ quench are consistent with a rescaling of energy scales for all quantities. This is expected because the random hopping model has only one energy scale $J_{2}$, and the analog of the reparameterization in Eq. (54) applies here. For analytical expressions of the spectral functions we refer to Appendix A. We further comment, that whenever an energy scale appears to be one, we measure all other physical quantities with respect to this scale.

III.2 $(0,0)\rightarrow(0,J_{4})$ quench: From bare Majoranas to the SYK model

In Fig. 4, we show the spectral function long after suddely switching on the quartic interaction in the SYK model, starting from bare, noninteracting Majorana fermions. For an analytical expression of the thermal spectral function we refer to Appendix B. When interpreting results of this quench protocol, it is important to keep in mind that the Hamiltonian before the quench is zero, so that any finite $J_{4,f}$ is an arbitrarily strong perturbation.

In this scenario we find that the energy of the system does not change during the quench. It is zero before and after the quench, i. e. the quench did not pump energy into the system. The Keldysh component of the Green’s function vanishes before and after the quench. This is to be expected for free Majorana fermions and is consistent with the temperature of the system being infinite, i.e. $\beta_{f}=0$ after the quench.

We interpret this as a result of the particle-hole symmetry of the system before the quench when there is only one energy level at $\omega=0$. The quartic interaction broadens this energy level, but does not break the particle-hole symmetry. Thus, the spectral weight of the fermions is distributed over the entire spectrum. Alternatively, one can argue that the Hamiltonian before the quench is zero, so that any added interaction is in fact infinitely strong.

III.3 $(J_{2,i},0)\rightarrow(0,J_{4,f})$: From a quadratic to a quartic model

The quench from the purely quadratic to the purely quartic model decouples the regions $t_{1}(t_{2})<0$ from $t_{1}(t_{2})>0$ as can be seen by inspecting the structure of the Kadanoff-Baym equations and of the self-energy, which reads

[TABLE]

for this quench protocol. Inserting this into the Kadanoff-Baym equations for $t_{1}>0$ and $t_{2}>0$, one can see that all time integrals are restricted to positive times and the initial condition does not matter: Although $G^{>}(t_{1},t_{2})$ shows some time evolution for $t_{1}\gtrless 0$ and $t_{2}\lessgtr 0$ due to integrals involving the region with $t_{1}<0$ and $t_{2}<0$, $G^{>}$ in this region does not influence the time evolution at positive $t_{1}$, $t_{2}$ because $\Sigma(t_{1},t_{2})=0$ when $t_{1}\gtrless 0$ and $t_{2}\lessgtr 0$. Thus, the relevant initial condition for the time evolution of $G^{>}(t_{1},t_{2})$ is $iG^{>}(t_{1}=0,t_{2}=0)=1/2$. Hence propagating $iG^{>}$ forward in time using the self-energy of the quartic model, we obtain the same time evolution at positive times as when starting from bare Majorana fermions.

III.4 $(J_{2,i},J_{4,i})\rightarrow(0,J_{4,f})$ quench: $J_{2}$ + $J_{4}$ model to the SYK model

In this section, we discuss quenches from a Hamiltonian with a quadratic and a quartic term to a Hamiltonian, where the quadratic term is switched off. Particularly, the time scale of thermalization is the main interest. We study the validity of the flucutation-dissipation relation and the behavior of the effective temperature $T_{\text{eff}}(\mathcal{T})=\beta_{\text{eff}}^{-1}(\mathcal{T})$ in order to detect thermalization. The effective temperature can either be obtained from the derivative of $iG^{K}(\mathcal{T},\omega)/A(\mathcal{T},\omega)$ at $\omega=0$ or as fit of $\operatorname{tanh}[\omega/(2T_{\text{eff}})]$ to that function to some frequency interval. As shown in Fig. 5, long after the quench the numerical results for $iG^{K}(\mathcal{T},\omega)/A(\mathcal{T},\omega)$ are well described by $\operatorname{tanh}[\omega/(2T_{\text{eff}})]$, indicating that the system has thermalized. At intermediate times, the low frequency behaviour is still well described by this function. At high frequencies, deviations are visible which could arise due to the fact that $iG^{K}(\mathcal{T},\omega)/A(\mathcal{T},\omega)$ is non-thermal.

In Fig. 6 we show the time-dependence of the effective temperature and a fit of $T_{\text{eff}}^{-1}$ to Eq. (8). After the quench, the effective temperature shows exponential behavior. In Fig. 7 we show the thermalization rate as obtained from such fits. It appears proportional to the final temperature $T_{f}$ at small final temperatures as noted in Eq. (9), and saturates at higher temperatures. It is difficult to determine $\Gamma$ numerically at large $T_{f}/J_{4}$, and this likely leads to the oscillatory behavior present in Fig. 7.

IV Kadanoff-Baym equations: Large $q$ limit

This section will consider a model with a $q$ fermion coupling $J(t)$, and a $pq$ fermion coupling $J_{p}(t)$. For now, we keep the time-dependence of both couplings arbitrary. We will find that the initial state is exactly solvable for $p=2$ and $p=1/2$.

The Kadanoff-Baym equations in Eqs. (40), (41) and (42) for $G^{>}(t_{1},t_{2})$ now become

[TABLE]

where we have defined

[TABLE]

Also, recall that Eq. (31) relates $G^{>}$ to $G^{<}$.

We note a property of the Kadanoff-Baym equations in Eq. (52), connected to a comment in Section I above Eq. (6). If we set $J_{p}(t)=0$, then all dependence of Eq. (52) on $J(t)$ can be scaled away by reparameterizing time via

[TABLE]

This implies that correlations remain in thermal equilibrium in the new time co-ordinate. However, when both $J_{p}(t)$ and $J(t)$ are non-zero, such a reparameterization is not sufficient, and there is non-trivial quench dynamics, as was shown by our numerical study in Section III. Below, we will see that the quench dynamics can also be trivial in the limit $q\rightarrow\infty$, even when both $J_{p}(t)$ and $J(t)$ are both non-zero. But this result arises from fairly non-trivial computations which are described below, and in particular from an SL(2,C) symmetry of the parameterization of the equations.

In the large $q$ limit, we assume a solution of the form in Eq. (10). Inserting (10) into (52), we obtain to leading order in $1/q$

[TABLE]

where

[TABLE]

It is a remarkable fact that these non-linear, partial, integro-differential equations are exactly solvable for our quench protocol and all initial conditions, as we will show in the remainder of this section. The final exact solution appears in Section IV.4, and surprisingly shows that the solution is instantaneously in thermal equilibrium at $t=0^{+}$.

Taking the derivatives of either equation in Eq. (IV) we obtain

[TABLE]

We will look at the case where $\mathcal{J}_{p}(t)$ is non-zero only for $t<0$, while $\mathcal{J}(t)$ is time independent:

[TABLE]

Then Eq. (57) is the Lorentzian Liouville equation in three of the four quadrants of the $t_{1}$-$t_{2}$ plane (labeled as in Fig. 8), and its most general solution in these quadrants is

[TABLE]

Eq. (59) will also apply for $\alpha=C$ when $\mathcal{J}_{p}=0$. Eq. (11) implies that the functions $h_{A1}(t)$ and $h_{A2}(t)$ obey

[TABLE]

and similarly

[TABLE]

Note that $g_{\alpha}(t_{1},t_{2})$ remains invariant under the SL(2,C) mapping in Eq. (14). Given the symmetries of Eq. (IV) and the thermal initial conditions, we look for solutions which obey

[TABLE]

This property can be related to the causality of the Kadanoff-Baym equations as is explained in Appendix D. We can satisfy Eq. (62) in quadrants A, B, D by the relations

[TABLE]

The system is in equilibrium in quadrant C, and so the Green’s function depends only upon time differences. We describe the nature of this equilibrium solution in Appendix E. In the following subsections, we describe the results of inserting the parameterization in Eq. (59) back into Eq. (IV) to obtain ordinary differential equations for $h_{\alpha 1}$ and $h_{\alpha 2}$.

IV.1 Quadrant B

From Eq. (IV), we obtain for $t_{1}<0$ and $t_{2}>0$

[TABLE]

Inserting Eq. (59) in Eq. (64) we obtain

[TABLE]

We also have the compatibility condition at the boundary between regions B and C, which is

[TABLE]

taking the derivative of this equation we obtain precisely the first equation in Eq. (65) after using Eq. (109). This reassures us that the system of equations are not overdetermined. We can easily integrate Eq. (67) to obtain

[TABLE]

IV.2 Region A

Now $t_{2}>t_{1}>0$. From the second Eq. (IV) we obtain

[TABLE]

Inserting Eq. (59) into Eq. (69), we obtain

[TABLE]

From the first Eq. (IV) we obtain

[TABLE]

Inserting Eq. (59) into Eq. (71), we obtain

[TABLE]

From Eqs. (IV.2) and (IV.2) we obtain

[TABLE]

This is precisely the logarithmic derivative of the Majorana condition in Eq. (60). The compatibility condition at the boundary of region A and region B is

[TABLE]

This can be integrated to

[TABLE]

where $c_{4}$ is a constant of integration.

IV.3 Combined equations

We adopt a simple choice to solve Eq. (75)

[TABLE]

Then we find that Eqs. (IV.1) and (IV.2) are consistent with each other. Collecting all equations, we need to solve

[TABLE]

It can be verified that all expressions above are invariant SL(2,C) transformations of the $h_{\alpha 1}$ and $h_{\alpha 2}$ fields. Given the values of $h_{B1}(-\infty)$, $h_{B2}(0)$, $h_{B2}^{\prime}(0)$, and $g_{C}(t)$, Eqs. (79-82) determine the values of $h_{A1}(0)$, $h_{A1}^{\prime}(0)$, and $h_{A2}(0)$. Then Eqs. (77,78) uniquely determine $h_{A1}(t)$ and $h_{A2}(t)$ for all $t\geq 0$. Because of the SL(2,C) invariance, the values chosen for $h_{B1}(-\infty)$, $h_{B2}(0)$, $h_{B2}^{\prime}(0)$ won’t matter for the final result for $g_{A}(t_{1},t_{2})$.

IV.4 Exact solution

A solution of the form in Eq. (123) does not apply to Eqs. (77-82) in region A because it does not have enough free parameters to satisfy the initial conditions. However, we can use the SL(2,C) invariance of Eqs. (77-79) to propose the ansatz similar to that in Eq. (15) (here, we redefine constants by factors of $e^{i\theta}$):

[TABLE]

We can now verify that Eq. (83) is an exact solution of Eqs. (77-82). This solution is characterized by 4 complex numbers $a$, $b$, $c$, $d$ and two real numbers $\theta$, $\sigma$. These are uniquely determined from the values of $h_{B1}(-\infty)$, $h_{A1}(0)$, $h_{A1}^{\prime}(0)$, and $h_{A2}(0)$ by the solution of the following 6 equations

[TABLE]

It is now easy to verify that $g_{A}$ takes the form in Eq. (17), and so all of quadrant A is also in thermal equilibrium. We also numerically integrated Eqs. (77-82), starting from generic initial conditions, and verified that the numerical solution obeyed the expressions in Eqs. (83) and (84).

The above solution determines the values of $\sigma$ and $\theta$ in the final state, and hence the value of the final temperature via Eq. (18). In general, this will be different from the value of the initial temperature in quadrant C.

Eqs. (83) and (84) also determine the solutions in the other quadrants via expressions specified earlier. The solutions in region $\widetilde{A}$ follow from the conjugacy property in Eq. (63). In quadrant B, we have $h_{B2}(t)=h_{A2}(t)$ in Eq. (76), while $h_{B1}(t)$ was specified in Eq. (68) using the initial state in quadrant C. Note that the solution in quadrant B is not of a thermal form, as $h_{B1}$ and $h_{B2}$ are not simply related as in Eq. (83). The solution in quadrant D follows from that in quadrant B via the conjugacy property in Eq. (63). And, finally, the initial state in quadrant C was described in Appendix E.

V Conclusions

Quantum many-body systems without quasiparticle excitations are expected to locally thermalize in the fastest possible times of order $\hbar/(k_{B}T)$ as $T\rightarrow 0$, where $T$ is the absolute temperature of the final state Sachdev (1999). This excludes e.g. the existence of systems in which the local thermalization rate, $\Gamma\sim T^{p}$ as $T\rightarrow 0$ with $p<1$, and no counterexamples have been found.

In this paper we examined SYK models, which are systems which saturate the more rigorous bounds on the Lyapunov time to reach quantum chaos Maldacena et al. (2016). Our numerical study of the model with a final Hamiltonian with $q=4$ showed that this system does thermalize rapidly, and the thermalization rate is consistent with $\Gamma=CT$ at low $T$ where $C$ is dimensionless constant, as indicated in Eq. (9) and Fig. 7.

We also studied a large $q$ limit of the SYK models, where an exact analytic solution of the non-equilibrium dynamics was possible. Here we found that thermalization of the fermion Green’s function was instantaneous. It will be necessary to study higher order corrections in $1/q$ to understand how this connects to the numerical $q=4$ numerical solution: does the constant $C\rightarrow\infty$ as $q\rightarrow\infty$, or (as pointed to us by Aavishkar Patel) is thermalization at large $q$ a two-step process. In two-step scenario, a very rapid pre-thermalization (which we have computed) is followed by a slower true thermalization of higher order corrections.

Finally, we comment on a remarkable feature of the large $q$ solution given by Eq. (13) and (15): its connection with the Schwarzian. The Schwarzian was proposed as an effective Lagrangian for the low energy limit of the equilibrium theory. Specifically, consider the Euler-Lagrange equation of motion of a Lagrangian, $\mathcal{L}$, which is the Schwarzian of $h(t)$

[TABLE]

The equation of motion is

[TABLE]

It can now be verified that the expressions for $h_{1,2}(t)$ in Eq. (15) (and Eq. (83)) both obey Eq. (86). We note, however, that we did not obtain Eq. (15) by the solution of Eq. (86): instead, Eq. (15) was obtained by the solution of the Schwinger-Keldysh equations of the large-$q$ SYK model in Eq. (IV). This connection with the Schwarzian indicates that gravitational models Sachdev (2010a, b); Almheiri and Polchinski (2015); Kitaev (2015); Maldacena and Stanford (2016); Jensen (2016); Engelsöy et al. (2016) of the quantum quench in AdS2, which map to a Schwarzian boundary theory, exhibit instant thermalization as in the large $q$ limit.

Indeed, the equation of motion of the metric in two-dimensional gravity Almheiri and Polchinski (2015) takes a form identical to that for the two-point fermion correlator in Eq. (12). And studies of black hole formation in AdS2 from a collapsing shell of matter show that the Hawking temperature of the black hole jumps instantaneously to a new equilibrium value after the passage of the shell Engelsöy et al. (2016); Engelsöy (2017). These features are strikingly similar to those obtained in our large $q$ analysis. There have been studies of quantum quenches in AdS2, either in the context of quantum impurity problems Erdmenger et al. (2017), or in the context of higher-dimensional black holes which have an AdS2 factor in the low energy limit Withers (2016); it would be useful to analytically extract the behavior of just AdS2 by extending such studies.

We thank J. Maldacena of informing us about another work Kourkoulou and Maldacena (2017) which studied aspects of thermalization of SYK models by very different methods, along with connections to gravity on AdS2.

Acknowledgments

We would like to thank Wenbo Fu, Daniel Jafferis, Michael Knapp, Juan Maldacena, Ipsita Mandal, Thomas Mertens, Robert Myers, Aavishkar Patel, Stephen Shenker, and Herman Verlinde for valuable discussions. JS was supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE1144152. AE acknowledges support from the German National Academy of Sciences Leopoldina through grant LPDS 2014-13. AE would like to thank the Erwin Schrödinger International Institute for Mathematics and Physics in Vienna, Austria, for hospitality and financial support during the workshop on ”Synergies between Mathematical and Computational Approaches to Quantum Many-Body Physics”. VK acknowledges support from the Alexander von Humboldt Foundation through a Feodor Lynen Fellowship. This research was supported by the NSF under Grant DMR-1360789 and MURI grant W911NF-14-1-0003 from ARO. Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Economic Development & Innovation. SS also acknowledges support from Cenovus Energy at Perimeter Institute.

Appendix A Spectral function and Green’s functions of random hopping model

The random hopping model of Majorana fermions, *i. e.*the SYK Hamiltonian for $q=2$, can be solved exactly Maldacena and Stanford (2016). We repeat the calculation within our notation in the following. The Matsubara Green’s function follows from the solution of

[TABLE]

This can be rewritten as a quadratic equation for \bigl{(}G(i\omega_{n})\bigr{)}^{-1} and matching the branches with the correct high-frequency asymptotics and symmetry yields the propagator

[TABLE]

Analytically continuing to the real frequency axis from the upper half plane yields the retarded Green’s function,

[TABLE]

where $\delta$ can be set to zero due to the presence of the finite imaginary part. Note that we exploited that $\operatorname{sgn}\omega_{n}\rightarrow 1$ in the analytic continuation from the upper half plane. We obtain

[TABLE]

which is the well-known semicircular density of states due to random hopping.

This yields

[TABLE]

where $J_{1}$ and $H_{1}$ are the Bessel function of the first kind (BesselJ in Mathematica) and the Struve function (StruveH in Mathematica), respectively.

Appendix B Spectral Function in the Conformal Limit

In the scaling limit at non-zero temperature the retarded Green’s function is given by the following expression

[TABLE]

where $\Delta=1/q$ is the fermion scaling dimension. At $q=4$ we obtain

[TABLE]

The Wigner transform of the retarded Green’s functions is given by

[TABLE]

where $b=({4\pi J^{2}_{4}})^{-{1}/{4}}$. The associated spectral function is

[TABLE]

Appendix C Details on the Numerical Solution of Kadanoff-Baym equation

The Green’s function is typically determined on two-dimensional grids in $(t,t^{\prime})$ space with $8000\times 8000$ or $12000\times 12000$ points, where the quench happens after half of the points in each direction. When starting from initial states in which only $J_{2}$ is finite, the Green’s function decays algebraically in time. This leads to significant finite size effects in Fourier transforms. In order to reduce the latter, most numerical results were obtained by starting from a thermal state in which $J_{2}$ and $J_{4}$, or only $J_{4}$, are finite. In these cases, the Green’s function decays exponentially as a function of the relative time dependence.

We checked the quality and consistency of the results by monitoring the conservation of energy, the normalization of the spectral function and the real and imaginary part of the retarded propagator are Kramers-Kronig consistent for long times after the quench.

In order to time-evolve the Kadanoff-Baym equations we have to determine $G^{<}(t_{1},t_{2})$ for $t_{1},t_{2}$. When quenching the ground state of the random hopping model, we can use the exact solution for $G^{>}$ as initial condition. When we do not have an analytical expression the Green’s function we solve the Dyson equation self-consistently according to the following scheme:

Prepare $iG^{>}$ with an initial guess, for example the propagator of the random hopping model. 2. 2.

Computation of retarded self-energy in time domain:

[TABLE] 3. 3.

Fourier transformation

[TABLE] 4. 4.

Dyson equation:

[TABLE] 5. 5.

Determine spectral function

[TABLE] 6. 6.

Determine $iG^{>}(\omega)$ from spectral function,

[TABLE]

Note that this is the only step in the self-consistency procedure where the temperature $\beta^{-1}$ enters through the Fermi function $n_{F}$. 7. 7.

Fourier transformation to time domain

[TABLE] 8. 8.

Continue with step 2 until convergence is reached.

Appendix D Conjugacy property of $g(t_{1},t_{2})$

In this appendix we illustrate that the causality structure of the Kadanoff-Baym equations eq. (IV) leads to the property

[TABLE]

in all four quadrants of Fig. 9. Since the system is thermal in quadrant $C$ the conjugacy property $[g(t_{1},t_{2})]^{\ast}=g(t_{2},t_{1})$ can be read off the thermal solution for $t_{1},t_{2}\leq 0$. Next, we consider the propagation from the line $\left\{(t_{1},t_{2})\in\mathbb{R}^{2}\,|t_{2}=0,t_{1}\leq 0\right\}$ for an infinitesimal time $\Delta t$ in the $t_{2}$ direction to the line $\left\{(t_{1},t_{2})\in\mathbb{R}^{2}\,|t_{2}=\Delta t,t_{1}\leq 0\right\}$. We discretize equation (IV) as follows

[TABLE]

and solve the resulting equation for $g(t_{1},\Delta t)$. Taking the complex conjugate leads to

[TABLE]

On the right hand side we are allowed to use the property $g(t_{1},t_{2})^{\ast}=g(t_{2},t_{1})$ since all $g$’s are still living in the $C$ quadrant. We obtain

[TABLE]

Next we let the system propagate from the line $\left\{(t_{1},t_{2})\in\mathbb{R}^{2}\,|t_{2}\leq 0,t_{1}=0\right\}$ for an infinitesimal time $\Delta t$ in the $t_{1}$ direction to the line $\left\{(t_{1},t_{2})\in\mathbb{R}^{2}\,|t_{1}=\Delta t,t_{2}\leq 0\right\}$. We discretize equation (52) as follows

[TABLE]

We set $t_{2}=t_{1}$ in the last equation and solve for $g(\Delta t,t_{1})$ leading to

[TABLE]

Comparing the right hand side of (104) with (107) we conclude $g(t_{1},\Delta t)^{\ast}=g(\Delta t,t_{1})$. Furthermore, the point $g(\Delta t,\Delta t)$ fulfills the conjugate property trivially. In total we propagated the conjugate property one time slice. Repeating this argument for every time slice of size $\Delta t$, the property will hold in all four quadrants of Fig. 9.

Appendix E Initial state in the large $q$ limit

In quadrant C in Fig. 8, all Green’s functions are dependent only on time differences, and and then we can write

[TABLE]

and Eq. (IV) as

[TABLE]

This implies the second order differential equation (also obtainable from Eq. (57))

[TABLE]

Eq. (110) turns out to be exactly solvable at $p=2$ (pointed out to us by Wenbo Fu, following Jian et al. (2017)). We write $g_{C}(t)=\ln(-1/f(t))$ and then Eq. (110) becomes

[TABLE]

The solution of this differential equation yields

[TABLE]

with

[TABLE]

By evaluating the Fourier transform of the Green’s function we find $G^{>}(-\omega)=e^{-\beta_{i}\omega}G^{>}(\omega)$, or by analytically continuing to imaginary time, we find that the initial inverse temperature is

[TABLE]

We will ultimately be interested in the scaling limit in which $\theta\rightarrow 0$ and $\sigma\ll\mathcal{J},\mathcal{J}_{2}$.

Eq. (110) is also exactly soluble at $p=1/2$ by

[TABLE]

where now

[TABLE]

The value of the inverse initial temperature remains as in Eq. (114).

Note that both solutions in Eqs. (112) and (115) obey

[TABLE]

It is also useful to recast the solution in quadrant C for the case $\mathcal{J}_{p}=0$ in the form of Eq. (59). We subdivide quadrant C into two subregions just as in quadrant A. From Eqs. (IV) and (59) we obtain for $t_{2}>t_{1}$

[TABLE]

Adding the equations in Eq. (118), we have

[TABLE]

which integrates to the expected

[TABLE]

So the final equations for the thermal equilibrium state are

[TABLE]

Unlike Eqs. (77,78), Eqs. (121,122) have to be integrated from $t=-\infty$. One solution of Eqs. (121,122) is

[TABLE]

with

[TABLE]

Note that the $g_{C}(t_{1},t_{2})$ obtained from this solution agrees with Eq. (112) at $\mathcal{J}_{2}=0$ and with Eq. (115) at $\mathcal{J}_{1/2}=0$.