Quantum chaos transition in a two-site SYK model dual to an eternal traversable wormhole

It has been recently proposed by Maldacena and Qi that an eternal traversable wormhole in a two dimensional Anti de Sitter space (${\rm AdS}_2$) is the gravity dual of the low temperature limit of two Sachdev-Ye-Kitaev (SYK) models coupled by a relevant interaction (which we will refer to as spin operator). In this paper, we study spectral and eigenstate properties of this coupled SYK model. We have found that level statistics in the tail of the spectrum, and for a sufficiently weak coupling, shows substantial deviations from random matrix theory which suggests that traversable wormholes are not quantum chaotic. By contrast, for sufficiently strong coupling, corresponding to the black hole phase, level statistics are well described by random matrix theory. This transition in level statistics coincides approximately with a previously reported Hawking-Page transition for weak coupling. We have shown explicitly that this thermodynamic transition turns into a sharp crossover as the coupling increases. Likewise, this critical coupling also corresponds with the one at which the overlap between the ground state and the thermofield double state (TFD) is smallest. In the range of sizes we can reach by exact diagonalization, the ground state is well approximated by the TFD state only in the strong coupling limit. This is due to the fact that the ground state is close to the eigenstate of the spin operator corresponding to the lowest eigenvalue which is an exact TFD state at infinite temperature. In this region, the spectral density is separated into blobs centered around the eigenvalues of the spin operator. For weaker couplings, the exponential decay of coefficients in a tensor product basis, typical of the TFD, becomes power law. Finally, we have also found that the total Hamiltonian has an additional discrete symmetry which has not been reported previously.


I. INTRODUCTION
The dynamics of quantum many-body systems is very rich and, in general, highly dependent of the details of both initial conditions and the Hamiltonian that governs its time evolution. However there are situations where the time evolution has universal features. A paradigmatic example is that of a quantum system whose classical counterpart is chaotic for time scales of the order of the inverse of the mean level spacing, the so called Heisenberg time. In that region, where the discreteness of the spectrum is relevant, and therefore quantum effects are always strong, quantum chaotic systems, assuming localization effects are negligible, relax to a fully ergodic state that only depends on the global symmetries of the system. Level statistics are a powerful tool to detect and classify universal features in this region, because it is expected that the predictions of random matrix theory will apply in this case. Indeed, agreement with random matrix theory results has been found in a broad variety of systems: highly excited states of nuclei [1], weakly disordered systems with [2] and without interactions [3], deterministic quantum chaotic systems [4] or the low energy limit of QCD in a box [5,6].
The field of quantum chaos received an important boost after it was claimed [7] that quantum black holes in the semiclassical limit are quantum chaotic. More specifically, it was proposed [7] that the growth of certain out-of-time-ordered correlation functions in the semiclassical limit, and for intermediate times of the order of the Ehrenfest time, is exponential as for quantum chaotic systems [8]. Moreover, the growth rate, given by the Lyapunov exponent, obeys a universal bound which is saturated for field theories with a gravity dual. Later this feature was also observed [9,10] in the strong coupling limit of the Sachdev-Ye-Kitaev (SYK) model [9,[11][12][13][14][15][16][17][18][19][20], a zero dimensional fermionic model with infinite-range random interactions. This simple quantum mechanical model has attracted a great deal of attention because, despite being quantum chaotic and strongly interacting, it can be tackled analytically [9,10,[21][22][23][24][25][26][27][28][29] and has all the features expected for a system with a gravity dual [9,10,27,30]: it does saturate the bound on chaos mentioned above, it also has a finite zero temperature entropy [18], a linear specific heat in the low temperature limit [9,10,21,31] and a density of low energy excitations [10,21,31] that increases exponentially. Level statistics of the SYK model was shown to be well described by random matrix theory [21,29,[31][32][33] in the spectral region where the SYK model is expected to have a gravity dual which confirms that it is also quantum chaotic for long times. Based on the same pattern of conformal symmetry breaking, [34,35], the gravity dual of the SYK model has been identified to be a Jackiw-Teitelboim [36][37][38] background.
A natural question to ask, assuming that the holographic duality applies for sufficiently long times, is whether these quantum chaotic features are exclusive to quantum black holes geometries, such as the Jackiw-Teitelboim background, or also apply to other backgrounds. Recently, Maldacena and Qi [39], see also [40,41], constructed a near AdS 2 background representing an eternal traversable wormhole. The quantum feasibility of a traversable wormhole was previously suggested [42] by showing that a double trace deformation could make the quantum matter stress tensor to have a negative average null energy without violating causality. An example of such construction for AdS 2 geometries was worked out in [43] and for higher dimensions, and other backgrounds, in [44][45][46][47] Out of equilibrium features for intermediate times were investigated in Ref. [48] and entanglement and other quantum information observables were studied in Ref. [49].
Returning to eternal traversable wormholes, it was conjectured [39] that its ground state is well approximated by a thermofield double state (TFD) and that its field theory dual is a two-site SYK model coupled by a relevant interaction in the limit of weak coupling and low temperature. It was found that for sufficiently high temperatures, or stronger coupling, the system undergoes a thermodynamic Hawking-Page transition from the wormhole phase to a two black hole phase. By increasing the strength of the coupling, the transition takes place at higher temperatures and finally, after reaching the critical coupling, this first order transition terminates. Intriguingly, the wormhole phase is still controlled by a generalized Schwarzian action when the residual SL(2) symmetry of the standard SYK model is broken [39]. The associated Liouville quantum mechanics problem [25] has graviton bound states which are interpreted as excited states of the quantum wormhole background. A non-random SYK model with the same large N limit was investigated previously in Ref. [50] with similar conclusions.
In this paper, we study different aspects of this coupled SYK model. Our motivation is twofold. First, we aim to investigate the long time dynamics of this two-site coupled SYK model by level statistics in order to clarify whether quantum chaos and random matrix theory are generic in quantum gravity backgrounds or are restricted to black holes. Second, we would like to gain a more detailed understanding of the phase diagram of the system. More specifically, we aim to determine the range of validity of the TFD as a ground state of the system and also the exact nature of the termination of the first order Hawking-Page transition at finite coupling. We have found that indeed these two main goals are closely related. We show that the first order transition suddenly becomes a sharp crossover above the critical coupling. Around the same coupling strength, the number of non-zero coefficients, in a tensor product basis of the two independent SYK's, that contribute substantially to the ground state of the system increases sharply. This is typical of chaotic-integrable transitions. More specifically, except for strong coupling and in the range of sizes we can explore numerically, we have observed that the energy dependence of these coefficients is power law rather than the exponential decay typical of the TFD state.
Level statistics are also sensitive to the value of the coupling constant. In the tail of the spectrum, and couplings weaker than the critical one, corresponding to the wormhole phase, level statistics are not described by random matrix theory. However, well above the critical coupling, spectral correlations in the full spectrum are well described by random matrix theory.
The organization of the paper is as follows. First, we introduce the model with a special emphasis on the description of its symmetries and its ground state. In particular, we identify a global spin symmetry in the combined system that has an important effect on spectral statistics and thermodynamic properties. In Section III, we investigate the ground state of the system. For the number of Majoranas we can simulate by exact diagonalization, the ground state is well approximated by TFD only in the limit of strong coupling between the two SYK's. Thermodynamic properties are investigated in Section IV. We have found that the first order transition at a finite coupling, reported in Ref. [39], is followed by a sharp crossover for stronger coupling. At the critical coupling the overlap between the ground state and TFD is smallest. In the strong coupling limit, the spectral density develops well separated blobs centered around the spin eigenvalues that control the free energy. Section V is devoted to the analysis of level statistics. For weak coupling, and in the infrared region of the spectrum, corresponding to the traversable wormhole phase, we did not find agreement with the random matrix theory, thus suggesting that the gravitational bound states in this phase are not quantum chaotic. However, for sufficiently large coupling, corresponding to the black hole phase, level statistics are well described by random matrix theory. Finally, in Section VI, devoted to outlook and conclusions, we speculate that the termination of the observed transition at strong coupling is reminiscent of a Gross-Witten transition [51], induced by the gradual reduction of the effect of interactions by increasing the coupling between the two SYK models.

II. THE MODEL AND ITS SYMMETRIES
We study N Majorana fermions in (0 + 1) dimensions where the first N/2 fermions, ψ L , are labeled by Left (L), and the remaining N/2, labeled by R, ψ R , will be called the Right (R) fermions. In each of these two subspaces, the dynamics is governed by the SYK model, where ψ L,i , ψ R,i are Majorana fermions {ψ A,i , ψ B, j } = δ AB δ i j (A, B = L, R) and J i jkl are Gaussian distributed random variables with J i jkl = 0 and standard deviation J 2 i jkl = 4 The two SYK model are coupled by the following interaction: so that the total Hamiltonian is simply given by, For the moment we set α = 1 in which case the Hamiltonian has an extra global spin symmetry that we will discuss in detail shortly. For the study of level statistics we will find it computationally advantageous to break this symmetry by setting α 1 to a slightly different value which does not affect qualitatively most of the key features of the model [39], like for example the observed gap between the ground state and the first excited state.
The anti-commutation relations of ψ L,i , ψ R,i can be realized by introducing N Gamma matrices Γ i as Since each term in the full Hamiltonian is an even power of the Gamma matrices, the full Hamiltonian preserves chirality, defined as the eigenvalue of Γ c = Γ 1 Γ 2 · · · Γ N . We have summarized the notation for the Gamma matrices and related symmetries in appendix A.

A. S mod 4 symmetry
Let us define the spin operator S as that is, H int = kS . Below we show that S has the discrete spectrum S = −N/4, −N/4 + 1, · · · , N/4 and that, if α = 1, the full Hamiltonian (3) preserves S mod 4 in addition to the chirality (S mod 2). We shall call S the spin operator.
First, let us deduce the spectrum of S . If we define Γ ± i as then S can be written as The operators Γ ± i also satisfy the following (anti)commutation relations from which we can derive the following spectrum: where |S = −N/4 is the lowest spin state satisfying Next we discuss the S mod 4 symmetry. If we define the full Hamiltonian (3) with α = 1 can be written as The P i j and Q i j operators satisfy the following commutation relations with S That is, P i j preserves the spin S . To interpret the commutation relation for Q i j it is convenient to act Eq. (13) on an eigenstate of S , with the spin eigenvalue s, which we shall call |φ . One easily derives This implies that an expansion of Q i j |φ contains only spin eigenstates with S = s ± 2. Hence S mod 4 is a symmetry of Q i j Q k . Putting the above computations together we find that the full Hamiltonian Eq. (3) preserves also the S mod 4 symmetry.
It is instructive to express the exponent of the spin operator as a product For θ = π we recover the chirality operator, while θ = π/2 gives the S mod 4 symmetry, This also shows that the chirality symmetry operator is the square of the S mod 4 operator. In appendix B,

III. COMPARISON OF THE TFD STATE WITH GROUND STATE
Here we consider the case with identical L and R systems, i.e. α = 1. For k = 0 (and N = 0 mod 4 so that N/2 is even), the system we are considering is H = H S Y K ⊗ 1 + 1 ⊗ H S Y K . That is, the system consists of two decoupled SYK systems with N/2 fermions, the L and R systems, which are identical to each other.
It was already pointed out in [39] that when the interaction between L-system and R-system is turned on the ground state has a large overlap with the Thermo Field Double (TFD) state of H S Y K ⊗ 1 + 1 ⊗ H S Y K . In this section, we shall study this similarity in more detail.
Because the left and right systems are identical, the left and right Γ-matrices can be represented as where U is a unitary matrix (we use the notation that the lower case gamma matrices are of size 2 N/4 × 2 N/4 and upper case gamma matrices are of size 2 N/4 × 2 N/4 ).
The factor γ c ensures that the left Γ matrices anti-commute with the right Γ matrices. Below, we will call this basis the tensor basis. If the eigenvalues and eigenstates of H L are given by then Let us define the TFD state, where we have included a phase factor, exp iφ n and CK is the charge conjugation operator. We will next show that is the ground state of the spin operator. To this end we calculate the expectation value where χ m is the chirality of the state |m , Choosing |m R = |m L (i.e. U = 1 in Eq. (17)) the sum over m can be eliminated by completeness.
which is the desired result. Therefore |I is the state with the lowest spin. A similar prove can be used to show that Γ − i |S = −N/4 2 = 0 (no sum over i). This shows that the state |I is annihilated by all lowering operators Γ − i . The details of the ground state which depend on N/2 mod 8 are worked out in Appendix D for N/2 mod 8 = 4, in Appendix C for N/2 mod 8 = 2 and in Appendix E for N/2 mod 8 = 0.
At finite temperature, we define the TFD state |TFD β as where Z(N, β) is a normalization factor such that TFD| β TFD β = 1. As confirmed by numerical calculations for N up to 32, the complex phase factor of the |n ⊗ |n components of the ground state at finite coupling k, does not depend on k. This is a nontrivial result which we only understand perturbatively.
A. |TFD β versus ground state of coupled SYK In [39] it has been argued that the ground state |gs of the coupled system at finite coupling k is close to the TFD state |TFD β by studying the overlap TFD| β gs . Here the inverse temperature β is a function of k determined by maximizing the overlap. However, it should be stressed that a large value of this overlap does not necessarily imply that the coefficients ( n| ⊗ n|)|gs , taken as a function of the energy levels of the single SYK system, behave like e −2βE n . The reason is that for small coupling most of the strength of the wave function can be localized in one or few components of the wave function. In this section we will show numerically that, at finite N ≤ 28, this is actually the case and that for sufficiently small values of the coupling k the coefficients display large deviations from the exponential behavior typical of the TFD state.
Before moving to a detailed analysis, let us first consider the results in some limiting cases. In the limit k → ∞, since H is dominated by H int = kS , the ground state |gs coincides with |TFD β=0 , see Eqs. (22), (24) and (25). In the limit k → 0, on the other hand, the tensor product states |m ± L ⊗ |n ± R are the exact eigenstates of H with the eigenvalues E m + E n , where the E n are the eigenvalues of a single SYK, H S Y K |n χ L,R = E n |n χ L,R , and the index χ = ± denotes the chirality of the state. The ground state is |1 χ L ⊗ |1 χ R which may be degenerate depending on the value of N. Hence the ground state in the k → 0 limit is given by the |TFD state for β = ∞. Note that this does not necessarily imply that the ground state is close to a TFD state for small but nonzero k.
Numerically, the phases of the c n do not depend on the coupling constant and are thus given by the phases of the ground state of S (perturbed by a (H L + H R ) with a small constant). This is not the case for the phases of the coefficients b mn . We can understand this fact to first order in perturbation theory in 1/k with S as the zeroth order Hamiltonian, with |S m an eigenstate of S with eigenvalue S m and S 0 the spin of the ground state. Inserting a complete set of eigenstates of H L + H R , m,n |m ⊗ |n m| ⊗ n|, and noticing that m| ⊗ n| (H L + H R )|I = 0 for m n, the matrix element can be written as Only those states with the same S mod 4 spin as the ground state contribute to the sum over n. It turns out that the complex phase of the matrix elements S m | |n ⊗ |n does not depend on m apart from an overall phase which cancels with the contribution to the matrix element which is therefore real. The complex phase of the first order correction is therefore due to the phase of n| ⊗ n| |S m which is the same as the phase of n| ⊗ n| |I (up to a minus sign). The same argument cannot be made for third and higher order perturbative corrections because off-diagonal states, |m ⊗ |n In the numerical studies below we will only consider the cases N = 20 and N = 28, where |TFD β is defined as, see Eq. (25), and Z(β) = 2 2 N/4−1 n=1 e −4βE n is the normalization factor. Since k is the parameter of the Hamiltonian that interpolates between a ground state |TFD β=∞ and |TFD β=0 which also corresponds to the ground state of the spin operator, we expect that the effective inverse temperature β(k) monotonically decreases as k increases. Furthermore, we have computed | S = −N/4|gs | 2 and found that this overlap is almost 1 already around k ∼ 1 (see Fig. 1). Hence we expect the effective temperature changes substantially only in the regime 0 < k < 1.
Below, to be self-contained, we first redo the comparison of the ground state and the TFD state by maximization of their overlap. This computation was already performed in [39]. We have checked numerically, for N = 20 and N = 28, that the maximized overlap with respect to β, between the TFD state and the ground state, is close to 1. For this purpose, we have considered the ansatz Eq. (30) for the ground state, in which the E n are the energy levels of a single SYK model and β is the fitting parameter determined by maximizing the overlap. We normalize this state as TFD| β TFD β = 1 by choosing Z(β) as 1. Generate an ensemble of 2 N/4 ×2 N/4 matrices H S Y K = 1 2 i< j<k< J i jk γ i γ j γ k γ with Gaussian random J i jk , and compute their eigenvalues {E n } and eigenvectors {|n ± } through exact numerical diagonalization. Then construct the tensor product |n ± L ⊗ |n ∓ R .
2. For the same ensemble of random couplings, generate 2 N/2 × 2 N/2 matrices representing the Hamiltonian Eq. (3) with α = 1 using the tensor basis (17) and compute the ground state |gs through the exact diagonalization.
3. Compute c ± n = ( n ± | L ⊗ n ∓ | R )|gs . We have observed that c + n = c − n and c + n > 0 for all n = 1, 2, · · · 2 N/4 /2, up to an overall phase provided that we define the phase of |n ± L,R such that Eq. (30) holds without any extra phase (see appendix D for detail). Hence we shall denote c + n = c − n simply by c n . 4. For each ensemble realization determine β by maximizing the overlap TFD| β gs .
5. Take average of both β and TFD| β gs over the ensemble.
Following this procedure we obtained both the effective inverse temperature β(k) and the overlap TFD| β gs .
The results are displayed in Fig. 2. As we can see, the inverse temperature is indeed growing for small k and it tends to zero at strong coupling as expected. In agreement with [39], we confirm that TFD| β gs is pretty close to 1 with relatively large deviations around k ∼ 0.1. More quantitatively, for small k, the inverse effective temperature depends on k roughly as 1/k 1/6 while at k ≈ 0.1, there is a transition to a 1/k dependence.
As was already mentioned, the large overlap between the ground state and the TFD, does not necessarily imply that the coefficients c n have an exponential behavior in energy. Indeed, it could be simply that one or a few coefficients in (26) are much larger than the others. In this case, the overlap TFD| β gs would be very close to be 1, even though all the other coefficients could behave in a completely different way.
Hence, to check whether the coefficients of the ground state have an exponential profile, we have studied their average on a logarithmic scale versus the ensemble average of each of the E n (for the case N = 28 only).
Results, depicted in Fig. 3, clearly show that an exponential form is not a good fit to the |c n | coefficients.
For comparison, we have also included the coefficients of the TFD state with β determined by maximizing the overlap (green curves) as described before, and fitting the logarithm of the ensemble average of |c n | to the TFD form also leaving the normalization constant as a free parameter (red curve). None of these fitting to the TFD state give good results. By contrast, a power law ansatz (black line),      provides an excellent fitting. In order to compare the two parameter fits, we fix p at p = 2.178. The fitted value of e 0 is somewhat smaller than the ground state energy E 1 for small k and decreasing for increasing values of k. The value of c is in principle determined by the normalization, but leaving it as a free parameter gives a better fit. For small coupling k < 0.5 the value of c is quite different from the one obtained from the normalization. The reason is that most of the strength of the wave function is in c 1 which does not contribute much to the χ 2 when you fit the logarithm of the |c n |. One could also take p as a k-dependent fitting parameter which gives slightly better fits. In any case, except in the regime of strong coupling, power law fits provide a much better description of the numerical results than exponential fits related to the TFD state. For k ≥ 1 the TFD ansatz is a good description of the ground state though. This is consistent with the expectation that at large coupling the system is dominated by the spin operator and the TFD state (with ponentials with a χ 2 that is comparable to the power law fit, but it has one more parameter. A Gaussian dependence, exp −βx − αx 2 , first introduced in the context of a nuclear many body system [52], also gives a reasonable fit in this parameter region.
Once again, we stress that the deviation from the TFD form is not in contradiction with the good overlap between the TFD and the ground state, since the dominant coefficients are well reproduced by the TFD at least for small k, as it can be seen from Fig. 3.
These numerical results suggest that the coefficients of the ground state are not well approximated by the TFD ansatz, and that instead they are very close to a power law ansatz. However, we should note that N = 28 may not be close to the large N limit of this system where analytical arguments [39], show that the ground state is well approximated by a TFD state. This indicates that the convergence to the large N limit is slow and most likely nonuniform. Interestingly, it has been observed in [53], that the Hamiltonian Eq. (3) is just an approximation of a Hamiltonian whose exact ground state is the TFD.
In Fig. 4 we show the ratio |c 2 /c 1 | (left) and the inverse participation ratio (IPR) (right) defined as We observe a crossover at k ≈ 0.1 where β(k) ≈ E 2 − E 1 . Around the same value of k, the inverse participation ratio increases dramatically. This is a feature typical in metal-insulator and integrable-chaotic transitions.
The results of this section suggest that the ground state of the model undergoes a qualitative change around k ∼ 0.1 − 0.2. We shall see in next section that precisely in this region the first order transition mentioned in the introduction turns into a sharp crossover. We turn now to the study of low energy excitations by the analysis of thermodynamic properties and spectral correlations. Ref. [39], in more detail.

A. Numerical Study
In order to investigate the thermodynamic properties of the Hamiltonian Eq. (3), we compute its spectrum by exact diagonalization techniques for up to N = 34 Majoranas. We are mostly interested in the low temperature limit where, according to the results of Ref. [39], a first order Hawking-Page transition occurs at a finite value of the coupling k. The transition, according to Ref. [39], seems to terminate for k ≥ k c though it is not clear whether it becomes second order or just a crossover.
We study the free energy as a function of k and N. For any N 1 we have found, see For larger k (right) the peak is broad which suggests a crossover. For small k (left), finite size effects seems to be stronger and we cannot reach any firm conclusion.
of the critical temperature can be obtained from the specific heat. For a first order phase transition, the specific heat diverges at the transition, for a second order transition it shows a finite jump and for higher order phase transitions, the specific heat is smooth.
In Fig. 6, we show results for the size dependence of the specific heat for different values of k. We

B. Solution of Schwinger-Dyson Equations
The above analysis, though illustrative, is not conclusive because of the limited size we can study by exact diagonalization. In order to gain a more quantitative understanding, we have computed the free energy and specific heat in the large N limit by solving the Schwinger-Dyson equations numerically. The details of this calculation have been already explained in great detail in [39], to which we refer. The purpose of this section is to perform a more precise analysis around the critical coupling k c . free energy on the temperature step: the first order phase transition is clearly visible by noticing that, up to k = k c , the free energy shows a hysteresis curve as function of the temperature, see Fig. 7.
The critical value of k c ∼ 0.175 is consistent with the one, obtained in the previous section, at which the ground state changes qualitatively. We have also observed, see Fig. 8, that for slightly larger values of k ≥ 0.177, this first order phase transition turns out to be a sharp but smooth crossover, which becomes broader with increasing k. To reach this conclusion, it has been necessary to study how the peak in the specific heat behaves when decreasing the temperature step: as shown in Fig. 8, by further decreasing the temperature step to the very small value of dT = 5 · 10 −7 , the peak turns out to be a sharp, but smooth, crossover. We therefore conclude that the window for a hypothetical second order phase transition must be very narrow, 0.175 ≤ k ≤ 0.177. We note that our results are in perfect agreement with those of [50], where a detailed analysis of the Schwinger-Dyson equations associated with a very similar model without disorder was carried out.
We have also computed, see Fig. 9, the energy gap E g between the ground state and the first excited state as function k for a fixed low temperature by a careful fitting of the exponential decay of the relevant Green's function. The same computation has been already performed and explained in [39]. Here, we just reproduce the result for completeness. The dependence of the gap on the coupling k is a further confirmation of the existence of a transition. Around k ∼ k c we observe a change of behavior from E g ∝ k 2/3 for k k c to E g ∝ k for k k c . The latter linear behavior is expected if the gap is due to direct hopping induced by the coupling term in the Hamiltonian which is one-body (two Majoranas) and therefore non-interacting.
Physically, by increasing the coupling k, we are effectively reducing the many-body interaction inside each SYK model in favor of a direct coupling between the two SYK models which is a diagonalizable one-body (two Majoranas) interaction. This is vaguely reminiscent of a Gross-Witten transition [54], or a BEC-BCS crossover [55], where the reduction of the interaction strength induces a higher order transition or just a crossover.
For larger k > 1, we have found, see Fig. 10, that the spectral density starts to split into different regions.
There is a simple explanation of this phenomenon. The interaction term in the Hamiltonian, which becomes dominant in the large k limit, is proportional to the spin operator which has a discrete spectrum with large degeneracies (see (9)). Therefore we expect that, in this region, the spectral properties are largely controlled by this coupling term and the interactions inside each of the SYK models spread out the degenerate states.
Since the eigenvalues of the coupling term range from −kN/4 to kN/4, we expect the full spectrum of the model to cluster around these few eigenvalues which leads to the appearance of the above mentioned blobs in the spectral density. This is exactly what we observe in the large k limit (see Fig. 10).
The free energy can be computed analytically in the limit k → ∞ as follows. If we ignore the SYK terms, we are left with the discrete spectrum just given by Eq. (9), for which the free energy is given by In Fig. 11, we compare the numerical free energy for k = 1, 2 and 5 with this analytical approximation and find that the agreement becomes excellent around k ∼ 5.
All these results indicate that the physical reason for the termination of the Hawking-Page transition is just the gradual reduction of the interaction strength. This is a strong indication that a gravity interpretation is restricted to low temperatures and weak coupling where the first order transition takes place.

V. LEVEL STATISTICS
We now turn to the study of level statistics, with the main aim to characterize the dynamics in the traversable wormhole phase observed in the low temperature limit of the free energy. For that purpose, we obtain the exact spectrum of the model by exact diagonalization techniques for up to N = 34 Majoranas.
We are especially interested in the spectral correlations of the smallest eigenvalues above the ground state, which are likely to be closely related to gravitational modes of the wormhole phase.
For a meaningful analysis of spectral correlations it is in general necessary to unfold the spectrum so that the mean level spacing is the same across the spectrum. For that purpose, we employ the splines method that fits locally consecutive subsets of many (> 10) eigenvalues with low order polynomials. Results are insensitive to the degree of the fitting polynomials.
We aim to study the evolution for long times, of the order of the Heisenberg time, in order to find the range of applicability of random matrix results, a feature of quantum chaotic systems, rather than finite size deviations from it. Hence, we will focus on short range spectral correlators, such as the level spacing distribution, P(s), and the adjacent gap ratio. The former is defined as the probability to find two consecutive average local level spacing). For a fully quantum chaotic system it is given by Wigner-Dyson statistics [56] which is well approximated by the so called Wigner surmise which depends on the universality class [57]. For the Gaussian Orthogonal Ensemble (GOE), corresponding to systems with time reversal invariance, it is given by: P W,GOE (s) ≈ π 2 s exp −πs 2 /4 . For an insulator, or a generic integrable system, it is given by Poisson statistics, P P (s) = e −s . The adjacent gap ratio is defined as [2,58,59], for the ordered spectrum For a Poisson distribution it is equal to r P ≈ 0.386 while for a random matrix ensemble it depends on the symmetry class, with r ≈ 0.530 for the GOE [60]. The advantage of r over P(s), is that it does not require to unfold the spectrum. The adjacent gap ratio was also studied in [61] for another model related to the traversable wormhole [62].
We note that for a correct analysis of spectral correlations it is necessary to consider only eigenvalues with the same (good) quantum number. This means that we have to consider eigenvalues only of a given spin and parity sector. Due to the block structure of the Hamiltonian, it is straightforward to select eigenvalues of the same parity. For the spin symmetry, we could not find a numerically inexpensive method to go beyond comparatively small sizes N = 28. We avoid this problem by introducing an asymmetry in the couplings (α = 1.15 in Eq. (3)) between left and right SYK models. It has been argued [39] that the main features of the model, including the existence of a gapped phase, though with a smaller gap, for low temperature which is a signature of the wormhole phase, are robust to a small interaction strength asymmetry which breaks the spin symmetry. We start our analysis of level correlations with the case of very small coupling. in the bulk are well described by Poisson statistics. This is the expected behavior for a system which is defined as the tensor product of two many-body quantum chaotic systems whose spectral correlations are described by random matrix theory. For very small k 2 −N/2 , the spectral degeneracy is barely lifted. As a consequence, the distance between eigenvalues that were degenerate for k = 0 is much smaller than that between neighboring eigenvalues which makes the analysis of level statistics difficult. For k 2 −N/4 , the degeneracy is already lifted, so the statistical analysis does not require any artificial pruning of the spectrum.
In this region, the level statistics in the tail and in the bulk of the spectrum are qualitatively different.
In Fig. 13, we plot the nearest neigbor spacing distribution in the bulk of the spectrum for small k and different values of N. Level repulsion is observed in all cases, but for N = 26 and N = 28 the decay for large distances is clearly exponential with an exponent close to the one corresponding to a Poisson distribution.
However for N = 32 and N = 34 we find good agreement with random matrix correlations. Although a careful finite size scaling analysis would be necessary to reach a final conclusion, the latter is a strong suggestion that in the bulk of the spectrum, corresponding to high temperatures, even a small coupling k, makes the system quantum chaotic. However, see Fig. 13, it seems that the lowest eigenvalues in the same region of small k, for N = 32 and N = 34 are correlated to Poisson statistics like for k = 0. Rather than a metal-insulator transition at a certain value of k, this behavior suggests that for sufficiently weak coupling Agreement with Poisson statistics is very good. This is a strong indication that the infrared limit of the spectrum for sufficiently small k is similar to that for k = 0 corresponding to two black holes. the two SYK's, each of them quantum chaotic, are effectively disconnected systems for low energies (tail) while for sufficiently high energies (bulk) the behavior is similar to a single quantum chaotic system whose levels are correlated according to random matrix theory. It is tempting to speculate that this indicates a transition from two black holes to one black hole for sufficiently high temperatures. However, the fact that we do not have a clear geometrical understanding on how this transition can occur, together with the limited range of N values which we could study numerically, prevents us to reach any firm conclusion.
Summarizing, our results show that, while in the bulk of the spectrum a very small value of k is sufficient to induce the transition from Poisson statistics to random matrix theory spectral correlations, in the tail of the spectrum, Poisson statistics is robust, and a larger coupling is necessary to induce the transition to level statistics described by random matrix theory. We investigate this case next.

B. Critical k
As k increases, the bulk of the spectrum is still correlated according to random matrix theory with no quantitative change from the region of small k investigated previously. For this reason, we will focus exclusively on the spectral correlations in the tail of the spectrum.
In Fig. 14  block that includes the ground state. In the tail of the spectrum we do observe clear deviations from the GOE prediction that become gradually smaller as k increases. Interestingly, for k ≈ 0.14, the ratio for the lowest eigenvalues becomes suddenly very close to the GOE value while the subsequent eigenvalues still deviate substantially from it. A further increase of the coupling k leads to more of the lowest eigenvalues becoming correlated according to the GOE while the subsequent ones still show deviations. Finally, for k ∼ 1, all eigenvalues in the spectrum are GOE correlated. We stress that the value of k for which the lowest eigenvalue becomes GOE correlated is very similar to the one at which the termination of the Hawking-Page transition, between the traversable wormhole and the black hole phase, takes place. This is a strong suggestion that deviations from random matrix theory in the tail, with values close to Poisson statistics, may be related to the wormhole phase and that the sudden appearance of the GOE correlated eigenvalues deep in the infrared may signal an instability of the wormhole phase. This instability might be expected close to the transition to a black hole phase.
As a further check that random matrix theory describes the level statistics of the lowest eigenvalues for k > 0.14, we compare the level spacing distribution P(s) of these eigenvalues with the Wigner surmise.
Results, depicted in Fig. 14 confirm that for k ∼ 0.5, the level spacing distribution for the lowest eigenvalues is in excellent agreement with the Wigner surmise P W,GOE (s).
The coupling between the left and right site is a relevant perturbation which we expect to become dominant in the tail of the spectrum. Interestingly, in Ref. [39] it was reported that, assuming N fixed and k sufficiently large, some of the lowest energy excitations of the effective Hamiltonian are of gravitational origin and not related to the breaking of conformal symmetry. Physically this corresponds to excited states of the wormhole geometry. It is tempting to speculate that there is a connection between these excitations of the wormhole geometry and the deviations from random matrix theory in the tail of the spectrum below the critical coupling. If this picture applies, these gravitational modes related to the wormhole phase are not quantum chaotic and the sudden transition to random matrix theory in level statistics is a dynamical aspect of the thermodynamic Hawking-Page transition. That would imply that random matrix theory, and therefore quantum ergodicity, is only a signature of quantum black holes but not of an AdS graviton gas. In other words, the Hawking-Page transition can be dynamically characterized as a chaotic-integrable transition.
In order to test this hypothesis, we compare explicitly the number of eigenvalues with random matrix like correlations with the spectral density (one point function) in that region. For that purpose, we employ again the adjacent gap ratio r , but with a larger N = 34, so that more eigenvalues show the anticipated anomalous behavior. Regarding the spectral density, for which the prediction of gravitational excitations of Ref. [39] applies, we employ the local spectral density and d ≡ E i+1 − E i . Results, depicted in Fig. 15, clearly show that the number of eigenvalues with random matrix correlations is in good agreement with the number of eigenvalues whose d shows a pronounced dependence on k for sufficiently large k. Although further research is required to confirm this point, this quantitative agreement is encouraging. It strongly suggests that full quantum ergodicity is typical only of black holes, while other geometries like quantum wormholes may be closer to Poisson statistics typical of integrable dynamics. We now move to study level statistics for α = 1 in order to confirm that these findings are not particular to the value of α ∼ 1 used in that as the asymmetry increases the separation between the low lying eigenvalues gradually disappears. Fig. 16, clearly show similar features as in the slightly asymmetric case: deviation from GOE for small k followed by the sudden appearance of the lowest eigenvalues correlated according to random matrix theory at around the same value k c ∼ 0.15. However, for N = 28, the crossover region is very narrow, only the two lowest eigenvalues become GOE correlated before the full spectrum becomes GOE correlated, and the distinction between the tail and the bulk of the spectrum becomes more difficult to identify.
The spin symmetry is broken by the coupling asymmetry. However, as mentioned previously, a small asymmetry does not change the physics of the model but makes it technically easier to study the level statistics. In order to confirm this point more explicitly, we investigate the dependence of level statistics on the asymmetry parameter α in Eq. (3). In Fig. 17, we depict results for the adjacent gap ratio and different values of α. As α increases, the differences between the tail and bulk of the spectrum becomes increasingly small. Eventually, the full spectrum becomes correlated according to random matrix theory even with a comparatively weak coupling. Physically this is an indication that the wormhole phase related to the graviton gas disappears if the asymmetry between left and right is sufficiently strong. This is consistent with the results of Ref. [39] where it was found that both the ground state energy gap and the left-right correlation decreases as the asymmetry increases. However further research is required to determine whether there is a sharp maximum value of the asymmetry for the wormhole phase to occur.

VI. OUTLOOK AND CONCLUSIONS
We have studied a coupled two-site SYK model whose gravity dual is conjectured to be [39] a traversable wormhole geometry in the low temperature, small coupling, limit and dual to a black hole geometry for sufficiently strong coupling or high temperature. We have analyzed the spectral and thermodynamic properties of this model and have found that it has a previously unknown discrete spin symmetry which among other things is important for the analysis of level correlations.
In Ref. [39], it was shown that the phase transition of this model is first order for weak coupling and seems to terminate for stronger coupling. An important point in the analysis of [39] is that the ground state in the wormhole phase is well approximated by a TFD state. In this paper we have found that except in the limit of strong coupling between the two SYK models, TFD is never a good approximation of the ground state for the number of Majoranas (N ≤ 34) that we can investigate by exact diagonalization. For weaker coupling, coefficients of the expansion of the ground state in a tensor product of the two separate SYK models have a slower power law decay, in contrast to the fast exponential decay with the energy expected in a TFD state. Physically, this suggests a non-thermal state with entanglement properties qualitatively different from those of the TFD state. As the coupling k increases, the overlap of the ground state with the TFD state decreases with a minimum at k c ∼ 0.1. For stronger coupling, the ground state becomes increasingly well approximated by the TFD state. In this large coupling limit, eigenstates of the associated spin operator with the lowest eigenvalues are both close to the ground state of the system and to a TFD state at infinite temperature. Note that the deviations from the TFD state at weak coupling mentioned above do not implicate that the overlap between the ground state and the TFD state is small. The reason is that for small coupling the strength of the wave function in a tensor basis of H S Y K ⊗ H S Y K is concentrated in only one or a few components of the eigenvectors which determine the overlap with the TFD state independently of the distribution in energy of the components.
We have also studied thermodynamic properties and level statistics that provide valuable information on the long time dynamics of the model. Regarding thermodynamic properties, we have fully confirmed [39] that the transition is of first order for weak coupling and it terminates abruptly at k c ∼ 0.2. A careful analysis of the free energy and specific heat shows that for stronger coupling, at fixed temperature, there is no second order transition but just a sharp crossover. To a good approximation, this is also the value of the coupling where the k dependence of the energy gap changes from a power law to a linear dependence which is typical of a non-interacting system. In the strong coupling limit, the spectral density develops well separated blobs centered around the eigenvalues of the spin operator.
Level statistics are also affected by the coupling strength. For k < k c , we observe in the tail of the spectrum important deviations from the random matrix prediction, that we speculate, is an indication that gravitational bound states related to the wormhole geometry are not quantum chaotic. Around the transition we have observed that, suddenly, the lowest eigenvalues become correlated according to random matrix theory. This quantum chaotic behavior is an expected feature [32] of quantum black hole backgrounds. For larger k, the number of eigenvalues in the tail of the spectrum, that is well described by random matrix theory, increases. Eventually, for k ∼ 1, the full spectrum is correlated according to the random matrix theory prediction. Physically, random matrix spectral correlations are an indication of quantum chaotic features and full quantum ergodicity which is believed to be a distinctive feature of quantum black holes.
Therefore, the sudden appearance of random matrix correlations deep in the infrared may be an indication of the instability of the wormhole phase towards the black hole phase. On the (quantum) gravity side, the exact meaning of a sharp crossover between the two geometries is unclear.
In the opposite limit of very small k, well below the transition k = k c , we observe Poisson statistics for very small coupling strength k c 0.03 (once the degeneracies are removed) which suggests that, at least in the range of finite N we can explore, a minimum coupling is necessary to access the wormhole phase. For larger coupling, the bulk of the spectrum becomes quickly correlated according to random matrix theory while deviations in the tail, very likely related to the wormhole phase, are still present until the transition occurs. We did not manage to quantitatively understand these deviations. The limited range of sizes we can explore numerically prevent us to perform a finite size scaling analysis. Additional research is required to clarify whether this limit of very small k or the region of sudden appearance of random matrix correlations mentioned above broadens the phase diagram of the system.
We note that the coupling term is one-body and therefore integrable. This leads us to speculate that for a fixed temperature, the effective interactions among left and right Majoranas decrease as the coupling k increases. It is tempting to speculate that the full phenomenology of the model is consistent with a strongweak coupling transition as k increases such as the Gross-Witten transition [54] or the BEC-BCS crossover [55] observed in cold atom physics.
Interesting problems for further research include, the universality of these results for non-minimal couplings between the two sites or to clarify whether similar physics is observed in supersymmetric analogues of this problem. Another still unresolved problem is the gravity dual interpretation, if any, of the sharp crossover that we observed around the transition. More specifically, it would be interesting to clarify whether both geometries can somehow coexist in this region or simply the system ceases to have a gravity dual interpretation. Another point which would be nice to investigate further is the study of the level statistics for the case of equal couplings: in particular, it would be important to push the numerical studies to larger N to see if the distinction between the bulk and the tail of the spectrum becomes more prominent. To reach this goal, we believe that the method of [63] could be useful. The non TFD nature of the ground state for weak coupling, characterized by a power law decay of the coefficients in a tensor product basis, such as entanglement properties, deserves further investigation as well. Does this power law dependence have a gravity dual interpretation? We plan to address some of these problems in the near future.
We also define which is a symmetry of the total Hamiltonian. The chiral projectors are given by If N/4 is even, the charge conjugation operator commutes with the corresponding projection operator, but this is not the case when N/4 is odd. Then {Γ R 5 , C R } = 0 and if we have an eigenstate of H R , H R |n = λ|n , with Γ 5 |n = g 5 |n (A10) then also To see that a product of the four factors 1 + Γ k+N/2 Γ k (B3) also commutes with terms from the Hamiltonian that contain the index k, i.e. terms of the from it is simplest to use the decomposition Γ k+N/2 Γ l+N/2 Γ m+N/2 Γ n+N/2 + Γ k Γ l Γ m Γ n = 1 2 (P kl P mn + Q kl Q mn ) (B5) with P kl = Γ k Γ l + Γ k+N/2 Γ l+N/2 , with F kl = −(1 + Γ k+N/2 Γ k )(1 + Γ l+N/2 Γ l )/2. We have that , so that the commutator vanishes. This proves the S mod 4 symmetry.
Appendix C: Ground state of the spin operator for N/2 mod 8 = 2 In this appendix we derive that the ground state of the spin operator N/2 mod 8 = 2 is given by for an appropriate choice of the phases of the states |n ± L,R .
Let us consider the case N = 20, 28, · · · (N/2 mod 8 = 2). In this case the charge conjugation matrix anti-commutes with the chiral γ matrix with γ c = i − N 4 γ 1 γ 2 · · · γ N/2 . Both γ c and CK commute with the Hamiltonion of a single SYK system (K is the complex conjugation operator) In a chiral basis with eigenvectors v n,± of H SYK we thus have that CK v n,± is also an eigenvector of H SYK but with opposite chirality because of the anti-commutation relations (C2).
With this convention we obtain m ± | R γ i |n ∓ R = ±i m ± | L γ i |n ∓ L , n ∓ | L γ i |m ± L = ( n ± | L γ i |m ∓ L ) * , where the first equation holds due to the prefactor exp ∓ πi 4 in |n ± R , and the second equation is a consequence of the relation v n,− = C( v n,+ ) * . Combining these two equations we find m ± | R γ i |n ∓ R = ±i n ± | L γ i |m ∓ L . (C6) To show that the state (C12) is the ground state of the spin operator we calculate the expectation value Since γ i flips the chirality, only terms with ± = ∓ contributes: where we have used γ c |n ± = ±|n ± and (C6). Now again noticing the fact that γ i flips chirality, we can safely replace n |n ∓ L n ∓ | L → n (|n + L n + | L + |n − L n + | L ) = 1 in the last line This is exactly the structure that was used in the general proof (20).
The phases of the right states will be chosen as For the last equality we have used that that the γ-matrices for N/2 mod 8 = 0 can be chosen real symmetric.