Traversable wormhole and Hawking-Page transition in coupled complex SYK models

Recent work has shown that coupling two identical Sachdev-Ye-Kitaev (SYK) models can realize a phase of matter that is holographically dual to an \emph{eternal traversable wormhole}. This phase supports revival oscillations between two quantum chaotic systems that can be interpreted as information traversing the wormhole. Here we generalize these ideas to a pair of coupled SYK models with {\em complex} fermions that respect a global U(1) charge symmetry. Such models show richer behavior than conventional SYK models with Majorana fermions and may be easier to realize experimentally. We consider two different couplings, namely tunneling and charge-conserving two-body interactions, and obtain the corresponding phase diagram using a combination of numerical and analytical techniques. At low temperature we find a charge-neutral gapped phase that supports revival oscillations, with a ground state close to the thermofield double, which we argue is dual to a traversable wormhole. We also find two different gapless non-Fermi liquid phases with tunable charge density which we interpret as dual to a `large' and `small' charged black hole. The gapped and gapless phases are separated by a first-order phase transition of the Hawking-Page type. Finally, we discuss an SU(2)-symmetric limit of our model that is closely related to proposed realizations of SYK physics with spinful fermions in graphene, and explain its relevance for future experiments on this system.


I. INTRODUCTION
The Sachdev-Ye-Kitaev (SYK) model [1][2][3][4] has emerged recently as a powerful toy model allowing to glean insight into the behavior of non-Fermi liquids, quantum chaos [5][6][7][8], holography [9][10][11] and strange metallic transport [12][13][14][15]. The model consists of N Majorana fermions coupled via all-to-all, random Gaussian interactions. It derives its predictive power from the fact that, once averaged over quenched disorder, the model can be solved exactly at low temperatures through a large-N saddle-point expansion. Remarkably, at large-N the model can be shown to be maximally chaotic: its out-of-time-order correlators (OTOCs) exhibit an exponentially growing regime with a Lyapunov exponent that saturates the bound on many-body quantum chaos [5]. Such growth indicates fast scrambling of quantum information in the system and underlies the holographic connection between the SYK model and black holes which are also fast scramblers. Motivated by these exciting predictions, a number of proposals have emerged for the physical realization of the SYK model and its variants in atomic [16], optical [17] and solid-state [18][19][20] platforms, or using quantum simulators [21,22].
Interesting physics also occurs when two identical SYK models are coupled by interactions [23] or tunneling [24]. At low temperature the coupling can drive phase transitions to symmetry-broken states [23] or, remarkably, to a phase holographically dual to an eternal traversable wormhole [24] with an AdS 2 throat (AdS 2 refers to the 1+1-dimensional anti-de Sitter spacetime). This phase * ssahoo@mail.ubc.ca † lantagne@phas.ubc.ca enables the transmission of information between two chaotic systems through 'revival dynamics' [25,26] corresponding, in the gravity interpretation, to sending particles through a wormhole [27][28][29][30][31][32][33]. Proposals for the physical realization of such coupled SYK models in condensed matter platforms have also been discussed [34].
A key concept at the heart of the traversable wormhole proposal is the thermofield double (TFD) state. Given two identical copies of a quantum mechanical system, a TFD state is defined as Zβ n e −βEn/2 |n 1 ⊗ |n 2 .
Here Zβ = n e −βEn is the partition function of a single system at inverse temperatureβ, and |n = |Θn where Θ is an anti-unitary symmetry. The TFD is a purification of a thermal ensemble at inverse temperatureβ: it is an entangled state of the two copies, such that tracing over either copy recovers the thermal density matrix for the other. As such, the TFD can be used as a resource state to study thermal properties in quantum simulators [35]. The TFD also obeys an inverted version of time-translation invariance which enables accessing OTOCs using ordinary time-ordered measurements [34] and teleporting states or operators between the two subsystems [36,37]. Recent work has shown how to construct model Hamiltonians with a ground state close to a TFD [24,38] by coupling identical systems, the Maldacena-Qi (MQ) wormhole model [24,39,40] being an example of such a construction.
Variants of the SYK model built from ordinary complex fermions, rather than real Majorana fermions, have also been studied and exhibit similar properties [2,41]. The main difference between the complex and canoni-1 -1 0  [24] to complex fermions, while the α = 1 line with SU(2) symmetric interactions maps to the graphene flake proposal of Ref. [42].
cal SYK model is that the former has a conserved U(1) charge. Importantly, these complex SYK variants, henceforth abbreviated as cSYK, might be easier to realize in condensed matter systems [42,43], where Majorana fermions are notoriously difficult to obtain and control. Further, a number of experimental probes such as spectroscopy [42,44], electrical conductance [45] and thermopower [46] have recently been proposed to identify signatures of complex SYK models.
In this work we investigate the physics of coupled complex SYK models using a combination of analytical arguments, exact diagonalization for small N and saddlepoint solutions for large N . We consider two types of couplings: a simple tunneling term with strength κ, and random two-body interactions with strength α that conserve charge in each system separately. The form of the couplings, illustrated in Fig. 1 (a), is partly motivated by their natural connection to the disordered graphene flake proposal of Ref. [42].
We first show that coupling identical cSYK models with the tunneling term leads to similar physics as the MQ wormhole model: at low temperature the system is gapped, with a charge-neutral ground state close to a TFD. The correct definition of the TFD state in the presence of a U(1) symmetry is however subtle and we explain it in detail. We obtain the finite-temperature phase diagram of the model, showing that the gapped phase with a TFD ground state is separated from a gapless cSYK phase at high temperature by a first-order transition line ending at a critical point. Further, we investigate the dynamics of the system and find that two-point correlations between the two subsystems decay as a power-law in the high-temperature SYK phase, but show periodic revivals in the gapped phase [25,26]. These results lead us to conjecture that the gapped phase of this model, similarly to Ref. [24], is holographically dual to a traversable wormhole. The first-order phase transition between the wormhole phase with a TFD ground state and the charged black hole (cSYK) phase can thus be interpreted as a Hawking-Page type transition [47].
We then consider two-body interactions with strength α that conserve charge in each system separately -leading to a U(1)⊗U(1) symmetry group. At charge neutrality we find two low-temperature phases depending on the interaction strength α. For α < 0 or α > 4, we observe a gapped phase with spontaneous symmetry-breaking from U(1)⊗U(1) down to global U(1), while for 0 ≤ α ≤ 4 we obtain a gapless non-Fermi liquid phase with tunable charge density, with properties similar to the cSYK phase. Combining both types of couplings we find the phase diagram illustrated in Fig. 1 (b) which contains a dome-shaped gapless cSYK phase, separated from the surrounding gapped phase by a zero-temperature firstorder phase transition. Surprisingly, away from charge neutrality we uncover another first-order phase transition to a different gapless non-Fermi liquid describing a smaller black hole with half of its degrees of freedom gapped out. We finally discuss in detail the α = 1 limit of our model, featuring SU(2)-invariant interactions, which is directly relevant to the proposed graphene flake realization of the cSYK model [42]. Using this new connection we revisit the results of Ref. [42] and argue that the irregularly-shaped graphene flake model for weak applied magnetic field admits a cSYK non-Fermi liquid phase.
The rest of this paper is organized as follows. In Sec. II, we review properties of the TFD and show how to define it for a pair of identical complex SYK models with U(1) symmetry. In Sec. III, we couple the two systems via the tunneling term and discuss the resulting gapped phase which we conjecture is dual to a wormhole in the gravity description. In Sec. IV, we add interaction terms between the two systems and explore the resulting phase diagram. In Sec. V we provide a connection of our model to the graphene flake proposal [42] with SU(2)-symmetric interactions, and discuss its relevance for future experiments on this system. More technical contributions, including details of the TFD construction and solutions of the saddle-point equations, in both imaginary and real time, are presented in the Appendices.

II. THERMOFIELD DOUBLE STATE CONSTRUCTION
Consider two identical copies of a quantum system described by Hamiltonians H 1 = H 2 , with eigenstates |n 1 and |n 2 of a common eigenenergy E n . A TFD is an entangled state of the two copies defined by Eq. (1). Such a definition represents a one-parameter family of states which can be generated from where |I = |TFD 0 . (We denote the effective inverse temperature of the TFD state asβ to differentiate it from the physical temperature T = 1/β of the system).
The advantage of the representation in Eq. (2) is that |I is a maximally entangled state between the two copies, and is thus usually simple to write down.
The TFD state has a few important properties which follow from its definition, Eq. (1). First, it is a purification of the thermal density matrix, such that for any operator V a (a = 1, 2) that acts only on one subsystem. Here the expectation value is taken in the state |TFDβ , and the trace is over the Hilbert space of a single subsystem. Although the TFD is not an eigenstate of the full system's Hamiltonian H 1 +H 2 , it is an eigenstate of the difference which implies that two-point correlation functions respect an inverted version of time-translation invariance, for any operators V 1 and W 2 . This property has an analog in the gravity context, where the TFD is used to describe traversable wormholes and time effectively flows in opposite directions on its two sides. Our goal in this section is to construct the TFD state when each subsystem is described by a cSYK model. As we shall see this construction entails a subtlety: Whereas in the canonical SYK model the state |n can be chosen as equal to |n (up to a phase), this is not the case for cSYK where |n carries a different charge quantum number.

A. Complex SYK model
The complex SYK model is the charge conserving variant of the canonical (Majorana) SYK model. It is described by the Hamiltonian where the c i , c † i are N fermionic operators satisfying {c i , c † j } = δ ij . The coefficients J ij;kl are complex Gaussian random numbers with and satisfy the symmetry constraints J ij;kl = −J ji;kl = −J ij;lk = J * lk;ji (8) imposed by the fermionic commutation relations. The model has a global U(1) symmetry, c j → c j e −iφ that expresses the conservation of the total charge A detailed discussion of the cSYK model and its physical properties can be found in Refs. [2, 41, and 48]. Additionally, at charge neutrality (µ = 0) the model has an anti-unitary particle-hole symmetry if we constrain the tensor of couplings J ij;kl to be fully antisymmetric. This anti-unitary symmetry is generated by with K the complex conjugation operator, and transforms c j ↔ c † j up to a sign that depends on N . One can check that P −1 HP = H and Such a symmetry has been discussed in Refs. [41,[48][49][50] and is useful to simplify calculations. Ref. [48] implements it by including additional hopping terms to cancel out the unwanted interaction terms with repeated indices. Here we simply consider fully anti-symmetric J ij;kl , and use the anti-unitary symmetry P to define the states |n that are required to construct the TFD state.
Leveraging the U(1) symmetry, the Hamiltonian Eq. (6) can be block diagonalized in symmetry sectors labeled by the eigenvalues q of the charge operator Q. Since P −1 QP = −Q, at charge neutrality the model has a two-fold spectral degeneracy guaranteed by the mapping between states in the q and −q sectors. For even N there is a special zero-charge (q = 0) sector which maps onto itself under P ; this enforces a two-fold Kramers degeneracy when P 2 = −1, for N mod 4 = 2 (see Table I).
B. Thermofield double state for complex SYK As discussed above, the eigenstates |n q of the cSYK model, in the presence of the anti-unitary symmetry P , form particle-hole pairs in charge sectors ±q. We thus define the TFD as a zero-charge state with where the state |n −q is equivalent to P |n q up to a phase which depends only on the symmetry label q. (An arbitrary phase θ n is not allowed as it cannot be absorbed by a gauge transformation in |n and |n .) This phase can be fixed most easily by choosing a specific infinitetemperature (β = 0) TFD state |I and identifying it with |TFD 0 . The definition in Eq. (2) then automatically fixes the same phases for the finite-temperature TFD states. The state |I must be maximally entangled, as it is the purification of the infinite-temperature density matrix. This motivates the choice of a product of N Bell pairs between the two systems [51], In Appendix A we explicitly show that the state in Eq. (13) is equivalent to with the anti-unitary symmetry defined as Here q and Γ = (−1) q+ N 2 are the charge and fermion parity, respectively, of the state |n q and η is a sign depending on the total number N of fermions, since P −1 c i1 P = ηc † i1 . The parity-dependent phase factor has been discussed for the Majorana version of the TFD state [39]. The charge-dependent phase is required to cancel the minus signs appearing when fermionic creation or annihilation operators are taken across eigenstates of subsystem 1 to act on subsystem 2. Note that the TFD for a bosonic Hamiltonian would not have these phases.
state of this Hamiltonian is the 'Bell pair' state |I φ = |TFD 0 , whereas for κ/J = 0 it is simply the product state |TFD ∞ . As we shall see below using numerical exact diagonalization, the model admits a ground state close to |TFDβ for all κ/J, with a parameterβ that is a monotonically decreasing function of κ/J. Note that a gauge transformation on the fermion operators in either subsystem can absorb the phase φ. Hence in the following we consider, without loss of generality, a purely imaginary tunneling term with φ = π/2, corresponding to the simplest case where the charge-dependent phases in Eq. (15) disappear.
The symmetries and the ground state degeneracy of the model, deduced by simple arguments and verified using exact diagonalization, are summarized in Table I. For κ = 0 we have two decoupled cSYK models, where the charge is separately conserved in systems 1 and 2 (U(1)⊗U(1) symmetry). There are two anti-unitary symmetries P 1 = P ⊗ 1 and P 2 = 1 ⊗ P , where P acts in one subsystem. The ground state degeneracy is the product of that of each SYK model, which is unique for N mod 4 = 0 and doubly degenerate otherwise. For even N , the ground states of the coupled system are always in the q = 0 sector, while for odd N they are distributed in the q = 0, ±1 sectors as shown in Table I. The tunneling term κ breaks the charge conservation in each system down to total charge conservation (global U (1) symmetry), and also breaks the anti-unitary symmetries P 1 and P 2 -thus the ground state for κ > 0 is unique. Finally, there is a discrete mirror symmetry that exchanges fermion operators between systems 1 and 2. For φ = π/2 this symmetry transforms c 1 → c 2 and c 2 → −c 1 , and constrains the two-point correlators of the system as discussed in Sec. III B.

A. Exact diagonalization: TFD ground state
We first perform an exact diagonalization study of the model, Eq. (17), to confirm that it admits a ground state close to a TFD. To do this we construct the family of TFD states with parameterβ from the eigenstates of a single SYK model, using the definition Eq. (12). We then compute the overlap of this family of TFD states with the numerical ground state of the coupled model, Eq. (17) and select the TFD with the parameterβ that maximizes the overlap. As shown in Fig. 2 this best-fit overlap is always close to 1, with a minimum of ∼ 0.96 at a value κ/J ∼ 0.1 which roughly corresponds to the end of the finite temperature first-order transition line seen in the large-N calculation (see Sec. III B). The parameter β characterizing the best-fit TFD is monotonically decreasing with κ. In the gravity interpretation of the MQ model the parameterβ is proportional to the length of the wormhole, or equivalently to the period of the revival oscillations between the two sides [24,25].

B. Large N saddle point solution
We now derive the large-N saddle point equations of the model Eq. (17) in imaginary (Euclidean) time τ , and solve them numerically to investigate its properties in the thermodynamic limit. The partition function of the system at inverse physical temperature β is given by Upon disorder averaging (keeping only replica-diagonal terms) and integrating out the fermion fields we arrive at the following effective action (see Appendix B), where G ab (τ, τ ) = 1 N T i c ia (τ )c † ib (τ ) is the averaged time-ordered correlator at the saddle point, and Σ ab (τ, τ ) are Lagrange multipliers which can be interpreted as fermion self-energies. Varying the effective action Eq. (19) with respect to G ab and Σ ab leads to large-N saddle-point equations. Using the time-translation symmetry G ab (τ, τ ) = G ab (τ − τ ) and the mirror symmetry which enforces G 11 (τ ) = G 22 (τ ) and G 12 (τ ) = −G 21 (τ ), the saddle-point equations can be reduced to the form where ω n = (2n + 1) π β are fermionic Matsubara frequencies, and the parameters J, µ and κ are real. When κ = 0 and thus G 12 = 0, these equations reduce to the canonical cSYK model, which can be solved in the low-frequency or long-time limit by appealing to an emergent conformal invariance. The zero-temperature result is [2] where b = 1/(4πJ 2 ) 1/4 . The 'twist' parameter E, which leads to a spectral asymmetry in the frequency domain, comes about because the chemical potential µ = 0 breaks the particle-hole symmetry of the problem. The resulting The asymmetry parameter E is also related to the thermodynamic quantity where S 0 is the residual zero-temperature entropy density of the cSYK model. In the holographic picture E describes the electric field near the charged AdS 2 black hole horizon, and respects Eq. (24) with S 0 the Bekenstein-Hawing entropy density of the horizon [2,53]. Let us first consider the model at charge neutrality, µ = 0. When κ (and thus G 12 ) is non-zero, an exact solution of the saddle-point equations (20) in the lowenergy limit cannot be obtained. Instead we solve them numerically by iterating until a self-consistent solution is found, choosing the initial seeds for the iteration to be the non-interacting solution G 11 (τ ) = 1 2 sgn(τ ), and G 12 = i with small . We find that only the real and imaginary part of G 11 and G 12 , respectively, are nonzero. For κ = 0 we recover the conformal result Eq. (22) for long times and low temperatures J −1 τ β. As shown in Fig. 3(a), when turning on a small coupling κ/J = 0.03, a gap opens at low temperature, as indicated by the exponential decay of the correlators G 11 and G 12 . For high temperatures the correlators instead decay as a power law. In Fig. 3(b) we show that the energy gap extracted from the exponential decay of G ab (τ ) at low temperature βJ = 10 4 scales as ∼ ( κ J ) 2/3 for κ/J 1. This is consistent with the scaling of the analogous MQ model [24].
The MQ model also exhibits a first-order phase transition at finite temperature for small values of κ. In the gravity context this transition was interpreted [24,54] as a Hawking-Page transition [47] because it separates a stable AdS 2 black hole at high temperature from a lowtemperature phase (the wormhole) which appears thermal for an observer having access to only one subsystem. This is manifest in the wormhole phase admitting a TFD ground state (see also Eq. 3). Such a transition can be identified from the thermodynamics of our complex fermion model. The free energy F = −T ln Z is obtained by substituting the saddle point solutions in the action, Here we regularized the free energy using its value for 2N non-interacting complex fermions, 2 n ln(iω n ) = 2 ln 2, to cancel out divergences at large frequencies in the numerical evaluation of D(iω n ). In Fig. 3c we present the numerically obtained free energy density when sweeping from high to low temperatures, and vice versa. For high temperatures the gapless cSYK solution is favored, whereas for low temperatures the gapped solution prevails. These two phases can be identified from the temperature dependence of the free energy: a constant negative slope S = −∂F/∂T at low temperatures indicates an SYK phase with residual entropy S 0 , while the gapped wormhole phase with a unique ground state shows zero slope. We obtain a clear hysteresis between the two solutions indicating a first-order phase transition. The transition temperature T c , identified with the crossing point of the gapped and gapless free energies in Fig. 3(c), increases monotonously with κ/J until the phase transition line terminates at a critical point, as shown in Fig. 3(d).
When introducing a non-zero chemical potential µ, the ground state is not necessarily charge neutral. For the canonical cSYK model this results in a spectral asymmetry or a 'twist' E in the conformal limit of the imaginarytime Green's functions [2,41], see Eq. (22). The corresponding U(1) charge density Q can be read off from the value of the imaginary-time Green's functions near τ = 0, In Fig. 4 we show the U(1) charge, obtained numerically from Eq. (26) as a function of parameters κ and µ. We find that the gapped wormhole phase is stable to the inclusion of a finite chemical potential µ, even though µ breaks the microscopic anti-unitary symmetry P used to define the TFD state. The charge density of the wormhole phase remains zero throughout. Increasing µ drives the system into a gapless phase with a tunable charge density Q ∈ [− 1 2 , 1 2 ] and then finally to a gapped, polarized state with the maximal charge Q = ±1/2. The transition to the polarized state is first order, where the extensive entropy of the charged black hole phase jumps to zero, similar to the transition seen in Ref. [55].  26) at low temperature βJ = 500. The gapped wormhole phase is stable to inclusion of a finite chemical potential µ and remains charge neutral throughout. For larger |µ| the system transitions to a gapless phase with tunable charge density, then finally to a gapped polarized phase with Q = ±1/2. The dashed lines denote first-order phase transitions .

C. Real-time dynamics
In order to probe the dynamical behavior of the model we now switch to real-time representation of the saddlepoint equations (see Appendix D for details). Following Ref. [25] we focus on the transmission amplitude which expresses the probability amplitude of recovering a fermion in system a at time t after inserting the corresponding fermion in system b at time 0, averaged over all fermionic modes j in the system. Fig. 5a shows the transmission amplitudes for small κ in the low temperature regime T /J = 0.0001 κ/J. They exhibit sharply peaked revival oscillations in both T 11 and T 12 that are out-of-phase, consistent with the propagation of the fermions back and forth between the two chaotic systems. As with the Majorana case [25] we find that the sharp revivals rely on a tower of equally-spaced states in the spectral function, see the inset in Fig. 5a, which occur at harmonics of the gap ∼ κ 2/3 . The overall decay of oscillations is due to the width of those spectral peaks, which increases when going to higher frequencies and/or temperatures [26]. By comparison, for temperatures above the first-order transition shown in Fig. 3d we observe a smooth, power-law decay ∼ |t| −1/2 of the transmission amplitude characteristic of the SYK non-Fermi liquid phase.
Based on the results of this Section we conjecture that two identical cSYKs models coupled with a weak tunneling term admit a low-temperature phase which is holo-graphically dual to a traversable wormhole. We rely on the following observations: (i) the presence of a TFD ground state with largeβ, (ii) a first-order phase transition separating the presumed (gapped) wormhole phase from a gapless cSYK phase at high temperature and (iii) revival dynamics showing the transmission of excitations between the two chaotic subsystems.

IV. COUPLING cSYK MODELS WITH INTERACTION TERMS
In this section we couple the two cSYK models with four-fermion interactions which conserve charge on each system. We consider the Hamiltonian introduced in the previous section, modified by an extra term where the coupling constants J ijkl are identical to those within each cSYK system. Related models have been studied before in the context of symmetry-broken ground states [23] and superconducting instabilities of SYK non-Fermi liquid phases [56]. An additional motivation to study such an interaction term, as explained in more details in Sec. V, is that it naturally arises in proposed physical realizations of the SYK model in graphene flakes [42]. The additional α term alters the symmetries of the model, resulting in different degeneracies. For α = 0 but κ = 0, P 1 and P 2 are no longer symmetries but the combined anti-unitary symmetry P 12 = P 1 P 2 remains, where P 2 12 = (−1) N (2N −1) . This symmetry guarantees a twofold degeneracy for odd N (see Table I). As before, degeneracies are lifted at κ = 0 as the tunneling term breaks the anti-unitary symmetry for odd N , P −1 12 KP 12 = −K. The large-N saddle-point equations are obtained in a similar way as Eqs. (20), with details delegated to App. B. The equations for the Green's functions G ab (iω n ) are unchanged while the expressions for the self-energies Σ ab (τ ) acquire additional terms The low-temperature phase diagram of this model, obtained from the self-consistent numerical solution of the above equations, is shown in Fig. 1 and analyzed in Fig. 6. We find two phases whose properties are discussed in the next subsection: a gapless phase with tunable charge density similar to the cSYK non-Fermi liquid, and a gapped  phase that is adiabatically connected to the α = 0 wormhole solution of Sec. III. The gapped phase persists down to κ = 0 for α < 0 or α > 4 through a symmetry-breaking mechanism, where a finite expectation value K is generated spontaneously.

A. Phase diagram at charge neutrality
The low-temperature phase diagram of the model, Eq. (28) consists of two phases near charge neutrality µ = 0. For α < 0 or α > 4 the system is in a gapped phase for all values of κ, as indicated by the exponential decay of the two-point correlators G ab (τ ) at late times 1 τ J βJ. In Fig. 6a,b we show the gap extracted from that exponential decay at low temperature βJ = 10 4 as a function of κ and α. For 0 ≤ α ≤ 4, we find that a gapless phase survives for a range of κ inside a dome-shaped region, where the correlator G 11 (τ ) ∼ τ −1/2 decays as a power-law with the same exponent as in the SYK phase. For α = 4 a gap opens for any κ = 0, similarly to the α = 0 case analyzed in Sec. III. This can be easily understood from the saddle-point equations (29,30). Using the symmetry of the two-point correlators at charge neutrality, G 11 (−τ ) = −G 11 (τ ) and G 12 (−τ ) = G 12 (τ ) we see that at α = 4 the saddle-point equations reduce to which are just the equations for two decoupled cSYK models (α = 0) but with a renormalized J eff = 3J. When turning on a finite κ we thus expect the same low temperature 'wormhole' physics as for α = 0, but with a renormalized gap E gap ∼ κ 2/3 J 1/3 eff . We verify this scaling from our numerical simulations, as shown in Fig. 6b. The mapping between α = 0 and α = 4 is reminiscent of the duality that exists in the analogous Majorana model in Ref. [23]. However, here the duality is only emergent We then compute the free energy of the model, using the same approach as in Sec. III, Eq. (25). The free energy shows hysteresis across the phase transition between the gapless and gapped phases for 0 ≤ α ≤ 4. The phase transition lines for various α are shown in Fig. 6c. As α increases the first-order transition lines move up the κ axis and have a non-zero intercept, such that the gapless phase extends down to zero temperature for non-zero κ. Thus for 0 ≤ α ≤ 4, the first-order Hawking-Page phase transition occurs at zero temperature upon varying the tunneling strength. This is also indicated by the discontinuous jump in the gap magnitude across the transition shown in Fig. 6b.
Using exact diagonalization we obtain the overlap between the ground state of the coupled model and the TFDs defined in Eq. (12). In the gapped phase the overlap with the TFD decreases continuously when moving away from the α = 0 'wormhole' line, and sharply drops to zero upon entering the gapless phase. Interestingly, at α = 4 the ground state is not well approximated by a TFD, a further indication of the duality between α = 0 and 4 being only valid in the large-N limit.

B. Spontaneous symmetry breaking
When κ = 0 and α < 0 or α > 4, the ground state develops a non-zero expectation value for the tunneling operator K, which can be read off numerically from K κN = 2iG 12 (τ = 0), as shown in Fig. 8. We therefore conclude that the gap opens at κ = 0 through a sponteaneous symmetry-breaking mechanism, from the U (1) ⊗ U (1) charge conjugation symmetry down to U (1). For finite N the spontaneous symmetry breaking can be analyzed using exact diagonalization, providing a simple explanation of the phase transition observed at κ = 0. For even N the ground state is unique and we always find that K = 0. For odd N however, we have two degenerate ground states which, in the gapped phase, are located in the q = 0 sector (see Table I). We can perform a basis rotation in this twofold degenerate space to obtain two eigenvectors of K with opposite eigenvalues, as shown in Fig. 8. The system can thus spontaneously choose a ground state which breaks U (1) ⊗ U (1) symmetry, as in the saddle-point result. In the process the anti-unitary symmetry is also spontaneously broken as P −1 12 KP 12 = −K for odd N . In contrast, in the gapless phase (for 0 ≤ α ≤ 4) the two ground states are in charge sectors q = ±1 and have K = 0 since K conserves charge. There is thus no possible symmetry breaking, in accordance with the large N result.

C. Revival dynamics
The transmission amplitudes [Eq. (27)] for non-zero α at low temperature βJ = 10 4 are shown in Fig. 5b,c. For α = −1 and small κ, deep inside the gapped phase, we again find revival oscillations. Those are notably less sharp than at the MQ point α = 0, consistent with the observation that the ground state is not well approximated by a TFD. The reason is that the gap remains large as κ → 0. Therefore, there is only a small number of states in the conformal tower (at harmonics of the gap) that can fit within the energy scale J which limits the conformal scaling behavior [25]. To this end compare the spectral functions at α = 0 (Fig. 5a inset), showing a large number of evenly spaced peaks, and α = −1 (Fig. 5b inset), showing only few and far-spaced spectral peaks. Thus at α = −1 the revivals are controlled by a few spectral peaks rather than an extensive tower of states. For α = 1 and small κ (deep inside the gapless dome) we find a power-law decay of the transmission down to the lowest accessible temperatures, closely tracking SYK behavior, but now with T 12 (t) = 0 also decaying as a power-law. The spectral function ρ 11 (Fig. 5c  inset) shows gapless power-law behavior as expected in the SYK phase, however ρ 12 = 0 indicates the presence of correlations between both subsystems.

D. Conformal solution and SU(2) symmetry
We now consider whether a low-energy conformal solution of the saddle-point equations (29,30) can be found to describe the gapless phase. When κ = 0, we have G 12 (τ ) = 0 and the saddle-point equations can be solved at low energies to show that G 11 (τ ) is a conformal cSYK correlator with J 2 replaced by J 2 α = J 2 (1 + 1 2 α 2 ). Note that this connects smoothly with the behavior at α = 0, 4 discussed in Eqs. (31). With κ > 0, G 12 (τ ) does not vanish which renders the analytical solution of the saddle-point equations more difficult. In the low energy limit, we find that a conformal solution for both G 11 and G 12 is in general not possible, except at the special point α = 1 discussed in Appendix C. This limit admits a power-law solution with the same power for both correlators, G ab (τ ) ∼ b ab |τ | −1/2 in the long-time limit 1 τ J βJ, and where the two coefficients are related by b 4 12 = b 4 11 − 1/6πJ 2 . This is demonstrated numerically in Fig. 9. Note that the saddle-point equations alone are not sufficient to fix the coefficients of the two power laws. Instead one must impose a constraint linking microscopic and conformal physics, in a way similar to the way U(1) charge enters the conformal solution in the cSYK model [2,52] (see Appendix C).
The α = 1 point is special because it has SU(2) symmetric interactions. An important consequence is that the tunneling term K is now a symmetry of the model, [H α=1 , K] = 0. Hence K is a conserved quantity which can be tuned by the tunneling parameter κ, in analogy with the U(1) charge density Q tuned by the chemical potential µ. However, the two symmetries have different signatures: introducing κ leads to a non-zero value of iG 12 but does not introduce a twist parameter leading to a spectral asymmetry, as occurs with non-zero µ. Another consequence of the SU(2) symmetry is that the non-interacting ground state |TFD 0 of K is an eigenstate of the full model for any κ. In fact, exact diagonalization shows that for κ > κ c 0.27J the model admits the |TFD 0 state as an exact ground state. Because the interaction and tunneling terms commute, the only way to change that ground state is through an energy level crossing, which occurs at the first-order transition at κ c , when the gapless phase becomes favored.

E. Moving away from charge neutrality: a tale of two black holes
We now discuss the physics of our model away from charge neutrality. The α = 0 limit was analyzed in Fig. 4 which showed the stability of the wormhole phase to a finite µ, but where no new phases (i.e. absent at µ = 0) appeared apart from the gapped, polarized phase with Q = ±1/2. As we show below the situation is more interesting when both α and µ are non-zero.
We focus on the SU(2) symmetric point α = 1 and investigate its low-temperature phase diagram in the κ, µ plane. To this end we show in Fig. 10 the U(1) charge obtained numerically from Eq. (26) as a function of parameters κ and µ, and the residual entropy density S 0 at low temperature βJ = 500. We recover the known gapped phases discussed above: for large κ, which admits the non-interacting charge-neutral |TFD 0 ground state, and for large |µ|, which corresponds to the polarized Q = ±1/2 states. The boundaries of the two gapped phases host first-order phase transitions to gapless non-Fermi liquids, as indicated by the discontinuous jump in entropy density at the phase boundaries. Surprisingly, we find not one but two such gapless phases. Near charge neutrality, we obtain a phase smoothly connected to the conformal solution discussed above. In this phase both correlators G 11 and G 12 show power-law decay which indicates strong correlations between the two subsystems. This phase can be thought of as a single cSYK phase with 2N fermions, dual to a 'large' black hole comprising all degrees of freedom in the combined system.
Farther from charge neutrality we find another firstorder phase transition to a different gapless phase with charge density Q 0.25 and about half of the residual entropy at charge neutrality. To understand this, note that at α = 1 where the interactions are SU (2)  we can rotate to a new basis c j± = 1 In this basis one can interpret the system as two cSYK models with different chemical potentials µ ± = µ∓κ and SU(2) invariant interactions between them. Let us first focus on the κ = µ line where the chemical potentials are simply 0 and 2µ. When µ increases, eventually one of the cSYK models undergoes a first-order phase transition to a gapped, polarized state with Q = 1/2. At low energies (below the gap), its degrees of freedom thus decouple and we are left with the other cSYK model at charge neutrality Q = 0. The combined system thus has exactly Q = 1/4 and S 0 = S cSYK , as observed in the saddle-point solutions. In Fig. 11 we show the spectral functions ρ ± for the rotated basis fermions c ± . We observe a power- law scaling at low frequency in the ρ + channel while the ρ − channel is gapped, confirming the argument above.
In imaginary time, the corresponding correlators G ± (τ ) show power-law and exponential decay, respectively, at long times.
In the vicinity of the µ = ±κ lines the residual entropy and charge density change smoothly, see Fig. 10, as one half of the system is in a compressible cSYK state while the other half remains gapped. We interpret this phase as dual to a 'small' black hole, comprising half of the degrees of freedom of the combined system, with the other half decoupled and frozen into a fully polarized state.

V. PHYSICAL REALIZATION IN GRAPHENE FLAKES
We now turn to potential physical realizations of the model introduced above. As originally described in Ref. [42], a promising platform for realizing cSYK physics is a mesoscopic graphene flake under a perpendicular magnetic field B, with the chemical potential µ lying within the zeroth-Landau level (LL 0 ). An Aharonov-Casher argument [57] implies that LL 0 remains sharp as long as the chiral (sublattice) symmetry of the model is unbroken, thus forbidding two-fermion terms that would destroy the non-Fermi liquid physics at low energies. Disorder that preserves this chiral symmetry (such as an irregular boundary) then imprints disorder on the LL 0 wavefunctions without lifting their degeneracy, leading to random and all-to-all interactions between them.
Due to the negligible spin-orbit coupling in clean graphene, it is reasonable to assume identical wavefunctions for the two spin components. This should still hold in the presence of non-magnetic disorder (such as an irregular boundary) which preserves the SU(2) symmetry of Coulomb interactions. The graphene flake setup thus naturally leads to two identical copies of the cSYK model, one for each spin component. In Ref. [42] the authors argued that the Zeeman splitting obtained by applying a magnetic field to the sample generates a large spin gap (augmented by exchange interactions), which effectively reduces the problem to a single cSYK model. In this Section we revisit this analysis by looking more carefully at the role of spin in the above proposal. Using a mapping to the model studied in Sec. IV we conclude that the graphene flake model with a weak Zeeman splitting is in a gapless cSYK phase with fermion scaling dimension 1/4 and tunable charge density. In contrast to expectations that a strong Zeeman splitting should give rise to a cSYK phase [42], we find that it instead leads to a gapped phase with an exact |TFD 0 ground state.
The Coulomb interactions between electrons in the graphene flake read where ρ r = ρ r↑ + ρ r↓ is the total charge density at a point r in space and V (r − r ) is the screened Coulomb potential. The electronic charge densities can be expressed in terms of the eigenfunctions φ j (r) of the noninteracting Hamiltonian which, neglecting spin-orbit coupling effects, are independent of spin σ =↑, ↓, Projecting to the LL 0 wavefunctions, Eqs. (34,35) lead to a (normal-ordered) interaction Hamiltonian with all-toall couplings H int = ijkl J ij;kl c † iσ c † jσ c kσ c lσ and spinindependent coupling constants Assuming spatially random wavefunctions and strong screening these become random Gaussian coupling constants [42] (see Ref. [58] for the influence of varying screening lengths in a related model).
Adding the Zeeman term the Hamiltonian describing the low-energy physics of the graphene flake becomes The Coulomb interactions in this model are invariant under SU(2) rotations in spin space. We can thus perform a basis change corresponding to a rotation of the spinor (c i,↑ c i,↓ ) by π 2 along the x-axis, which gives This has the same form as the model in Eq. (28) with α = 1, if we identify gµ B B with the tunneling amplitude κ, and up/down spin projections with subsystems 1 and 2. As mentioned in the previous section, the SU(2) symmetric interactions commute with the tunneling term. In the original basis of Eq. (37), this is easily seen by noting that the Zeeman term is proportional to the total spin projection S tot z = i σ z i . Using this mapping we thus expect that the graphene flake remains in the gapless cSYK non-Fermi liquid phase up to a threshold value κ c = gµ B B c ∼ 0.27J. Above the critical field strength B c the system becomes gapped and the ground state is the infinite-temperature TFD state |I , which is a product state of spin singlets. In Ref. [42] the authors estimate that for B ∼ 20 T the Zeeman splitting gµ B B ∼ 2.4 meV while the characteristic Coulomb interaction strength J ∼ 25 meV, thus placing the system within the gapless phase. Given various uncertainties in these estimates J could well be smaller in which case the system would realize the gapped phase indicated in Fig. 10. Further, the power-law scaling characteristic of the conformal regime is expected for temperatures T J eff ∼ 3 2 J at the SU(2) symmetric point α = 1, corresponding to T 430K which should enable exploration of the low temperature regime.
As shown in Sec. IV D the cSYK non-Fermi liquid phase for B = 0 is stable against inclusion of a chemical potential up to a threshold value µ c which corresponds to adding charge density Q = 1/2, filling all the states in LL 0 . For B = 0 the system first transitions to an intermediate non-Fermi liquid phase corresponding to a small black hole. This transition could be explored by tuning the chemical potential in the graphene flake by external gates. Experimentally, the non-Fermi liquids could be distinguished from each other, and from the gapped phases at large κ or large µ, by measuring their spectral function through scanning probe techniques or their charge compressibility ∂Q/∂µ [10,41].
An important caveat is that in Eq. (38) the coupling constants J ij;kl are only restricted to be antisymmetric under exchanging (iσ, jσ) and (kσ, lσ) by fermionic commutation relations. This is a weaker requirement than the full antisymmetry present for the Majorana SYK model and assumed in this work (the same assumption was made in Ref. [42]). For example, interaction terms with two pairs of matching indices, corresponding to either direct (density) interactions or to exchange (spin) interactions [59] are excluded from our model. A detailed analysis of such effects is left for future work.

VI. CONCLUSION AND OUTLOOK
In this work we generalized the 'eternal traversable wormhole' construction of Maldacena and Qi [24] to a system of coupled complex SYK models with a global U(1) charge symmetry. We explained how to define the TFD state in the presence of a U(1) symmetry, and showed that the model admits a gapped phase with a ground state close to a TFD for all values of tunneling κ between the two subsystems. Whether the weaktunneling and low temperature limit of the model admits a gravitational dual similar to the wormhole of Ref. [24] remains an intriguing open question for the high-energy community. The presence of a gapped ground state close to a largeβ thermofield double, a first-order Hawking-Page phase transition to a gapless cSYK non-Fermi liquid at high temperature and sharp revivals in fermion transmissions are however highly suggestive.
Further, we considered the effect of four-fermion interactions between the two cSYK models that are disordered identically to the interactions within each system. We explored the phase diagram of the system as a function of tunneling, interactions and chemical potential in Figs. 4, 6, and 10. At low temperature we obtain three non-trivial phases: a gapped, charge-neutral phase which is adiabatically connected to the conjectured wormhole and two gapless, compressible non-Fermi liquid phases which describe either a 'large' or a 'small' (i.e. with half of its degrees of freedom gapped out) charged black hole with an AdS 2 horizon. All phases are separated by firstorder phase transitions exhibiting extensive residual entropy jumps. The transition out of the gapped wormhole phase can be understood as a Hawking-Page transition, as it separates a black hole from a phase appearing locally thermal, a consequence of its TFD ground state. The phase diagram also contains the special case α = 1 with SU(2)-symmetric interactions which admits a conformally invariant solution at low energies, and is directly relevant to the graphene flake proposal of Ref. [42].
We conclude by highlighting a few caveats of our analysis and point out interesting directions for further work. First, in order to define the anti-unitary particle-hole symmetry P which enables the TFD state construction, we restricted the model to only contain interactions that are completely antisymmetric in the indices i, j, k, l. In the large-N limit these should dominate as their num-ber scales as N 4 (in contrast, the number of terms with one or two pairs of identical indices, of the form J ij;ik or J ij;ij , respectively scales as N 3 and N 2 ). However, for mesoscopic realization of SYK physics with finite N these terms could be important; understanding their effect will be an important step towards connecting our results with ongoing experimental efforts.
A different approach to define a TFD for cSYK models could be to rely on an anti-unitary time-reversal (rather than particle-hole) symmetry. Time-reversal symmetry is obviously broken for a single cSYK model (as manifest by the coupling constants J ijkl being complex), but can be restored globally by considering a pair of timereversed cSYK models. We anticipate that the TFD construction, saddle-point physics and physical realizations will be different in this case, and leave its detailed study for future work. Another interesting topic concerns the quantum chaotic properties of our model. It would be of great interest to explore how scrambling, as captured through out-of-time-ordered correlators, behaves across the Hawking-Page transition which separates the wormhole and black hole phases.
Finally, in light of the interesting phases of our model for generic α = 1, it is interesting to ask about their potential physical realization. In the graphene flake proposal, α = 1 is enforced by the SU(2) symmetry of Coulomb interactions. However, in other model systems where the two subsystems are realized by two surfaces, such as multilayer graphene [34] or a topological insulator flake, one could imagine changing the distance between the two surfaces as a way to tune the ratio of inter-system to intra-system interactions α. In this way one could potentially explore the Hawking-Page phase transition between the wormhole and black hole phases discussed in this work.
is a unique ground state of the tunneling Hamiltonian term, K = κ i c † i1 c i2 + H.c. with κ a complex coefficient |κ|e iφ . Showing that the overlap I|TFD 0 = 1, where |TFD 0 is defined in Eq. (14), is equivalent to showing that the expectation value of K in the TFD state has minimum value i.e. −|κ|N . Consider the general definition of TFD in Eq. (1) atβ = 0 and evaluate the expectation value We insert an identity operator using complete set of basis states with appropriate normalization, I = q ,q m ,n |m q 1 ⊗ |n q 2 n q | 2 ⊗ m q | 1 to separate the operators acting on system 1 and 2, The phase factor e +iπ(q+ N 2 ) arises when we commute c i2 across states in system 1 to act on site i of system 2. Summing over Kronecker delta symbols we find (A4) We now observe that matrix elements in Eq. (A4) can be evaluated separately for system 1 and 2 which, importantly, are identical. The matrix elements can thus differ between the systems at most by a phase. We choose this phase so that it cancels the phase factors present in Eq. (A4), namely This corresponds to the following definition of the anti-unitary symmetry where Γ = (−1) q+ N 2 is the fermion parity of the SYK eigen-state |n q and P = i (c † i1 + c i1 ) and η = (−1) is a sign that depends on total number of fermion in SYK model such that P −1 c i1 P = ηc † i1 . To check this note that the expectation in Eq (A5) is only non-zero when q = q − 1 and a similar argument holds also for κ * term when q = q + 1 s.t. n −q | 2 c i,2 |m −q 2 = e −iφ e i π 2 e −iη π 4 (Γ −Γ) n q | 1 P −1 c i1 P |m q 1 δ q,q −1 (A7) where phases are complex conjugated. Noticing (Γ − Γ) is 2 or −2 when q + N 2 is odd or even, we can replace the phase simply by −ηe −iφ e −iπ(q+ N 2 ) . Using P −1 c i1 P = ηc † i1 , the above expression becomes Eq. (A6) gives the TFD definition quoted in the main text, Eq. (14), with the required phase factors. The remaining task is to show that such an expression for TFD gives the proper expectation value for the tunneling operator. Substituting Eq. (A5) in Eq. (A2), the expectation value becomes where both the states and operators now refer to system 1. Summing over |m q and recalling that κ = |κ|e iφ , we have We used the fact that the number of states in charge sector q with fermion number n q = q + N 2 is N nq and that N nq=0 n q N n q = N 2 2 N .
The reason is that, for the saddle-point solution of the SYK model, the replica off-diagonal terms can be ignored as they do not contribute to zeroth order in 1/N . Performing the replica trick with only replica-diagonal terms is formally equivalent to the quenched disorder average. Combining with the free part of the action, we finally obtain the averaged partition function for the fermions Z = D[c † , c]e −S with the effective action We now integrate out fermions by introducing the averaged Green's functions G ba (τ , τ ) = 1 where the Lagrange multipliers Σ ab (τ, τ ) play the role of fermionic self-energies. After integrating out fermions, the effective (G, Σ) action for the averaged Green's functions and self-energies becomes where the dots represent terms not explicitly containing Σ(ω) and Finally using δ ln Det(M (ω))/δΣ mn (ω ) = 2(M −1 ) nm δ(ω − ω ) one obtains the saddle point equations (20) and (29) quoted in the main text.
This solution represents decoupled SYK models with no correlations (at the saddle-point level) between the two sides, where the only effect of α was to renormalize the constant b 11 . This solution can thus only represent the κ = 0 limit of our model. This is indeed the numerically obtained solution for κ = 0 and 0 < α < 4, in the gapless phase. On the other hand, we do not find a conformal solution for α < 0 or α > 4 as the system develops a gap through the symmetry breaking mechanism discussed in the main text. For non-zero κ we always obtain a non-vanishing correlator G 12 which is inconsistent with the solution above. Thus a conformal solution cannot be found for generic points inside the gapless dome. However, at the SU(2) symmetric point α = 1 the equations above have additional structure. Using the relation s 11 b 12 = −s 12 b 11 , it is clear that the two equations (C8) are now equivalent. We thus have an under-constrained system, which we can solve for b 4 12 = b 4 11 − 1 6πJ 2 . (C10) As shown in Fig. 9, this relation appears to be satisfied numerically for α = 1 and small κ. To fix the value of b 11 , we notice that there is another constraint coming from the relation where K is the tunneling operator. For α = 1 the value of K is a good quantum number of the system, because [K, H α=1 ] = 0. This is similar to the case of finite chemical potential µ, for which [2] G 11 (τ → 0 + ) = 1 2 + Q, where Q is the conserved U(1) charge of the system, and is related to the asymmetry parameter E appearing in the low-energy Green's function, Eq. (22) in imaginary time [1,2]. It is possible to directly compute the value of Q from the microscopic theory [52], and relate it to the conformal scaling parameter E through Eq. (23). Similarly, here we have at zero temperature which provides the second constraint allowing to fix b 11 and b 12 . It is not clear if an analytical result can be obtained for this constraint, as the full form of G 12 (ω) at all energies is needed.