Traversable wormhole in coupled SYK models with imbalanced interactions

A pair of identical Sachdev-Ye-Kitaev (SYK) models with bilinear coupling forms a quantum dual of a traversable wormhole, with a ground state close to the so called thermofield double state. We use adiabaticity arguments and numerical simulations to show that the character of the ground state remains unchanged when interaction strengths in the two SYK models become imbalanced. Analysis of thermodynamic and dynamical quantities highlights the key signatures of the wormhole phase in the imbalanced case. Further adiabatic evolution naturally leads to the 'maximally imbalanced' limit where fermions in a single SYK model are each coupled to a free Majorana zero mode. This limit is interesting because it could be more easily realized using various setups proposed to implement the SYK model in a laboratory. We find that this limiting case is marginal in that it retains some of the characteristic signatures of the wormhole physics, such as the spectral gap and revival dynamics, but lacks others so that it likely does not represent the full-fledged wormhole dual. We discuss how this scenario could be implemented in the proposed realization of the SYK physics in quantum wires of finite length coupled to a disordered quantum dot.


I. INTRODUCTION
The proposal of the SYK model [1][2][3] as dual holographic description for a quantum black hole in (1+1) dimensional anti-de Sitter space AdS 2 has motivated a flurry of research in scrambling of quantum information and many-body quantum chaos. The model consists of N Majorana fermions interacting through random all-to-all interactions and is defined by the Hamiltonian Here, χ i are Majorana fermion operators that satisfy the anti-commutation relations {χ i , χ j } = δ ij and J ijkl are real constants drawn from a Gaussian ensemble obeying J 2 ijkl = 6J 2 /N 3 . Two identical SYK models coupled via one-to-one bilinear coupling have been shown to realize a holographic dual of an eternal traversable wormhole [4]. The quantum ground state corresponding to the wormhole is a thermofield double (TFD) state that can be regarded as a particular purification of the thermal Gibbs ensemble at inverse temperatureβ and is given by Zβ n e −βEn/2 |n L ⊗ |Θn R .
Here, |n L and |n R are eigenstates of two copies of the original system, denoted as 'left' and 'right'. They are related by an anti-unitary symmetry Θ such that |n and Θ |n are degenerate with energy E n , and Zβ = n e −βEn is the partition function of a single system. Preparing the TFD state for an arbitraryβ is a nontrivial task, which was made possible in the Ising model through a variational method [5,6]. In a parallel series of studies it was shown that TFD state can occur as an exact or approximate ground state of two identical copies of a system, provided they are coupled in a specific way [4,7,8]. Here, the coupling strength sets the temperaturẽ β of the TFD. For example, in Ref. 4

by Maldacena and
Qi (MQ), the ground state of two identical SYK Hamiltonians (1) coupled by has been shown to satisfy with high fidelity. In the gravity context, such coupling leads to negative null energy making the wormhole traversable, i.e., particles can be sent through it [9][10][11][12][13][14]. The MQ model exhibits sharp revivals of scrambled local perturbations [15] and undergoes a finite temperature Hawking-Page transition [16] between wormhole and black hole phases [17]. These effects have been identified as characteristic signatures of the wormhole phase. Significant efforts have been made to physically realize the SYK model in atomic [18], optical [19], and solid-state [20][21][22] platforms, or using quantum simulators [23,24]. Proposals for implementing and observing the wormhole physics in the MQ model and its complex fermion variants have also been made [8,25,26]. The key challenge in these proposals is that one must realize two SYK systems with random couplings J ijkl that are identical at the microscopic level. In solid-state realizations such randomness inevitably comes from disorder which cannot be controlled, hence the challenge.
Here, we present an approach to this challenge based on a simple adiabatic argument. We consider a variant  21. Majorana zero modes localized at the two ends of proximitized quantum wires are weakly coupled to two disordered quantum dots, L and R. We consider relatively short wires such that χ L j and χ R j are weakly hybridized due to the overlap between their wavefunctions inside each wire. Because microscopic realizations of disorder will inevitably be different in L and R quantum dots the corresponding J L ijkl and J R ijkl will be uncorrelated. However, when the R dot is removed altogether then the setup trivially realizes the maximally imbalanced MQ model defined by Eq. (7) with η = 1. of the MQ model with identical interactions in L and R subsystems except for an overall scale factor, leading to an MQ model with imbalanced interaction strengths. We show that wormhole behavior persists in this imbalanced model all the way to the limit when the interaction strength in one of the systems vanishes. This maximally imbalanced limit represents an interesting marginal case which, as we show, still exhibits some vestiges of the wormhole physics. Crucially, this case could be amenable to experimental study because it only requires realization of a single SYK system whose constituent fermions are each coupled to a free Majorana mode.
The idea relies on the identity which can be verified directly from the definition of the TFD state, Eq. (2). Combining this with Eq. (4) it follows that |TFD also satisfies with the Hamiltonian which defines the imbalanced MQ model. If |TFD was an exact ground state of the original MQ Hamiltonian H η=0 then by virtue of identity (5) it would remain an exact eigenstate of H η for any value of η. Because, however, Eq. (4) holds only approximately, it remains to be seen how well Eq. (6) holds away from η = 0. On general grounds we expect that when η is slowly increased from 0, |TFD remains the ground state as long as the gap in the spectrum does not close. We also note that, as pointed out in Ref. 27, for noninteracting fermions (i.e. systems described by H L,R quadratic in fermion operators) the ground state of H η=1 is exactly the same as the ground state of H η=0 .
In the rest of this work we show, through a combination of exact diagonalization (ED) and large-N saddle point solutions, that the system defined by Hamiltonian (7) indeed remains gapped for η between 0 and η c at any finite µ. The critical imbalance parameter η c > 1 depends on µ and approaches 1 from above as µ → 0. At the special point η = 1, the TFD state can be realized approximately as the ground state of 2H L SY K + H LR . As mentioned above this construction requires only a single realization of the SYK model and therefore circumvents the enormous challenge of fabricating two SYK models with identical couplings J ijkl . We will discuss below how a small modification to the setup proposed in Ref. 21 to realize SYK physics in Majorana wires coupled to a quantum dot could be used to implement and probe the wormhole physics in the maximally imbalanced η = 1 limit. The proposed setup is sketched in Fig. 1.
Previously, three key signatures of the wormhole phase have been identified in the literature: (i) the presence of the thermofield-double (TFD) ground state, (ii) a temperature-driven first order phase transition of Hawking-Page type, and (iii) revival dynamics showing transmission of excitations between two maximally chaotic subsystems. We will investigate these indicators and find that all three persist in the imbalanced MQ model for η < 1. At η = 1 the Hawking-Page transition disappears and a weak form of the revival dynamics is observed. These observations suggest that the imbalanced MQ model is dual to a traversable wormhole for all η < 1. The maximally imbalanced case η = 1 shows some tantalizing signatures of the wormhole physics but is likely not a full-fledged wormhole dual.
In Sections II A and II B we present an analysis of the ground state and the phase diagram for general values of the tunneling strength µ and the imbalance parameter η. We discuss the revival dynamics of the ground state for the imbalanced case in Sec. II C. In Sec. III, we analyze the revival dynamics for the η = 1 limit and discuss how it differs from the characteristics revival dynamics of a canonical MQ wormhole. The realization of such a limit in the proposed SYK model implementation of Ref. 21 is also discussed here. We conclude with a summary and open questions in Sec. IV.

II. THE MODEL AND ITS PROPERTIES
The Hamiltonian of the imbalanced MQ model can be written more compactly as where index α assumes values (+1, −1) for (L, R) in numerical expressions. A constraint on the coupling, J ijkl = −J jikl = −J ijlk = J jilk , is imposed since fermion operators anticommute. Unlike the canonical MQ model Hamiltonian (8) lacks the discrete mirror symmetry defined by χ L j → χ R j , χ R j → −χ L j . The balanced limit (η = 0) also has an additional global Z 4 "spin symmetry" which the imbalanced model lacks. If the spin operator is defined as Therefore S mod 2, i.e. chirality, and S mod 4 are both symmetries of the MQ Hamiltonian. When η = 0 we are left with only chirality as a symmetry. These differences have implications for the level statistics analysis [17] and for the ground-state degeneracy. For any η the model has antiunitary symmetry Θ that transforms χ L → χ L , χ R → −χ R along with complex conjugation. This symmetry is crucial for our purposes as it enters the definition of the TFD state in Eq. (2).

A. TFD ground state
We use a spinor representation (2 N dimensional matrix) of Majorana fermion χ α i operators satisfying the anticommutation algebra χ α i , χ β j = δ ij δ αβ to perform an ED study. Here N denotes the number of Majorana fermions in one subsystem for the total of 2N Majoranas. The degeneracy of the SYK ground state depends on N with a unique SYK ground state guaranteed when N is a multiple of 8. Here we focus on systems with 2N = 16 and 32 to numerically find the lowest 50 eigenvalues E n and corresponding eigenvectors |E n . We determine the TFD inverse temperatureβ by minimizing the ground-state fidelity f TFD (E 0 ) = | E 0 |TFDβ | 2 with respect toβ. For this TFD temperature we then calculate f TFD (E n ) for the excited states. While the TFD inverse temperatureβ of the ground state depends on the coupling strength µ, we findβ to be nearly independent of η. Fig. 2(a-c) shows the evolution of the 50 lowest energy levels E n as a function of 0 ≤ η ≤ 2 for several different values of coupling strength µ. The colorscale indicates the overlap f TFD (E n ). In all three panels we clearly observe a single state (highlighted in red) with the overlap close to 1. As expected from our adiabatic argument, the energy of that state is nearly unaffected by η, at least for the small number of Majoranas 2N = 16 and 32 considered in our ED study. In the first panel (µ = 0) the two SYK models are uncoupled and hence the eigenstates are of the form |n ⊗ |m with energies E n,m = E n + E m + η(E n − E m ). In the large-N limit, the model is gapless, however in the ED calculation, the model has a finite-size gap that is exponentially small in N . This allows for observation of the TFD state even at µ = 0 which would otherwise be hidden in the degenerate ground state manifold.
For µ > 0 the energy level dependence on η becomes nonlinear. The gap between the ground state and the excited states decreases as η grows and eventually closes for critical strength of η. In the limit of µ → 0, the gap closes exactly at the special point η = 1. For finite µ, the critical point shifts to larger η > 1. Below we compare this with the phase diagram in the parameter space of µ and η. as the ground state is only close to a TFD.

B. Phase diagram
We supplement the above discussion based on ED for small number N of Majorana fermions with a saddle-point solution that becomes exact in the large N limit, βJ N . Standard procedure (outlined in Appendix A) leads to a system of saddle-point equations for the time-ordered fermion correlator G αβ (τ, τ ) = and the associated self energy Σ αβ (τ, τ ) which, for the imbalanced MQ model, take the form whereᾱ = −α and We used the anti-symmetry property of the fermion fields to impose G αβ (τ ) = −G βα (−τ ) and thus reduce the number of equations. For µ = 0 these equations can be exactly solved in the low frequency limit (βJ ω), giving a conformal solution for the diagonal Green's functions, G αα (τ ) = bτ − 1 2 , with b = 1/(2πJ 2 (1 + αη)) 1 4 , while the off-diagonal Green's function, G LR (τ ) vanishes. From the numerical solution of Eq. (10) [see Fig. 3(a) and discussion below], we find that for η ≤ 1, an infinitesimal µ gaps out the model, while for η > 1, the model remains gapless until a critical value for µ. In the parameter range where the model is gapless, the low frequency conformal solution has a scaling symmetry in η. Given a solution by rescaling the correlators according tõ We find that the conformal solutions at different points in the gapless region of the (η, µ) phase diagram are related by an overall scale factor and decay with the SYK power law G ∝ |τ | −1/2 . The boundary separating the gapless phase from the gapped phase is shown in the (η, µ) phase diagram Fig. 3(b). For small values of µ the gapped phase exhibits an approximate conformal symmetry where G αβ (τ ) can be approximated by reparametrization of the conformal correlator and the gap scales as E gap ∼ µ 2/3 J 1/3 for η < 1. For general µ, the numerically obtained exponential decay of the imaginary time Green's function gives the gap dependence on η as shown in Fig. 3(a). The transition between the gapless and the gapped phase along the η axis is continuous; the order parameter vanishes smoothly and the free energy smoothly changes across the transition (not shown). This motivates us to fit the gap to a BCS-type gap function. The functional form a 1 tanh(a 2 η c /η − 1) is used to fit to the dependence of the gap on η and to find the phase boundary in (µ, η)parameter space close to zero temperature. The system free energy F is the large-N action given as Eq. (A6) in Appendix A evaluated for the saddle point solutions and reads The free energy difference ∆F between gapped and gapless phases, as obtained by sweeping the temperature T up and down during the self-consistent solution of the saddle point equations (10), is displayed in Fig. 4. For η < 1 we observe pronounced hysteresis indicating a first-order Hawking-Page type transition [16] between the low temperature SYK phase and the gapped wormhole phase. The region in the µ-T phase diagram where the two phases can coexist shrinks as η increases only to finally disappear when η → 1.

C. Revival dynamics
A prominent signature of the wormhole behavior is the revival dynamics [15], manifest as an oscillatory behavior of the transmission amplitude Here, the greater Green's function G > αβ (t) = 1 N θ(t) i χ α i (t)χ β i (0) may be interpreted as the amplitude of recovering the excitation χ β i |GS , inserted on side β = L/R and probing on side α at a later time t, averaged over all modes i. Inside the black hole phase, any excitation is rapidly scrambled and the transmission amplitude decays according to the power law t −1/2 . In the wormhole phase, periodic revivals are expected where T LR and T LL oscillate out of phase [15].
To compute the real-time Green's functions, we analytically continue the imaginary-time saddle-point equations (10). The detailed procedure is given in Appendix B, yielding the following set of equations for retarded propagators and self energies: Σ ret αβ (ω) = −iJ 2 (1 + αη)(1 + βη) Here, Equations (14) are solved by numerical iteration using the fast Fourier transform algorithm to pass between time and frequency domains. The spectral function is obtained as the imaginary part of the mode-averaged retarded Green's function, ρ αβ (ω) = − 1 π ImG ret αβ (ω). Using the relation one can easily extract the greater Green's function, and thus the transmission amplitude T αβ (t), from the knowledge of the spectral function. At low temperatures, for 0 ≤ η < 1 the ground state is gapped as shown before, and additionally the spectral function ρ(ω) develops a series of equidistant peaks shown in Fig. 5(a). For µ/J small, they represent the conformal tower of states in the MQ model, E n = E gap (4n + 1) with n integer [4,15]. Here we confirm that this also holds in the imbalanced case as shown in Fig. 5(b) for η = 0.5, where we plot the peak energies as a function of µ. Fourier transform of the spectral function autocorrelation is equal to the absolute square of the transmission amplitude [15]. Hence, the equidistant peak structure of ρ LL and ρ LR leads to revival oscillations in the transmission amplitude of the excitation. Since the peaks in ρ LR have alternating sign structure, there is a π/2 phase difference between T LL and T LR as shown in Fig. 5(c). The large-N approach can be used to indirectly test for the proximity of the ground state to the |TFD by considering the expectation value H L SYK − H R SYK which, according to Eq. (5), should vanish in the TFD state. As we show in the Appendix C, the gapped ground state obtained by solving the large-N equations (14) indeed has this expectation value very close to zero for all 0 ≤ η < 1 and non-zero µ. This suggests that the ground state of imbalanced MQ model remains close to |TFD for all N .

III. THE MAXIMALLY IMBALANCED LIMIT
A. Wormhole or black hole?
The saddle point solution at η = 1 shows that the ground state is gapped for finite µ but the Hawking-Page transition does not occur as temperature increases. This already indicates that the physics for η = 1 is different compared to the ordinary MQ model. Additional differences are observed in the revival dynamics as indicated in Fig. 5(d-f). The spectral function still shows a series of peaks but their positions no longer approach the conformal tower values E gap (4n + 1) for small µ. The revival oscillations for the transmission amplitudes, T LL and T LR , are now in phase. A related observation is that unlike what we observe for η < 1, the signs of ρ LR peaks do not alternate.
Saddle point equations can be used to further analyze the peak structure in the special case η = 1. For µ J, the MQ model (η = 0) has approximate conformal symmetry, and hence the Green's function in the wormhole phase have the approximate functional form [4] G αα (t) = c |2J sin(w re t/2)| 1/2 , G LR (t) = ic |2J cos(w re t/2)| 1/2 , where w re = 4E gap ∼ µ 2/3 , and the distance between the consecutive peaks w re (t − t ) is 2π. This conformal form of the correlator is obtained using time-reparametrization symmetry of the one sided SYK correlator, G(|t where the time on the left side transforms as t → − cot(w re t/2) and on the right side as t → tan(w re t/2). Note that this results in a π/2 phase difference between T LL (t) = G LL (t > 0) and T LR (t) = G LR (t > 0). For other η, we speculate that the conformal invariance is broken when µ/J (1 − η 2 ) no longer holds, rendering reparametrization invalid.
In Fig. (6) we show, for η slightly less than 1, a transition from out-of-phase to in-phase oscillation as µ increases. At η = 1, we find that the oscillations are in phase for all µ, which can be understood by the following argument. The only non-zero self energy in this case is Σ LL , thus we can relate diagonal Green's function with off-diagonal one by One can check that taking a numerical derivative of T LR with appropriate multiplicative constant gives us T LL with peaks in phase with T LR . For η = 1, although the peaks in T LR do not decay sharply, they suggest that the ground state could be a TFD, as also indicated by our ED results. Whether or not it describes a full-fledged wormhole in the gravity context remains an interesting open question which we leave for future study.

B. Proposed physical realization
We consider a setup proposed by Chew, Essin and Alicea [21] for the physical realization of the Majorana SYK model. The setup is depicted in Fig. 1 and consists of an array of proximitized quantum wires weakly coupled to a disordered quantum dot. Wave functions of Majorana zero modes localized on the left end of wires delocalize into the quantum dot and, in the presence of interactions, are described by the SYK model in Eq. (1). In Ref. 21 very long wires are assumed such that zero modes localized at the right end of wires do not play any role. Here we consider a case with short wires such that χ L j and χ R j are weakly coupled thus furnishing the H LR term required to construct the imbalanced MQ model. One expects that µ ∼ e −2L/ξ cos(k F L), with ξ the effective coherence length [28]. Intuitively one also expects that when the R Majorana modes are not coupled to a quantum dot this setup will realize the maximally imbalanced MQ model. Below we examine the mapping of Ref. 21 to the SYK model more closely and confirm the validity of this intuition in the limit of short wires.
We first consider the non-interacting Hamiltonian H 0 describing the setup in Fig. 1 with the R quantum dot absent. According to Ref. 21 the relevant model belongs to the BDI symmetry class [29,30], where time-reversal is preserved, and is described by the Hamiltonian Here where χ L = χ L 1 , · · · χ L N are the new N left zero modes. The LR hybridization term in this basis becomes Since the dot energy level separations are much larger than the tunneling strength µ, we make a reasonable approximation that the first term that hybridizes dot fermions with the right zero modes can be ignored.
If the matrix O χ,χ was orthogonal, then by defining χ R = O T χ,χ χ R , we would have the same one-to-one tunneling between the primed Majorana modes. This does not apply trivially because the sub-matrix O χ,χ should be first appropriately normalized to be a N dimensional orthogonal matrix, and even then it is only approximately orthogonal. This can be quantified by writing O χ,χ as a sum of its nearest orthogonal matrix O o and a nonorthogonal correction, such that where O o = U V T and Σ = U T O χ,χ V is the singular value decomposition of the sub-matrix. After an orthogonal rotation of the right Majorana zero modes, with µ eff = µ N/(N dot + N ), which implies a small all-to-all random tunneling term between the L and R primed fermions in addition to the desired one-to-one tunneling.
Two questions now naturally arise: (i) how large is the undesirable coupling t, and (ii) will it spoil the MQ physics that we would like to implement. To answer the first question we assume M to be a random matrix chosen from the Gaussian orthogonal ensemble [31]. We then find t to also be a random matrix with zero mean. The variance of the distribution is inferred by evaluating the Frobenius norm of the matrix which is invariant under any orthogonal transformation where {i,j}=1 t 2 ij and 0.2N is the typical value obtained from the random matrix numerical simulation.
The effect of all-to-all random tunneling term on the imbalanced MQ model is tractable within the replica approach and we analyze it in Appendix A. We find that the gap in the energy spectrum survives in the presence of this additional perturbation and only reduces in magnitude, thus indicating that the ground state for this setup is adiabatically connected to the η = 1 ground state found previously in the absence of the perturbation. These considerations lead us to conclude that the setup depicted in Fig. 1 with the R quantum dot completely absent could plausibly be used to realize the maximally imbalanced version of the MQ wormhole model.

IV. DISCUSSION AND OUTLOOK
Our exact diagonalization and large-N saddle-point results demonstrate that the thermofield double state defined through the SYK model eigenstates occurs as the ground state of the imbalanced Maldacena-Qi model. The latter consists of a pair of SYK models with identical interactions up to an overall scale factor. We showed that the gap, the structure of the excitations, and the revival dynamics of two-sided retarded Green's functions all persist as the interaction strengths in two subsystems become different. These results strongly suggest that the imbalanced MQ model remains dual to a traversable wormhole. It would be interesting to find the corresponding Schwarzian action and the reparametrization symmetry when left and right conformal correlators are different. Physical consequences of the time evolution asymmetry of the ground state in the imbalanced case are also worth exploring.
In the maximally imbalanced limit, i.e. when one SYK Hamiltonian completely vanishes, the solutions indicate a mixed picture with some wormhole and some black hole characteristic features. For a small coupling between leftand right-side fermions, the early time dynamics of the retarded Green's function shows SYK power-law decay, and a weak revival oscillations are observed at later times. The LL and LR oscillations are in phase, which is unlike the canonical MQ model where they occur π/2 out of phase, consistent with a gravity interpretation of a particle bouncing back and forth between two ends of the wormhole. The ground state is shown to be gapped and, importantly, our results indicate that it is close to the TFD state. This could be useful since the TFD state can be prepared as a ground state with just a single system emulating the SYK model, thus obviating the difficult challenge of preparing two microscopically identical SYK systems. Finally, we have shown how this maximally imbalanced limit could be physically realized in the previously proposed setup, thus making our analysis relevant for future experiments seeking to produce the elusive TFD state in a solid-state setting.
We conclude by noting that the adiabaticity argument we employed here is completely general and applies to any model where the ground state of two coupled subsystems is close to TFD.

ACKNOWLEDGMENTS
We would like to thank Stephan Plugge andÉtienne Lantagne-Hurtubise for helpful discussions and GoogleX for hosting the workshop "Quantum gravity in the lab" where the basic idea behind this work was originally conceived. This work was supported by NSERC, the Max Planck-UBC-UTokyo Centre for Quantum Materials and the Canada First Research Excellence Fund, Quantum Materials and Future Technologies Program. TH's research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Industry Canada and by the Province of Ontario through the Ministry of Colleges and Universities.
where H is given in Eq. (8). We perform a quenched average over random couplings J ijkl using the replica trick. The variance for the random distribution is σ 2 = 6J 2 /N 3 with integration measure DJ ijkl = i<j<k<l dJ ijkl , where the product is over all independent couplings. The partition function in Eq.(A1) can be written using Gaussian integral (dχ L ia dχ R ia ) and Z 0 the non-interacting part of the action. We made the usual assumption that the solution of interest is diagonal in the replica indices. The averaged partition function is thus Introducing averaged Green's function G αβ (τ, τ ) = 1 N T i χ α i (τ )χ β i (τ ) using self energies Σ αβ (τ, τ ) as Lagrange multipliers, we have following form for the action S, Note that this is exactly the Maldacena-Qi action [11] except for the η-dependent prefactor on the second line. Time-translation symmetry, G(τ, τ ) = G(τ − τ ), can be used to further simplify the action. The quadratic functional integral over the Fermi fields is performed in the frequency space where the non-interacting part of the action becomes diagonal. Using the following convention for Fourier transforms with Matsubara frequencies ω n = (2n + 1)π/β, the effective action becomes The saddle-point equations in Eq. (10) of the main text are obtained by minimizing the action, taking δS/δG = 0, and δS/δΣ = 0.

Large-N equations for disordered tunneling
As discussed in Sec. III, when Majorana zero modes hybridize within the wire, after projecting onto states in the disordered quantum dot, the partition function has an extra all-to-all tunneling. The non-interacting part of the Hamiltonian is modified with a random tunneling part and reads where (µt) 2 ij = s 2 /N (0.2µ 2 eff )/N and µ eff = µ N/(N dot + N ). The t ij are independent Gaussian random variables and we average over them in the same way we did for J ijkl to get get the modified partition function Deriving the new set of saddle point equations we find that self energies get modified as follows, Σ LL = Σ LL + s 2 G RR , and Σ RR = Σ RR + s 2 G LL . The solution of the modified equations is plotted in Fig. 7 and demonstrates that the gap does not close. This suggests that the ground state of the Chew-Essin-Alicea setup with intra-wire Majorana hybridization is adiabatically connected to the TFD state. The first two of Eqs. (10) are transformed to real-time representation by the usual prescription that connects retarded and Matsubara propagators, namely G ret (ω) = G(iω n → ω + iδ) with δ a positive infinitesimal. The equation for the self energy requires more effort as we want to recast it in a form that will be suitable for iterative numerical solution. To this end we employ the procedure originally devised by Banerjee and Altman [32] and subsequently used in various related contexts [8,15,25].
For the sake of simplicity we consider the self energy equation Σ(τ ) = G(τ ) 3 where, compared to Eq. (10), we have temporarily suppressed all prefactors and focused on a single component of the 2 × 2 matrix equation. After Fourier transforming we obtain We express each Green's function in its spectral representation G(iω k ) = ∞ −∞ dωρ(ω)/(iω k − ω), where ρ(ω) = − 1 π ImG ret (ω) is the spectral function and write We next perform the two Matsubara sums using standard techniques. Details of this calculation can be found e.g. in Appendix D of Ref. 8 and the result is We analytically continue this expression, iω n → ω + iδ, to obtain the retarded self energy where P (ω 1 , ω 2 , ω 3 ) denotes the numerator in Eq. (B2). Finally, defining with s = ±, we can rewrite the retarded self energy in a more compact form as Restoring the prefactors and matrix indices leads to the last expression in Eq. (14). The full set of equations (14) can be used to solve for real-time Green's functions in a way that requires only algebraic manipulations and Fourier transforms which can be performed efficiently using the fast Fourier transform algorithm. Let us finally discuss the dependence of the ground state wave function |GS on η. The actual ground state of the imbalanced MQ model can be expressed in terms of the uncoupled SYK basis |GS = nm ψ nm |n ⊗ |m . Deviations of |GS from the TFD state in Eq. (2) can occur in two ways: (i) either the diagonal coefficients ψ nn do not exactly follow the Gibbs distribution e −βEn/2 , or (ii) off-diagonal coefficients become non-zero, i.e. ψ nm = δ nm ψ nn . In the latter case we find that (H L − H R ) |GS = 0 which would result in increasing deviations of the ground state from |TFD as η grows. Our ED calculations reveal that for 2N = 16 off-diagonal terms are exactly zero up to machine precision. Here, the fidelity is less than 1 only for reason (i). For 2N = 32, off-diagonal amplitudes become finite as shown in Fig. 8. The exact TFD state satisfies (H L SYK − H R SYK )|TFD = 0. In Fig. 8, the numerical ground state for η = 0 shows deviations from the TFD state for 2N = 32, indicating that the above condition is only approximately satisfied.
A weaker condition that verifies if the ground state is a TFD in the large-N limit is to check if the ground state expectation value H L SYK − H R SYK = 1 (1+η) H L − 1 (1−η) H R vanishes, as expected on the basis of Eq. (5). To evaluate this expression we take advantage of the fact that the expectation value of the system Hamiltonian can be computed from the Green's function. Specifically, for the MQ model one can show that the following expression holds [25] for α = L, R. Passing to the Matsubara frequency domain we find We next use the spectral representation, G αβ (iω n ) = ∞ −∞ dω ρ αβ (ω) iωn−ω , such that ρ LL = − 1 π Im [G ret LL (ω)] and ρ LR = i π Im [iG ret LR (ω)]. Summing over Matsubara frequencies we have where we used In Fig. 9, we plot a ratio, H L SYK − H R SYK / H L SYK + H R SYK , that should vanish for an exact TFD. Note that for η = 0, this is exactly zero by the symmetry of Green's function where G LL = G RR , but this symmetry argument does not work for non-zero η. If the ground state is close to a TFD we expect the ratio to be close to zero. This is indeed what Fig. 9 shows for all values of µ and η, except for smallest µ in the close vicinity of η = 1, in accord with our finding that the system becomes gapless in this limit.