Quantum Quenches and Thermalization in SYK models

We study non-equilibrium dynamics in SYK models using quantum quench. We consider models with two, four, and higher fermion interactions ($q=2, 4$, and higher) and use two different types of quench protocol, which we call step and bump quenches. We analyse evolution of fermion two-point functions without long time averaging. We observe that in $q=2$ theory the two-point functions do not thermalize. We find thermalization in $q=4$ and higher theories without long time averaging. We calculate two different exponents of which one is equal to the coupling and the other is proportional to the final temperature. This result is more robust than thermalization obtained from long time averaging as proposed by the eigenstate thermalization hypothesis(ETH). Thermalization achieved without long time averaging is more akin to mixing than ergodicity.


Introduction and Summary
The study of non-equilibrium dynamics is becoming important both in condensed matter physics [1][2][3][4][5][6][7][8][9] as well as in string theory [10][11][12][13][14]. One of the most interesting question in this field is to understand patterns of thermalization in the systems which are out of equilibrium. For example, it is important to know under what conditions a closed quantum system thermalizes, i.e., for a system prepared in a pure excited state, and undergoes unitary evolution, determine how the late time limit of the expectation values of certain observables are effectively described by a thermal ensemble 1 . Interest in the non-equilibrium dynamics from string theory point-of-view stems from black hole physics. The AdS/CFT correspondence(or the holographic principle, in general) says that a black hole corresponds to thermal ensemble in the boundary quantum theory, and the thermalization process in the quantum system is conjectured to be dual to black hole formation in the bulk gravitation theory.
On the bulk gravity side it has been conjectured that black holes are fast scramblers [15]. This proposal led to another conjecture [16] that the chaotic behaviour, that leads to scrambling, which is parametrized by the Lyapunov exponent λ L has an upper bound, and that upper bound is saturated by black holes. This naturally gave additional impetus to the study of non-equilibrium dynamics in systems which exhibit chaos, especially if the Lyapunov exponent of the theory saturates the upper bound.
The eigenstate thermalization hypothesis (ETH) is an attempt to explain how closed unitary quantum systems in pure excited states can thermalize [17,18]. Thermalization with ETH crucially involves long time averaging of the observables under consideration. It is, however, not clear what is the precise relation between chaos and ETH. In many studies of quantum systems, thermalization is observed even without long-time averaging [3]. Thermalization has also been seen in the integrable systems without long time averaging. The late time behaviour of integrable models is described by the generalized Gibbs ensembles [19,20]. These ensembles have fugacities turned on for several conserved charges of the integrable system. The integrable model, by definition, is not chaotic on its own.
The Sachdev-Ye-Kitaev(SYK) model which is a (0+1) dimensional model of Majorana fermions with all to all q-body random interactions. The q = 4 and higher models were studied by Kitaev[4], and by Maldacena and Stanford [21]. They showed that the out of time ordered four point correlators in these models saturate the upper bound on the Lyapunov exponent, in addition, these models also satisfy ETH. There has been a lot of work on this model, its variants and their bulk duals [4-6, 12, 21-48]. The q = 2 SYK model is not chaotic and also does not satisfy ETH. However, unlike integrable local quantum systems, it does not have infinite number of conserved charges in spite of them being exactly solvable. In fact, it has only one conserved charge which is the total energy of the system. With this background in mind, we study non-equilibrium dynamics of excited states in q = 2, 4, and higher SYK models.
The most convenient method for studying non-equilibrium dynamics, both theoretically [5-9, 11, 14, 49-51] and experimentally [52,53], turns out to be quantum quench. In other words, quantum quenches are the most convenient way of generating non-trivial excited states of the theory. In quantum quench one abruptly changes parameters of the Hamiltonian of the system starting from an equilibrium configuration(generally a thermal state or the ground state) of the system. The change in the coupling generally excites the system and the system evolves non-trivially with the final Hamiltonian. The evolution of the system is examined by calculating the expectation values of some of the observables of the system. If the expectation values of those observables approach the expectation values in a thermal ensemble, the system is said to have thermalized.
Certain aspects of quantum quenches in SYK models have been studied in [5]. In this paper, using similar numerical techniques, we will study quantum quenches in q = 2, 4, and higher SYK models. We will consider one particular observable which is the greater Green's function G > (t 1 , t 2 ).
For majorana fermions, all other two-point functions can be calculated from G > (t 1 , t 2 ). The non-trivial time evolution of G > (t 1 , t 2 ) can be examined by exactly solving its equations of motion which are the Kadanoff-Baym(KB) equations. Our analysis will involve changing various parameters with two different kinds of time dependence. The usual quench protocol in condensed-matter literature is changing, suddenly 2 or smoothly but rapidly, the parameters from one value to another different value. We will consider sudden change from one value to another, which we call step quench. In addition to this, we will also study bump quench, in which the coupling changes for a finite time interval before returning back to the original value 3 . We follow the convention q = k quench when the final hamiltonian of the system has k fermion interaction and the couplings J q undergo quench with q = k. We will also consider only sudden limit for both step and bump quenches.
The quenches which are relevant for our main results are: • Quenches in q = 2 theory: We use four, six and eight fermion interactions (J 4 , J 6 , and J 8 couplings) separately to quench the system for both step and bump protocols.
• Quenches in q = 4 theory: For this theory, we use J 2 (two fermion interaction coupling), J 6 , and J 8 with both step and bump protocols.
For quenches in q = 4 theory, we start from finite temperature thermal states which reduce finite-size effect drastically and ensure good numerical accuracy. For quenches in q = 2 theory, we start from the ground state as well as from finite temperature states. We observe that if we take the initial thermal states to be of sufficiently low temperatures, the effect of the initial temperature becomes insignificant. We consider only the greater Green's function G > (t 1 , t 2 ), because for the Majorana fermions, all other two point functions can be expressed in term of G > (t 1 , t 2 ) alone.
The main technical results of our analysis are as follows: • In q = 2 theory, the two point functions do not thermalize in all the quench scenario. But an interesting observation is that the two point functions equilibrate instantaneously as soon as both the time arguments are outside the quench region.
• In q = 4 theory, the two point functions thermalize for all the quench scenario. G > (t 1 , t 2 ) converges exponentially towards its equilibrium expectation value. This exponential behaviour is observed as soon as both the time arguments are outside the quench region.
• In q = 4 theory, we also identify two exponents, of which, one is equal to the coupling and the other is proportional to the final temperature. The first one is the exponent of G > (t − t a , t) as a function of t with t a fixed, while the other is the exponent of G > (t, t b ) as a function of t with t b fixed.
• We compute the thermalization rate of the effective temperature in both step and bump quench. We show that the thermalization of the effective temperature fits an exponential ansatz. We also find that the thermalization rate is independent of the coupling constant.
An important aspect of the present work is to check if step quenches produce special fine-tuned pure states which looks exactly thermal. These pure states are inspired by the Euclidean evolved boundary states of Calabrese and Cardy [56]. These states, which we will refer to as Kourkoulou-Maldacena (KM) states below, have interesting bulk duals [30]. The details of these pure states can be found in section 2.3. We observed that the final states of quantum quenches using disordered couplings are not KM states. But one can use mass like terms to perform the sudden step quenches for which the final states are the KM states.
The thermalization we observe in q = 4 theory without long time averaging, is much more robust than what one expects from the ETH. We therefore believe that thermalization in a chaotic system is more akin to mixing in classical systems which is a stronger condition than ergodicity.
The outline of this paper is as follows: In section 2, we will briefly recall the SYK model. This will also be used to fix our notation. We will write down the Schwinger-Dyson equation for a model with both q = 2 and 4 interactions. The couplings for these terms will have arbitrary time dependence to start with. We will then set up the Kadanoff-Baym equations for this system which can be easily generalized for higher q models. Finally we will briefly discuss the eigenstate thermalization hypothesis(ETH). In section 2.3, we discuss Kourkoulou-Maldacena states with an eye on possible relation between our results and these excited states. In section 3, we discuss various quench protocols that we study in the SYK model and present results of our numerical computations. Section 4 contains conclusion and discussion where we wrap up our results and discuss about ways to prepare Kourkoulou-Maldacena states and the implications of thermalization without long-time averaging.

The SYK models
We begin with a review the model studied by Sachdev et al. [5]. This will help set up notation for subsequent sections. Our starting point is the SYK model with the hamiltonian H that contains q-point (q even) interaction between N Majorana fermions, The coupling J i 1 ,i 2 ,..,iq is random with gaussian distribution, with vanishing mean value and the width of the gaussian is given by To compute correlators at finite temperature the Schwinger-Keldysh formalism is employed in which, the observables are computed by integrating along the closed-time contour C. The initial state is evolved along this contour both forward and backwards in time. The contour-ordered Greens function is defined as [5], The correlation function in the path integral formalism is computed by inserting the components of fields on the forward and return path of the contour. The components of the matrix Green's functions that we will be interested in are called greater (lesser) Green's functions, denoted as G >(<) (t 1 , t 2 ), and are defined in the following manner 4 (2.4) 4 We use the commutation relation {ψ i , ψ j } = δ ij . So, G > (t, t) = −i/2 and G < (t, t) = i/2.
where by t + i we mean t i on the upper contour and t − i denotes t i on the lower contour, and the contracted index i simply denotes a sum over i. The relative minus sign above is due to swapping of the position of two Majorana fermions under contour ordering. From the above definitions, for Majorana fermions, This relation holds even for non-equilibrium dynamics [5,57]. This model exhibits conformal symmetry in the infrared which is spontaneously broken by the h = 2 mode, where h is the quantum number of the SL(2) subgroup of the conformal symmetry. This h = 2 mode has chaotic behaviour for q ≥ 4. It turns out that the h = 2 mode saturates the chaos bound λ L = 2π/β [16]. The model with only q = 2 term, however, does not have chaotic behaviour. This is clearly due to the quadratic nature of the action and as a result the model is integrable. We are interested in studying the SYK model with time dependent coupling which can exhibit different behaviour by virtue of having the coupling as a function of time.
Our main object of interest is the Kadanoff-Baym equations which we will use to analyse the non-equilibrium dynamics of the SYK model. Before we set up the Kadanoff-Baym equations, let us consider the Schwinger-Dyson equation.

The Schwinger-Dyson(SD) equations
We will consider the time dependent Hamiltonian which describes different quench protocols depending on the kind of time dependence we allow for the couplings of the theory. To simplify the matter we will extract the time dependence of the couplings and write it in terms of separate functions of time. For example, up to the quartic fermion interaction i.e., q = 4, the Hamiltonian is where, f 2 (t) and f 4 (t) contain the time dependence of the couplings. The partition function of this model is written in terms of the action functional, (2.7) All the interaction terms in the SYK model couple all fermions to each other and have random couplings. The randomness of the coupling is meant to mimic the disorder in the system. We will average the partition function over the gaussian distributed random couplings, where the gaussian weight functions, P 1 (J 2,ij ) for the quadratic coupling and P 2 (J 4,ijkl ) for the quartic coupling have width 2J 2 2 /N and 12J 2 4 /N respectively. Usually in the quenched disorder the integration over the random variables is carried out at the end of the computation, however, in the large N limit we can reverse the order. Carrying out the gaussian integral over the quadratic and quartic couplings gives us the effective action (2.9) In this effective action the sum runs over all values of i, j, k, l and the combinatoric factors take care of the ordering of fermions in each term. Following [5], we will write this effective action in terms of auxiliary fields and convert it into a quadratic action in terms of the fermions. The path integral in terms of the auxiliary functions, suggestively named as G(t) and Σ(t), where, The auxiliary field Σ is introduced so that we can implement the constraint (2.11) as an equation of motion of Σ. This is done by implementing the constraint through the δ-function. This procedure reduces the action (2.10) to quadratic form in terms of the fermions. We can now integrate out the Majorana fermions and write the effective action S[G, Σ] purely in terms of G and Σ, (2.12) An advantage of this form of the effective action is that the Schwinger-Dyson equations can be derived as equations of motion of this action, A similar analysis can be carried out for the six and higher fermion interactions in an analogous manner. Let us now consider the eq.(2.14) and take the convolution product with G(t 1 , t 2 ) from both right and left, this procedure gives us two equation, To study the Kadanoff-Baym equations besides eq. (2.15), (2.16) we will need the retarded, the advanced and the Keldysh Green's functions which are defined as Along these lines define the retarded, advanced self-energy in the following manner.
In the next subsection we will use these ingredients to derive the Kadanoff-Baym equations.

The Kadanoff-Baym (KB) equations
Equations (2.15) and (2.16) can be manipulated using the real space representation of G −1 0 on the left hand side and contour deformation on the right hand side to write Note that the contour starts from some time t 0 and the operators are inserted in the correct order for different values of t 1 and t 2 and then comes back to t 0 . For quenches starting from a thermal state, the contour further goes down in the imaginary time direction for an interval of length β i which is the inverse temperature of the initial thermal state ( Figure 1).
If one takes the limit t 0 → −∞ then for all observables at finite time, the contribution from the imaginary time interval can be neglected which follows from the Bogoliubov principle of weakening correlations [58]. 5 We will briefly explain derivation of (2.22) using the Langreth rules below. Derivation of (2.23) follows in an analogous manner. The left hand side of (2.22) can be derived starting from the equation(2.15), and choosing the Green's function G(t 3 , t 2 ) to be the greater Green's function G > (t 3 , t 2 ), and integrating by parts to get where we have used the fact that G −1 0 is given by the derivative of the δ-function. The right hand side of (2.15) is Using the contour deformation we can rewrite (2.25) as The first term in (2.26) can be written as where,τ = t 1 − τ . Inserting Heaviside Θ(τ ) function in the term involving Σ < we can extend the integration limit from (0, ∞) to (−∞, ∞). After substitutingτ = t 1 − τ , the integral remains invariant. So we get, (2.28) Similar manipulations can be carried out for the second term in (2.26) to get,

Eigenstate Thermalization Hypothesis
It has been shown that the q = 4 SYK model with Majorana fermions [13,25] and complex fermions [24] with large but finite N satisfy the eigenstate thermalization hypothesis (ETH). Although it has been claimed [59] that q = 2 SYK model with complex fermions satisfies ETH, it was later found that the finite N scaling in q = 2 SYK model with Majorana fermions does not scale correctly with the system size [25]. It has therefore been suggested that q = 4 SYK model should thermalize while the q = 2 model should not. Our results do not conflict with this suggestion, however, note that ETH necessarily involves long-time averaging of the observables [17,18,60]. Long time averaging is not necessary for thermalization or equilibration in many scenario of quantum quenches [3], even in free theories [14]. In fact, it is not even clear what is the relation of ETH with such thermalization or equilibration processes which do not involve long-time averaging after quantum quenches. Also note that in black hole collapse geometries [54,55,61], there is no long-time averaging invloved. These geometries are the bulk duals of thermalization in the corresponding boundary CFT.

Kourkoulou-Maldacena states and Instantaneous thermalization
In this section we will introduce certain pure excited states in SYK models. The motivation for constructing these states comes from the boundary state ansatz of quantum quenches in 1D systems in the thermodynamic limit [56]. The ansatz by Calabrese and Cardy corresponds to starting from the ground state of a gapped theory and quenching it to a gapless theory (1+1D CFT), the final state obtained after the quench has the generic form where κ > 0 is a parameter fixed by the quench process, H CF T is the Hamiltonian of the final gapless theory and |B is a conformally invariant boundary state (B state) of the CFT. We will refer to these states as Calabrese-Cardy(CC) states. Determination of the particular B state that is relevant for the description of the post quench state of the system for a specific quantum quench is a non-trivial problem [62]. Nevertheless, using conformal symmetry of the final theory, it can be shown that expectation values of one-point and two-point functions effectively thermalize, where the expectation values in the long-time limit are described by a thermal ensemble with inverse temperature β = 4κ. In fact, it has been shown that finite subsystems thermalize where again the long-time limit is described by a thermal ensemble with inverse temperature β = 4κ [7,10]. Since the quench process started from the ground state, the system always remains in a pure state. An interesting aspect of this process of thermalization of subsystems is that correlation functions of holomorphic operators of the final CFT thermalize instantaneously [14,63].
We will now consider certain pure excited states in SYK models. These states were first constructed by Kourkoulou and Maldacena in [30]. Considering N majorana fermions, the analogous B states are defined as Hence, there are 2 N/2 number of such B states. These are high energy states. One can produce lower energy states by evolving these B states for a finite euclidean time κ.
We will refer to these low energy states as KM states.
An interesting feature of KM states is that, in the large N limit, "diagonal" twopoint functions ψ i (t 1 )ψ i (t 2 ) are "instantaneously thermalized"(using the 1+1D CFT terminology used above) where the effective inverse temperature β = 2κ. The "off-diagonal" two-point functions ψ 2k−1 (t 1 )ψ 2k (t 2 ) have non-trivial time dependence and decay to zero in the long-time limit. These "off-diagonal" two-point functions are zero in a thermal ensemble. The KM states also have interesting bulk duals in AdS 2 . Unlike in 2D CFT quenches, we could not find any quench scenario with disordered couplings where the final state is the KM state. This work was initially inspired by our curiosity about the possibility of the KM states being the final states of step quenches but not for bump quenches in SYK models. The negative result that the final states in quenches in SYK models are not KM states leads to deeper understanding of the thermalization process in chaotic theories. We will comment further on this issue in the concluding section 4.

Quantum Quenches in SYK models
The KB equations are solved numerically after discretizing the two time arguments t 1 and t 2 . For quenches in q = 2 theory, we could start from the ground state, since the Green's function oscillates and decays fast with time. For all other cases, we start with a thermal state which gives an exponential decay of the initial data as a function of the relative time difference. Moreover, since we start from a stationary state, all the initial data in the third quadrant are shifted functions of the data on (t 1 < 0, t 2 = 0) line and (t 1 = 0, t 2 < 0) line. We use a grid of the kind bounded by red coloured lines in figure 3. Since the terms far away from the diagonal fall of exponentially fast, the grid points in the second and fourth quadrant lying outside the red coloured lines are ignored in our numerical code. We used grids of three different sizes 2001 × 1001, 3001 × 1501 and 4001 × 2001 points. The computation time grows very fast with increasing grid size. We also used a fixed time step size dt = 0.05. 6 In the rest of the paper, we will suppress factors of this time step size dt. So, unless it is explicitly mentioned all the times are measured in units of dt. In step protocols, the quenches happen at t 1 = 0 and t 2 = 0. For all the cases with bump protocol, the perturbations 7 are turned on between t 1 = 1 and t 1 = 10, similarly between t 2 = 1 and t 2 = 10 for the other direction. The KB equations are solved selfconsistently in this grid using the Predictor-Corrector method. The predicted values on line A are calculated causally from the data on line B as shown in figure 3. The predicted values are then corrected until the desired accuracy is obtained.
For most of quenches we are considering here, the initial data is obtained by solving the SD equation numerically for finite inverse temperature β [5]. For step quenches in q = 2 theory in which J 2 interaction is dominant, we can start from the ground state. The initial data are obtained by solving the SD equation in the ground state (β → ∞) numerically. In this case we use In case of the bump quench in q = 2 theory, for cases in which we start from the ground state, the initial data is calculated using the analytic expression for G > (t 1 , t 2 ). The greater Green's function in ground state for q = 2 theory is Calculation of final temperature: The temperature in the long time limit is calculated using the relation [5] iG K (ω) which is a function of only t 1 − t 2 in a thermal ensemble and A(ω) is and it holds for all fermionic theories. We can therefore conclude that the system under consideration has thermalized only if the quantity on the LHS of (3.3) has tanh profile as a function of the frequency ω. Note that for the determination of the final temperature we also have to use the relation between greater and lesser Green's functions (2.5).
Check for energy conservation: We also check for energy conservation to ensure that our numerical results are correct. From (2.10), the total energy as a function of time t 1 is given by In the second line, the first term arises from the upper half of the contour and the second term arises from the lower half of the contour. We have also used (2.5) for the second term.
The quench processes we are considering, merely satisfying (3.3) in the long time limit is not sufficient to guarantee thermalization. This is because, as we mentioned above, all fermionic theories at finite temperature satisfy the relation (3.3). So, to check thermalization, we first calculate the final temperature using the above relation. The SD equation of the final theory is then solved at the calculated final temperature and in the end we check if the generated real time two-point functions agree with the two-point functions obtained from the quench process.

Quenches in q = 2 SYK model
In this subsection we will study quantum quenches in which the final theory is the q = 2 SYK model, that is the model which only has 1-body (quadratic, J 2 ) interaction. These quenches are special cases because the two-point functions equilibrate instanteneously. From (2.22, 2.23), for q = 2 final theory, This is observed in our numerical solutions of the KB equations below. However, note that the instanteneously equilibrated configuration is not a thermal ensemble, so the final state cannot be a KM state. Since, the initial theory is J 2 dominant(for step quench) or a q = 2 theory, we can start the quench from the corresponding ground state. We will present here only cases in which J 4 interaction is used to perform both step and bump quenches. We also found similar results for quenches using J 6 and J 8 interactions, as we expect from (3.7). The results are qualitatively similar for quenches starting from thermal state.
The value of the J 2 coupling is always fixed at 1. We will present results for step quench with initial J 4 = 2 which is suddenly turned off at time t = 0. For bump quench, we turn on J 4 = 5 for a time duration of 9 × dt = 9 × 0.05 = 0.45 from time step t = 1 to t = 10. This same quench parameters are used for all quenches starting from different initial temperatures including the ones starting from ground state.
The step quench happens at t = 0, the two time arguments of G > (t − 100, t) are outside the quench region if t > 100. The bump quench happens between t = 0 and t = 11 so the two time arguments are outside the quench region if t ≥ 111. Figure  (4) are plots of the real and imaginary parts of G > (t − 100, t) as a function of time t for step and bump quenches starting from ground states. One can see that the Green's function freezes or equilibrates instantaneously once the two time arguments are outside the quench regions. But the equilibrated value is different from the thermal expectation value. Figure (5) compares iG K (ω)/A(ω) with tanh(β f ω/2) for step and bump quenches starting from initial inverse temperature β = 10.

Quenches in q = 4 SYK model
In this subsection we will consider quantum quenches in which the final theory is q = 4 SYK model which only has 2-body (quartic, J 4 ) interaction. We will present results for which the interaction terms used for the quench process is J 2 . We also found similar results for quenches with J 6 and J 8 interactions. For the initial thermal states, we considered three different inverse temperatures β i = 10, 20, and 30. We find that increasing the inverse temperature from 20 to 30 does not affect the results much. This is expected since for a fairly large β, the fermion distribution function is well represented by the step function (3.1). So, we expect that the quench starting from β = 20 and 30 should also be qualitatively similar and quantitatively close to the quenches starting from ground states.
Three different values of J 4 are used, namely, 0.5, 1 and 1.5. For step quenches, we start from a theory with J 4 and J 2 . At t = 0, the J 2 coupling is suddenly changed to 0. For the bump quenches, starting from a theory with only J 4 , J 2 is turned on for a time duration of 9 × dt = 9 × 0.05 = 0.45 from time step t = 1 to t = 10. As mentioned above, we will use this time interval for all bump quench protocol. Changing this time interval does not affect our main results. Longer time interval only injects more energy into the system resulting in higher final temperature. (a) Once both the time arguments are outside the quench region, we find that the greater Green's function thermalizes rapidly but not instantaneously, as can be seen in Figure(6). Figure (7a, 7b) are two resolved plots of G > (t − 100, t) for different initial inverse temperatures as a function of t for step quenches. Since the step quench happens at t = 0, both the time arguments are outside the quench region if t > 100. Immediately after time t crosses 100, G > (t − 100, t) changes rapidly and exponentially towards its equilibrium thermal value. The evolutions for t > 100, both real and imaginary parts, fit exponential functions very well. The two exponents of the two exponential fits for real and imaginary parts are roughly equal. This behaviour is not a numerical artifact. The exponents do not change with change in time step size. We have checked for different time step sizes dt = 0.05 and dt = 0.025. Moreover, we have also checked energy conservation using (3.6).  Similarly, for bump quenches in Figure (7c, 7d), once the two time arguments are outside the quench region, the Green's function thermalizes rapidly and its real and imaginary parts fit exponential functions very well. Below, we will consider only the exponent for the imaginary part which we will denote by γ Itt .
The bump quench happens between time steps t = 0 and t = 11, so the two time arguments of G > (t − 100, t) are outside the quench region if t ≥ 111. One of the most interesting numerical result of this work is that we find that This can be seen from Fig. (8) and Table 1.
We also check if the final stationary limit is described by a thermal ensemble. For which we compare iG K (ω)/A(ω) with tanh(β f ω/2) for some final temperature β f . Figure (9a, 9b) are two such comparisons. Figure ( Since we observe thermalization, another observable of interest is G > (t, t 2 ) where t 2 is fixed. In the hydrodynamics limit [54] of large t, both the real and the imaginary parts of the expectation value of this observable are again exponential functions with both the exponents equal. We will consider the exponent of the imaginary part which we denote by γ It . This exponent is equal to the exponent of the retarded Green's function in a thermal ensemble with temperature equal to the temperature of the final thermalized limit of the quench process. We will denote the exponent of the retarded Green's function by γ ret .
At low temperature, γ It is proportional to the final temperature.
In a thermal ensemble, the retarded Green's function is a function of the relative time difference. In the conformal limit of SYK model, the retarded Green's function in a thermal ensemble of inverse temperature β is where ∆ = 1/q = 1/4 and b = (4πJ 2 4 ) −1/4 . In the conformal limit, the exponent is .  Figure 10: The exponent γ It as a function of γ conf = π/(2β f ). γ ret is exactly equal to γ It as we can see from Table 1 so γ ret is not plotted here.
At high temperatures, we find that the exponent of G R (t) gets significant correction compared to its value at the conformal limit. The corrected value of the exponent, which we have denoted by γ ret above, is calculated by solving the SD equation numerically. Important numerical results for the step and bump quenches with J 2 , starting from different initial temperatures, are summarized in Table 1. We also calculate the exponent γ Itt for G > (t, t − 100), G > (t − 300, t), G > (t, t − 300), G > (t − 500, t) and G > (t, t − 500). The numerical values do not change significantly compared to the values given in Table 1 for G > (t − 100, t) hence, we can conclude that G > (t − t a , t) and G > (t, t − t a ) thermalize exponentially with the same exponent for arbitrary t a .  Let us now look at the thermalization rate of the effective temperature in the case of step as well as bump quench 8 . Once the system thermalizes it acquires the final equilibrium temperature, which we denote as 1/β f . Following [5], we assume that the relaxation behaviour of the effective temperature is given by 14) The β ef f settles down to the β f in the long time limit. Therefore one needs to determine α and Γ in (3.14). To do that we make a change of variables from (t 1 , t 2 ) to T = (t 1 +t 2 ), and t = t 1 − t 2 and analyse the ratio of the Keldysh Green's function with the spectral function for different values of t but holding T fixed. We repeat this for different values of T . It is convenient to work in the frequency space, therefore we carry out a partial Fourier transform with respect to t and compute We then calculate the effective temperature by taking small ω limit of the above quantity by fitting it to a hyperbolic tangent function as given in (3.3). The effective temperature obtained in this procedure is then fitted to the exponential ansatz given in (3.14). We summarise our results in the figure (11a) and (11b).

Bump
Step 10  It is clear from Fig. 11b that like in the case of the step quench as was observed in [5], the thermalization rate Γ for the bump quench is also proportional to the final temperature. It is also interesting to note that the constant C is independent of the coupling J 4 .

Conclusion and Discussion
We studied quench in the SYK model with different quench protocols. While we have presented results for q = 2 theory, and q = 4 theory with step and bump quench protocols, we have carried out this analysis for q = 6 as well as for q = 8 models. We find that the qualitative features of the results are similar to the q = 4 cases.
We observed that the q = 2 theory does not thermalize for any of the quench scenario we considered. We considered quenching of J 4 , J 6 , and J 8 using step and bump protocol. The initial states that we considered are thermal states of inverse temperature β i = 10, 20, and 30 as well as the ground states. An interesting aspect of all the quenches is that the greater Green's function G > (t 1 , t 2 ) equilibrates instantaneously as shown in (3.7). Its expectation value freezes once both the time arguments are outside the quench region. Although in the final states G > (t 1 , t 2 ) equilibrates instantaneously, its equilibrium value is not the same as the thermal ensemble expectation value.
The instanteneous equilibation or freezing that we observed is like a glassy state. It can be shown that if the final theory have both J 2 and J 4 couplings, then two point functions always thermalize. This is true even for arbitrarily small J 4 coupling in the large N limit that we are considering. We expect that this would change if we consider effects subleading in N, where J 2 and J 4 couplings would truly start competing [65].
It would be interesting to identify the final state after each of these quenches. It is, however, beyond the scope of the present work since we are working only with the equations of motion of the G > (t 1 , t 2 ) and solving them as an initial value problem. The q = 2 theory is not chaotic and does not satisfy the ETH, nevertheless thermalization in this theory is possible if the final state were a KM state (2.32). This, for example, happens quite often with step quenches in 1+1 dimensional theories (even in integrable theories) where the analog of KM states are the CC states (2.30).
In q = 4 theory, we find that thermalization happens in all the quench scenario we considered. We considered quenching with J 2 , J 6 , and J 8 using step and bump protocols. The initial states are thermal states of inverse temperature β i = 10, 20, and 30. We examined two kinds of greater Green's functions, G > (t − t a , t) and G > (t, t b ) as a function of time t with fixed t a and t b .
When both the time arguments t − t a and t are outside the quench region, both the real and imaginary parts of G > (t − t a , t) are exponential functions with the same exponent. This exponent γ Itt is equal to the value of coupling J 4 of the system. It would be interesting to compare this result with the bulk calculation. The scalar Feynman propagator is known to thermalize instantaneously in the AdS 2 -Vaidya spacetime [66,67]. This instantaneous thermalization is seen only for AdS 2 and is not observed, say, in AdS 3 . There is no analogous calculation for fermions is available in the literature.
The long time limit of both the real and imaginary parts of G > (t, t b ) are exponential functions with the same exponent. This exponent γ It is equal to the exponent γ ret of the retarded Green's function G R (t 1 , t 2 ) in a thermal ensemble (3.12) with temperature equal to the final temperature of the quench process. This is obvious at least for the imaginary part of the G > (t, t b ) since the system thermalizes. G R (t 1 , t 2 ) is a simple multiple of the imaginary part of G > (t 1 , t 2 ). As one can see in Figure 12, the long time limit of G > (t, t b ) is calculated in a subset of the large (t 1 − t 2 ) of G R (t 1 , t 2 ).
We also studied the thermalization rate for both step and bump protocols. In both Figure 12: The large t limit of greater Green's function G > (t, t b ) is calculated in the large (t 1 − t 2 ) region of the retarded Green's function G R (t 1 , t 2 ). Moreover, in this region, the system has more or less thermalized. Hence, γ It = γ ret .
cases we find that the thermalization rate is proportional to the final temperature 1/β f and the proportionality constant C is independent of the coupling constant of the final theory.
One clear and important observation that we can make is that the thermalization in q = 4 theory is not because the final state is a KM state. If the final state had been a KM state, G > (t − t a , t) would have thermalized instantaneously once both its time arguments are outside the quench region. 9 This is because the 'diagonal' two-point function that we are considering are already thermalized in a KM state.

How to prepare KM states
The KM state, in principle, can be prepared by performing a sudden quantum quench starting from the ground state using the extra term where s k 's specifies the particular |B s defined in (2.31). This new term has been used in a different but related context in [30]. The argument behind this assertion is similar to the argument provided in [68] for the preparation of thermofield double state by performing a sudden quench. We will consider small µ limit. The full Hamiltonian before the quench at t = 0 is H + H µ = (i) q/2 1≤i 1 <i 2 <...<iq≤N J i 1 ,i 2 ,..,iq ψ i 1 ψ i 2 ....ψ iq + iµ N/2 k=1 s k ψ k ψ k+1 (4. 2) The ground state of the above Hamiltonian is the state which minimizes the second term. But minimizing the second term corresponds to strong positive or negative correlation of ψ k and ψ k+1 depending on the value of s k . Strong correlation of ψ k and ψ k+1 is the basis of the definition of |B s in (2.31). In hindsight, it is in some sense obvious why the KM states were not obtained from the step quench using the disordered couplings like j 2,ij or J 2 . This is because not just two fermions, but all the fermions were randomly and strongly correlated in the ground states of the initial Hamiltonians.

Ergodicity versus Mixing
In this work, we don't consider long time averaging. q = 4 theory satisfies eigenstate thermalization hypothesis(ETH). But thermalization from ETH crucially requires long time averaging. Thermalization without long time averaging has been observed in many other works but in most of the cases it is because the final state turns out to be a very particular state like CC states. So, in this sense, the thermalization that we observed is much more robust than what one expects from ETH. Thermalization with ETH follows from quantum ergodicity. But what we observe is more akin to a quantum version of mixing.
In classical theories, mixing is a much stronger phenomenon compare to ergodicity. Figure (13a) shows ergodic evolution in the classical phase space. The initial state is described by an ensemble concentrated in the deformed rectangle in the phase space. The volume is conserved under time evolution due to the Liouville theorem for a closed system, but the shape can change. For ergodic systems, the shape of the initial sample hardly changes but it sweeps out the entire allowed space under time evolution. So, a long time averaging gives the expectation value in the micro-canonical ensemble. In mixing, as shown in Figure (13b), the initial sample spreads out and reaches infinitesimally close to all the points in the allowed region of the phase space. So, without time averaging, mixing gives the expectation value in the microcanonical ensemble.
Using this classical analogy, we believe that even in quantum systems, chaos is a much stronger condition for thermalization than the eigenstate thermalization hypothesis(ETH). Our results on thermalization in the quenched SYK model seem to suggest that quench without long time average is a quantum analog of mixing. It would be interesting to make this more concrete. We hope to return to this soon.  Figure 13: (a) Ergodicity: the shape of the initial sample only changes slightly but sweeps out the entire allowed region under time evolution, (b) Mixing: the initial sample spreads out and reaches infinitesimally close to all the points in the allowed region of the phase space. Figure adopted from [50,69].