Keldysh Wormholes and Anomalous Relaxation in the Dissipative Sachdev-Ye-Kitaev Model

We study the out-of-equilibrium dynamics of a Sachdev-Ye-Kitaev (SYK) model, $N$ fermions with a $q$-body interaction of infinite range, coupled to a Markovian environment. Close to the infinite-temperature steady state, the real-time Lindbladian dynamics of this system is identical to the near-zero-temperature dynamics in Euclidean time of a two-site non-Hermitian SYK with intersite coupling whose gravity dual has been recently related to wormhole configurations. We show that the saddle-point equations in the real-time formulation are identical to those in Euclidean time. Indeed, an explicit calculation of Green's functions at low temperature, numerical for $q = 4$ and analytical for $q = 2$ and large $q$, illustrates this equivalence. Only for very strong coupling does the decay rate approach the linear dependence on the coupling characteristic of a dissipation-driven approach to the steady state. For $q>2$, we identify a potential gravity dual of the real-time dissipative SYK model: a double-trumpet configuration in a near-de Sitter space in two dimensions with matter. This configuration, which we term a Keldysh wormhole, is responsible for a finite decay rate even in the absence of coupling to the environment.

The existence of an environment, either thermal, quantum, or the very measurement process, makes the time evolution of a quantum system not strictly unitary, complicating the description of its time evolution. An elegant approach to this problem, advocated by Caldeira and Leggett [1], is to describe the system plus environment by an Hermitian Hamiltonian. The integration of the environmental degrees of freedom then results in nonunitary dynamics. Especially for many-body systems, it is technically hard to tackle the resulting problem, either numerically or analytically, because the tracing out of the environment leads to nonlocal interactions in time.
In the limit where the environment has a sufficiently short correlation time, the Liouvillian generator that governs the quantum time evolution of the density matrix of the reduced system can be written in the so-called Lindblad form [2][3][4][5][6]. Compared to the usual von Neumann equation of motion, the Lindblad equation contains additional terms modeled by the so-called jump or Lindblad operators [3] that describe the coupling of the system to the environment. The advantage of this approach is that the quantum master equation is expressed only in terms of the system degrees of freedom, which facilitates its solution. In recent times, there has been a renewed interest in different aspects of the Liouvillian dynamics, especially in random systems [7][8][9][10][11][12][13][14][15][16], including the study of the spectrum [7,10,11,14], with the aim to extending the Bohigas-Giannoni-Schmit conjecture [17] to dissipative quantum systems [18][19][20], the characterization of the late-time dynamics [8][9][10]15] towards a steady state [10,16,21,22] or the robustness of dissipative quantum chaotic features [23].
The Keldysh path integral [24,25] is the standard procedure used to investigate the nonequilibrium time evolution of the density matrix. A key difference to the standard path integral, which will be very important in the following, is that, in order to set up the calculation for the density matrix, it is necessary [24,25] to consider the time evolution of a state matrix, not a state vector. This results in two time branches, one forward and one backward in time, and therefore to the doubling of degrees of freedom of the path integral (vectorization). In this representation, the density matrix is a state in the doubled Hilbert space and the Liouvillian has two parts: the first is anti-Hermitian, corresponding to the Hermitian Hamiltonian, before any coupling to the environment but multiplied by the imaginary unit, and its conjugate such that the whole system is anti-Hermitian; the second contains Lindblad jump operators depending on fields living in both copies and includes (i) an explicit coupling between the degrees of freedom of the two copies and (ii) intrasite interactions effectively modifying the (anti-)Hermitian part. Using standard field theory techniques, it is then possible to investigate the out-of-equilibrium time evolution of a strongly interacting system either coupled to a bath or undergoing a continuous measurement process. If we are interested in the steady state, and in sufficiently long timescales, it is plausible to assume that the initial state becomes irrelevant, leading to further simplifications, namely, a time-translational invariant relaxation.
Intriguingly, operators similar to the vectorized Liouvillian are being employed in a completely different context: as non-Hermitian two-site Sachdev-Ye-Kitaev (SYK) Hamiltonians [26][27][28][29][30][31][32][33][34], N Majoranas with infinite range q-body interactions in zero dimensions, whose gravity dual is a wormhole [35][36][37]. Earlier, a single-site SYK model was proposed as an example of the holographic duality based on the fact that features such as the saturation of a bound on chaos [38], spectral correlations described by random matrix theory [39,40], and an exponential growth of the density of low-energy excitations [41][42][43] are shared by certain near-AdS 2 backgrounds [44][45][46][47] in the low-temperature limit. It has also been noticed [35] that the gravity dual of two identical two-site Hermitian SYKs with a weak intersite coupling is a traversable wormhole [36]. These traversable wormholes are the dominant configurations only for sufficiently low temperature. At higher temperature, a transition to a black hole configuration occurs, which can also be characterized by an analysis of level statistics [48].
By contrast, a pair of non-Hermitian PT-symmetric SYKs with no explicit intersite coupling in the low-temperature limit is dual to a Euclidean wormhole [37]. The presence of an explicit coupling triggers [49] a transition from Euclidean to traversable wormholes. Both traversable and Euclidean wormholes are characterized by a gapped ground state and a first-order phase transition in the free energy separating the wormhole and black hole phases. Differences include a qualitatively different dependence of the gap on the parameters of the model and a different pattern of oscillations of Green's functions in real time directly related to different symmetries of the solutions: U (1) for Euclidean and SL(2, R) for traversable. For earlier work on wormholes, mostly in higher dimensions, see [50][51][52][53][54].
In this paper, we show that the near-zero-temperature dynamics of a two-site SYK Hamiltonian in Euclidean time, dual to gravitational wormholes in certain limits, and the real-time Liouvillian describing a single-site SYK coupled to a Markovian bath close to the infinite-temperature steady state coincide. More specifically, we show that the equations of motion in the Euclidean problem in the low-temperature limit are identical to the equations of motion in the real-time Liouvillian evolution close to the steady state. This implies that the retarded Green's function of the real-time problem is identical to the equivalent one in the Euclidean problem. As a consequence, the so-called gap in the Euclidean formulation that characterizes the wormhole phase is equal to the dominant decay rate in the real-time Liouvillian evolution. This equivalence only holds when neglecting the boundary conditions, which are generally different for the two approaches. In practical terms, the equivalence works for sufficiently low temperatures in the Euclidean problem and for times such that the system is sufficiently close to the steady state in the real-time problem.
Finally, we note that the gravitational dual of a dissipative field theory with decoherence was investigated in Ref. [55]. We also stress that there are already several works that have studied different aspects of the non-Hermitian SYK model: from the use of the Lindblad formalism [56,57] and the investigation of the entanglement entropy growth [58,59] and decoherence effects [60][61][62], to its relation to Euclidean wormholes [37,63], replica symmetry breaking [49,64], and a symmetry classification [65]. However, these studies have not investigated the mentioned intriguing relation between wormholes in Euclidean time and strongly interacting dissipative systems in real time. We shall see that this finding has a strong impact on important features of the dynamics, such as an enhancement of the decay rate in the limit of very weak coupling to the environment.
The remainder of the paper is organized as follows. In the next section, we introduce the open SYK model, which will be our main object of study, and the equations of motion in the ΣG formulation [34] in both real and Euclidean time. In Sec. III, we calculate the long-time behavior of the Green's function by numerically solving the Schwinger-Dyson (SD) equations for q = 4 and find an anomalously large decay rate (or gap) for small intersite coupling. For q = 2, the Green's functions can be evaluated analytically and, because the model is not quantum chaotic, no anomalous relaxation is found, which is discussed in Sec. IV. The large-q limit is worked out analytically in Sec. V. The gravity dual of the SYK Hamiltonian coupled to an environment, dubbed Keldysh wormhole, is discussed in Sec. VI and concluding remarks are made in Sec. VII. Additional technical details are worked out in six appendices.

II. LINDBLAD EQUATION AND SYK HAMILTONIANS
We study a Hermitian SYK Hamiltonian with a q-body interaction of infinite range coupled to a Markovian environment. The SYK model is defined by the Hamiltonian where J i 1 ···iq are random numbers extracted from a Gaussian distribution of zero average and variance J 2 i i ···iq = (q − 1)!J 2 /N q−1 , with i 1 , . . . , i q = 1, . . . , N . We will mostly focus on the cases q = 2 (integrable), q = 4, and large q → ∞. The last two are quantum chaotic. The ψ i are Majorana operators defined by the commutation relation {ψ i , ψ j } = δ ij .

A. The vectorized Lindblad equation
The real-time evolution of the density matrix ρ [5] of this open SYK is simply dρ/dt = L(ρ) where L is the Liouvillian which can be conveniently expressed in the Lindblad form: Note the −i difference in dρ/dt = L(ρ) with respect to the time evolution of a state in the Schrödinger picture.
The Markovian environment is described by Lindblad jump operators of the form with i = 1, . . . , N and µ a positive real number. We believe that most of our main conclusions apply to more general jump operators but we will see that the relation to wormholes is better understood in this case.
For sufficiently long times, and taking into account that our jump operators are Hermitian, the density matrix will relax to a maximally entangled state at infinite temperature, characterized by L(ρ ∞ ) = 0. We note that for jump operators more general than those of Eq. (3), the system's density matrix may decay either to a maximally entangled state at finite temperature or to a nonequilibrium steady state, depending on the details of the Hamiltonian. Our main interest is in the approach to the steady state (4). For that purpose, we study the retarded Green's function where Θ(t) is the Heaviside function.
In order to compute this quantity, it is useful to employ the Keldysh path integral that involves the doubling of the degrees of freedom and the vectorization of the Liouvillian (see Appendix A for details). Here we state the final result for the partition function, where L, R stand for the two copies of the original degrees of freedom, i.e., the Hilbert space has been doubled H = H L ⊗ H R . The appropriate action is given by with the vectorized Liouvillian Ref. [56].

B. The Keldysh approach
The Keldysh approach relies on the closed-contour representation of the time integration, see also Refs. [22,24,25,56]. We introduce fields ψ i (t ± ) living on the closed-time contour C = C + ∪ C − , with real time running from −∞ to +∞ (branch C + ) and then back again to −∞ (branch ) corresponding to Majoranas acting on the left (right) of the density matrix. After disorder averaging, the action (7) can be written in terms of the Green's functions with a, b = +, − corresponding to the respective contours, and its self-energy Σ ab (t 1 , t 2 ) (see Appendix B for details). The saddle-point equations in terms of G ++ (t 1 , t 2 ) (both times on C + ), G +− (t 1 , t 2 ) (time t 1 ∈ C + and time t 2 ∈ C − ), and the corresponding self-energies are given by Here, we have assumed translational invariance of the solutions, so that they only depend on the relative time t = t 1 − t 2 , which makes it possible to write the first line of equations in frequency space. The equations can be decoupled in terms of the combinations This results in with self-energies given by At the saddle point, we have that the time-ordered and lesser Green's functions are, respectively This results in and we recognize G R/A as the retarded and advanced Green's functions. Moreover, at the saddle point, there is a single independent Green's function. From now on, it is understood that we are always at the saddle point and drop the brackets · · · throughout.
Numerically, we solve the saddle-point equations in terms of the spectral function see Eqs. (B30) and (B31) in Appendix B. In general, we expect the system to relax exponentially to its steady state at a rate Γ (also known as the gap). In that case, we have In general, the rate Γ can be obtained from the numerical solution of the saddle-point equations (see Sec. III for q = 4), but also analytically in special cases (see Sec. IV for q = 2 and Sec. V for the large-q limit).

C. The Euclidean approach
The Liouvillian of Eq. (8) can also be interpreted as a two-site non-Hermitian but PT-symmetric SYK Hamiltonian with intersite coupling [35,37,63,64]. Ignoring the irrelevant constant term, in this section we study the Hamiltonian with the sum running over the complete basis of eigenstates of H SYK L,R , K the complex conjugation operator, and the unitary operator U determined by the conditions (see Appendix A for an explicit construction of U ): Because we thus have that and |0 is also the ground state of the coupling term. Naturally, the β = 0 TFD is the vectorization of the (unnormalized) steady-state density matrix (5), which, as mentioned before, is the maximally The only difference with the real-time action is that the integral is over [0, β] instead of the real axis. The so-called ΣG action is obtained in two steps [34]. First, the integration over the Gaussian disorder results in a bilocal fermionic action. The original fermionic fields are then integrated by introducing two auxiliary fields: the Green's function, and the self-energy Σ ab (τ 1 , τ 2 ), a Lagrange multiplier that fixes the above definition of G ab in the path integral, where a, b = L, R. The resulting Euclidean action is given by with a, b ∈ {L, R}, s LL = s RR = 1, s LR = s RL = (−1) q/2 , t LL = t RR = 1, t LR = t RL = −1.
In the large-N limit, assuming time translational invariance, a saddle-point analysis results in the following SD equations: where the first two equations are in Fourier space and τ = τ 1 − τ 2 . The details of the derivation of the Euclidean SD equations are given in Appendix C. If the system is gapped, naively, we expect that G LR (τ ) decays exponentially at a rate given by the spectral gap ∆, namely, the gap between the ground state energy and the first excited state. However, we will see that, in the weak-coupling limit µ 1, the decay rate E g , also referred to as the gap, that governs the exponential decay of G LR is much larger than the spectral gap E g ∆.
In the time domain, at τ = 0, the saddle-point equations in the low-temperature limit admit solutions satisfying G LL = αG LR , which upon substitution gives α 2 = −1. In order to reproduce the delta function we thus find the solution where we have used that, in the limit of zero temperature, Using a suggestive notation, the saddle-point equations can be simplified by introducing the This results in the saddle-point equations (14) we found for the real-time approach: and For q = 4, these equations simplify to [66] Using the relation (32) we find that and we recognize G R and G A as the retarded and advanced Green's functions, introduced in the Keldysh approach of Sec. II B.
Importantly, we note that in the limit of zero temperature, and assuming no memory of the initial state, both actions, the Euclidean one for the Hamiltonian and the Lorentzian one for the Liouvillian are identical, upon setting τ = t. We also have seen that the Green's functions satisfy the same Schwinger-Dyson equations, and at low temperature the solutions will be the same. More discussion of this equivalence can be found in Appendix D.
We have thus established a complete equivalence between the two problems: the real-time path integral of an SYK coupled to a bath close to the steady state and the Euclidean path integral of a two-site non-Hermitian SYK related to wormhole configurations in the low-temperature limit.
For instance, the value of E g in the Euclidean problem will be the same as the decay rate Γ that describes the typical relaxation time to a steady state in the real-time problem. We have shown that the retarded Green's function in the real-time problem is equal to −iG LL − G LR in the Euclidean problem. We would like to remind the reader that, although we use ψ L , ψ R in both Liouvillian and Euclidean discussions, they are not exactly the same: the former is for real time and L, R represent different branches of the Keldysh contour, while the latter is for imaginary time and L, R denotes two sites of the system. Rather than using a Keldysh contour, we could also have solved the real-time problem as a two-site system.
In the remainder of the paper, we discuss the solutions of the saddle-point equations which illustrate the equivalence of the Keldysh approach and the Euclidean approach at low temperature.
We will start with a study of the q = 4 SYK model, whose dynamics is quantum chaotic, by solving numerically the saddle-point Schwinger-Dyson equations. We will then continue with analytical studies both for q = 2, whose dynamics is integrable, and in the large-q limit, whose dynamics is quantum chaotic.
Finally, we will explore the possible existence of a gravity dual and its relevance in the late-time dynamics of our real-time SYK model. We are motivated by the recent result that, in the lowtemperature limit, a similar non-Hermitian SYK model, with additional real couplings, is dual to a Euclidean wormhole [37]. This possibility was originally noticed in Ref. [35] and later formulated in detail in Ref. [37], see also Refs. [49,63,64,67]. We shall show that qualitatively similar results are found in our model and we will dub the possible gravity dual a Keldysh wormhole. Interestingly, this identification between strongly correlated dissipative quantum matter and quantum gravity strongly suggests that the observation of a gap or decay rate, even if there is no coupling to the bath, is closely related to the existence of a wormhole gravity dual.

REAL-TIME OPEN SYK FOR q = 4
We now proceed with the explicit calculation of Γ, E g that control the exponential decay of Green's functions in the real-time and Euclidean problem respectively, by solving the Schwinger-Dyson equations of both the Euclidean problem related to the non-Hermitian q = 4 SYK model, or the Lorentzian problem related to Liouvillian dynamics, which, as we have seen in the previous section, give the same results. For that purpose, we will compute the relevant Green's functions by solving the SD equations (30) and (31) at β = 500 and the SD equations (B30) and (B31), respectively.
Similar calculations have been carried out recently. In Ref. [57], the real-time calculation was carried out but no results were presented in the small-µ region of interest. In this limit, we shall see that the large-q calculation shows some deviations from the result obtained from the explicit solution of the SD equations for q = 4. The Euclidean calculation for a similar model, where the SYK couplings are not purely imaginary but also have a real part, was carried out in Ref. [67].
However, we shall see that E g and other observables are qualitatively different.
In Fig. 1, we depict the Euclidean, G LR (τ ), and retarded Lorentzian, G R (t)/2, Green's functions for different values of the coupling to the environment µ. As argued above, the two computations coincide. We stress that we are not making any analytical continuation but just compare Euclidean and Lorentzian times. For small µ, we observe oscillations in the Green's functions and a fast decay for long times. In both cases, as shown in the figure, we find an excellent agreement with a fit to The exponential decay, controlled by E g , persists for any coupling strength µ provided that temperature is sufficiently low. Following previous literature [35,37,48,67], we have denoted the decay rate as E g in reference to a gap. However, we stress that, especially for zero or small µ, this is not related to a gap in the spectrum between the ground state and the first excited state. It is rather the order parameter that, in the gravity dual, characterizes the wormhole phase. Therefore, it is a feature of the model resulting after ensemble averaging. A heuristic explanation of its origin is that the low-energy excitations in the partition function are nullified because the spectrum is complex and PT-symmetric. The ground state is not affected because it is real. As µ increases, E g approaches the spectral gap ∆. In Fig. 2, we show the gap E g (left), or decay rate Γ for the real-time calculation, and the frequency of the oscillations (right) of the Green's function as a function of the coupling µ. The curves for E g (µ), Γ(µ) and Ω(µ) are barely distinguishable which further supports that, in the T → 0 limit, the Hamiltonian (22) is identical to the Liouvillian (8) that represents an SYK coupled to a Markovian bath. Both E g and Ω have an intriguing µ dependence that we now analyze in detail.
For µ = 0, the gap is finite and the frequency of the oscillations in the Green's function is largest.
Physically, this corresponds to a complex order parameter whose real part controls the exponential decay and the imaginary part the frequency of oscillation. As µ increases, the imaginary part of the excited eigenvalues becomes smaller. Therefore, a larger µ is expected to reduce the value of Ω(µ). Indeed, we find that the frequency decreases monotonically with µ. Interestingly, Ω(µ) shows excellent agreement with a simple ansatz Green's function is purely exponential, which implies the vanishing of the frequency of oscillations.
This is consistent with the existence of a second-order phase transition at µ = µ c .
The gap or decay rate in the small-µ limit, which has not been investigated before, has also a rather unexpected behavior. The decay rate (or gap) is not zero even in the limit µ → 0 of the coupling to the bath. We stress that, as mentioned earlier, E g at µ = 0 is not related to the spectral gap. On the Euclidean side, an analogous feature has already been predicted [35] by noticing that the existence of a finite E g , which is a signature of the dual wormhole configurations, only requires that the variance of the left-right coupling is stronger than the left-left or right-right ones. This can be easily achieved with complex couplings in the SYK model. Indeed, the purely imaginary couplings of our model are sufficient since the variance of the left-left, or right-right, couplings are negative while the left-right ones are positive [35]. Strictly speaking, the gravity dual cannot be a Euclidean wormhole because the variance of real couplings (zero in our case) is not larger than the imaginary ones [35]. However, we shall see that our model still has wormhole-like features similar to those of a Euclidean wormhole. Its counterpart in the real-time calculation is the existence of a finite decay rate even if there is no coupling with the environment, which is expected provided that the SYK model is quantum chaotic (q > 2) [33,34,39,41]. We term the wormhole configuration in real time and µ 1, responsible for the finite decay rate, a Keldysh wormhole. We shall give a more precise characterization of this possible gravity dual in Sec. VI.

For a small but finite
This makes sense physically since a stronger coupling to the bath results in a larger decay rate. We think this linear behavior is related to the fact that the energy of the first excited state is always real and equal to 4µ. However, the overall effect on E g is relatively small compared with the µ = 0 contribution, namely, the spectral gap ∆ ≡ E 1 −E 0 , where E 0 and E 1 are the energies of the ground and first excited state, is only a small contribution to E g . The linearity suggests that the slope can be obtained by a perturbative treatment in µ. This is in contrast with the µ 2/3 behavior observed for a two-site SYK with µ 1 whose gravity dual is a traversable wormhole [35]. Therefore, we can rule out any quantitative relation between our setting and this gravity dual.
At a critical value µ c ≈ 0.15, corresponding to the vanishing of the oscillation frequency, the linear increase of E g stops abruptly. A further increase in the coupling µ leads to a sharp, likely continuous, drop in E g or Γ. This is rather unexpected as an increase in the coupling to the bath should always lead to an increase in the decay rate as the approach to equilibrium should occur faster. Similarly, a stronger explicit left-right coupling in the Euclidean problem is expected to give a larger gap because the hopping probability increases. However, we clearly observe a drop as µ 0.15 increases. A possible explanation of this nonmonotonic behavior is that the oscillations become imaginary, Ω → iΩ, for µ > µ c , so that they effectively become a negative contribution to the gap. More generally, we believe that this sharp decrease is directly related to the effective weakening of the strength of the imaginary couplings in each SYK due to the left-right term which is Hermitian in our model. For a sufficiently large µ, in agreement with Ref. [57], the decay rate Γ, or E g , approaches a linear dependence E g ∼ µ, which is the expected result [25] for a strong coupling to the environment. In this region, the spectrum is mostly real and E g is close to the spectral gap ∆.
In order to gain a more complete understanding of the physical meaning of E g and Ω, we computed the finite-N Green's functions, see Appendix E, by exact diagonalization (ED) techniques.
Most of our finite-N results have been obtained for N = 9, and at this point we are satisfied with a qualitative agreement between the finite-N results and the SD results. For µ < 0.15, we find a large difference between the average of the Green's function and the average of the absolute value of the Green's function. In both cases, we find an initial fast exponential decrease followed by a slower exponential decrease at longer times. However, the range of the initial exponential decay increases substantially upon ensemble averaging because of cancellations in the average of the Green's function. We expect that in the large-N limit, and in the limit of a very large number of realizations, where these cancellations become manifest, the initial exponential decrease will become dominant and will converge to the result obtained from the solution of the Schwinger-Dyson equations. Indeed, we find qualitative agreement with the results of this section, see Fig.   9 of Appendix E. In particular, we find a finite gap at µ = 0 which increases linearly for small µ.
The long-time tail remains dominant for the average of the absolute value of the Green's function.
For small µ its decay rate converges to the result given by the spectral gap, which is equal to , while for large µ the decay rates of the initial exponential decay and the longtime exponential decay become equal. At finite N , we observe oscillations realization by realization but, because their period fluctuates, they are mostly averaged out after ensemble averaging.
The detailed nature of the cancellation mechanism of the long-time tail is not clear at this point.
For µ > 0.175 no cancellations happen, but the initial decay rate is still different from the long-time decay rate. A better understanding of this difference requires a detailed finite-size scaling analysis.
Results on the finite-N scaling of the Green's function are relatively scarce in the SYK literature.
They have been studied to some extent for out-of-time-order correlators (OTOC) [68] but in that case the problem is that the exponential behavior is barely visible at finite N , while in the present problem the exponential decay is very clear for the long-time tail. Also the free energy, which is studied in Sec. VI, only shows a weak N dependence at low temperatures, much weaker than the one observed in a different non-Hermitian SYK model [64]. One of the reasons is that, in the present case, the finite-N ground state is reproduced identically by ED techniques. The free energy also has finite-N issues that are not well understood: for small µ, it has oscillations that are not observed in the solutions of the SD equations. Another potential source for the absence of the long-time tail in the SD equations is the averaging procedure. In the exact diagonalization, the average is done directly on the individual Green's function while in the SD analysis the average over disorder is taken at earlier stages of the calculation. At present, we cannot rule out completely that this is in part related to other solutions of the saddle-point equations. In any case, progress in this problem requires substantially larger values of N than those discussed in Appendix E and a more precise understanding of the relation between the spectral gap ∆ and the decay rate E g .
Finally, a natural question to ask is to what extent the anomalous behavior we have found in the weak-coupling limit is universal or whether it is a particularity of the SYK model. More specifically, is it related to the fact that the system is strongly interacting and quantum chaotic? In order to answer this question we turn to the study of a nonchaotic q = 2 SYK coupled to a Markovian bath described by the same Lindblad operator as for q = 4.

IV. ANALYTICAL COMPARISON BETWEEN THE EUCLIDEAN TWO-SITE SYK
AND THE REAL-TIME OPEN SYK FOR q = 2 The q = 2 SYK model is not quantum chaotic and it is not dual to any wormhole background.
Therefore, it is an ideal candidate to investigate whether the anomalous behavior observed for weak coupling to the environment is related to quantum chaos and field theories with gravity duals.
We start with the Euclidean problem at µ = 0, for which it is possible to find analytical solutions of the SD equations in Fourier space given by In terms of the G R/A variables introduced in Eqs. (34) and (35), the equations simplify to solution that vanishes at ω → ±∞ is given by with → 0 + , and the advanced Green's function is equal to Taking linear combinations of G R and G A , we have In the time domain, the retarded and advanced Green's functions are given by where J 1 is a Bessel function of the first kind. Note that these solutions are also well defined for finite > 0. Then, letting → 0, we find in agreement with the fundamental identity (32). The same result is obtained for real time. We observe, see Fig. 3, that oscillations persist as in the q = 4 case, but, in stark contrast to it, the decay of Green's function is power law, not exponential. Therefore, the gap vanishes for µ = 0.
We now turn to the region of nonzero µ. In fact, we already solved this problem because µ plays the role of at µ > 0. In this case, the saddle-point equations for G R/A are given by This results in the Green's functions Fourier-transforming to the time domain, we find  The analytical result for the real-time problem is identical with τ = t. Bottom left: −iG LL (ω) and −iG LR (ω) in Fourier space. The Green's function iG LR (ω)/π coincides with the spectral function of the real-time problem, given by Eq. (54). Bottom right: the gap E g as a function of µ. We find E g = µ for all values of µ which, in the real-time problem, is an indication that the process of relaxation is always dissipation-driven.
For completeness, we also write down the single-particle spectral function of the model: As before, the saddle-point equations for the real-time calculation are the same, resulting in the same solutions provided that boundary/initial conditions are not important (i.e., at low temperature in the Euclidean calculation and for times such that the system is close to the steady state in the real-time calculation).
In Fig. 3, we compare the analytical solution discussed above with the exact numerical solution of the SD equations for both the real-time and Euclidean-time problems. We find that both are indistinguishable from the analytical result without using any fitting parameters which confirms that also for q = 2 the real-time calculation agrees with the Euclidean calculation in the lowtemperature limit. However, unlike the q = 4 case, we do not observe an anomalously large value for the gap or the decay rate for small values of µ. For q = 2, E g = Γ = µ for all values of µ, so the approach to equilibrium is dissipation-driven. The reason is the integrable nature of the q = 2 SYK, which requires a coupling to the environment to reach a steady state. Therefore, in the absence of any coupling to the environment, (many-body) quantum chaos and strong interactions seem necessary conditions to reach a steady state.
Finally, we note that, for q = 2, the analytic prediction for E g is in agreement with the spectral gap. In this case, the finite-N spectrum is given by (for k = 0, 1, . . . , N ) with e k,n real, and it is realization dependent except for the single real eigenvalue for each value of k. The gap between the first excited state energy with zero imaginary part and the ground state energy is thus equal to µ. So, unlike q = 4, the resulting spectral gap is equal to the decay rate of the Green's functions, E g = ∆ = µ. The observed oscillations in the Green's functions are determined by the imaginary parts of λ k,n .

V. ANALYTICAL COMPARISON BETWEEN THE EUCLIDEAN TWO-SITE SYK AND THE REAL-TIME OPEN SYK IN THE LARGE-q LIMIT
We now show that the equivalence between Euclidean and Lorentzian problems also holds in the large-q limit [34,35]. In order to carry out the large-q analysis for the two-site non-Hermitian SYK model, we follow closely the method introduced in Ref. [34] for the one-site SYK model and applied in Ref. [35] to a two-site SYK model dual to a traversable wormhole.
The Green's function in the large-q limit can be written as To leading order in 1/q, the SD equations simplify to sign(τ )(∂ 2 τ g LL (τ ) + 2J 2 e g LL (τ ) ) = 0, where J 2 = 2 q−1 J 2 /q and µ =μ/q. To satisfy the causality relation (32) at low temperature we need that which requires that g LL (τ ) = g LR (τ ). In Eq. (58), theμδ(τ ) term does not contribute because of the sign(τ ) factor. The boundary conditions are [35] g LL (0) = 0, ∂ τ g LR (τ = 0) = −μ, and g LL (τ ) − g LR (τ ) → 0 as τ → ∞. We note that the above SD equations and boundary conditions coincide with those of the real-time problem of the open single-site SYK model [57]. In the lowtemperature limit, the solution is given by with α 2 /J 2 cosh 2 γ = 1, 2α tanh γ =μ. Using Eqs. (56) and (61), the gap is given by which is also the result for the decay rate of the real-time one-site SYK coupled to a bath [57]. The µ = 0 limit of the decay rate was derived earlier in Ref. [34].
We know that it is expected [35] that this large-q limit is qualitatively different from the q = 4 case, which is also quantum chaotic. For instance, Green's functions have no oscillations and the gap is monotonic with µ. Still, some similarities remain, and the gap tends to a constant in the µ → 0 limit, which indicates that for weak coupling the gap and the decay rate are fully controlled by internal dynamics and not by the coupling to the bath. Likewise, in the µ → ∞ limit, the decay rate is fully controlled by the environment.
A similar calculation, see Appendix F, leads to the free energy in the large-q limit at low temperature T : which is the value that is obtained from the ground-state energy of the Hamiltonian discussed in the next section, see Fig. 4.

VI. GRAVITY DUAL OF A DISSIPATIVE SYK: KELDYSH WORMHOLE
We now turn to the tentative gravitational interpretation of the results obtained above. As was mentioned earlier, it is well established [35,37] that, in the low-temperature limit, the dual field theory of two-dimensional anti-de Sitter (AdS 2 ) wormholes is a two-site SYK. Depending on the setting, the gravity dual can be a traversable or a Euclidean wormhole. Since the path integral describing the dissipative SYK is identical to the path integral of the two coupled non-Hermitian SYK models, its gravity dual will be the same. Therefore, the findings of this section could shed light on the gravity dual of a single-site SYK coupled to an environment in real time.
Our two-site SYK model is qualitatively similar to that recently investigated in Ref. [67], which reported a Euclidean-to-traversable wormhole transition by tuning the intersite coupling. However, there is an important difference: the random couplings in our model are purely imaginary while those of the model of Ref. [67] are complex with the variance of the real part larger than the variance of the imaginary part. Therefore, the low-energy effective action of our model cannot be easily identified with a Euclidean or a traversable wormhole. The reason is that in our model the coupling J in the Schwarzian term of the low-energy effective action becomes complex, J L → iJ and J R → −iJ, with respect to that in the Euclidean and traversable cases. This is a direct consequence of performing the ensemble average over a two-site SYK Hamiltonian that is purely imaginary after vectorization. Intriguingly, this complexification of couplings is precisely the necessary prescription [69] to obtain the correct Schwarzian action describing a double-trumpet configuration in a nearde Sitter background in two dimensions (dS 2 ) [69] in Lorentzian time from the equivalent in a near-AdS 2 in Euclidean time [44,45]. If this is the case, our setting would correspond to that of a wormhole in this near-dS 2 geometry. Indeed, the so-called bra-ket wormholes [70,71] 1 with a similar double-trumpet geometry have been identified in near-dS 2 backgrounds coupled to conformal matter. These wormholes are connected geometries between bra (ket) states prepared forward (backward) in time following the Keldysh protocol. Finally, it has been speculated that the SYK dual of a one-boundary near-dS 2 background [72] in Lorentzian time, not related to wormholes, would require the same coupling complexification J → iJ. Beyond the explicit comparison of the low-energy effective action, a potential way to identify the gravity duality more explicitly is to  For the latter, we use a different convention than for the SD calculation: {ψ i , ψ j } = 2δ ij and J 2 ijkl = 6J 2 /(2N ) 3 . For small µ, we identify a first-order phase transition where the low-temperature phase is characterized by a flat free energy. For large µ, in line with the result observed for Euclidean [37] and traversable [35] wormholes, the transition becomes a crossover. For intermediate temperatures above the first-order transition, we observe oscillations only in the exact diagonalization calculation which are suppressed as µ increases. The change in slope observed for intermediate µ at a temperature above the first-order transition, which likely signals a higher-order transition, is reproduced by both calculations.
compare the gravitational quasinormal modes in a near-dS 2 setting and the decay rate and frequency that characterizes the approach to a steady state in the dissipative SYK. For an interpolating AdS 2to-dS 2 geometry, the study of quasinormal modes was carried out in Ref. [73,74]. More recently [75,76], quasinormal modes leading to a purely exponential decay-with no oscillations-of the retarded Green's function have been computed in a dS 2 static patch. It has been conjectured [75,76] that the field theory dual may be a certain Hermitian single-site SYK in the double scaling limit [77]. We note that this is qualitatively different from our non-Hermitian two-site SYK model whose retarded Green's function in the small µ limit also has oscillations and an exponential decay that does not require large q even for µ = 0. Therefore, it is expected that the gravity dual of both models is qualitatively different as well.
We now turn to confirm the wormhole features of the low-temperature limit of the two-site non-Hermitian model by the study of the free energy, where wormhole configurations are characterized by an almost flat free energy at low temperature that ends abruptly at a certain temperature where a first-order phase transition takes place. In the first two rows of Fig. 4, we present results for the free energy of our model by solving the SD equations for different values of the coupling µ. We note that the free energy is obtained [34] by inserting the solutions of the SD equations back in the action. The free energy has two branches from different solutions of the SD equations. For very low temperatures, and finite but not very strong coupling, we observe a flat part typical of wormhole configurations. This branch eventually crosses a fast decreasing branch, signaling the presence of the first-order transition consistent with a wormhole-black hole transition [35,48]. For a stronger coupling µ, the transition becomes a crossover, a feature that has been observed in the free energy corresponding to a two-site SYK model dual to wormhole configurations [48,67]. For intermediate temperatures above the first-order transition, where the two branches coexist, the free energy undergoes a close to linear decrease. At a higher temperature, close to the end of the wormhole branch, the slope of F (T ) in the high-temperature phase changes rather abruptly. This is consistent with the existence of a higher-order phase transition separating this linear temperature behavior from the black hole phase. We note that around the same temperature at which this change in slope occurs, there is also a qualitative change in the Green's functions: G LL (τ ) and G LR (τ ) have finite real and imaginary parts in this range of temperatures, while they are, respectively, purely real and imaginary elsewhere. These results for the free energy and the ones of Appendix F strongly overlap with those of Ref. [78], which was posted on the arXiv a few days after this paper. corresponding to the region of coexistence of two solutions in the large-N calculation above, we observe oscillations that are suppressed as µ increases. This is an indication of a pathological behavior characterized by a negative specific heat and a negative entropy. For µ 0.125, there are no oscillations and the free energy of the large-N calculation is reproduced qualitatively including the change in slope in the high-temperature limit mentioned above. Further research will be needed to find out both the physical mechanism causing the change of slope and the origin of the observed oscillations in the finite-N free energy. Overall, the first-order transition is rather similar to that reported in other SYK settings [67] where the wormhole gravity dual is explicitly known.

are oscillations of Euclidean
Green's functions after analytic continuation to real time. The frequency of the oscillations is related to E g . For Euclidean wormholes, G LR and G LL are in phase [67], while for traversable wormholes, they are out of phase [66,79]. Since we have that G LL (τ ) = −i sign(τ )G LR (τ ), the Green's functions of our Euclidean model in real time will oscillate in phase.
Although we cannot identify precisely the dual gravity background, we summarize, based on results in this and previous sections, its main features: a gapped ground state, with a gap that has a nontrivial dependence on the couplings but which is different from that of Euclidean or traversable wormholes; oscillating Green's functions in real time with |G LL (t)| = |G LR (t)|, as observed in Euclidean wormholes; and a flat free energy in the low-temperature limit that terminates in a firstorder phase transition (for higher temperatures, we observe a change of trend in the free energy consistent with a higher-order phase transition). Moreover, as mentioned above, the Schwarzian action of our model should correspond to that of bra-ket wormholes [70] in a near-dS 2 background.
Since these features are similar to those corresponding to Euclidean or traversable wormholes, but with important qualitative differences, we have termed Keldysh wormholes the field configurations with wormhole-like features that govern the path integral describing the time evolution of an SYK model coupled to an environment.

VII. CONCLUSIONS
We have shown that in the near-zero-temperature limit, the dynamics of a two-site non-Hermitian SYK model in Euclidean time are equal to the real-time evolution of the Liouvillian describing a one-site SYK model coupled to an environment close to the infinite-temperature steady state. The two-site model has features, like Green's functions that decay exponentially with a decay rate that is finite even if there is no coupling to the environment and a free energy with a first-order phase transition, consistent with those found earlier in SYK models with a wormhole gravity dual.
Indeed, the low-energy effective action is that of a Euclidean wormhole but with the coupling J in the Schwarzian term complexified, J → iJ. This complexification of the coupling has been related [69,80] to wormholes in a near de-Sitter background and we have termed it a Keldysh wormhole.
However, we stress that a precise identification of the gravity dual requires further research.
We have also shown that the decay rate of the open SYK is finite even if the coupling to the environment is zero provided that the dynamics is quantum chaotic (for q = 4). Only in the limit of strong coupling to the bath, the approach to the steady state is controlled by the environment.
Interestingly, only in the region of weak or vanishing coupling do Keldysh wormholes control the path integral. An analytical study of the large-q limit, which is also quantum chaotic, confirms that the decay rate is finite even if there is no coupling to the environment. Another distinct feature of these wormhole configurations is that Green's functions show damped oscillations whose amplitude decreases as the coupling to the bath increases. For a sufficiently strong coupling, the oscillations vanish and, later, the decay rate becomes a monotonic function of the coupling. In contrast, for a dissipative integrable SYK model (q = 2), the decay rate, which we have obtained analytically, is always dominated by the coupling to the environment. Indeed, contrarily to the quantum chaotic case, it vanishes if there is no coupling to the bath, showing no anomalous relaxation. Also unlike the quantum chaotic case, the pattern of oscillations is observed for any coupling.
Natural extensions of this work are a more detailed analysis of the low-energy effective action in Euclidean signature in order to establish a precise relation to wormhole configurations and de Sitter holography [69,80]. We also plan to study the relation between the decay rate and the spectrum The idea of the Choi-Jamiolkowski (CJ) transformation is to write an operator as a vector. In

bra-ket notation
If the operator A acts on the space V , then, after the CJ transformation, A is represented as a vector in V ⊗ V . The first space will be labeled by L and the second one by R.
The Lindblad operator L acts on the density matrix ρ. After the CJ transformation, ρ becomes a vector in V ⊗ V , so that L maps to an operator acting on V ⊗ V . The basis of V will be denoted by |k ⊗ |l or, in short, |k |l . After the CJ transformation, we thus have To work out the Lindblad operator after the CJ transformation, we first need to derive what is known as operator reflection.

Operator reflection
The idea of operator reflection is to find, for a given A L , the operator A R satisfying where K is the complex conjugation operator and U is unitary. For the system of interest A L = ψ L and A R = αψ R with α a constant, and we determine U and α such that (A3) is satisfied. The representation of γ matrices is taken to be with γ k the gamma matrices for the left SYK model, and γ c the chirality matrix. Note that with this choice, A R in (A3) also acts on the first |j in (A3). This leads to the requirement (repeated indices summed) for all l and m. This can be simplified to Therefore, we must have that Since the charge conjugation matrix C satisfies we choose This gives the condition αγ c e iπγc/2 = 1, We conclude that The overall sign from C T = ±C is irrelevant and has been ignored.

The vectorized Lindblad operator
a. The MQ term in the Lindblad operator. Let us consider the L n ρL † n term in the Lindblad evolution. We take L n = √ µψ n with ψ n a Majorana operator. The CJ transformation in this case is given by (repeated indices are summed) µ i|ψ n ρψ n |j |i j| → µ i|ψ n L ρψ n L |j |i |j , where the factor U K has been absorbed in the second |j of |j |j . Using completeness and the reflection property (A12), we obtain µ i|ψ n L ρψ n L |j |i |j = −iµψ n L ρ|j ψ n R |j = −iµ k|ψ n L |l l|ρ|j m|ψ n R |j |k |m = −iµ m|ψ n L |l k|ψ n R |j |m |k l|ρ|j .
Comparing to Eq. (A2), we can read off the operator acting on V ⊗ V : After the CJ transformation, the −iH SYK ρ term becomes After the CJ transformation, the iρH SYK term becomes (for even q) This gives the term c. The constant term. For − 1 2 µ ψ n † ψ n ρ + ρψ n † ψ n we use that (note the summation convention) with N the number of L (or R) Majorana fermions. Then, we obtain This gives −(N µ/2)1 L ⊗ 1 R in the SYK operator.

Appendix B: Details on the real-time evolution of the open SYK model
The partition function is with action and Liouvillian As discussed in the main text, we now move to the closed-time contour representation (C = C + ∪C − ), On the Keldysh contour the action reads as (note that dz = dt on C + and dz = −dt on C − ) with dissipative kernel Here, we have dropped any dependence on the initial conditions, since we are interested in the long-time relaxation, and the immaterial constant −N µ/2. Performing the disorder average over the couplings J i 1 ···iq we obtain Introducing the collective field and its Lagrange multiplier Σ(z, z ), the action reads as The Schwinger-Dyson (saddle-point) equations on the contour are given by We now return to real times (t 1 , t 2 ). Different components of the real-time Green's function are obtained by restricting (z, z ) to the two branches C + and C − : with a, b = +, − corresponding to the respective contours. At the saddle point, we introduce the greater, lesser, time-ordered, anti-time-ordered, retarded, advanced, and Keldysh components, respectively, as The components of the contour self-energy Σ(z, z ) satisfy the same relations as the Green's functions. In the long-time limit T = (t 1 + t 2 )/2 → ∞, the system loses any information about the initial condition as it relaxes to the steady state. Time-translational invariance in T emerges and the Green's functions depend only on the relative time t = t 1 − t 2 . Moving to Fourier space with continuous frequencies ω using the convention we define the real quantities and their Hilbert transforms The spectral function ρ − (ω) is symmetric and normalized as dω ρ − (ω) = 1. Because the jump operators are Hermitian, the system relaxes to the maximally mixed state. Close to the steady state, ρ + (ω) and ρ − (ω) are related by a fluctuation-dissipation-like relation, ρ + (ω) = tanh(β TFD ω/2) ρ − (ω), where β TFD =0 is the parameter of the infinite-temperature TFD steady state. Hence, ρ + (ω) vanishes identically. In terms of the two-site non-Hermitian SYK model, the vanishing of ρ + can be traced back to the PT symmetry of the model [49]. Note that for more general jump operators that lead to finite-temperature equilibrium or nonequilibrium steady states we have, in general, ρ + (ω) = 0.
In terms of ρ ± , the components of the Green's function read as [in the last equality of each line we use the special property ρ + (ω) = 0 of our model] Exactly the same relations hold for the self-energies σ ± .
In terms of the quantities ρ − (ω) and σ − (ω), the Schwinger-Dyson equations (10) and (11) can be rewritten as [56] ρ − (ω) = σ − (ω) where (ρ − ) * n (ω) denotes the n-fold convolution of the spectral function with itself. In real time t, our main quantity of interest is G R (t) which can be obtained from the spectral function by Appendix C: Details on the Euclidean evolution of the two-site SYK model The partition function can be expressed as a path integral Integrating out the random couplings J i 1 ,··· ,iq satisfy the normal distribution introducing G ab (τ ) and Σ ab (τ ) with identities and finally integrating out the fermions, we obtain the partition function in terms of G ab and Σ ab variables, and the effective action where a, b ∈ {L, R} and s LL = s RR = 1, s LR = s RL = (−1) q/2 , t LL = t RR = 1, t LR = t RL = −1.
Changing variables according to the time derivative in this action appears as i∂ τ and it directly maps to the real-time action (B8).
This will be discussed in more detail in the next appendix.
The equations of motion, the Schwinger-Dyson equations, follow from the variation with respect to G ab and Σ ab and are where * denotes convolution.
Appendix D: Relation between the Euclidean and real-time calculations As established above, the relation between the real-time Majoranas on the Keldysh contour and the Euclidean-time Majoranas is given by The real-time and the Euclidean Green's functions are defined by respectively, where a, b = −, + in real time and a, b = L, R in Euclidean time. Using the identification (D1), we find the following relations between the components (with time-translation invariance) These identifications are consistent with the symmetries of the Green's functions in both real and Euclidean time: G < (t) and G LR (τ ) are imaginary and even (about t = 0 and τ = 0, respectively); G T (t) is imaginary and odd and G LL (τ ) is real and odd. In order to keep the term ΣG in the action invariant, the self-energy components must be identified as where the minus sign in the first relation results from the backward time direction in the integration.
Let us now compare the Schwinger-Dyson equations. In real time they are given by where the equations in the first line are in frequency space and the ones in the second line are in the time domain.
In Euclidean time we have Using the identifications of Eq. (D3) in Eqs. (D6) we recover the saddle-point equations (D5).

Appendix E: Exact finite-N results
In this appendix, we discuss exact diagonalization results obtained for a finite number N of Majorana particles. We expect that for large N these results will approach the results obtained from the SD equations. Before discussing the numerical finite-N results, we first discuss the exact analytical results we have used as checks of the numerical finite-N calculations, i.e., the eigenfunctions and eigenvalues of the ground state and the first excited state, and the large-ω behavior of G LR (ω) in the low-temperature limit.

Analytical finite-N results
The Hamiltonian is given by It has the ground state (see Appendix A) which satisfies and The excited states of iµ k ψ k L ψ k R are obtained by acting with raising operators on the ground state. The lowest excited states are given by which can be shown by separation of the k = m and k = m terms. The next excited state is obtained by acting twice with the raising operator on the ground state, etc. Since we can take as a basis for the eigenstates of the Hamiltonian.
Numerically we find that, for small µ and sufficiently large N (N > 8), the energy of the first excited state is given by This excited state |1 is the next eigenstate of iµ k ψ k L ψ k R which can be similarly annihilated by suggesting that |1 is an entangled state composed of states |j |j , or alternatively can be expressed as a linear combination of Indeed, this has been observed numerically.
Next we discuss analytical checks for the Green's functions. As discussed in the main text, in the low-temperature limit we must have that G LL (τ ) = −iG LR (τ ) for 0 < τ β/2 and β → ∞.
The symmetry of G LL (τ ) and G LR (τ ) gives the relation for −β/2 τ < 0. Combining the two we The values at τ = 0 ± follow from the properties of the gamma matrices. The result for G LL is given by which are the Green's function in free limit, as is expected. The time-ordering operator has been denoted by T , which for fermionic operators satisfies T ψ(0)ψ(t) = −ψ(t)ψ(0) for t < 0. For the calculation of G LR (0) in the low-temperature limit we use the property (ψ k L + iψ k R )|0 = 0 to obtain These values for G LL and G LR agree with what we observe in finite-N calculations and for solutions of the SD equations.
As a last check we discuss the asymptotic large-ω behavior of G LR (ω), the Fourier transform of is equal to the Green's function for the ith particle G ii LR (τ ), which reads In the low-temperature limit we have contributions from τ /β 1 and (β − τ )/β 1. Using the approximation, we obtain for τ /β 1 Because G LR (τ ) = −G LR (β − τ ) the contribution for (β − τ )/β 1 is given by −G LR (β − τ ).
The Fourier transform G ii LR (ω) of both contributions adds up to We have used that for fermionic Matsubara frequencies, ω n β = 2π(n+ 1 2 ), so we have exp(−iβω n ) = −1. Therefore, We expand this result to second order in 1/ω. The first-order term cancels and from the secondorder term we obtain To evaluate this expectation value we use that as well as To derive this result we have used the commutation rules of the fermion operators and properties of the ground state, Our final result for G ii LR (ω) is given by 2. Numerical finite-N results for q = 4 In this section, we discuss the low-temperature finite-N results for G LR (τ ) obtained by exact diagonalization. The Green's functions can be evaluated using Eq. (E14) and the completeness relation of the eigenstates of H. Care must be taken that the right eigenvectors are only orthogonal to the left eigenvectors. The right eigenvectors satisfy where D is a diagonal matrix with the eigenvalues of H. Inserting V T V T −1 after the exponential into Eq. (E14) within the region τ /β 1, we obtain The solutions of the SD equations (blue curve) show the same large-ω behavior but deviate from the finite-N results at low and intermediate frequency.
Next, we study the large-time behavior of G LR (τ ). In Fig. 6 we show results for N = 10 and two values of µ (see legend) and compare them to the result from the SD equations. We observe a discrepancy between the SD results and the finite-N behavior of the Green's function at long times, especially when µ is small. At this moment we do not have a complete understanding of this discrepancy. It might have its origin in the large-N limit of the spectrum, but we do not see this from the finite-N calculations. As an example, in Fig. 7 we give scatter plots of spectra for N = 12 and various values of µ. Next, we will show that the discrepancy is mainly due to cancellations . Spectrum for µ = 0.02, µ = 0.2, µ = 0.5, and µ = 0.7, N = 12, and q = 4. In the small-µ limit, most eigenvalues are complex and are located close to the imaginary axis, but low-lying eigenvalues remain real. When µ gets larger, the number of real eigenvalues increases, but the gap in units of µ between the first excited state and the bulk of the spectrum decreases. For large µ, we observe vertical stripes corresponding to the eigenvalues of iµ k ψ k L ψ k R .
that become manifest only for a large number of realizations.
In Fig. 8, we show results for N = 9 and various values of µ (see legends). For small µ, we can distinguish two exponents, an early rapid decay and a much slower decay at large times. The exponent of the early rapid decay is in qualitative agreement with the results of the Schwinger-Dyson equations, but the exponent of the long-time tail is definitely not in agreement with the SD equations for µ < 0.15. We have performed the same analysis for N = 11 and N = 13, which shows a similar behavior of G LR (τ ).
To obtain a better understanding of the discrepancy between the long-time tail of the Green's function obtained from the SD equations and from exact diagonalization, we plot in Fig. 8 Figure 9. The gap E g as a function of µ for N = 9 (red curve) obtained from the initial exponential decay of the Green's function. The green line shows the analytical result for the spectral gap for small µ. The results obtained from the SD equations are represented by the black curve. The blue curve represents E g obtained by fitting the long-time exponential decay of the Green's function. For µ = 0.075 the average is over 25,000 realizations while for other values of µ in the range 0.066 ≤ µ ≤ 0.1125 we used 10,000 realizations. For all other points we used 4000 realizations.
by E g , and the long-time exponent (green curve), denoted by E g , for N = 9 are shown in Fig. 9.
The gap is obtained by fitting exp(−E g τ ) to the initial decay (green line in Fig. 8). For the decay rate of the long-time tail, we fit exp(−E g τ ) (red line in Fig. 8) to the numerically obtained long-time tail of G LR (τ ). Also shown is the result obtained from the solution of the Schwinger-Dyson equations (black curve). We observe that the result for the Schwinger-Dyson equations is in qualitative agreement with the exponent of the initial decay rate. For small µ, the gap in the spectrum is well defined and is given by (see green line in Fig. 9) The finite-N numerical results for the long-time tail are in good agreement with this analytical prediction, but we did not find solutions of the SD equations with this behavior.
Our finite N results show that the cancellations that take place for µ < 0.15 increase the range of the initial exponential behavior, but it requires a very large ensemble to realize this. We expect that, in the large-N limit, this initial exponential behavior will become dominant. A detailed finite-size scaling analysis required to show this will be deferred to future work.
Appendix F: Large-q limit in the Euclidean problem Following the method in Ref. [35], the SD equations can be solved analytically in the large-q limit at low temperature. The main idea is to solve the equations in small-τ and large-τ region with proper boundary conditions, then match these two solutions to fix the undetermined parameters.
Recalling the Schwinger-Dyson equations, where J 2 = qJ 2 /2 q−1 , we first notice that they give rise to the free solutions in the τ → 0 which leads to Σ LL = − J 2 q e q−1 q g LL ≈ − J 2 q e g LL , Substituting the expressions for G ab and Σ ab into SD equations, and applying an additional partial derivative, the leading order terms cancel and at order 1/q we obtain sign(τ ) ∂ 2 τ g LL + 2J 2 e g LL (τ ) = 0, ∂ 2 τ g LR + 2J 2 e g LR (τ ) + 2μδ(τ ) = 0, withμ = µq.
The boundary conditions for g LL and g LR at τ = 0 are which follow from the free limit of G LL at τ → 0, and the integral over the second saddle-point equation of Eq. (F4) over an infinitesimal region τ ∈ [0, ], respectively. In addition, when T → 0, g LL and g LR satisfy the boundary condition which will select the exponentially decaying solutions. At finite but low temperatures, g LL and g LR are not equal, and we need to solve the SD equations to find the relation between g LL and g LR .