Thermalization in different phases of charged SYK model

We study thermalization of charged SYK model in two different phases. We show that both the highly chaotic liquid phase and the dilute gas phase thermalize. Surprisingly the dilute gas state thermalizes instantaneously. We argue that this phenomenon arises because the system in this phase consists of only long-lived quasi-particles at very low density. The liquid state thermalizes exponentially fast. We also show that the additional introduction of random mass deformation (q = 2 SYK term) slows down thermalization but the system thermalizes exponentially fast. This is observed despite the fact that the addition of large q = 2 SYK interaction forces spectral statistics to obey Poisson statistics. An interesting new observation is that the effective temperature is non-monotonic during thermalization in the liquid state. It has a bump at relatively long time before settling down to the final value. With non-zero chemical potential, the effective temperature oscillates noticeably before settling down to the final value.


Introduction and summary
Non-equilibrium dynamics of interacting quantum systems has been a subject of great interest for a long time [1,2]. It is a subject of interest in various fields of physics, e.g. various aspects of condensed matter physics, heavy ion collisions, AdS/CFT correspondence and black hole dynamics, quantum cosmology, etc. It also has wide applications in engineering sciences (e.g. working principle of everyday semiconductor devices), biophysics (e.g. application in protein folding), etc. It is also the underlying principle for the impending rise of quantum computers in the next few decades [3].
In this paper we will examine non-equilibrium dynamics of chaotic quantum systems. We will be considering closed quantum systems. It is generally expected that these systems thermalize. We will be considering thermalization in the sense that, starting from a thermal state, after some finite perturbation the system comes back to thermal equilibrium. 1 The most interesting aspect of this work is that we find certain special states which thermalize instantaneously.
The original SYK model consists of Majorana fermions but we will be considering complex fermions with which the model has a conserved charge other than the Hamiltonian.

JHEP04(2021)157
For the rest of the paper, we will refer to a single interaction term as, for example, q = 2 SYK term. This interaction term can be a part of a Hamiltonian consisting of multiple SYK interaction terms. We will refer to the full model as, for example, (q = 2, 4) SYK model. This model has both q = 2 and q = 4 SYK terms. Another example is (q = 2, 4, 6) SYK model which will consists of q = 2, q = 4, and q = 6 SYK interaction terms.
In [8], the different phases of charged SYK model have been studied in the presence of chemical potential. It has also been examined if (q = 2, 4) SYK model undergoes a similar phase transition (in the absence of chemical potential). In this work, we will study nonequilibrium dynamics in different phases of charged SYK model and (q = 2, 4) SYK model. We do this by taking the systems out of equilibrium by performing quantum quenches. We will solve the non-equilibrium Green's functions using their equations of motion. We will be using lesser Green's function G < (t 1 , t 2 ), greater Green's function G > (t 1 , t 2 ), retarded Green's function G R (t 1 , t 2 ) and advanced Green's function G A (t 1 , t 2 ). The definition of these different Green's functions are as follows Non-equilibrium dynamics of SYK models and other related models have been studied in [9][10][11][12][13]. Pure excited states of SYK models have also been studied in [14][15][16].
We will briefly clarify on what we meant by a chaotic quantum system. There are various diagnostics of quantum chaos. The two most popular and well studied diagnostics are comparison of spectral statistics with random matrix theory (BGS conjecture, [17]) and exponential decay of Out-of-Time-Ordered correlators (OTOC) [18]. For a system like SYK model with widely separate time scales of dissipation and scrambling, OTOC decays exponentially. Considering the operators to be the microscopic fermionic degrees of freedom, the OTOC is C(t 1 , t 2 ) = Tr e −βH SYK /4 Ψ † i (t 1 )e −βH SYK /4 Ψ † j (0)e −βH SYK /4 Ψ i (t 2 )e −βH SYK /4 Ψ j (0) where H SYK is the Hamiltonian. λ L is the Lyapunov exponent. We have taken the regularized OTOC as in [18]. We will be working with two-body q = 2, four-body q = 4, and six-body q = 6 all-to-all random interactions. So the general Hamiltonian is

JHEP04(2021)157
Since we will be performing quantum quenches using q = 2 and q = 6 SYK terms, we have considered these interactions to be time-dependent. The mass term introduces an effective chemical potential η = µ. We can also consider thermal states with explicitly chemical potential η turned on in which case the total effective chemical potential is η + µ.
It has been shown that in the presence of chemical potential there is a first order liquidgas phase transition [8,19,20]. The phase transition is between the highly chaotic state of the SYK model which we will call the liquid phase and the dilute gas phase. 2 Unlike the liquid phase, the dilute gas phase has a unique well-separated ground state so the zero temperature (T → 0) entropy is log 1 = 0. This phase is also highly compressible as the name suggests. Again, as the name suggests, the charge or occupation number is also very low in this phase. For example, Q = 0.00035 at β = 30 and µ = 0. 27 (1.9) The phase diagram is produced in [20]. The phase transition happens in a finite range of the chemical potential. There is also a range of temperature in which the system can be in either of the two phases. So, one of the phases is metastable in this temperature range. In the liquid phase, the system is strongly interacting and there is no quasi-particle dynamics. While in the dilute gas phase, the system consists of long-lived quasi-particles at very low density. This is evident from the plot of the spectral function in the dilute gas phase. The spectral function is defined as the imaginary part of the retarded Green's function For the sake of simplicity the phase transition is always studied without the two-fermion random interaction q = 2 term. The Lyapunov exponent of the liquid phase is suppressed exponentially when the chemical potential is turned on [8,21]. This implies that chaos is suppressed. The Lyapunov exponent in the dilute gas phase has also been calculated in [8].
It is very small at very low temperature. But it is large at relatively high temperature especially in the temperature range where the system can exist in either of the two phases. Figure 1 is a reprint of the plot of the normalized Lyapunov exponent for the two different phases with varying temperature and a fixed chemical potential. It is conventional wisdom that quantum chaos implies thermalization. So the behaviour of the Lyapunov exponent suggests that turning on the chemical potential will slow down the thermalization process in the liquid phase. We will indeed find this to be the case. With the same reasoning, we expect that the dilute gas state would thermalize very slowly (or not at all) at very low temperature. At higher temperature in the dilute gas phase, we still expect that the system would still thermalize very slowly because the system still consists of long-lived quasi-particles at very low density. Surprisingly, we will find that the system in the dilute gas phase thermalizes instantaneously. The final state is a thermal state with a different temperature and a different chemical potential as compared to the initial state. We also like to point out that in purely (q=2) SYK model, the system equilibrates instantaneously but the final state is not a thermal state [10].  Figure 1. Normalized Lyapunov exponent with varying inverse temperature in the two different phases for η = 0.27 (reprinted from [8]). The sharp changes in the blue (orange) curve represents the transition from the highly chaotic liquid (dilute gas) state to the dilute gas (liquid) state. So, the blue (orange) curve mostly constitutes the liquid (dilute gas) phase.
Highly chaotic systems thermalize exponentially fast [9]. So, we believe that the phenomenon of instantaneous thermalization of the dilute gas state has an underlying physics different from the chaotic dynamics. It is interesting that a simple harmonic oscillator (SHO) thermalizes instantaneously [8]. This is true for both fermionic and bosonic SHOs.

Consider the fermionic SHO Hamiltonian
Consider we start from a thermal state of this system with the density matrix e −βH fSHO . We can perform a quantum quench by making µ time dependent. For simplicity, let us consider that we take µ to zero suddenly. 3 So, µ(t) is a step function. The result is that as soon as µ goes to zero, the system is described by the thermal ensemble e −βηQ , where η = µ. So, in summary It is a special kind of instantaneous thermalization where temperature remains constant but the chemical potential changes. This still happens even when the initial Hamiltonian has other terms as long as Q is a conserved charge, meaning Q commutes with the other terms. In our case here, the other terms will be the SYK interaction terms. So far this does not explain the instantaneous thermalization in the dilute gas phase when we perform quantum quench using time dependent SYK term. The connection is that the dilute gas phase consists of only long-lived quasi-particle excitations at very low density. The quantum quench slightly changes the fuzzy dressing of the quasiparticles. It also leads to the change in the temperature during the quench process unlike the case of fermionic SHO above. This 3 But note that the following results apply to any functional form of µ as a function of time t. It can be slow quench, a zig-zag quench, . . . etc., or the final value of µ can be any real number.

JHEP04(2021)157
argument is supported by the fact that the spectral function changes only slightly during our quench process. On the other hand, the quench process in the liquid phase changes the spectral function significantly. Detailed comparison can be found in section 3.2.
The smallest time scale in our theory is 1/J 4 = 1. We used discrete time step-size dt smaller than this value. So, the instantaneous thermalization is not due to use of inadequate time scale in our numerics. Moreover, it has been found that the dilute gas phase of twosites coupled SYK model does not thermalize instantaneously [22]. We believe that this is due to the fact that the site hopping term does not commute with the SYK terms so the change in the dressing of the hopping term leads to a thermalization process which is not instantaneous. It is unlikely that dilute gas phase in higher dimensional theories would thermalize instantaneously. Mass quenches in higher dimensional free relativistic theories have non-trivial time dependence [23][24][25].
We will also study the non-equilibrium dynamics of (q = 2, 4) SYK model. The q = 2 SYK term alone is integrable. The j 2,ij couplings can be diagonalized resulting in a theory of free N fermions with random masses. There is no sharp phase transition if we consider the (q = 2, 4) SYK model in the absence of chemical potential [8]. The Lyapunov exponent is suppressed by the integrable interaction. We calculate the Lyapunov exponent of this system to high precision. But it does not sharply go to zero even when the (q = 2) interaction strength is very strong and at very low temperatures. We will also show that this system thermalizes exponentially fast. These observations imply that this system is always in the highly chaotic liquid phase without quasi-particle excitation. For a chaotic system the spectral statistics is described by Wigner-Dyson (WD) statistics. For a generic integrable system, the spectral statistics is described by Poisson distribution. It has been shown that (q = 2) SYK interaction forces the spectral statistics towards Poisson distribution [26]. 4 This exercise was done considering finite size systems and numerically diagonalizing the Hamiltonian. But as we have mentioned above, the system is always chaotic and always thermalizes exponentially fast. So, this is a shortcoming of the BGS conjecture. It is hard to quantify how much chaotic a system is or if the system is completely integrable. There has been other works in similar line where the spectral statistics is not WD statistics but the system is nevertheless highly chaotic [28].
Technically, we show the exponentially fast thermalization process in the liquid phase by calculating the effective temperature during the non-equilibrium dynamics. We find that, without chemical potential, the effective temperature has a single bump before settling down to the final value. With non-zero chemical potential, besides the bump the effective temperature oscillates noticeably before settling down to the final value. The oscillation frequency depends on the chemical potential and the frequency cutoff used to calculate the temperature. In case of the dilute gas phase, it is not possible to calculate the temperature. Even out of equilibrium, it is not possible to calculate the effective temperature. So, we have to employ a different technique to show the instantaneous thermalization. The details can be found in section 3.2.

JHEP04(2021)157
The conserved charge is It does not change during our quench processes using time-dependent SYK interaction terms. This is true in both the phases. This is because the time-dependent SYK terms commute with the charge operator. This implies that spectral asymmetry frequency (SAF), in the liquid phase, remains unchanged when we perform quantum quenches using timedependent SYK interaction terms. SAF is the shift in the peak of the product of the retarded Green's function and advanced Green's function [29] is an even function of ω. ω s is also the position of the peak of the spectral function A(ω) = −2 ImG R (ω). The relation between SAF and the conserved charge at low temperature is [29] It is interesting that the conserved charge Q depends only on SAF even at high temperature.
Other quantities like chemical potential and temperature changed when we compare the final state and the initial state after such quantum quenches. In summary, the technical results of this work are: 1. We show that both the highly chaotic liquid phase and the dilute gas phase of charged SYK model thermalize. A system in the dilute gas phase thermalizes instantaneously.
2. (q = 2, 4) SYK model is always in the highly chaotic liquid phase without any quasiparticle excitation and the system always thermalizes exponentially fast. We also calculate the Lyapunov exponent of this system with high precision.
3. In quantum quenches starting from the liquid state, the effective temperature is nonmonotonic.
Our initial aim of this work was to show that a system in the dilute gas state thermalizes very slowly (or fails to thermalize completely) but instead we find that the system thermalizes instantaneously. Failure of thermalization in highly interacting systems is a topic of intense research interest in experimental as well as theoretical physics. There are two popular paradigms in which a chaotic system fails to thermalize. The first one is the existence of quantum scars [30][31][32][33][34][35][36]. Quantum scars are eigenstates which violates eigenstate thermalization hypothesis (ETH). They also have finite energy density but anomalously low entanglement [37]. So states (pure or mixed) formed out of quantum scars do not thermalize. The other paradigm is many body localization (MBL) [38][39][40][41]. MBL is the suppression of chaos and slowing down of thermalization of an otherwise chaotic system due to the introduction of random disorder.
It would be interesting to examine closely the nature of the dilute gas phase. We believe that an analytical treatment might be possible especially for this phase. The dynamics of JHEP04(2021)157 weakly interacting classical systems are well understood. The most famous example being the Fermi-Pasta-Ulam problem where the system fails to thermalize for exponentially long times even when there is a small but finite non-linear term. Similar phenomenon in weakly interacting quantum systems called prethermalization has been a topic of great interest in recent times [42][43][44].
Instantaneous thermalization has previously been shown in quantum quenches starting from the ground state of a gapped theory to a gapless theory in two spacetime dimensions [23,24]. This applied only to correlators consisting of holomorphic operators of the 2D conformal field theory. Other correlation functions, say, of scalar operators, do not thermalize instantaneously in general. In holographic set-up, it has also been shown that one-point functions thermalize instantaneously [45] but two-point functions does not [46]. It is worth mentioning here that [47] claims that Hamadard functions thermalize instantaneously. But again, this has been explained in [46] that this particular two-point function does not probe inside the apparent horizon during black hole formation, which is why they thermalize instantaneously. Other two-point functions on the boundary theory do not thermalize instantaneously. Recently, it has also been shown that (q → ∞) SYK model thermalizes instantaneously [9].
We would like to point out the differences of our work from [11] which considers a modified system consisting of SYK model coupled to quadratic peripheral fermions. This modified system has the familiar highly chaotic non-Fermi liquid (NFL) phase and a Fermi liquid (FL) phase. Interestingly the above paper considers the modified model without any mass or chemical potential. In this setting, the system is either in the NFL phase or in the FL phase depending on a parameter p which is the ratio of the number of the quadratic peripheral fermions and the number of SYK fermions. So depending on the value of p, the Hamiltonian of the modified system dictates if the system is in the NFL phase or in the FL phase. This is different from our present case where for the same Hamiltonian (at a given temperature and a given chemical potential) the system can be either in the highly chaotic liquid phase or the dilute gas phase. Another important difference is that in the above paper the modified system in the FL phase thermalizes slowly but in our present work we observe that a system in the dilute gas phase thermalizes instantaneously.
The outline of this paper is as follows: we will work out the details of SYK model with complex fermions in section 2. The Kadanoff-Baym (KB) equations are derived in 2.1. We also explain the numerical methods involved in solving the equations. In section 3, we explain the phase transition and the details of the two phases. We examine in details the non-equilibrium dynamics of the liquid state in section 3.1. We show that a system in the dilute gas phase thermalizes instantaneously in section 3.2. In section 4, we examine the non-equilibrium dynamics of (q = 2, 4) SYK model. Section 5 consists of conclusions from this work.

SYK model with complex fermions
The general Hamiltonian is given in (1.8). To make the derivation simpler we will first consider only (q = 4) SYK model with explicit time-dependence for the moment. The

JHEP04(2021)157
generalization to the full Hamiltonian of (q = 2, 4, 6) SYK model is straightforward. The action on the Schwinger-Keldysh contour C is We work with quenched averaging of the coupling in the large N limit [9,10]. The contour ordered Green's function is defined as The partition function becomes Integrating out the fermions, the action in terms of the bilocal fields is The equations of motion of G and Σ are These are the Schwinger-Dyson (SD) equations in the Schwinger-Keldysh contour. Generalizing these equations for the general Hamiltonian (1.8) of (q = 2, 4, 6) SYK model, the SD equations are

JHEP04(2021)157
We obtain different Green's functions when we go from the Schwinger-Keldysh contour to the real time axis where t ± = t±i . Operator at t + comes before operator at t − . With this, the SD equations are where G R is defined in (1.5) and Σ R is also similarly defined. The equilibrium solution are solved numerically. The connection between the above two equations are the fluctuationdissipation relations which gives the expression of G > (ω) and G < (ω) in terms of the spectral function A(ω).
In the absence of µ and η, the SD equations are same as that of Majorana SYK models [9,10]. This is because in thermal equilibrium, This relation holds true even out-of-equilibrium during time evolution starting from thermal equilibrium. So in the absence of mass or chemical potential, quantum quenches using timedependent J 2 , J 4 or J 6 in SYK model with complex fermions are same as quenches with the same time-dependent couplings in SYK model with Majorana fermions. Also note that in thermal equilibrium, Again this relation also holds true out-of-equilibrium starting from thermal equilibrium. This relation halves the computer time for time evolution because we can solve G >(<) (t 1 , t 2 ) for either only t 1 ≥ t 2 or t 1 ≤ t 2 .
We also verify conservation of energy during the quench processes. The expression for energy is

JHEP04(2021)157
The temperature and chemical potential can be calculated from G > (t 1 , t 2 ) and G < (t 1 , t 2 ) using the relation The effective temperature during the non-equilibrium time evolution is calculated by using the method in [9,10]. For this, we perform a coordinate transformation from (t 1 , t 2 ) to Fourier transform w.r.t. t − and using (2.18) at small ω region gives the effective temperature as a function of t + . Note that t + increases or decreases in time step of 2×dt. Highly chaotic system without quasiparticles are expected to thermalize exponentially as a function of the final temperature where the thermalization rate is directly proportional to the final temperature. This indeed has been shown to be true for SYK model without chemical potential. The effective temperature is given by where T f is the final temperature.

Kadanoff-Baym equations and quantum quenches
We perform the quenches by solving the equations of motion of the Green's functions in integro-differential form which are well-known as Kadanoff-Baym (KB) equations. The details of the derivation of the KB equations from the SD equations can be found in [10].
After performing a convolution in eq. (2.8) with G we get The real time KB equations are obtained after contour deformation (Langreth rule). The imaginary leg in the contour is removed using Bogoliubov principle with weakening initial correlations.
Pure (q = 2) SYK model does not thermalize [10]. We can see from the above equations that (q = 2) SYK model even in the presence of chemical potential or a mass term does not thermalize at all. The system freezes conpletely as soon as the time-dependent perturbation stopped. This is because the expressions on the right side of (2.22) and (2.23) (or (2.24) and (2.25)) are same. So,

JHEP04(2021)157
Moreover, the final configuration is not thermal. So, (q = 2) SYK model equilibrates instantaneously but does not thermalize. Unlike the (q = 2) SYK model, quantum quenches in the presence of (q = 4) interaction are non-trivial. Using relation (2.16), we will only use (2.23) to evolve G > (t 1 , t 2 ) and (2.25) for G < (t 1 , t 2 ). The most convenient forms of the equations for numerical applications are Note that the second integrals are same in the two equations. These equations can be solved incrementally/causally using a Predictor-Corrector scheme. We use forward difference for the prediction and we use backward difference for the correction. In the presence of nonzero µ, the forward difference and the backward difference must be strictly taken otherwise the numerical scheme fails to converge. The predicted value is consists of the sum of the two integrals in (2.27). The correction is performed using the backward difference formula and simple mixing to achieve convergence.
An important check for our numerical codes is to perform equilibrium time evolution without quantum quenches. It serves three purposes. First it verifies that our codes are correct. It also verifies that our initial data are correct and accurate. The initial data are generated by solving the SD equations. Lastly, it verifies that numerical errors are under control.

Charged SYK model
In this section we will consider the system which has non-zero chemical potential or mass term in the Hamiltonian. We will consider only (q = 4) interaction. The Hamiltonian is 5 As we have mentioned in section 1, the system can undergo a liquid-gas phase transition in the presence of effective chemical potential (explicit chemical potential η or mass µ or both). In the highly chaotic liquid phase, the system does not have a quasiparticle picture.
In the dilute gas phase, the system consists of long-lived quasi-particles at very low density. Figure 2 are plots of spectral functions in the two different phases. The spectral function in the dilute phase has a sharp single peak. The two phases are separated by a large energy gap. Figure 3 are plots of energy in the different phases with the mass term or the chemical potential.
The dilute gas phase also has occupation number close to 0 or 1 depending on the sign of the effective chemical potential. It can be easily seen from the plot of the spectral function. The area under the plot is equal to 1, this is fixed by fermionic commutation relation. The occupation number is given by In case of the liquid phase, the spectral function is spreaded and the fermionic distribution function more effectively suppresses the occupation number. 5 Note that later we will perform quantum quenches in this system. Keeping the q = 4 SYK term unchanged, we will use time-dependent q = 2 or q = 6 SYK term to take the (q = 4) SYK system out of equilibrium. So in that context, the full time-dependent Hamiltonian will consist of q = 4 SYK term and a time-dependent q = 2 or q = 4 SYK term.  Figure 1 is the plot of normalized Lyapunov exponent for the two different phases with varying temperature and a fixed chemical potential. Lyapunov exponent in the liquid phase is large. The normalized Lyapunov exponent increases as we decrease temperature for a fixed effective chemical potential. At the transition point from the liquid phase to the dilute gas phase, the Lyapunov exponent decreases sharply. But note that Lyapunov exponent in the dilute gas phase is non-zero and it is large at higher temperature.

Thermalization in the liquid phase
In this subsection and the next subsection, we will consider quantum quenches in (q = 4) SYK model in the presence of chemical potential. We will study time evolution of the system after the system is taken out of equilibrium using time-dependent q = 2 or q = 6 SYK interaction terms. We will consider only bump quenches where the time-dependent term is turned on for a short time duration and completely turned off again. The timedependent Hamiltonians are or Starting from the liquid state, we find that the system evolves non-trivially even after the time-dependent perturbation has been turned off. The system thermalizes exponentially fast.   . Plots of the effective temperature as a function of t + starting from different states and for different µ. The quenches have been tuned so that the final inverse temperature is β ∼ 18. All the quenches are performed with time dependent q = 4 SYK interaction except for the red curve which was with time dependent q = 6 SYK interaction. The blue curve is for µ = η = 0. The inset shows the non-monotonicity of the effective temperature for this case. The orange curve is for initial state with η = 0.24. The green curve is for initial state η = 0 but with µ = −0.263. The red curve is for the same initial state as the green curve but the time dependent perturbation was (q = 6) interaction.
The initial equilibrium state is prepared by solving the SD equations (2.11) and (2.12) without the q = 2 and q = 6 SYK terms. To perform the quantum quench, the KB equations (2.27) and (2.28) are solved numerically. The time dependent terms are turned on from t 1 , t 2 = 1 × dt to t 1 , t 2 = s × dt where dt is the discrete time step-size. Figure 4 is the plot of effective temperatures for quenches in different settings but with the final temperature T f ∼ 0.55. As we can see the effective temperature is non-monotonic in all the cases. Without chemical potential, the effective temperature settles down to the final value. But with chemical potential, the effective temperature oscillates before settling down to the final value. As we have mentioned earlier, we use (2.18) at small ω region to calculate the temperature. The oscillation of the effective temperature actually depends on the frequency cutoff that we used for calculating the temperature. Figure 5a are the plots of effective temperature calculated using different frequency cutoffs.
Note that the value of the chemical potential changes during quenches. As mentioned in section 1, The charge operator commutes with the time-dependent SYK interaction terms. So, the conserved charge does not change during our quench processes. Accordingly, the spectral asymmetry frequency (SAF) also does not change during the quench processes. As defined in (1.14), SAF is the position of the maximum of G R (ω)G A (ω). Figure 5b are the plots showing the position of SAF before and after a quantum quench.
In the liquid state, the Green's functions also thermalize exponentially fast. The Green's functions change non-trivially even after both the time arguments have crossed the quench region. The Green's functions converge towards their final values exponentially [10]. Figure 6 is the plot of the real part of the greater Green's function G > (t−t a , t) as a function of t. We will find that in case of quenches starting from the dilute gas state, the Green's functions will reach their final values abruptly.

Thermalization in the dilute gas phase
We will now consider quantum quenches where the initial state is a dilute gas state. With µ = 0, the Green's functions oscillate since the effective theory is a weakly interacting massive theory. So, we have to consider very small time-step size dt which results in very large number of discretization points. In case of µ = 0 and η = 0, one can take large dt but since the effective theory is a weakly interacting (almost) massless theory, the Green's function do not decay fast so one has to consider again a very large number of discretization points although dt can be large. Another technical difficulty while dealing with the dilute state is that the temperature cannot be calculated using (2.18) even for an equilibrium state solution calculated directly from the SD equations. The numerical precision is not good enough to cancel the spectral function and reproduce the tanh function. Consider the case of µ = 0 and η = 0.27, as JHEP04 (2021)157 we can see in figure 2 the spectral function is peaked around ω = 0. From (2.18), the expression on the left hand side is zero and rapidly varies only around ω = −η. But around this region of frequency ω, G K (ω), G R (ω) and A(ω) are numerically very small. So during the quench process, we will compare the energy of the final state with the energy of dilute gas states generated using the SD equations with different temperature and chemical potential. After finding a match in the energy, we compare G >(<) (t 1 , t 2 )'s of the final state from the quantum quench and the thermal state.
As we have mentioned at the end of section 2.1, equilibrium time evolution without quantum quench is an important check. This check is very important for the dilute gas states due to the technical difficulties mentioned above.
We will present three cases: 3. We also considered a case in which the system is perturbed with time-dependent J 6 term. For this, we used the same initial state as in the second case, i.e., µ = 0, η = 0.27 and β = 50. For the quantum quench, we perturb the system with J 6 = 100000 from time t = 1×dt to t = 9×dt. The system thermalizes to the dilute gas state with η = 0.428934 and β = 31.471. The mass µ = 0 is not changed during the quench.
We find that the systems stop evolving instantaneously when the time-dependent perturbation is turned off. The Green's functions freeze as soon as the two time arguments t 1 and t 2 cross the quench region. Figure 7 are plots of the real parts of the greater Green's function G > (t − t a , t) for different cases that we are considering. The imaginary part of G > (t − t a , t) as well as the lesser Green's functions also freeze as soon as the two time arguments cross the quench region.
We could not calculate the effective temperature and the effective chemical potential of the final state. 6 So, we search for thermal states with energy equal to that of the final state from the quantum quench. This is done by solving the SD equations for different values of temperature and chemical potential. In this two dimensional parameter space, JHEP04(2021)157

Re[G > (t-160,t)]
(b) we could find a line for which the energy matches the energy of the final state. Then we compare the Green's functions of these thermal states with the Green's functions of the final state. Again we could find a unique thermal state for which the Green's functions match those of the final state. So, in conclusion, the system in the dilute gas state thermalizes instantaneously. The instantaneous thermalization happens irrespective of the perturbing term. We used time-dependent J 2 or J 6 . Figure 8 are the plots of the real parts of the greater Green's functions G > (t−t a , t) obtained after quantum quenches and thermal states generated using SD equations. It is evident that the final states are thermal states.
As we have pointed out in section 1, the spectral function changes only slightly during the quench process in the dilute gas phase. Figure 9 are plots of the spectral functions of the initial state and the final state in the two different phases. As, we can see the spectral function changes significantly during the quench process in the chaotic phase. Now we will explain why the temperature of the dilute gas phase also changes during the quench process unlike in the case of fermionic SHO in (1.12). Since, the spectral function has a sharp peak near µ in the dilute gas phase, the fermionic distribution function JHEP04(2021)157

Re[G > (t 1 ,t 2 )]
(b)  (e) JHEP04(2021)157 is close to 1 at the energy scale where the non-trivial physics in happening. For example, This is the same reason why we could not calculate the temperature in the dilute gas phase using (2.18). Hence, the temperature in this phase is almost fixed by the spectral function itself. During the quench process when the spectral function changes (even though slightly) it leads to a change in the temperature of the system.

(q = 2, 4) SYK model
In this section, we will consider (q = 2, 4) SYK model without chemical potential. The Hamiltonian is 7 The (q = 2, 4) SYK model is always in the liquid state. The system does not undergo liquid-gas transition. But the (q = 2) interaction suppresses chaotic nature of the q = 4 SYK term. This can be seen from figure 10. The Lyapunov exponent is greatly suppressed due to the presence of non-zero J 2 coupling but it does not sharply drop to a negligible value which is expected for a phase transition as in figure 1. The suppression of chaos has also been shown from spectral correlation calculation in [26]. The introduction of large J 2 coupling forces the spectral statistics towards Poisson statistics. The spectral statistics obeys Poisson statistics for a generic integrable system while it obeys Wigner-Dyson(WD) statistics for a chaotic system. But this analysis are performed in finite systems (fixed N of the order of 10) so it is rather hard to precisely identify if the system is fully integrable or chaotic. In case of (q = 2, 4) SYK model, the system is always in the highly chaotic liquid phase. In this work, we study the non-equilibrium dynamics of this system and show that JHEP04(2021)157 the system thermalizes exponentially fast which is expected for a chaotic system. We also calculated the Lyapunov exponent of this system with high precision.
To calculate the Lyapunov exponent we followed the steps in [26]. Note that in the absence of chemical potential, SYK model with N complex fermions is same as SYK model with 2N Majorana fermions in every respect at large N limit. So, we performed the Lyapunov exponent calculations with Majorana fermions. F 1 (t 1 , t 2 ) from (1.7) satisfies F 1 (t 1 , t 2 ) = dt 3 dt 4 K R (t 1 , t 2 , t 3 , t 4 )F 1 (t 3 , t 4 ) (4.2) where G R (t) is the retarded Green's function and G lr (t) is G < (t + iβ/2). Using the ansatz F 1 (t 1 , t 2 ) = e λ L (t 1 +t 2 )/2 f (t 1 − t 2 ) and going to the frequency domain using Fourier transforms, we get Tuning λ L numerically to satisfy (4.4) gives the Lyapunov exponent. It is well-known that the Lyapunov exponent is bounded by 2π/β. So, in the figures we have plotted normalized Lyapunov exponent instead which is defined as λ * L = λ L /(2π/β) (4.5) Figure 10 is the plot of the normalized Lyapunov exponent as a function of the inverse temperature for different values of J 2 . It is interesting that with increasing inverse temperature the normalized Lyapunov exponent decreases gradually after a peak while in the presence of chemical potential the normalized Lyapunov exponent increases gradually in the liquid phase and sharply drops to a negligible value for the dilute gas phase as shown in figure 1. To perform the quantum quenches, we consider two sets of the parameters which were considered in [8]. We consider J 2 = 0.5, J 4 = 1 and show thermalization below β = 55. Both the initial inverse temperature and the final inverse temperature are below β = 55. We also consider J 2 = 2, J 4 = 1 and show thermalization below β = 15. For these sets of parameters, it has been shown that the spectral statistics is almost completely Poisson statistics [26].
Here also we will perform bump quenches. We work with dt = 0.02. With J 2 = 0.5, J 4 = 1, the initial state is at inverse temperature β i = 70. We took the time range t 1 −t 2 ∈ {−5000×dt, 5000×dt}. The quantum quench is performed by turning on J 6 = 0.4 for the time duration 9×dt from time t = 1×dt to t = 9×dt. The final inverse temperature is β f = 60.7. Figure 11a is the plot of the effective temperature as a function of t + . With J 2 = 2, J 4 = 1, the initial state is at inverse temperature β i = 30. The time range was {−3000 × dt, 3000 × dt}. The quantum quench is performed by turning on J 6 = 0.7 for the time duration 9 × dt from time t = 1 × dt to t = 9 × dt. The final inverse temperature is β f = 17.1. Figure 11b is the plot of the effective temperature as a function of t + . The system thermalizes exponentially fast.

Conclusions
We show that the liquid state in SYK model with complex fermions thermalizes. The presence of chemical potential suppresses the Lyapunov exponent. The effective temperature equilibrates exponentially fast. Closer examination reveals that the effective temperature is non-monotonic. Without chemical potential, there is a single bump and the effective temperature settles down to its final value. In the presence of the chemical potential, there are damped oscillations during thermalization. The frequency of the oscillations depends on the frequency cut-off used to calculate the effective temperatures.

JHEP04(2021)157
We also show that the dilute gas phase in SYK model with complex fermions thermalizes instantaneously. The process of instantaneous thermalization suggests that the underlying physics must be different from the well-known chaotic dynamics. We argued that this is because of the long-lived quasi-particle nature of the excitations at very low density in this phase. Moreover, in case of pure (q = 2) SYK model, we have also shown that the system equilibrates instantaneously although the final state is not a thermal state.
On the other hand, the (q = 2, 4) SYK model always thermalizes expoenentially fast. This happens even when the normalized Lyapunov exponent is extremely small λ * L ∼ 0.003 for J 2 = 2, J 4 = 1 at inverse temperature β = 17.1. With these interaction strengths, the spectral statistics of the model is almost completely Poisson statistics. So in conclusion, despite the nature of the spectral statistics, the system is always in a highly chaotic liquid phase without quasi-particle excitation.