Euclidean-to-Lorentzian wormhole transition and gravitational symmetry breaking in the Sachdev-Ye-Kitaev model

We study a two-site Sachdev-Ye-Kitaev model with complex couplings and a weak inter-site interaction. At low temperatures, the system is dual to a Euclidean wormhole in Jackiw-Teitelboim gravity plus matter. Interestingly, the energy spectrum becomes real for sufficiently strong inter-site coupling despite the Hamiltonian being non-Hermitian. In gravity, this complex-to-real transition corresponds to a Euclidean-to-Lorentzian transition: a dynamical restoration of the gravitational SL(2,R) symmetry of the Lorentzian wormhole, broken to U(1) in the Euclidean wormhole. We show this by identifying an order parameter for the symmetry breaking and by matching the oscillating patterns of the Green's functions. Above the transition, the system can be continued to Lorentzian signature and is dual to an eternal traversable wormhole. Additionally, we observe a thermal phase transition from the wormhole to two black holes and provide a detailed matching of the associated physical quantities. The analysis of level statistics reveals that in a broad range of parameters the dynamics is quantum chaotic in the universality class of systems with time reversal invariance.

A theme that has emerged in recent years is that the semiclassical gravitational path integral appears to be dual to an ensemble average of theories [22,23]. The meaning of this average remains to be understood [48,49], see [50][51][52][53][54][55][56][57][58][59][60][61] for a sample of recent discussions on this point. The inclusion of complex metrics appears necessary but it is unclear how to specify the contour of integration [62][63][64][65]. In this paper, we study the gravitational path integral in a setup where we have analytic control on the gravity side and a dual microscopic description in terms of SYK, in order to shed some light on the rules governing the gravity path integral.
A notable example of nearly AdS 2 holography is the eternal traversable wormhole [26] dual to a two-site SYK model with a weak inter-site coupling. It was shown that both the wormhole solution in JT gravity and the two-site SYK model were described by the same Schwarzian effective action. The ground state was argued to be gapped and close to a thermofield double state. The gap induced by the weak inter-site coupling is enhanced by the strong interactions in each site, as can be observed in the real time probability of tunneling between the two sites [27,66] . As the temperature is increased, the system experiences a first order phase transition to a phase with two black holes. This transition affects the quantum dynamics [67] that is quantum chaotic only in the high temperature phase. Extensions of these results include replacing Majoranas with Dirac fermions [68][69][70][71] or the use of a sparse [72,73] two-site SYK [74] while potential applications in condensed matter were addressed in [75].
A purely Euclidean version of this story was obtained in [76] by studying a two-site non-Hermitian SYK model with complex couplings but without inter-site interaction. In this case, the Hamiltonian is non-Hermitian and does not define a Lorentzian system with unitary evolution. Nonetheless, the system can be studied as a purely Euclidean system, from the point of view of statistical mechanics. The low temperature phase was shown to be dominated by replica symmetry breaking configurations [77,78] corresponding to a Euclidean wormhole in JT gravity [76]. Here, the imaginary part of the Hamiltonian gives imaginary sources in the wormhole. As there is no coupling between the two sites, the Euclidean wormhole has to be the result of the SYK average. This leads to a factorization puzzle which can be resolved by finding half-wormhole solutions [54] which realize the proposal of [50] in nearly AdS 2 holography. Recent studies on the non-Hermitian SYK model include a symmetry classification of quantum chaotic dynamics, [79], a generalization with additional U(1) charge [80], measurement-induced transitions [81,82] and Lindbladian approach to dissipative quantum dynamics [83,84].
In this paper, we study a two-site SYK model with both complex couplings and a weak intersite interaction. This system combines the two effects discussed above and is dual to a Euclidean wormhole in JT gravity. Our main result is the observation of a dynamical complex-to-real transition where the spectrum of the Hamiltonian (obtained by exact diagonalization) becomes real for sufficiently strong inter-site coupling, despite the Hamiltonian being non-Hermitian. This shows that a purely Euclidean system can experience a transition above which the system has real energy spectrum and hence defines a Lorentzian system with unitary evolution.
In JT gravity, we show that this complex-to-real transition corresponds to a Euclidean-to-Lorentzian transition: the dynamical restoration of the SL(2, R) symmetry of the Lorentzian wormhole, broken to U(1) in the Euclidean wormhole. To be more precise, there are two distinct SL(2, R) symmetries which come from the isometries of the Lorentzian wormhole (the global AdS 2 geometry). There is an SL(2, R) gauge symmetry, which is part of the gravitational constraints of the theory, and an SL(2, R) global symmetry, corresponding to the physical action of the isometries on the asymptotic boundary [85,86]. These are gravitational symmetries because they are part of the diffeomorphism group of the theory. Both symmetries are broken to U(1) in the Euclidean wormhole. The Euclidean-to-Lorentzian transition corresponds to a restoration of the gravitational SL(2, R) symmetry (both gauge and global) in the purely Euclidean system (defined by imposing only the U(1) gauge constraint).
We study this complex-to-real transition by identifying, both in SYK and JT gravity, an order parameter for the gravitational symmetry breaking. We propose a mechanism for the transition involving a change of contour in the gravity path integral. The transition can also be diagnosed in a phase shift measured by the Green's functions. Above the transition, the Hamiltonian is pseudo-Hermitian, i.e. it has real eigenvalues despite being non-Hermitian, and defines a Lorentzian system with unitary evolution in the framework of pseudo-Hermitian quantum mechanics [87]. This shows that a pseudo-Hermitian Hamiltonian can be holographic, and is dual here to an eternal traversable wormhole.
In addition to the complex-to-real transition, we also observe a thermal phase transition from the wormhole phase to a phase with two black holes. This transition has been discussed in [26,76] for limiting cases of our system. Here, we study the parameters of this transition, such as the energy gap and critical temperature, as functions of the two parameters λ and κ. We find an excellent match between the numerical SYK results and the analytical JT analysis. Note that this transition constitutes another example of gravitational symmetry breaking where the U(1) × U(1) symmetry of the two black holes is broken to the diagonal U(1) of the Euclidean wormhole.
In this work we view JT gravity as an effective theory, an approximation of the (still elusive, see e.g. [88]) exact holographic dual of SYK, in the same way that semi-classical gravity is viewed in higher-dimensional examples of AdS/CFT. This has to be contrasted with the exact path integral computations of [23] where JT gravity is viewed as an exact theory, dual to a random matrix model (RMT). The imaginary sources allow us to evaluate the path integral using saddle-points which is closer to what we can expect in higher dimensions. The gravitational symmetry breaking/restoration described in this paper, interpreted as a Euclidean-to-Lorentzian transition, might also be possible more generally. In higher dimensions, Euclidean wormholes with similar properties were constructed in [89] and static traversable wormholes were obtained in [90,91]. Finally, we study the level statistics in the SYK model and find that the dynamics is quantum chaotic in a broad range of parameters as it is well described by the random matrix predictions for either real or complex eigenvalues.

II. THE TWO-SITE NON-HERMITIAN SYK MODEL
The system we consider consists of two SYK models with complex couplings, labeled left (L) and right (R), and with a weak inter-site interaction. The Hamiltonian is with where ψ L i , ψ R i , i = 1, . . . , N are Majorana fermions. After ensemble average, and using the standard decoupling procedure with we obtain the following action, where a = L, R, b = L, R and t LL = t RR = 1, t LR = t RL = −1.
In the large N limit, the saddle-point Schwinger-Dyson (SD) equations are given by, where the first two equations are expressed in the frequency domain, while the last two are in the imaginary time one.
A. Thermodynamic properties In our system, there are two parameters λ and κ. The parameter λ controls the strength of direct hopping between the two systems. As mentioned earlier, the wormhole phase is characterized [26] by a non-trivial dependence of the gap on λ. By contrast, κ controls the strength of the imaginary part of each separate SYK Hamiltonian [76,77] and it is not directly related to the coupling between the two systems.
In this section, we study the combined effect of the two couplings in the thermodynamic properties of the system. More specifically, we compute the free energy and the spectral gap E g as a function of temperature and the parameters λ and κ.
The free energy is obtained from the on-shell action where the Green's functions in the action are given by the solutions of the saddle-point SD equations (4), We know that for λ = 0 and finite κ [76,77], the system will stay in the wormhole phase, whose low-temperature free energy is lower than that of black hole phase. In the wormhole phase, the free energy is independent of the temperature.
The reason of the existence of the wormhole phase even without an explicit coupling term (λ = 0) can be understood from the SD equations by rewriting them as where Solutions of these equations are Green's functions with an exponential τ dependence. Specifically, if we assume G LR is purely imaginary and exponential around τ = 0, G LR will be non-zero only in the neighborhood of τ = 0. Therefore, we can make the approximation where λ ef f can be identified as the effective coupling constant, and its sign depends on that of Im[G LR ]. Comparing with the SD equations in [26], such approximation will lead to the exponential solutions of G LL , iG LR ∈ R characteristic of the wormhole phase. Indeed, the gap E g that separates the wormhole ground state from excited states can be extracted from the exponential decay of G LR , E g = − lim τ →∞ log |G LR |/τ .    Figure 2. G ab for κ = 0.75, λ = 0.02 and T = 0.005 for the two branches shown in Fig. 1. Different sign of G LR will enhance/reduce the effective coupling, leading to the same behavior for F for different signs of λ. When λ is positive, the left solution is chosen, and vice versa.    Figure 4. (a) Free energy F as a function of temperature T for various κ's and λ = 0.03; (b) Free energy F as a function of temperature T for various λ's and κ = 0.5. In all cases, we observe a first order phase transition with a critical temperature that is an increasing function of λ and κ.
We now proceed to investigate the case when λ is also turned on. In principle, from the form of the SD equations, one may think that the effect of a finite λ on the free energy may depend on its sign and it may not always enhance the gap that characterizes the the wormhole phase.
However, our results clearly indicate that this is not the case. Finite λ will always lower the free energy for any given κ. The reason is that, for given λ, κ, T , there are three solutions and only the one with the lowest free energy will be chosen since it is the dominant saddle point solution when evaluating partition function. Two among these three solutions, denoted by G LR and G LR , give approximately constant free energy and represent the wormhole phase, the third one corresponding to the black hole phase. When λ = 0, two wormhole solutions satisfy G LR and they give rise to the same free energy. However, if λ = 0, this equality becomes approximately correct, LR , and the free energy will split, so one of the solutions for G LR increases the free energy while the other decreases it for any temperature or κ (see results depicted in Fig. 1).
Specifically, when λ > 0, the solution with positive imaginary part Im[G LR (τ )] > 0 in τ ∈ (0, ), as illustrated in the left diagram of Fig. 2, is preferred for lower free energy, while λ < 0 we have opposite selection for the same reason. Therefore, finite λ can always make the free energy lower by choosing proper solutions. We will see in the next section that the same mechanism takes place in JT gravity: two wormhole saddle-points are exchanged under a change of sign of λ, ensuring that thermodynamical quantities are even functions of λ.
In the large N limit, the path integral, and therefore the free energy, is dominated by the saddlepoint configurations given by the solutions of the SD equations. The physical solution is the one with the lowest free energy. Therefore, as illustrated in Fig. 3, the free energy will always decrease as λ is increased from zero, so the wormhole phase becomes more thermodynamically stable by combining the effect of λ and κ.
Having provided a qualitative description of the impact of a finite λ, we now proceed to a more systematic analysis of the free energy F (T ) for various κ and λ. In all cases, see Fig. 4, there always exists a first order phase transition separating a nearly flat free energy in the low temperature limit from a high temperature phase which is approximately linear. The critical temperature is a increasing function of κ and λ reinforcing the idea that the two effects are additive regarding the stability of the wormhole low temperature phase.
Likewise, the dependence of E g on κ and λ in the low temperature limit is summarized in Fig. 5 and Fig .6. Inspired by the analytical dependence of E g on λ and κ in the κ = 0 [26] and λ = 0 [76,77] limit, we first test whether the combined effect of both couplings is E g = Aλ However, this relation seems to work only in the small κ region. From the study of the gravity dual, we shall derive an analytical expression for E g which is in agreement with the ansatz above for small κ and also describe well E g for other values of λ and κ, as shown in Fig .23.
We also study the critical temperature T c (κ, λ) of the first order phase transition. In Fig. 7, we observe that for sufficiently small λ and large κ we have T c ∼ κ 4 . This is consistent with a recent analytical prediction [77] for λ = 0. In the gravity section, we will derive an analytic expression of T c given in (51) which is in agreement with these numerical results even at finite λ.  B. Spectrum, spin symmetry and complex-to-real transition We now study qualitative features, and the impact of symmetries of the spectrum of the Hamiltonian (1). In the λ = 0 limit, the Hamiltonian has an approximate spin-like symmetry represented byŜ = N i ψ L i ψ R i , which is exact for κ = 0. For a certain κ = 0, eigenvalues of the Hamiltonian tend to cluster around those of the operatorŜ for sufficiently large λ. This complicates the calculation of spectral correlations, for ifŜ strictly commutes with the Hamiltonian, the spectral analysis must be restricted to eigenvalues within each cluster. In Ref. [67], this problem was solved by breaking this spin symmetry completely by considering couplings in each of the SYK's that differ by an overall constant α = 1. If α 1, the wormhole phase still survives [26] but the gap becomes smaller.
Another solution is to choose a basis in which the Hamiltonian is block-diagonalized where each block corresponds to an eigenvalue ofŜ. Even if the symmetry is broken for κ = 0, there still exists the parity symmetry which corresponds to a spin-like operator e iπS [67]. For numerical convenience, we choose the following Hamiltonian, which is equivalent to (2): Here, E is the identity matrix with the size 2 The parity operator ψ c is conserved([ψ c , H] = 0) with eigenvalues ±1. So we only need to reorder diagonal elements in descending order, then use exactly the same ordering to reorganize the Hamiltonian into block-matrix form of two parity sectors. Thus we can carry out the level statistic analysis on one of the two blocks separately. In case we are interested to study low temperature properties, we must choose the block that includes the ground state.
In order to assess the importance of the spin-symmetry mentioned above in our Hamiltonian (1), we represent in Fig. 8 the spectral density for κ = 1 and different values of the explicit coupling λ. We observe a rather symmetric distribution for λ = 0 which is in agreement with the results of Ref. [76,77,79]. As λ increases, gaps around the real axis start to form. The spectral density for larger values of λ shows that a growing part of the spectrum becomes real but we still observe complex eigenvalues in certain regions. The existence of the latter is directly related to the spin operatorŜ which is an approximate symmetry of the Hamiltonian for sufficiently large λ. Complex eigenvalues are restricted to the area between nearby eigenvalues ofŜ on the real axis, and hence their number depends on N . Real eigenvalues starts to cluster around the eigenvalues of the spin operator (which we do not plot explicitly). By contrast, the maximum density of complex spectrum is located around double cones with tips on the real axis between the nearby eigenvalues ofŜ.
Interestingly, a further increase in λ leads to an unexpected result. Even though the Hamiltonian is non-Hermitian, the whole spectrum becomes real for λ > λ c (κ), see Fig. 9. Upon a further increase in λ, the spectral support of the already real spectrum is split in separate intervals. For λ large enough, we observe that these intervals are centered around the eigenvalues of the spin operator.
As was shown in Ref. [67], this clustering can be shown analytically by taking λ 1 so that the other terms in the Hamiltonian are a small perturbation.
In Fig. 10(a), we present the ratio of complex and real eigenvalues for different κ and λ.
We also find that though the percentage of the real spectrum increases with λ monotonically, this increase slows down when the percentage of complex eigenvalues is small (< 10%). In Fig. 10 we depict λ c (κ), the minimum λ for which all eigenvalues are real for a given κ. As was expected, it shares similarities with the contours in Fig. 10(a).
We also note that the existence of a critical λ = λ c is not an approximate result: the imaginary parts of the eigenvalues are strictly zero within the numerical precision 10 −15 . We shall see in the following section that this transition has an observable impact on the oscillation patterns of real time Green's functions, which are related to quantum tunneling for κ = 0. Moreover, the transition does not require large N or disorder average, although λ c is sensitive to the disorder realization.
It is instructive to consider the case N = 4 for which the complex-to-real transition can be seen explicitly. In this case we have simply H L = −(J + iK)ψ 1 ψ 2 ψ 3 ψ 4 where J and K are arbitrary real numbers. The Hamiltonian is easily diagonalized and the eigenvalues are where the first eigenvalue is threefold degenerate and the third one is fourfold degenerate. The complex-to-real transition in this case corresponds to the fact that the spectrum becomes real for λ > λ c = K.  An order parameter for the complex-to-real transition is the thermal expectation value of the imaginary part of the Hamiltonian, which is defined as where Z = Tr e −βH . This is an order parameter for the transition since Im H β vanishes at the transition point, because the spectrum of the Hamiltonian becomes real. This quantity can be derived from the partition function so it will be possible to compute it in the next section using the gravity path integral.
We plot the order parameter in Fig. 11. It initially increases, then exhibits a rather sharp decrease for some small λ from the effect of the spin-operator. It increases again within the range of λ ∼ (0.06, 0.12), then abruptly decreases to nearly zero and finally vanishes at λ = λ c . The value of λ c is consistent with the one obtained from the level statistics or the analysis of the real time Green's functions. It is also possible to study the higher-moments of Im H by taking additional derivatives with respect to κ. The fact that all the moments vanish at λ = λ c then implies that the spectrum becomes real.

C. Tunneling amplitude and real time evolution
In this section we study the evolution in real time of the SYK model (1) by solving the SD equation in real time. In this context, a similar study was first carried out in Ref. [69] for the case of a two-site Hermitian SYK model dual to traversable wormholes. We first present results for G > ab (t) when κ = 0 or λ = 0 respectively. We then provide a heuristic description of its main features. Finally, we carry out a detailed numerical computation of the combined effect of a finite λ and κ in real time Green's functions. Note that although the system only makes sense in Euclidean signature, it is a well-defined procedure to analytically continue the Green's function to real time.
Before proceeding, we firstly review the results for κ = 0 and finite λ as discussed in Ref. [69].
The quantum dynamics in this case is controlled by the SD equations in real time resulting from the analytical continuation of (5), We define whose Fourier trasnform is with a, b = L, R. Results for κ = 0, the case already studied in Ref. [69], are depicted in Fig. 12.
As is observed, |G > LL | and |G > LR | are out of phase, with maxima (minima) of one the functions corresponding to minima (maxima) of the other. The maxima of G > ab appear in different t, which can be understood, when the system is dual to an eternal traversable wormhole, as the propagating time through the bulk. These real time results in Fig. 12 are directly related to the imaginary solutions in Ref. [26] by applying a Wick rotation.
We have found out, see Fig. 12(b), that the largest two peaks of ρ ab (ω), have the lowest frequency and are symmetric with respect to ω = 0. Therefore, ρ ab is well approximated by with A = 0.5. By employing the definition for the retarded Green's function, we obtain Therefore, the imaginary time Green's functions are given by, As is illustrated in Fig. 13, G ab (τ ), calculated from (19), is real for a = b = L or a = b = R and purely imaginary for a = L, b = R. Moreover, |G ab | decays exponentially in the small τ region, which is consistent with previous numerical results [26]. For finite κ, we expect [77] that G ab (τ ), and therefore ρ ab (ω), has also similar properties.
The results for the case λ = 0, κ = 0, depicted in Fig. 14, also display an oscillatory behavior but with a crucial difference: both functions are now in-phase. Physically, it is an indication that there is no real tunneling between left and time particles but rather a synchronization of the dynamics of both sites despite the fact that they are not directly coupled. The only coupling arise after ensemble average. This in-phase pattern will be explained in the next section as a localization in the gravity path integral.
Since |G > ab (t)| ∼ | ψ a (t)ψ b (0) | is the probability amplitude of observing a particle again at time t if this particle appears at t = 0 once, the overlap of |G > ab (t)| with a, b = L, R indicates the synchronization of the dynamics of the left and right sites even though they are not physically coupled.
Another interesting feature in the λ = 0 limit, see Fig. 14, is that ρ LL (ω) has negative peaks.
Superficially, this result looks surprising since for > 0, all the peaks by definition must be positive from the very definition of ρ ab (ω). However, this applies only to κ = 0 where the Hamiltonian is Hermitian and the eigenvalues are real. For κ = 0, the eigenvalues are in general complex and therefore, peaks can be either positive or negative.
We now turn to the study of the combined effect of a finite κ and λ. In Fig. 15, we can clearly observe that, for a fixed κ, the behavior of |G > ab (t)| are different for small and large λ. For small λ, |G > ab (t)| are similar to that in the case of λ = 0, namely, both Green's functions are in-phase. Likewise, ρ ab (ω) have a simple structure: two leading peaks and two subleading peaks whose sign depends one whether a = b or a = b are dominant. However, when λ is sufficiently large, the oscillations of |G > ab (t)| become qualitatively different, as for κ = 0, |G > LL (t)|, |G > LR (t)| are out-ofphase though the oscillating pattern becomes more intricate. This is fully consistent with ρ ab (ω) where more subleading peaks are observed. More specifically, a finite λ (κ) is responsible for positive (negative) subleading peaks in ρ LL (ω). These peaks are suppressed when λ gets larger. We could not study the nature of the transition between these two regimes because, for a fixed κ, we could not find solutions in this critical region λ ∼ 0.1. A possible reason is that enhanced oscillations around the transition makes difficult to find |G > ab (t)| numerically. Finally, we check the consistency of our results by comparing the temperature dependence of the real time Green's function with the thermodynamic properties investigated previously. Results for F (T ) and E g (T ), depicted again in Fig. 16(e) indicate that the thermodynamic phase transition temperature occurs at T c ∼ 0.055. Previous real time Green's functions were considered in the low temperature region T → 0 where the system is in the wormhole phase and the pattern of oscillations of |G > ab (t)| is quite rich. As temperature is increased (T > 0.02), we observe a gradual suppression of oscillations. This suppression is directly related to a broadening of the peaks in ρ ab (ω). This behavior has already been predicted in Ref. [92] for the κ = 0 case. Oscillations disappears already for T = 0.0666 which suggests a transition to the black hole phase.
Finally, we study the effect of sign switch λ → −λ in the real time evolution. In Fig. 17, we depict |G > ab (t)| and ρ ab (ω) for positive and negative values of λ. As is observed, the sign flip in λ just changes the sign of ρ LR and leaves |G > ab (t)| invariant. This is reasonable since we expect that λ → −λ induces and overall sign difference in G R LR (t) → −G R LR (t). One interesting question is whether the observed transition between different oscillation patterns is accompanied with a change in the free energy F . The answer to this question is negative. The reason is that the exponential behavior of the imaginary time Green's function G ab (τ ) is related to the leading peak of the corresponding ρ ab (ω). However the oscillating patterns depend instead on the superposition regarding the leading and all the subleading peaks. Since F is a function of G ab (τ ), we do not expect F to experience any significant change during the complex-to-real transition. Indeed, the thermal phase transition is different from the complex-to-real transition. As we will see next in gravity, the former is a transition between the wormhole and two black holes while the latter can be interpreted as a Euclidean-to-Lorentzian transition in the wormhole phase.

III. GRAVITY DUAL
In this section, we study the gravity dual of the SYK system described by the Hamiltonian (1).
At low temperatures, it is dual to a Euclidean wormhole in JT gravity with two parameters η and k, respectively dual to λ and κ. This is a generalization of the eternal traversable wormhole [26] (corresponding to k = 0) and the Euclidean wormhole without interaction [76] (corresponding to η = 0). The combined effect of η and k studied in this paper leads to a wormhole with similar thermodynamical properties and we will find an excellent match between the SYK and JT results.
Our system is a purely Euclidean system studied from the point of view of statistical mechanics.
It doesn't have a Lorentzian interpretation (with unitary evolution) as the energy spectrum is generally complex. Note that the Euclidean quantities can still be viewed as suitably analyticallycontinued versions of Lorentzian observables, which is akin to studying a partition function at imaginary value of the chemical potential, see e.g. [93,94].
The main result of the previous section is the observation of a complex-to-real transition, where the energy spectrum becomes real for sufficiently strong inter-site coupling, despite the Hamil-tonian being non-Hermitian. The Euclidean wormhole has a gravitational U(1) symmetry while the Lorentzian wormhole has a gravitational SL(2, R) symmetry [85,86]. We will show that the complex-to-real transition corresponds to the dynamical restoration of the SL(2, R) symmetry of the Lorentzian wormhole, and thus can be interpreted as a Euclidean-to-Lorentzian transition.

A. Wormhole solutions
The theory we consider is Jackiw-Teitelboim gravity with a massless scalar field χ and a static Gao-Jafferis-Wall interaction [26,95] involving N fields of dimension ∆. To compare with SYK, we should take ∆ = 1/4. The theory is described by the action where The solution we consider is the Euclidean wormhole Following [76], we deform the theory with boundary sources taken to be imaginary: The imaginary sources model the imaginary part of the SYK couplings as they correspond to a deformation of the Hamiltonian where O L , O R are the marginal operators dual to χ on each boundary. We see that this is a good model of the imaginary part of the SYK Hamiltonian (1) as we can identify the marginal operators with the SYK operators The boundary conditions for JT gravity are and we study the theory in Euclidean signature with the periodicity conditioñ u ∼ũ + β .

Schwarzian effective action
Nearly AdS 2 holography is a theory of a boundary graviton, or reparametrization mode, which can be described by an effective Schwarzian action [12][13][14]. The wormhole has two boundaries so it is described by an action for two reparametrization modes τ L (u) and τ R (u) after integrating out the matter degrees of freedom.
Without boundary sources (k = 0), the system is the eternal traversable wormhole whose action was derived in [26] in both SYK and JT gravity. The boundary sources give an additional contribution that can be computed by evaluating the action for a general solution for χ in the wormhole as a function of boundary sources χ L and χ R : where the bulk-to-boundary propagators in AdS 2 are The value of the sources chosen here are The contribution of the scalar field is then obtained by evaluating the on-shell action after acting with the diffeomorphisms corresponding to the two Schwarzian modes [14], see Appendix A of [54] for additional details.
At the end, the effective Schwarzian action of the system takes the form which we study in Euclidean signature. We see the effect of the boundary sources is to add a constant term in the action proportional to the wormhole size b. We use the physical time u which is related to the coordinateũ via and which is periodically identified u ∼ u + β where β =βN is the physical temperature. The coupling constants are which are taken fixed in the large N limit. As in [26], the validity of the action requires that η 1 and that the system develops an approximate conformal symmetry close to the ground state, which will be assumed here. The gravitational regime corresponds to large N . The overall factor of N ensures that the path integral localizes on its saddle-points in the large N limit.
We have not been able to derive the Schwarzian action directly from the SYK model because, unlike in JT gravity, it is harder to split the real and imaginary couplings which enter in a rather symmetric way. Nonetheless, we expect, due to the strikingly similar properties of both systems, that the Schwarzian effective action will be the same in SYK, with k proportional to κ. This is also expected from universality if we view the Schwarzian action as a type of effective hydrodynamics, the contribution from the imaginary sources corresponding to a marginal deformation.

Wormhole solutions
The Euclidean wormhole corresponds to the solution which gives the Euclidean action This action needs to be minimized with respect to the wormhole size b. For this purpose, it is useful to introduce the variable X defined from the relation so that the action becomes where we have set ∆ = 1 4 . The action is a quartic polynomial in X. The saddle-point in X gives a cubic equation This equation is also equivalent to the vanishing of U(1) charge which is required as the U(1) symmetry is a gauge symmetry. In fact, (40) is the integrated Hamiltonian constraint of JT gravity and implies the gravitational equations of motion. For JT gravity with matter, it takes the form [85,86] on a Cauchy slice Σ. Here, the functional measures the energy at each boundary in terms of the boundary graviton. As an aside, we note that this form is similar to the integrated Hamiltonian constraint in higher-dimensional AdS. In [96], this was used to prove a perturbative version of the holography of information [97]. This suggests that a similar statement should be possible in JT gravity with matter for excitations of the eternal traversable wormhole, i.e. on the solution corresponding to the global AdS 2 geometry.
The cubic equation can be solved analytically using Cardano's method [98]. The three roots can be written as where j = e 2iπ/3 and The discriminant of the equation vanishes when η = η c with We will argue that η c is the counterpart of λ c in SYK, the critical value for the complex-to-real transition.
These three solutions give rise to three wormhole solutions which, for lack of a better terminology, we will refer to as the first, second and third saddle-points. We can see from Fig. 18 that the size b remains real for |η| < η c but can become complex above the transition, even though the dominant solution (in the canonical ensemble) always have real b. The appearance of similar complex saddlepoints was observed in [54] and will be related here to the complex-to-real transition.

Free energy
The free energy of the wormhole is We see that the free energy is independent of the temperature which reflects that the phase is gapped. The real and imaginary parts of the free energies for the three wormhole saddle-points are plotted in Fig. 19. For |η| > η c , the free energy of the two subleading wormholes become complex, although the total free energy remains real. This reflects the fact that these subleading wormholes become complex geometries as b acquires an imaginary part.
For η > 0, we use here the dominant solution X = X 1 . For η < 0, we should use X = X 2 as the two saddle-points get exchanged. This mechanism was also observed in SYK and implies that the thermodynamic quantities will be symmetric under η → −η as illustrated in Fig. 21. For this reason, it is enough to focus on the region η > 0.

Energy gap
The energy gap is given, for η > 0, by [26] where we take ∆ = 1/4.  Explicitly, the energy gap takes the form The energy gap for various values of η and k is plotted in Fig. 20. It is similar to previous SYK results. A more quantitative comparison by using k = Aκ, η = Bλ between the JT and SYK parameters where the coefficients can be fixed in various ways. In Fig. 23, we compare E g in gravity (48) with the numerical SYK predictions for E g using A, B as fitting parameters for λ ∈]0, 0.09] and different λ. For κ sufficiently small, we find an excellent agreement between the gravity and the SYK predictions for A ≈ 1.585 and B ≈ 27.056. Details of the fitting are given in Appendix A.
The fact that k = Aκ should be taken small is a consequence of the scaling regime (34). Ask should be fixed in the gravitational theory, k = Aκ scales as N −1/2 and must be small in the large N limit. For larger κ, we observe deviations in the SYK model, for example due to the fact that the Schwarzian terms are renormalized by a factor 1 − κ 2 . This effect is subleading in the large N limit. It should be possible to interpret it as a subleading (e.g. one-loop) effect in gravity but in this work we focus on the leading large N behavior.

Transition to two black holes
The system with the chosen boundary conditions has another saddle-point corresponding to two black holes with free energy In the canonical ensemble, the dominant solution is the one with the smallest free energy. We observe a phase transition at low temperatures where the two black holes become a wormhole. This transition was already studied in limiting cases in [26,76].
The critical temperature is the value T c for which the free energy of the dominant wormhole is equal to the free energy of the two black holes where we assume S 0 1. For η > 0, the dominant saddle-point is X 1 in (43). Explicitly, the critical temperature takes the form For negative η, the dominant saddle-point is the second one, so X 1 should be replaced by X 2 in the expression of T c .
Results for T c (k), for different values of η are depicted in Fig. 22. The critical temperature obtained from SYK in Fig. 7 has a similar behavior. A quantitive comparison is performed in  B. Gravitational symmetry breaking

Lorentzian phase
In JT gravity, the SL(2, R) symmetry corresponding to the isometries of AdS 2 has to viewed as a gauge symmetry. In the path integral formulation, this is because configurations related by an SL(2, R) symmetry should be viewed as equivalent [14]. From a Lorentzian point of view, this SL(2, R) symmetry is part of the diffeomorphism group. As a result, the associated charges must vanish Q 0 = 0, This is part of the gauge constraints of the theory and can be derived, for example, by integrating the gravitational constraints (or equations of motion) on a Cauchy slice.
Imposing Q ± = 0 implies that τ L (u) = τ R (u) ≡ τ (u) so the left and right boundary modes get identified [26]. The last constraint gives in terms of the Liouville variable ϕ(u) = log τ (u) and where the potential is In this case, the Euclidean action can be written as a Liouville action so we see that indeed the vanishing of Q 0 is equivalent to the equations of motion.

Euclidean broken phase
In the Euclidean wormhole, we don't have to impose the constraints Q ± = 0 since these generators are not isometries of the geometry. This is the identification τ ∼ τ + b breaks the SL(2, R) symmetry to U(1). We should still impose Q 0 = 0 which gives explicitly where E is the energy functional defined in (42). We can see that this corresponds to a system of two interacting boundary modes τ L (u) and τ R (u).
For η = 0, the system has a simpler description in terms of Liouville variables ϕ = log τ given by the equation which corresponds to two correlated (but non-interacting) Liouville particles.

Phase shift
In the broken phase, additional classical solutions are possible because we don't impose τ L = τ R .
In particular, we can have a shift by an arbitrary constant α This corresponds to a global U(1) symmetry, which can be interpreted as relative shifts between the two sides. This symmetry is present in the Euclidean wormhole at η = 0 but is explicitly broken by the inter-site coupling.
In the path integral, α enters as another moduli on which we should integrate. After evaluating the action on the classical solutions (58), we should perform the path integral over b and α. We can consider doing first the integral over α. This is where we highlighted the α-dependence. This integral has a saddle-point at α = 0 and this will be the dominant solution in the Euclidean path integral. As a result, the effect of non-zero α is subleading in the thermodynamics and cannot be easily measured in Euclidean signature.
Rather, we will see that the variable α can be measured in the Wick-rotated two-point functions.
We first define and the classical solution (58) gives the contribution The effect of α can be studied by continuing to real time. This is a well-defined procedure, even though we don't expect the theory to always have a Lorentzian interpretation. The Wick rotation can be achieved by using u 1 − u 2 = iv and we obtain where we used ∆ = 1 4 and defined α = iα. This shows that α measures the phase-shift between G LL and G LR in real time. Here, the frequency of the oscillations is given by The Wick rotation also acts on α since it is defined as the difference of two times. In other words, the real time Green's functions should be computed by integrating over the contour α ∈ iR. For α = 0, the Green's function are in-phase while forα = π they are out-of-phase.
In fact, since the integrand is periodic, we should only integrate on the circle α ∼ α + 2iπ. In terms of the variableα = −iα, the integral becomes where the three saddle-points in the b-integral corresponding to the three roots (43).
Note that for η = 0, the three saddle-points reduce to a single one. In this case, the different choices of α are exactly degenerate which corresponds to a global U(1) axial symmetry. To explain this, note that translations on τ L and τ R give a U(1) × U(1) symmetry. The inter-site interaction only preserves the diagonal U(1) diag symmetry (which correspond to the U(1) isometry of the Euclidean wormhole) and explicitly breaks the U(1) axial symmetry (which is the one that shifts α).
In the path integral, this soft breaking corresponds to the fact that the inter-site interaction gives a different action for different values of α.
For η > 0, the dominant saddle-point is the first saddle corresponding to X 1 . The effect of α can be accounted by replacing η with η/ |cos(α/2)| starting from theα = 0 configuration. We can see that the action is strongly localized aroundα = π as we have The divergence atα = π implies that we cannot use a saddle-point approximation here. Rather, the factor e −S 1 inserts a delta function in the path integral which localizes it on the valueα = π. This value implies that G LL and G LR should be in-phase. For η < 0, the first and second saddle-point are exchanged, it is then S 2 that dominates and inserts a similar delta function in the path integral.
This localization mechanism explains that the Green's functions are in-phase for small η. For η > η c , the restoration of SL(2, R) symmetry described below imposes τ L = τ R . As a result, we must haveα = 0 so that G LL and G LR must be out-of-phase.
These two regimes for the Green's functions are plotted in Fig. 25

Order parameter
The deformation by boundary sources ±ik corresponds to adding to the Euclidean action where we have defined the operator Writing Z as a path integral shows that its expectation value can be obtained by taking a derivative with respect to k: where we used that k = 2 πNk .
We will see that the expectation value of this operator is an order parameter for the SL(2, R) symmetry breaking. To study the SL(2, R) transformation, we formally continue to Lorentzian signature using τ = it and use that under an infinitesimal SL(2, R) transformation This leads to so this operator is invariant under U(1) but transforms non-trivially under the other generators of The path integral computation of O LR requires a suitable choice of contour for the integral over b. The general procedure to selects the contour is not well understood, see for example [62][63][64][65]. In principle, the contour could be determined by the appropriate analytic continuation from Lorentzian signature if the system has a Lorentzian interpretation. In the theory defined by the Euclidean path integral, Picard-Lefschetz theory could be used to understand the correct contour prescription. Here, we will give a natural contour prescription that reproduces the SYK results, and leave a better understanding of the choice of contour for future work.
For |η| < η c , we integrate over b ∈ [0, +∞) which is the contour in Fig. 26a. This gives a non-zero expectation value where b * is the size of the dominant wormhole. This non-zero expectation value spontaneously breaks the SL(2, R) symmetry to U(1).
For η > η c , we propose that we should choose the vertical contour in Fig. 26b. This leads to The result is depicted in Fig. 27. We see that at η > η c , the order parameter vanishes. The transition also happens in the negative η region, with the first (blue) and second (orange) saddlepoints exchanged. We can compare this to Fig. 11 in SYK and we see a good qualitative match.
This shows that the change of contour appears to be the right gravity mechanism to account for the complex-to-real transition observed in SYK. Note that the different choices of contour should reflect the ambiguity in the derivative of the partition function with respect to k due to the branch cuts appearing in (44).
Note that there is a larger value of η max = 3 3 2 k 3 above which the classical solution for the small wormholes cannot be trusted. This is the value at which the complex conjugated wormholes have Re b = 0. When Re b β, one-loop effect cannot be ignored and will modify the analysis.
Such one-loop effects were studied in [76] for η = 0. In this paper, we stay in the classical limit which is sufficient for the comparison with SYK results in the large N limit.
This mechanism explains the surprising observation in SYK that the Hamiltonian becomes real for |η| > |η c |. Indeed, we have from (11), which shows that this is, up to the factor β, the same order parameter that was previously identified in SYK. The higher moments of Im H can be obtained by taking more derivatives with respect to  For |η| > η c = k 3 , we observe a transition where O LR = 0 and the SL(2, R) symmetry is restored. This is in good agreement with the SYK result Fig. 11. Note that the additional peak in the latter is due to spin symmetry which is a non-universal featured related to the chosen minimal left-right coupling, k. The same mechanism suggests that these higher moments vanish above the transition, and the vanishing of all the moments implies that the energy spectrum must be real. This shows that the preservation of SL(2, R) symmetry implies that the energy spectrum must be real.
To summarize, we have proposed a mechanism for the restoration of SL(2, R) gauge symmetry in terms of a change of integration contour. We have shown that the restoration of SL(2, R) symmetry is equivalent to the energy spectrum becoming real by identifying an order parameter for the transition. At the moment, we cannot fully justify the change of contour in gravity, as the precise rules governing the path integral are not well understood. We leave a better understanding of this mechanism for future work. In terms of the variable X = b/β, the path integral takes the form The Picard-Lefschetz theory of similar integrals was considered in [99] and one might hope that this could shed some light on this transition.

C. Euclidean-to-Lorentzian transition
In this section, we explain in what sense the complex-to-real transition can be understood as a Euclidean-to-Lorentzian transition.
The Lorentzian wormhole (eternal traversable wormhole) corresponds to the global AdS 2 geometry. The SL(2, R) isometry group of the background translates into an SL(2, R) gauge constraint Q ξ = 0 for any Killing vector ξ. In general, we can write where Σ is a Cauchy slice with volume form ε ν . Here T µν contains a gravity and matter part obtained by varying the action with respect to the metric. In JT gravity, T grav µν is the stress-tensor of the JT dilaton (it is proportional to the Einstein tensor in higher dimensional gravity). The equations of motion imply that T µν = 0. This implies that the SL(2, R) charges have to vanish: The Euclidean wormhole is obtained by doing the Wick rotation from global AdS 2 and identifying periodically the time coordinate. This identification breaks the SL(2, R) symmetry to U(1) so we should only impose the vanishing of U(1) charge: In Euclidean signature, we must view this as a constraint on the configurations entering in the path integral. We see that the purely Euclidean system, where we only impose Q 0 = 0, has more configurations than the Lorentzian system. In the Schwarzian theory, this corresponds to configurations with τ L = τ R while we must have τ L = τ R in Lorentzian signature.
The system we study is viewed as a purely Euclidean system dual to a non-Hermitian Hamiltonian. The Euclidean-to-Lorentzian transition is the restoration of SL(2, R) symmetry of the Lorentzian geometry. Using the order parameter, we have seen that this implies that the energy spectrum is real. Conversely, a real spectrum implies the existence of a Lorentzian continuation with unitary evolution, so the SL(2, R) symmetry must be restored. This shows that the complex-to-real transition in SYK corresponds to a Euclidean-to-Lorentzian transition in JT gravity.
The SL(2, R) symmetry discussed above is a gauge symmetry in JT gravity that is part of the gravitational constraints of the Lorentzian wormhole. There is also a global SL(2, R) symmetry which corresponds to the isometry group of AdS 2 acting on the matter sector, appropriately dressed to commute with the SL(2, R) gauge symmetry. This global symmetry exists because the isometries are large diffeomorphisms at the asymptotic boundaries, see [85,86]. In the Euclidean wormhole, we similarly have a U(1) gauge symmetry accompanied by a U(1) global symmetry. The SL(2, R) → U(1) breaking/restoration involves both the gauge and global symmetry. These symmetries are gravitational in origin become they come from the isometries of the semi-classical background.
Although the system can always be viewed as a Euclidean wormhole, the restoration of SL(2, R) symmetry implies that above the transition, it can be continued to Lorentzian signature and is dual to an eternal traversable wormhole. 1 The Hamiltonian is non-Hermitian but its eigenvalues are real so it can be used to define unitary evolution with a suitably modified inner product [87].
In other words, the eigenstates of the Hamiltonian are not orthogonal with respect to the usual scalar product but they become orthogonal with a new inner product. Intriguingly, this indicates that pseudo-Hermitian Hamiltonians can appear in holography.
A final comment is that this Euclidean-to-Lorentzian transition is only interesting in the wormhole regime. For T > T c , the system is dual to two black holes. In this case, the Euclidean and Lorentzian symmetries are both equal to U(1) × U(1) so there can be no symmetry breaking/restoration.

IV. LEVEL STATISTICS AND QUANTUM CHAOS
We end the paper by investigating the nature of the quantum dynamics for long time scales of the order of the Heisenberg time.
The real level statistics in the κ = 0 limit, was addressed in Ref. [67]. Sufficiently far from the wormhole ground state, the long time dynamics is quantum chaotic as level statistics agrees well with the random matrix prediction [100]. More specifically, it agrees with the Gaussian Orthogonal Ensemble which corresponds to systems with time reversal invariance. Unlike single-site SYK models whose global symmetries depend in general on the number of Majoranas, for an SYK model with two identical sites, time reversal invariance is always present. However, the fact its low energy excitations deviate strongly from this universal result suggests that, the low temperature phase transition between the wormhole and the two black holes is accompanied by a qualitative change in the quantum dynamics.
Level statistics in the limit λ = 0 has also been investigated recently [79]. The level statistics of the combined system is trivially Poisson because both SYK are not explicitly correlated. However, the spectral correlations of each SYK separately agrees well with the predictions of non-Hermitian random matrix theory. This SYK model, depending on the number of Majoranas N and the q−body interacting Hamiltonian, reproduces many of the different universality classes predicted [101] in non-Hermitian quantum chaotic systems.
We now study the combined effect of a finite λ and κ in the spectral correlations. More specifically, we aim to clarify whether a weak explicit coupling λ 1 is enough to make the dynamics of the combined non-Hermitian system quantum chaotic, at least for sufficiently high energies.
Therefore the level statistics will be well described by random matrix theory. This is important as a further confirmation that even in a non-Hermitian setting, quantum black holes are related to quantum chaotic motion [102]. Likewise, we would like to explore whether the observed deviations from random matrix theory that characterize the quantum motion in the real case are also present in our model. It is also of interest to investigate whether in the region λ > λ c , where the spectrum is real, a finite κ is of any relevance in the description of the level statistics.
In order to avoid the spin-symmetry mentioned previously, we will perform the similarity transformation mentioned above and diagonalize numerically each parity block separately. We employ the complex spacing ratio as a spectral observable [103], that does not require the unfolding of the spectrum. This is especially important for a two-dimensional spectra where the unfolding procedure suffers in some cases from ambiguities. The complex spacing ratio is a short-range spectral observable that probes the quantum dynamics for times scales longer than the Heisenberg time, which is originally introduced to study correlations of real spectra [104][105][106]. In the complex case We do observe the halfeaten doughnut shape typical of random matrix behavior [103]. However, sizable deviations are observed as λ increases because part of the spectrum becomes real and the complex spacing ratio is no longer applicable.
it is defined as where E k is the complex spectrum for a given disorder realization, E NN k is the nearest eigenvalue to E k and E NNN k is the next to nearest eigenvalue to E k . In order to eliminate statistical fluctuations, we carried out ensemble-averaging until we obtain at least 10 6 eigenvalues for each choice of parameters N, λ, κ. In Fig. 28 and Fig. 30(a), we depict the distribution of the complex spacing ratio z k for κ = 1 and different values of λ. For larger values of λ, it is necessary to use a better ensemble averaging to obtain similar results, since from Fig. 10, the real eigenvalues become less with λ larger. As was expected, for λ = 0, we do not see any trace of revel repulsion for small distances.
However, even for very small values of λ = 0.015, the characteristic [103] half-eaten donuts shape that indicates level repulsion and potential quantum chaotic behavior is clearly visible. 2 In Fig. 29, similar results are obtained for different values of N and κ > 0.2. There exists obvious deviations from RMT for sufficiently small κ. This is expected because the RMT results assume a spectral density more or less locally symmetric in the complex plane. However, in the κ → 0 limit, the spectral correlations are greatly enhanced along the real line so they cannot be described by complex gap ratio. As a result, it is necessary to use a minimum κ min so that for κ > κ min , the complex part of the eigenvalues is much larger than the mean level spacing and non-Hermitian RMT results apply.
The distribution of the spacings does not allow a quantitative comparison with random matrix predictions. For that purpose, it is more convenient to use the angular ρ(θ) and radial ρ(r) distributions [79,103] of the complex spacing ratios (θ and r corresponding to the angular and radial variable in polar coordinates).
We observe in Fig. 30(b,c) that, for λ = 0 and κ = 1, both distributions are very close to the prediction of Poisson statistics typical of an integrable or localized systems. By contrast, see Fig. 31 and Fig. 32, even for small values of λ = 0.015, the agreement with the random matrix prediction is excellent for all κ and N considered. Unlike other SYK models, the universality class is always that of systems with time reversal invariance corresponding to the universality class AI † [101], related to transposition symmetry and not to the Ginibre Orthogonal Ensemble [107] related to complex conjugation symmetry. Although, strictly speaking, we do not have the equivalent of a Bohigas-Giannoni-Schmit conjecture [100] for non-Hermitian systems, we believe that this agreement with random matrix theory still provides evidence of quantum chaotic motion triggered by a weak explicit coupling λ. For larger values of λ ≥ 0.04, we start to observe growing deviations from the random matrix results, which is likely due to the fact that a growing number of eigenvalues become strictly real and the chosen spectral observables are intended for the analysis of complex eigenvalues. For instance, for λ = 0.1 in Fig. 10, about 85% percent of the spectrum is real. In any case, this intermediate region is not universal and therefore is of less interest.
We now investigate whether anomalies, even for small λ, are observed in the infrared part of the spectrum corresponding with the eigenvalues with the largest negative real part. This is precisely the region related to the Euclidean wormhole phase where the spectrum of the ensemble average system has a gap. An immediate problem is that it is not yet clear how to order the complex spectrum and therefore to define precisely the part of the spectrum related to the Euclidean wormhole. Strictly speaking, the Euclidean wormhole is associated with the eigenvalue with the largest negative real part. This eigenvalue E 1 is always fully real for all range of parameters considered. We thus computed the variance R 1 of the probability distribution of |E 1 |, and normalize it by ensemble averaging. The results show its excellent agreement with the RMT prediction. For instance, for λ = 0.015 and κ = 1, R 1 equals to 1.00332 while the random matrix prediction is 1.00315. This agreement extends to other values of λ where the spectrum is still complex. What's more, the agreement with the random matrix prediction goes beyond R 1 . In Fig. 33, we compare the full distribution of |E 1 | with the the random matrix prediction, namely, the Tracy-Widom distribution [108] for systems with time reversal symmetry β = 1. After the preceptive rescaling and shifting, we obtain a good agreement with the Tracy-Widom distribution for the distribution of the eigenvalue with the largest real negative part. However, substantial differences are observed even for the distribution of the eigenvalue with the third largest real negative part. We note that the distribution of the eigenvalue with the largest real negative part corresponding to a random matrix belonging to the AI † universality class is qualitatively different. As was expected, the agreement is worse if we fit it to a Gaussian distribution. To some extent, the complex level statistics results are expected as there is no a visible gap in the spectrum for λ = 0 before any ensemble average. This is another indication that the wormhole phase, which is still characterized by a gap, requires dominance of off-diagonal replica configurations [77] and therefore ensemble average.
We now turn to the analysis of spectral correlation for λ > λ c . The spectrum becomes real even if the Hamiltonian is non-Hermitian when κ > 0. We again employ the adjacent gap ratio which for a real spectrum [104][105][106][109][110][111][112] is given by where δ i = E i − E i−1 and the spectrum is assumed to be ordered.
For random matrices, and for uncorrelated eigenvalues, it is possible [105] to find explicit analytic expressions for both its average and the full distribution function. For instance, for a quan- tum chaotic system with no translational symmetry , the averaged gap ratio r ≈ 0.530 while r P ≈ 0.386 for Poisson distribution corresponding to uncorrelated eigenvalues. We shall also study the level spacing distribution P (s), namely, the probability to find two consecutive eigenvalues E i , E i+1 at a distance s = (E i+1 −E i )/∆, where ∆ is the mean level spacing in that region of the spectrum. For a fully quantum chaotic system, P (s) is given by the random matrix theory results which depends on the global symmetries of the system. In the case of time reversal invariance, it is well approximated by the so-called Wigner surmise P W,GOE (s) ≈ π 2 s exp(−πs 2 /4) for the Gaussian Orthogonal Ensemble (GOE), while for uncorrelated eigenvalues, corresponding to integrable non-degenerate or Anderson localized systems, it is given by Poisson statistics (P P (s) = e −s ). Technically, the calculation of P (s) requires unfolding the spectrum so the average local level spacing is one. We carried it out by employing a low order, in most cases six order, polynomial to fit the average spectral density. The level spacing distribution, which is complementary to the adjacent gap ratio, provides information, specially its tail, about the dynamic of the system for time scales of the order of the Heisenberg time. The adjacent gap ratio, in the other hand, probes the dynamics to even longer scales.
As it can be seen in Fig. 9, the approximate spin-symmetry has an important effect: the spectrum is concentrated around the eigenvalues ofŜ. The results are indeed very similar to that of the κ = 0  Figure 35. Level spacing distribution P (s) for λ > λ c where the spectrum is real. We find excellent agreement with GOE level statistics which is obtained by exact diagonalization of random matrices of size 1000 × 1000, even at the tail of the spectrum. Unlike the non-hermitian case, almost no size dependence is observed.   Figure 36. Distribution of the gap ratios P (r) for λ > λ c where the spectrum is real. We observe an excellent agreement with the GOE prediction even in the tail of the P (r), see insets in log scale.
case and therefore quite insensitive to κ. The same conclusion largely applies to level statistics. In such that the spectrum is always real. We find agreement for most parts of the spectrum with the random matrix prediction for systems with time reversal symmetry r RM T ≈ 0.529 [105]. The observed deviations occur for gap ratios corresponding to eigenvalues of the Hamiltonian located at the edges of sectors belonging to different eigenvalues ofŜ. The gap ratio r is very small since max{δ i , δ i+1 } is very large in this case. Therefore, these deviations with respect to the RMT prediction do not have a dynamical significance. We note that the infrared part of the spectrum, related to the wormhole phase, also fits well with the RMT prediction. This is in contrast with the κ = 0 case [67] where the wormhole phase is characterized by strong deviations from the RMT result.
Similarly, in Figs. 35 and 36, both the P (s) and the distribution function of the gap ratio P (r) of our model, agree well with the random matrix prediction. We note the agreement extends to even the tail of P (s) which probes the dynamics at time scales of not only the order, but longer, than the Heisenberg time. For sufficiently large λ, the eigenvalues of the Hamiltonian cluster around the eigenvalues of the spin-like operatorŜ. Thus it is required that each of these clusters is considered separately for level statistics analysis.

V. CONCLUSION AND OUTLOOK
We have investigated a two-site non-Hermitian SYK model with a weak inter-site coupling and its dual in JT gravity. In the SYK model we have employed exact diagonalization techniques and the numerical solution of the Schwinger-Dyson equations describing the large N saddle-points of the action. On the gravity side, we have derived an effective Schwarzian action and used it to compute the gravity path integral in the saddle-point approximation.
In both SYK and JT gravity, we have studied the thermodynamic properties and observed a thermal phase transition between the wormhole and two black holes. We have obtained an excellent match for the thermodynamic observables such as the free energy, the energy gap and the critical temperature, see Fig. 23 and Fig. 24. One of the motivation to introduce imaginary sources was to be able to study the gravitational path integral using saddle-points. We view here JT gravity as a low-energy approximation of the exact gravity dual of SYK, in a spirit similar to the higher-dimensional version of AdS/CFT.
The transition observed here could also happen in higher dimensions. The Euclidean version of a Lorentzian geometry can have smaller isometry group because the periodic Euclidean time identification may only preserve a subset of the isometries. It is interesting to note that this requires the geometry to be similar to a wormhole as black hole spacetimes don't have this property. In the Euclidean system, defined by the gravity path integral, we impose a smaller number of gauge constraints than in Lorentzian, which defines a purely Euclidean system without sensible Lorentzian continuation. It would be interesting to see whether a similar Euclidean-to-Lorentzian transition, i.e. the dynamical restoration of the Lorentzian symmetries in a purely Euclidean system, can be observed in higher dimensional AdS/CFT. Higher-dimensional Euclidean wormholes supported by boundary sources were described in [89] and have similar thermodynamical properties as our wormhole. Eternal traversable wormholes in higher dimensions are harder to construct [113,114] but can be obtained using the appropriate setup [90,91].
The Euclidean wormhole studied in this paper can also be used to address the factorization puzzle [48,49], which comes from the fact that the gravity path integral appears to compute an ensemble average, as recently discussed in [50][51][52][53][54][55][56][57][59][60][61]. For η = 0, there is no interaction between the boundary and the Euclidean wormhole is purely a result of the average. It was proposed in [50] that factorization could be restored by half-wormhole saddle-points. In our setup without inter-site coupling (η = 0), half-wormhole saddle-points were constructed in [54] and their behavior was matched with the SYK model at a single realization of the couplings. These half-wormhole solutions should be generalizable to non-zero inter-site coupling η. This would correspond to a system of two coupled half-wormholes which, as suggested by our results, could possibly transition to an eternal traversable wormhole. It would be interesting to make this more precise and obtain a gravity picture of the complex-to-real transition for a single realization of the SYK couplings.
The long time dynamics has been explored by the study of level statistics in the SYK model.
Both, for a real and a complex spectrum, spectral correlations are consistent with quantum chaotic motion as we find good agreement with the random matrix prediction in each case. Importantly, we have also found that the distribution of the energy level with the largest negative real part fit well the Tracy-Widom distribution which is an indication of quantum chaotic behavior in the wormhole phase at large N . In all cases, the universality class was that of systems with time reversal invariance. It would be interesting to investigate whether some aspects of the level statistics can be understood using the gravity path integral, which would require finer gravity observables than studied in this paper. We provide additional details about the method we have employed to compare the gap E g in JT gravity (48) with that in the SYK model.
The JT and SYK parameters should be proportional k = Aκ, η = Bλ where A, B are some constants. In order to determine these constants, we compare the numerical results for E g in SYK with the analytical expression E g in JT gravity (48) using A, B as fitting parameter. More specifically, we note that X = C 1 + C 2 , related to E g by (47), is the only real-valued solution of (39). For the sake of convenience, we define, we can then rewrite (39) as z(κ, λ) = A 2 x(κ, λ) + By(κ, λ). The functions x(κ, λ), y(κ, λ), z(κ, λ) can be calculated from the above expressions where E g is taken to be the numerical SYK result.
We see that A, B are just the coefficients of the linear equation z(x, y) = A 2 x + y which are easily calculated by numerical fitting. We find that for κ = 0, 0.1, · · · , 0.7, λ = 0, 0.01, · · · , 0.09 the best fit corresponds to A = 1.58497 and B = 27.05755. In Fig. 23 of the main text, we compare explicitly the numerical SYK results with those from the fitting above calculation and find that they are fully consistent especially when κ is small. We also compare T c in Fig. 24 using the same values of A and B and obtain an excellent match.

Appendix B: Real time calculation
We consider the Hamiltonian (2). In the large N limit, we obtain the effective action from the use of the replica trick after ensemble average. Here t ab is defined as The real time dynamics is studied after performing a Wick rotation −iω → ω + i to the Schwinger-Dyson equations for the imaginary time Green's function. We closely follow the method of Ref. [69,115].
With these definitions, we have and the relations G LL = G RR , G LR = −G RL . The retarded Green's function, after the Wick rotation, takes the form One of the main technical difficulties is the calculation of Σ r + . Following Refs. [66,115], we first calculate Σ LL (ω n ) andΣ LR (ω n ), which are simply the self-energy Σ LL (τ ) = J 2 G 3 LL (τ ) and Σ LR (τ ) = J 2 G 3 LR (τ ) but in the frequency domain. We then apply the Wick rotation to obtain, Σ + (ω) = Σ LL (ω) + iΣ LR (ω), with More details of the calculation can be found in Appendix F of [115] and Appendix D of [69].
Here we focus on the technical details required to solve numerically these real time saddle-point equations. The main difficulty lies in the choice of a proper cutoff L in t, and also a small but finite and the number of discrete points N in the sums, for given parameters T , λ, κ, J.
It can be seen from G(ω) = ρ(ω )dω ω−ω +i that the first relation that must hold is dω = 2π L . To make sure that (ω−ω )+ 2 ≈ πδ(ω − ω ), should be small enough so that we obtain a sharp peak, mimicking a delta function. We also need dω so that the profile is still smooth. For a given set of parameters, it is important to find the right balance between the necessary suppression of discretization artifacts that may obscure real physical effects and the optimization of computational resources in terms for instance to computation times are RAM usage.
Moreover, the choices of physical parameters and numerical parameters are not independent. From a rewriting of G r + (ω), we observe that we must impose that λ so that the effect of λ is not obscured by a too large .
Finally, the choice of the parameter N is driven by the required accuracy of the numerical integral over t. Since dt = L N , N should be small enough so that the accuracy of the integral meets a certain minimum. In addition, the range of ω is approximately ∼ 2π L N . Therefore a small enough L/N will also guarantee the irrelevance of the truncation for large ω. Taking all this into consideration we set L = 5 × 10 6 and N = 2 25 ≈ 3 × 10 7 .