SYK wormhole formation in real time

We study the real time formation of the ground state of two coupled SYK models. This is a highly entangled state which is close to the thermofield double state and can be viewed as a wormhole. We start from a high temperature state, we let it cool by coupling to a cold bath. We numerically solve for the large N dynamics. Our main result is that the system forms a wormhole by going through a region with negative specific heat, taking time that is independent of N. The dynamics is smooth everywhere and it seems to follow the equilibrium thermodynamic configurations of the microcanonical ensemble. Also we comment on the relation between this coupled SYK model and Jackiw-Teitelboim gravity theory with bulk fields.

1 Introduction and Summary 1

.1 Motivation
The Sachdev-Ye-Kitaev (SYK) model [1,2,3,4] is a strongly interacting but yet solvable model in the large N limit. At low energies, it displays an approximate conformal symmetry. In this region, the model has many features in common with nearly AdS 2 gravity, or Jackiw-Teitelboim (JT) gravity [5,6,7], coupled to matter fields. This is a simple two dimensional theory of gravity which describes some aspects of nearly extremal black holes in various dimensions.
An interesting variant is to consider a pair of identical SYK models coupled through a simple bilinear interaction [8], see also [9,10,11]. The ground state of this model has a gap, but its excitation spectrum also displays an approximate conformal symmetry. Furthermore, this ground state is close to the thermofield double state of two decoupled models. For reasons we explain below, we call the ground state of this coupled model "the SYK wormhole".
A conceptually similar state also arises when one considers two nearly extremal black holes that are relatively close, so that they are coupled. In this case, a traversable wormhole can connect the near extremal throats [12]. This can be effectively modeled by a nearly AdS 2 gravity theory where we have direct interactions between the values of the bulk fields near the two boundaries [13,8]. In other words, thinking of the Penrose diagram of AdS 2 as a strip, we put boundary conditions for the bulk fields that connect the two boundaries. The two boundaries are causally connected through the bulk, so that this spacetime describes a wormhole. This wormhole is the lowest energy configuration of the system and it also displays the approximate SL(2, R) isometries of nearly AdS 2 .
Given that this is a remarkable state, we are interested in knowing whether it is easy to get to it. In other words, if we start out from a general excited state of the coupled model, can we easily get to the ground state by cooling the system down? Or will the system get stuck in some other state? At first sight the answer seems straightforward, if it is the ground state, the system will surely find it if it can shed its excess energy to the bath. On the other hand, from the gravity perspective, the process involves a topology change. Such topology change might happen via a tunneling solution, but it would be exponentially suppressed in N (or the entropy of each separate black hole).

Wormhole formation in SYK
With this motivation in mind, we study this problem for the two coupled SYK models. We start with a relatively high temperature state of the coupled model which looks like two thermal density matrices, one for each SYK factor. Then we couple the system to a bath and study the evolution in real time by solving the large N Schwinger-Dyson equations. We find that the system indeed finds the "SYK wormhole" ground state in a time that is independent of N . In particular, there is no exponential suppression. Notice that the ability to efficiently find this ground state also makes it possible to prepare the thermofield double (TFD) state of the decoupled model, by simply switching off the interaction between the two sides [8], after we have found the ground state.
The approach we used is the following. The large N Dyson-Schwinger equations form a closed system for the two-point function [1,2,3,4]. In the out-of-equilibrium situation that we are considering, these equations are commonly referred to as the Kadanoff-Baym (KB) equations. We couple the system of two interacting SYK models to a large bath and find the real time dynamics using KB equations. The problem of coupling SYK to a bath was recently studied in [14] and we borrow some results from there. Also, the KB equations for a single SYK was recently studied numerically in [15,16]. Here we study this problem solving the dynamical equations at q = 4. The problem has many time scales and due to numerical limitations we could not separate them all by large amounts. However, our numerical results seem to confirm the picture where the system follows the microcanonical equilibrium curve. We now briefly review this equilibrium thermodynamics.

Equilibrium thermodynamics
In the canonical ensemble the system has two phases: the low temperature one corresponding to the ground state, the SYK wormhole, and its excitations; and a higher temperature phase which is closer to two separate thermal SYK systems. The two phases are separated by a first order phase transition. In the large q limit, the black hole phase and the wormhole phase are smoothly connected by a canonically unstable phase with negative specific heat [8]. However, in the microcanonical ensemble, we expect that the system smoothly interpolates between these two phases. In other words, in the microcanonical ensemble we expect no phase transition as we lower the energy.  coupled SYK models with J = 0.5, µ = 0.0053. Blue dots correspond to the "two black holes phase", whereas green dots correspond to the "cold wormhole phase". Red dashed line: curve for the "hot wormhole" phase expected from a low energy analytic analysis. The question mark "?" reminds us that we were not able to find it as a solution of the euclidean DS equations. Figure 1 shows energy vs inverse temperature β for q = 4. We use energy instead of free energy or entropy because we will be dealing mostly with Lorentzian non-equilibrium correlators numerically and it is easier to find the temperature and the energy from them. There are three different regions. At high temperatures T > T 2BH we have the phase we name the "two black holes phase". At low temperatures T < T WH we have the phase we call "cold wormhole" phase, which can be viewed a as a wormhole with few thermal excitations. The two phases overlap, since T 2BH < T WH . In the intermediate temperatures regime T 2BH < T < T WH we also expect a canonically unstable, but microcanonically stable, phase that we call "the hot wormhole phase". As we mentioned above, this phase can be found analytically in the large q limit. It has positive entropy but negative specific heat. However, at finite q we do not known much about this phase, since we have been unable to find it by solving the Euclidean Schwinger Dyson equations. We interpret this failure as resulting from its canonical instability. This is why we put a question mark in Figure 1. We will find evidence for this phase through the real time evolution, since we will find that the temperature goes up as the energy monotonically goes down. There is also analytic evidence from a low energy analysis, as we will review later. The names "hot" or "cold" wormhole refer to how these would feel to an observer who is inside the wormhole, at its center, in a gravity picture for these configurations. As is clear from Figure 1, there are outside temperatures where we can have both a "cold" and "hot" wormhole. When we talk about temperatures in this paper, we are always referring to the physical temperature as seen from the outside. Figure 1 also displays the critical temperature, T c , when two stable phases switch dominance in the canonical ensemble. For T 2BH < T < T c two black hole phase is thermodynamically metastable but is not a global minimum for the free energy. Similarly for T c < T < T WH and the cold wormhole phase. For our problem the microcanonical ensemble is more relevant. Notice that the different "phases" are continuously connected in the microcanonical ensemble, so they are not really sharply separated phases.
For small values of µ, the coupling between the two SYK models, we can make different analytic approximations for the different parts of the curve. For the two black hole region we can start with two separate thermal SYK models and use perturbation theory in µ. In this regime, the left-right correlator G LR is small and of order µ. The gravity picture is that we have two separate hyperbolic disks with a boundary perturbation that connects the bulk fields on the two disks. We find that T 2BH is in the region where this perturbation theory breaks down. For low temperatures the left-right correlator is of order one. We can access this regime by assuming that the system is close to the thermofield double state. The relevant part of the dynamics is captured by the Schwarzian mode. This aspect of the dynamics is the same for the SYK model and the nearly-AdS 2 gravity theory [17,18,19]. This describes both the cold wormhole and hot wormhole phases. In particular, we can see the existence of the hot wormhole phase in this approximation [8]. In particular, the temperature T WH can be found within this approximation. We review this description in Section 2.3.

Gravity picture
One of our motivations was to understand whether a similar wormhole formation process occurs in more general theories of nearly-AdS 2 gravity with matter.
With this goal in mind we will present a gravity picture for the transition we have in the SYK model. We do not know the precise gravity dual of the SYK model. But we consider a nearly-AdS 2 gravity theory that has some of the same features. For questions that mainly involve the Schwarzian mode, the SYK and nearly-AdS 2 answers match precisely [17,18,19]. However, wormhole formation goes slightly beyond this approximation, and we need to incorporate one important feature that is related to the origin of the ground state entropy, or "zero temperature" entropy, S 0 , of the SYK model. If we start from a phase consisting of two thermal states, then the entropy will have a large contribution of size 2S 0 (plus thermal corrections). As the wormhole forms, the system should be able to shed this large entropy into the bath. In gravity this involves topology change, which would naively be exponentially suppressed. On the other hand, as we discussed above, this happens without any such suppression in SYK.
We can reproduce this from a gravity picture as follows. First we view the two coupled systems as a nearly-AdS 2 gravity theory with N bulk fermion fields with Neumann boundary conditions. The two black hole phase consists of two hyperbolic disks with an interaction between the boundary values of the bulk fields. As we lower the temperature, this interaction effectively becomes strong and the theory flows to a new IR fixed point. The new fixed point is simply a theory with different boundary conditions, namely Dirichlet boundary conditions. This change in boundary conditions is similar to the one discussed in [20]. The two disks decouple again but the boundary conditions are effectively changed from Neumann to Dirichlet. Now we use the observation in [21], that e S0 is equal to the ratio of the Neumann vs Dirichlet disk partition functions for N fermions. This means that the effective theory in the IR, with Dirichlet boundary condition has now S ef f 0 = 0. This means that topology change "costs us nothing", and we can easily transition into the wormhole phase. In fact, by a similar argument we can say that the end of the hot wormhole phase also corresponds to the region where the interactions between the two sides of the global AdS 2 strip produce a flow that change the boundary conditions of the fermions from Neumann to Dirichlet.
In summary, we provide a qualitative gravity mechanism for the formation of the SYK wormhole. The purpose of this explanation was to contrast SYK with what we expect in a generic gravity theory. A generic gravity theory can have a number of fields much smaller than S 0 . In this case, the change in boundary conditions would not significantly change S 0 and it would still be difficult to change the topology. For this reason we could not answer the question of whether there is an "easy way" of forming the wormhole for more general gravity theories, such as the case of four dimensional magnetically charged wormholes in the Standard Model [12].
The paper is organized as follows. In Section 2 we review the two coupled SYK model [8]. We describe the perturbative approach at high temperatures, for the "two black hole phase". We also review the Schwarzian description of the low energy dynamics that describes the hot and cold wormhole phases. Section 3 contains our real time analysis of the formation of the wormhole. We set up the coupling to the bath, we write the Kadanoff-Baym equations (the real time Schwinger Dyson equations), and we present the result of a numerical analysis for some particular values of the parameters. In section 4, we discuss differences and similarities between SYK and nearly-AdS 2 (or JT) gravity and also provide the qualitative picture of the transition. Various computational details are discussed in the appendices. 2 The two coupled SYK model and its thermodynamics

Definition and properties of the ground state
Here we review properties of the two coupled SYK models introduced in [8]. The Hamiltonian of the model consists of two SYK terms coupled by an interaction where where the couplings are the same for both factors. They are Gaussian random variables with variance There is also a generalization where we consider a q fermion interaction term, instead of four. The interaction term has the form In the large N limit, µ and J stay fixed. We will mostly consider the case µ/J 1 and also consider temperatures T /J 1. This will be true even for what we call "high" temperatures. As an aside, let us mention that we can couple the two systems by an operator of dimension where F is the fermion number of O ∆ . We will mostly consider the case of (4) which corresponds to ∆ = 1/4 at low energies. However, we will give certain estimates for generic ∆. Like a single SYK, this model, (1), is solvable in the large N limit. We have four types of correlators: G LL , G RR , G LR , G RL , each defined in Euclidean space as Since we are dealing with Majorana fermions we have We have a closed system of Dyson-Schwinger equations for the two point functions [8] where the convolution * is taken along the Euclidean circle, a, b = L, R, and µ ab is given by The system has a Z 4 symmetry ψ L → −ψ R , ψ R → ψ L . Throughout our numerical computation we do not assume that this symmetry in unbroken. We find that it is unbroken, since the (Euclidean) correlators we obtained obey the following relations: A convenient expression for the energy is The interaction term (4) is a relevant perturbation, since for a single SYK model the fermion ψ has dimension ∆ = 1/4. Therefore at relatively high temperatures we expect that we have two weakly coupled SYK models, whereas at low temperatures the system flows into a gapped phase with a gap that scales as [8] E gap ∝ µ 2/3 J 1/3 , for µ J 1, and q = 4, ∆ = 1 4 Moreover the ground state is close to the TFD of the two models with effective (inverse) temperatureβ: The energy of the ground state, relative to the energy of the two decoupled SYK models, scales as And for general ∆, Since there is a gap and the ground state is unique, the entropy is small in the "cold wormhole phase". Whereas in the two black hole phase, we have a big entropy 2S 0 , where S 0 ≈ N × 0.23 is the "zero temperature" entropy of a single SYK model. The transition temperature T c is estimated by For arbitrary temperatures the Dyson-Schwinger equations can be solved numerically by starting from G LR = G RL = 0, G LL = G RR = 1 2 and then using an iteration procedure similar to one described in [3]. After obtaining the solution for some value of temperature, we can use it as a seed for the iteration procedure at higher/lower temperature. Figure 1 shows energy as a function of beta for particular values of parameters J = 0.5, µ = 0.0053.

Perturbation theory at high temperature
Here we use the term "high" temperatures for temperatures for the two black hole branch of the diagram T T 2BH , but still T /J 1. For µ = 0 we have two copies of the conventional SYK correlators [1,2] with b 4 = 1/(4π). Now we turn on a small value of µ (4). If we are at sufficiently high temperature then the coupled system is still in the phase with two separate black holes [8]. Nonetheless, the correlation between left-and right-SYK is not zero. We can try to use the conformal perturbation theory to study the system.
To linear order in µ, only G LR Green functions receive a correction: This integral is computed analytically in Appendix A. We can compare this leading order approximation against exact numerical solution of the Dyson-Schwinger equation for J = 0.5, µ = 0.05 and different βs, see Figure 2.   Also in Appendix A we computed the µ 2 correction to the energy. So that at low temperatures we have the following expression for the energy: where E SYK is the low-temperature result for a single SYK [3] and ∆E is the leading µ 2 correction derived in Appendix A. The comparison between (19) and the numerics is presented on Figure  6(b).
This approximation works better for high temperatures and then deteriorates close to T 2BH , where the phase is supposed to connect with the hot wormhole phase.
Let us find out until what temperature we can trust the perturbation series in µ, (4). The first point is that only even terms are non-zero. The term of order µ 2n contains a 2n-point function on the left and 2n-point function on the right, each of these now computed in a single SYK model. We are only interested in connected correlators for computing the corrections to the free energy. Higher point functions in SYK have two contributions: a purely conformal piece which is independent of β, up to an overall factor of 1/ √ Jβ 4n as in (17), plus contributions from the Schwarzian which are enhanced by an extra factor of βJ. We claim that the Schwarzian contributions are in fact zero, see Appendix B. The reason behind this cancellation is the following.
When we work at large N we are solving the classical equations. The reparametrization mode has a solution that is time translation invariant. The two sides are coupled by convolutions of Green's functions but this translation symmetry remains unbroken. This means that there is no source for higher Fourier components of the reparametrization mode, so that the standard thermal solution continues to be a solution. The integrals over time give β 2n . In total, we have µ 2n β 2n (Jβ) n . So the expansion parameter is So we can trust the above perturbative answer until temperature For ∆ = 1/4 this scales as µ 2 , whereas T c ∝ µ 4/3 is much larger.
We conjecture that the transition temperature T 2BH , when the two black holes phase cease to exist, in fact coincides with T pert , when the perturbation theory in µ breaks down We check this prediction for ∆ = 1/4 against the numerical phase diagram obtained in [8]. See Figures 3,5. One last comment on the leading result (18). In Appendix B, we studied the gravitational dressing of this term, searched for instabilities that would spontaneously break the U (1) time translation symmetry, but did not find any.

Low temperature thermodynamics using the Schwarzian
In this section we review the results of [8] on the Schwarzian description of the wormhole. We will see that the Schwarzian indeed admits a wormhole solution at low temperatures. Moreover, by including the matter contribution to the partition function one is able to see two phases which join at (inverse) temperature β WH . One phase has positive heat capacity and almost zero entropy. This is a cold and mostly empty wormhole. The second phase has higher energy longer throat and negative heat capacity. This is wormhole with extra matter excitations inside. We will see shortly that at small enough µ this approximation predicts β WH with good accuracy.
As we have mentioned before, the ground state of the system is close to the TFD state. Since we have global AdS 2 we have the following G LR correlator in Euclidean signature: where t L , t R are times on left/right.    [8] in Figure 3, using only data points with µ < 0.03. The fit is consistent with the analytical prediction µ ∼ √ T .
The action in the Schwarzian approximation now includes two kinetic terms a and an interaction with α S in (21). The wormhole solution is simply t L = t R = t u, where t is a constant. The effective temperature,β of the TFD state is given bỹ Inserting into the action we get the free energy We should also include the contribution from matter fields in the wormhole throat. If the temperature is low enough we excite only the lightest excitation in the bulk, which is the elementary fermion with mass 1/4. Its contribution to the free energy is: Extremizing the free energy with respect to t , which is the same as extremizing the full effective action, we have the following equation to determine t and correspondinglyβ: This equation has a solution with very small t which we can not trust, because we were assuming that the temperatures are low. For large enough β two additional solutions emerge. One of them correspond to what we call the cold wormhole and the other to the hot wormhole. The latter has negative heat capacity and can be viewed as a wormhole with more excitations in the throat. Figure 6 shows the two branches for two different values of µ/J and their comparison to the numerical solution to the DS equations.
We can calculate T WH from this equation. We simply need to find when these two solutions merge. To keep the discussion general, we consider general ∆, which corresponds to the case when the two sides are coupled through the product of two operators of dimension ∆. The equation for t now reads as: The cold wormhole branch can be approximately found [8] by neglecting the thermal excitations exponent in the above equation, so that t is equal Whereas the unstable branch with excitations can be approximated by neglecting the Schwarzian kinetic term ∼ α S t : Plugging the t from the first solution (34) into the above equation we find the T WH : Of course, both (34) and (35) are good for T T WH . Here we presented just estimates, but it is straightforward to solve (33) numerically, see Figure 6.
We can also compute the energy using b   (19) for the energy. Green curve: wormhole branch of eq. (32). Red curve: unstable branch of eq. (32). The uncertainties, represented as shaded regions, come from the uncertainties in E 0,SYK and α S . We see that in (b) the agreement is very good. However, in (a) the agreement is not so good, but the qualitative form of the curve is similar, if we joint the two end points of the dotted lines.
where E 0,SYK is the ground state energy of a single SYK, see (21).
In Figure 6 we have compared the results computed using eq. (37) with the numerical solution of the DS equation. For the value µ = 0.05, which is the one we will use for the real time numerical computation, the agreement is not very good, but the qualitative form of the curve is similar, see Figure 6(a). This means that that µ is not low enough for an accurate Schwarzian description. Indeed if we lower µ we get very good agreement. See Figure 6(b) for µ = 0.0053.
To summarize, for small µ, we have a hierarchy of temperatures where the rightmost term corresponds to ∆ = 1/4.

Coupling to a bath
In order to study the real-time formation of the wormhole, we need to cool down the system. Before considering real-time dynamics, first we need to understand how to couple our system to a thermal bath. Ideally we want the bath to be a large system in order to avoid back-reaction.
b When computing the derivatives one has to keep in mind that t is a function of β. And use (32).
Generally, we can couple a system's operator O S to a bath operator O B : where F is the fermionic number of O S . If V is small and the bath is large we can study this interaction in the Lindblandian approximation by considering the 1-loop result and assuming that there is no back reaction on the bath, such that we can substitute the product For our problem we have a varying temperature that sets an energy scale for the model. Specially for our numerical analysis, it is convenient to choose an interaction that is scale invariant (at least approximately), so that the effects of coupling to the bath are independent of the temperature. Otherwise the bath might be effectively decoupling in some temperature range and the system would take long to cool down.
A natural model for the bath is another SYK, possibly with larger number of fermions. We can consider the bath to be another single-SYK model with q = q B with large number of fermions M , much larger than the number of fermions in our system N . Recently this problem was studied in [14], we refer to this paper for details. If we denote the bath fermions by χ α , then the coupling can involve an arbitrary number s of system fermions and p bath fermions.
Again, to warm-up, let us first consider the case when the system consists of a single-SYK model. The coupling between the bath and the system has the form: where V is a random Gaussian-distibuted coupling. By choosing its variance appropriatly, the back reaction can be made of order N/M 1. The Euclidean Dyson-Schwinger equations for the system stay the same except for a correction to the self-energy: the self-energy acquires an additional term, where G B is bath two-point function and η is proportional to the variance of V . We can get a marginal interaction with s = 1, p = 3, when the bath consists of q = 4 SYK models (with, say J B = J S ). This is the bath we will use. More precisely, we introduce two separate baths, one for each SYK factor where V andṼ are independent Gaussian-distributed variables. This interaction leaves Σ LR unperturbed, but the other two self-energies have additional terms now: The above equations are written in Euclidean signature. We now turn to Lorentzian equations.

Kadanoff-Baym equations
We now write down the Lorentzian time version of the DS equations. For details see Appendix C.
Here we will discuss the non-equilibrum situation following [3] and [15]. It is convenient to work right away on the Keldysh time contour, see [22] for a comprehensive introduction. We will need Lorentzian time correlators which are not time ordered. This can be achieved by introducing a Keldysh time contour which runs from −∞ to +∞ and then back from +∞ to +∞. First, one introduces a Wightman function where t 1 , t 2 lie on different sides of the contour. This simply means that ψ b (t 2 ) is always located to the right of ψ a (t 1 ) in the correlator, regardless of time-ordering. This is why the Keldysh contour consists of two parts. Note the overall i in front of the correlator. Since we are dealing with simple Majorana fermions the "lesser" function G < ab is given by: Also we will need retarded and advanced Green functions: Dyson-Schwinger equations written on the Keldysh contour are known as Kadanoff-Baym equations, and are useful for non-equilibrium situations. Let us write them down explicitly for a single SYK: where the self-energy Σ > is given by These equations for the complete system of two interacting SYKs and a bath are derived using the path integral in Appendix D.
Remember that the "greater" Green function G > (t − 1 , t + 2 ) has time arguments lying on different sides of the Keldysh contour, this is why we do not have a delta-function on the right hand side of (48). The integral in the right hand side of (48), which involves different Green functions, is just a simple convolution Σ * G along the Keldysh contour [23]. We can show it by writing the anti-time ordered Σ = Σ > − Σ R and the time ordered G = G A + G > . One can easily see that equations (48) are casual.
Let us mention one subtlety. Strictly speaking, if one starts from a thermal state, then the precise Keldysh contour involves imaginary time stip at the end of the lower branch, at t = 0. This time strip has length β and prepares the thermal state. One can bypass this as follows. First we set the coupling to the bath to zero. Then we find the real time correlators at thermal equilibrium by solving the real time equations imposing the appropriate relations between the Green's functions, see e.g. (100). We then evolve the system for some time t β and then we turn on the coupling to the bath. For more details see Appendix E.
For two coupled SYK models one has to be very careful with the µ term. To understand its form on the Keldysh contour we can go back to G, Σ effective action derived in Appendix D We see that the µ contribution to Σ(t 1 , t 2 ) is In this expression t 1 , t 2 can be on either side of the Keldysh contour, this is why Σ does not have an additional index, like >, R, A. Notice that the delta-function δ C is defined on the Keldysh contour as well. It yields non-zero answer if and only if t 1 = t 2 and t 1 , t 2 are on the same side of the contour. Let us compute the contribitution of this term to the convolution Σ * G along the Keldysh contour: The final form of the Kadanoff-Baym equations, including the bath, is: where µ ab is defined in eq. (10) and the self-energy is

Forming the wormhole
Our numerical setup for solving KB equation (53),(54) is presented in Appendix E. We prepare initial Lorentz Green's functions using Lorentz-time Dyson-Schwinger equations described in Appendix C. Our initial Green's functions correspond to two interacting SYK models with nonzero µ at thermal equilibrium. In order to save computational time, the initial temperature is chosen close to (and slightly higher than) the transition temperature T 2BH . We extract the temperature using the Fluctuation-Dissipation Theorem(FDT) for the twopoint functions. Precisely, our numerical setup gives us the G > and G < Green's functions. At thermal equilibrium a certain combination of these Green's functions (eq. (119)) must be equal to tanh(βω/2) in the frequency domain. So we take the numerically obtained G > and G < , make a discrete Fourier transform and fit eq. (119) with the tanh. We refer to Appendix E for details about the precise choice of the Green's functions and the frequency domain for fit.
We can perform another check, this time taking η much smaller, namely η = 0.01. The result is shown on Figure 8. A few comments are in order. As is explained in detail in [14] the marginal system-bath interaction we will renormalize J, making it bigger. This is why expect that the actual transition will happen at higher β. This is indeed what we see. Moreover, the interaction with the bath will shift the ground state energy. To compensate for this we have shifted the energy vs beta curve in Figure 8 to match the final energy.
The red curve in Figure 8 has wild oscillations in temperature near the equilibrium for small η, see also Figures 7(b). The reason is the following. Because of the numerical error there is an additional flux of energy which pushes the system out of equilibrium. From the phase diagram (blue points) we see that the derivative dβ/dE is very large. If η is not big enough, the relaxation time is not small enough to smooth out these fluctuations.
The crucial question is whether we indeed have reached the wormhole phase or not. In principle, we might have ended up in some other phase. To verify that we have reached the wormhole we can make a precise check of the system's energy.
In the real-simulation the final value of the energy is(in units where J = 0.5): where the uncertainty comes from changing the size of the diagonal strip and changing the timestep. Also we can ask ourselves how carefully is the initial state prepared. (a) The energy as a function of time. The initial rise is due to the fact that we are coupling the bath to the system, and this changes the energy [14]. We then see the energy decreasing monotonicallly. (b) The inverse temperature as a function of time. We also see an initial sharp increase due to the coupling of the bath, then we see a decrease. Then a slight increase of the temperature that signals the phase with negative specific heat. Finally the temperature decreases again.
two-point function we know that G aa (0) = −0.5i. However, the iterations of the real-time DS equations have G aa (0) = −0.5002i. The error coming from this is esentially the same as in the above number.
How do we compare this result with the equilibrium phase diagram? In fact, we can solve Euclidean Dyson-Schwinger equation for the coupled system+bath and compare the equilibrium energy. We start from the Euclidean correlators in the wormhole phase, add coupling to the bath and solve the DS equations again. The value of the energy we obtained this way is(again in units where J = 0.5): It is represented as green dot in Figure 8. The uncertainty is estimated by changing the number of discretization points and imposing different cutoffs for the iteration procedure. We see an agreement with (55) within one standard deviation. This suggests that we indeed reached the wormhole.
To understand what happens near the transition we can look at the maximum value of G LR , see Figure 9 (lower part). We can notice that, during the transition through the unstable phase, the imaginary value of the correlator (which is proportional to the anticommutator) rapidly grows, indicating the growth in the information transmission rate.
Another thing we can see is that coupling to a bath generically thwarts the information transmission between the sides; for larger bath coupling η the ratio Im G LR / Re G LR is smaller. Figure 8 shows that the system is more or less following the thermodynamic curve. We see that the temperature and energy are smooth everywhere and the transition goes through a phase with negative heat capacity, where the energy decreases and the temperature increases.
To check whether the system remains thermal at all times we performed an additional check. Using the fluctuation dissipation theorem (FDT) we can find the temperature two ways: from LL correlator and the LR correlator, using (119). The result is shown on Figure 9 (upper part). We see that LL correlator is very close to thermal and the curve shows a clear period of temperature increase. In contrast, the temperature extracted from the LR correlator has big errorbars. This means that the LR correlator has larger deviations from precise thermality, and larger violations to the relation (119). We think that this happens for the following reason. We couple each SYK system to its own bath, so we have two uncorrelated baths. This introduces incoherence to the system, which can be seen on Figure 9: G LR is decreases when we increase η. Moreover we are in the regime of small µ, which means that the two sides interact weakly, so the equilibration time for G LR is much bigger than for G LL . Therefore we expect that the error bars for LR temperature are big because our system-bath coupling is too big. And, if we managed to cool the system down more slowly, then the LR correlators would remain thermal. Unfortunately, with our limited computer resources we could not go below η = 0.01. So we do not have a clear evidence for this interpretation.

Time to form the wormhole
In this subsection we provide analytic estimates for the time it takes to form the wormhole. We will first estimate the time it takes to reach T 2BH and then the time it takes to reach from there to T WH . At this point we basically have a cold wormhole, so we will consider it to be already formed. We could also consider it formed once we reach T 2BH and we start moving on the hot wormhole region.
In order to estimate these times we need an expression for the rate of energy emission into the bath. For a general coupling between a system and a large bath the energy loss rate can be written as (see [14] for more discussion) where factor of 2 comes from having both sides of coupled to a bath. In the above integral, we can replace the system two point function by the thermal one at the instantaneous temperature, assuming that the temperature varies slowly. Moreover, for SYK at low temperatures we can simply use conformal approximations for two-point functions.
If the bath temperature is much lower, than the system's temperature we can approximate bath Green's function by the zero temperature one

Reaching T 2BH
In the two black hole holes phase we assume that G LL is approximated by a single SYK thermal two-point function Plugging these into (57) we see that the answer is determined by dimensional analysis up to an irrelevant numerical coefficient c : Differentiating the energy expression (20) with respect to the time, and solving a simple differential equation for β(t), we find that it grows exponentially So that the time ∆t 2BH needed to go through the two black hole phase and reach T 2BH ∼ µ 2 /J depends only logarithmically on µ:

Reaching T WH
Now let us calculate time ∆t WH which is needed to go through the hot wormhole phase and reach T WH . To this end we will employ some results from the end of Section 2.3 about Schwarzian. This Schwarzian approximation breaks for very hot wormholes with temperature of order T 2BH , but holds for lower temperatures. Suppressing the numerical coefficients, G LL in this regime is given by: where t (t) is determined by the solution of eq. (32). The hot wormhole(unstable branch) is characterized by having t 1, so eq. (32) can be simplified by neglecting Schwarzian kinetic term (the first term in eq. (32)) exp This approximation breaks down near T ∼ T WH , so we further assume that we use this approximation for temperatures which are slightly below or of order T WH . From this equation, up to a logarithmic term in µ, t and β are related by t ∼ 1/β. Using (64), the energy (37) can be written as We see that the heat capacity is negative c We refer to [14] for the numerical coefficient.
We can compare the absolute value of this expression with the heat capacity of a regular SYK This ratio is much bigger than 1 for β 1/(J 1/3 µ 2/3 ) ∼ β WH . So apart from the region close to T WH the hot wormhole has a large negative heat capacity, compared to a single SYK model at the same temperature.
The energy flux can be computed using the expression (63) for G LL . The result is again determined by scale symmetry and it is again proportional to t 2 ∼ 1/β 2 as in (60). However, because of the big negative heat capacity, the time it takes to go through this region is much longer than (62). Solving for β(t) we get In our case we start from β 0 ∼ β 2BH and end with β 1 ∼ β WH . Since β 2BH β WH the overall time length is mostly determined by the region near T 2BH : This timescale is much larger than ∆t WH , (62), which scaled only logarithmically in µ. Moreover, it is mostly determined by the region near T 2BH , which is where the approximation is breaking down. So (69) should only be viewed as an order of magnitude estimate.
Our numerical results for β(t) on Figures 7 and 9 seem to qualitatively support these conclusions. Notice that, as expected, the times are inversely proportional to the coupling to the bath η.

Two coupled black holes in gravity
The low energy description of the SYK model has some features in common with certain two dimensional theories of gravity. In this section, we study a similar problem in a gravitational theory in order to compare to the answers we found above.
We consider a Jackiw-Teitelboim theory of gravity coupled to matter, see [17,18,19] for details. This gravity theory describes a two dimensional black hole with an AdS 2 geometry. The AdS 2 space has a boundary. We consider a system containing two such black hole exteriors and we introduce a coupling for the two dimensional matter fields propagating in the bulk. We assume that have N such matter fields. Let us say that χ is a matter field with a certain mass m in the bulk and quantized with Neumann boundary conditions so that its dimension is ∆, with ∆ < 1/2. We couple their boundary values through a term, for each field, were we imagine that J −1 is related to a cutoff in the radial AdS 2 direction d and u is the physical boundary time.

High temperature phase
We now consider the high temperature phase where in Euclidean space we have two separate disks that are connected through the interaction (70). Concentrating on the matter system, this interaction is easy to analyze because the full matter theory is just quadratic. In principle, we d With the AdS metric ds 2 = (dx 2 + dz 2 )/z 2 , this is the cutoff at z = , and we are defining J = 1/ . also need to consider the effects of gravity, and we will discuss them later. This interaction, (70), leads to the Feynman diagrams in Figure 10(a), which can be easily summed, as we explain below. Since the interaction is relevant, it becomes important at low temperatures. For sufficiently low temperatures, the net effect is to change the boundary conditions for the bulk fermions χ from Neumann to Dirichlet. Namely, at low temperatures we get two decoupled disks with Dirichlet boundary conditions for bulk fermions. We now discuss this more explicitly. Whenμ = 0 we have two separate disks and the matter partition function is just given by Z 2 N , namely the square of the partition function of a fermion with Neumann boundary conditions. Starting from this state we can now sum the diagrams in Figure (10)(a). For each fermion field, we get ] −2∆ as a matrix with indices u 1 , u 2 . We have set β = 2π for simplicity and we will restore it later.
For largeμ we find that the partition function gets an additional factor of the determinant of G ∆ . It turns out that this produces the Dirichlet partition function [21,24] In the last equality we neglected the energy contribution, since we will be focusing on the ground state entropy contributions. The conclusion is that if we start out with two disks with Neumann boundary conditions, after we turn on the relevant perturbation (70), for very low temperatures we get two decoupled disks again but with Dirichlet boundary conditions. This implies that for very low temperatures, the bulk fermion would be dual to an operator of dimension 1 − ∆.
Restoring the factors of β, this transition happens at β 2BH given by For q = 4 this reproduces (24). This is not surprising because we were summing the same type of diagrams. However, in the gravity case these are all the diagrams, so we can study the whole flow. The new IR fixed point simply corresponds to flipping the boundary conditions to Dirichlet. So nothing too dramatic happens in the gravity solution when we go to temperatures lower than the temperature T 2BH in (73), except that the change in the boundary conditions will change the value of the ground state entropy.
In this discussion, we have ignored the dynamics of gravity. In principle, we could wonder whether we should consider non-trivial solutions of the Schwarzian theory. If we assume that the solution is invariant under translations for each disk, then, up to gauge symmetries of the Schwarzian theory, the only solution is the usual one. In Appendix B, we examine whether non-constant Schwarzian modes could lower the action. We find that they do not, at least in the approximation we considered. In our analysis we assumed that gravity is classical, which is correct if φ r /β 1. Here φ r is the JT gravity analog of the coefficient of the Schwarzian, the analog of N α S /J . We have also assumed that we have a relatively low number of bulk quantum fields so that the effects of integrating them out does not significantly change the value of S 0 , the ground state entropy. This is the regime where the gravity theory is simplest. As we will discuss below, the SYK model is different in this respect.

Low temperature phase
At very low temperatures the coupling (70) leads to the formation of a wormhole [8]. This is identical to the small µ coupled SYK model analysis of section 2.3, since the effects of gravity can also be described in terms of the Schwarzian mode.
When we decrease the temperature along the negative specific heat region (the hot wormhole phase), the wormhole is getting longer and longer. Or, more precisely, there is a larger redshift factor between the boundary and the center of the wormhole. Then, the interaction, which is a relevant deformation, becomes stronger. When we considered the problem for the disks, we found that for strong interactions we get an effective change in boundary conditions from Neumann to Dirichlet. Here we expect the same phenomenon when t is becomes where t is the variable in (32), which is proportional to the value of the redshift factor at the center of the wormhole. In other words, t becomes of the order of the temperature T 2BH in (73). We refer to Appendix F for details. At this value of t the wormhole is so long that the approximations used in deriving (33) are no longer valid. Interestingly, due to (35), this happens also at a temperature of the order of T 2BH , which is the temperature where the two disk solutions starts being corrected. This might appear as a coincidence, but it is not. In the hot wormhole phase we find that the temperature sets the value of t and thus the amount of RG flow that the relevant left-right interaction undergoes. Therefore, this interaction becomes relevant at the same place.
This statement can be further verified by checking whether the hot wormhole thermodynamic curve (red curve in figure 6) will join with the two black hole phase (blue curve in figure 6) at T = T 2BH . In the hot wormhole phase the energy is given by (65) and in the two black hole phase by (19). Indeed, the two curves join at T ∼ T 2BH .
In a gravity theory with a relatively low number of fields, we expect that after T 2BH the wormhole phase might not exist any longer.
One conclusion is that, in a general JT gravity theory plus matter, we do not seem to be able to easily join the high temperature phase and the hot wormhole phase. This is mysterious in the gravity theory because it involves a topology change. Of course, the low temperature phase and the hot wormhole phase are connected smoothly at T WH in a region where we can trust the wormhole analysis in the Schwarzian approximation, as discussed near (33).

Comparison with the SYK model
In the SYK model, the addition of the interaction corrects the original diagrams by inserting µ terms in the propagators. If we insert them outside the self energy correction (the Σ bubble), then we get diagrams which are identical to the ones we discussed in gravity, see Figure 10(b).
However, in the SYK model we can also insert µ inside the self energy corrections, these are new diagrams that are not present in the gravity discussion, see Figure (10)(c).
As we mentioned above the region of the phase diagram near temperatures T ∼ T 2BH is different in a generic JT gravity theory plus matter than in SYK. However, we can consider the folloing gravitations model that would look more qualitatively similar to the SYK model.
First we note that the SYK ground state entropy is given by S 0 = N s 0 , where s 0 is given by [2,21] where the first equality follows from the usual G, Σ action at low energies. The second equality was mentioned in (72). This implies that if we want to describe the SYK model in terms of JT gravity, we should think that when the fermions have Dirichlet boundary conditions, the net entropy, or value of the topological terms in the action should be zero, φ 0 = 0. Then the actual value of the ground state entropy of the usual, single boundary SYK model, (75), is simply given by the contribution of changing the boundary condition for the bulk fields from Dirichlet to Neumann [21]. Returning now to the coupled model and starting from the high temperature phase, we see that when we reach the temperature T 2BH we are changing to a Dirichlet boundary condition. This means that the total S 0 now becomes zero, which implies that the topology change is easy. Similarly, if we start from the canonically unstable wormhole phase and approach T 2BH , we also see a change in the boundary conditions so that S 0 again becomes zero and topology change is easy. So we can join the two phases with a change in topology at T 2BH . In this way we can qualitatively understand the transition. We have given evidence that this is a smooth transition in the coupled SYK model. What we are discussing here is just a cartoon for a gravity picture of what is happening.
We also see why SYK is different than a generic JT gravity theory with a smaller number of fields. In such gravity theories the flow from Neuman to Dirichlet would not change S 0 by too much and the topology change remains suppressed. For this reason we have not been able to see a general mechanism for the transition that would also work in more general gravity theories, such as the Standard Model in the presence of magnetically charged black holes as discussed in [12].

Conclusion
In this paper we studied the approach to the ground state of the two coupled SYK models [8]. We first discussed the equilibrium thermodynamics picture. In the microcannonical ensemble we expect a continuous picture with no phase transition. As the energy decreases, the temperature decreases up to a value T 2BH where the system looks like two separate thermal SYK models with a small coupling. At T 2BH this coupling becomes strong and the system transitions to a "hot wormhole" phase with negative specific heat. Now the temperature increases up to T WH and then the wormhole becomes cooler and shorter, and the specific heat becomes positive again. This whole picture can be understood using simple analytic approximations, except for the transition region at T ∼ T 2BH .
We found that the real time evolution looked as if the system is following the above equilibrium phase diagram. Unfortunately, for the parameters we could use in our numerical computation, we could not trust quantitatively the simple analytic approximations. However, these gave a qualitatively correct answer. The system remained near thermal equilibrium as it cooled down, except for some deviations in the G LR correlator, which we think should disappear if we were to cool more slowly.
The conclusion is that, starting with a generic state of the two coupled SYK model, we can find the ground state by coupling the system to the bath and cooling it down. In particular, the system does not get stuck in a metastable state. This provides a feasible way to produce a state close to the TFD.
We also computed the time to form the wormhole. Most of the time is spent near the region with T ∼ T 2BH .
One of our goals was to extract some general lessons for wormhole formation in gravity. Unfortunately, the SYK model seems to be special, and its special features becomes manifest in the ease by which we can connect the two black hole phase with the hot wormhole phase near T ∼ T 2BH . These two phases do not seem to be so easy to connect in more general theories of gravity. We qualitatively explained why topology is simpler in a gravity theory that is similar to the SYK model, but harder in a more general theory of gravity.
Nevertheless we cannot say how hard forming a wormhole would be in a more general theory of gravity, such as the one describing the wormholes in [12]. It seems hard, but maybe there is an "easy" pathway to form it. It would be interesting to answer this question.
We expect that this article would be relevant for efforts that try to do it using the SYK model, see the proposal in [25], for example. Now, let us compute the correction to energy. From the path integral the leading correction to the free energy is We can recover the integral by taking τ → 0 limit in the conformal answer (78) for G LR . Unfortunatelt it produces a logarithmic UV divergence which we cut at τ = 1/J: where c 1 is the cut-off dependent constant. We can not find it from the conformal perturbation theory, because it is an effective low-energy theory with a build-in UV cutoff of order J. From the above expression we read off the energy correction: Also notice that the constant −c 1 − 2 is not simply a correction to the ground state energy. As we mentioned in the main text, perturbation theory in µ breaks down at large β 2BH ∼ J/µ 2 , so we can not take the limit β → ∞ in this expression. This is signalled by the presence of the logarithm. This divergence has IR nature, and it is not caused by using the conformal answer in the integral (80).
To extract c 1 we can compute integral (80) using the numerically obtained G LL . We find that c 1 = 1.66(1). This agrees very well with the actual numerical result for the energy - Figure  6.

B Checking whether Schwarzian fluctuations are stable
Let us start from doing a 1-loop calculation for Schwarzian first. We again assume that we in the phase with two separate black holes. Then the action will involve two Schwarzian terms plus the interaction piece. For a moment we can imagine that instead of a simple interaction µψ L ψ R we have a term involving composite operators i In the perturbation theory the leading contribution is of order µ 2 : Explicitly the action is: where ∆ is the dimension of O L(R) , for ψ L it is ∆ = 1 4 and {f, u} denotes Schwarzian derivative: Finite temperature solution without interaction reads as: We can perturb it by L , R : and to see whether the two black hole system has a perturbative instability at some temperature. For simplicity we can put β = π and expand in Fourier modes: Before doing an actual computation, let us stop and explain why Schwarzian does not contribute classicaly here and at higher loops. By classically we mean that its contribution is suppressed by 1/N . Since we are in a thermal state (87) is the a translation symmetry along the Euclidean time u. This is why after expanding in Fourier modes (90) we will not have terms linear in n . It means that the thermal solution (87) is still a classical solution of Schwarzian equations of motion even with complicated non-local interaction induced by loops. Since we have an overall N in front of the action, integrating out n will lead to a subleading correction.
Let us return to the actual 1-loop calculation. The only subtlety is that one has to be careful with the time ordering, since the denominator involves The result for the marginal deformation ∆ = 1/2: Now we want to switch to Lorentzian time. We define the Wightman function with an extra −i: As is well-known, upon the analytical continuation in time domain, the time-ordered Euclidean two-point function becomes the Wightman function, therefore e : The other DS equation is written in the frequency space, this is why after the analytic continuation it will involve the retarded components: So far we have not used any information about the state we are considering. This information is needed to connect G > and G R . In thermal state we can use Fluctuation-Dissipation Theorem(FDT): An example of how the Wightman's function look is presented on Figure 11.
e The minus sign is subtle: one can recover it either from the effective action (50) on the Keldysh contour or doing a careful analytic continuation through the frequency space as was done in [3] and the equation (99) should be understood as a matrix equation. For the diagonal Green's functions the FDT has the same form: However the imaginary part of the off-diagonal components is skew-symmetric in time, so we have: This system of equations can be solved numerically by the iteration procedure used for a single SYK model [3]. To ensure that the iteration procedure converges to an actual solution we monitor the discrepancies of eqs. (8): and make sure that d ab < 10 −10 . The typical number of discretization points is N points ∼ 2 17 . Before the transition the diagonal Green's function look similar to single SYK ones - Figure  12

D Derivation of the effective action
Let us write down explicitly the total action for the system on the Keldysh contour C. We will suppress the bath action. Bath fermions χ α andχ α , α = 1, . . . , M belong to independend q = 4 SYK systems with coupling J B . We denote their two-point functions by G B : The total action consists of four terms: S tot = S kin + S J + S µ + S bath (106) • S kin is a standard kinetic term for fermions: • S J is SYK interaction: • S µ is Maldacena-Qi interaction term: • Finally S bath is interaction with the bath: As usual, we can integrate out the disorders leading to bi-local expression in terms of ψ, χ and χ. Couplings V andṼ are Gaussing with the variance [14]: The action can be made quadratic in fermions by introducing the largrangian multiplier Σ, which is integrated over along the imaginary axis: Note that we have an overall minus if front of the action. It is important for the equation connecting the self-energies Σ and Green's functions. Integrating out the fermions produces we effective action (50): Tr log (ω − Σ ab (ω)) − ab C dt 1 dt 2 J 2 8 G ab (t 1 , t 2 ) 4 + 1 2 Σ ab (t 1 , t 2 )G ab (t 1 , t 2 ) + Variation of this action with respect to Σ ab and G ab yield the KB equations (53) and (54)

E Numerical method
Now let us describe the numerical method for solving Kadanoff-Baym equations. Numerical solution of Kadanoff-Baym equations for SYK model was described previously in [15,16] and our approach is essentially the same. We will use two-dimensional grid with uniform timestep to discretize (t 1 , t 2 ) plane. The timestep dt should be much smaller than the characteristic time-scales in SYK 1/J, 1/µ. Since µ J, the 1/J constraint is much stricter. We will work with J = 0.5 this is just a choice to fix energy units. Our time steps will be 0.2, 0.1, 0.05. The main numerical limitation comes from the fact that the Green functions have spread ∼ β, so we can not go to very big β, since we will have to use a huge grid. At finite temperature the Green functions decay exponentially, so to greatly speed up the computation we will concentrate on the strip |t 1 − t 2 | cβ max on the (t 1 , t 2 ) plane - Figure 13. We will assume that outside this strip all the Green functions are zero. The constant β max is the maximal β in the problem at hand. In our case β max = β B -bath's beta. One can verify that one c is big enough the result of the computation does not change.
With the computation power avaliable to us, in order to keep the computation time to be of order of dozens of hours, β should be less than 100. This limits us to µ 0.05. For µ = 0.05 the transition beta is ∼ 61 -see Figure 6.
The system-bath coupling η should be much smaller than J 2 = 0.25 so that the system remain thermal. We will use η in the range 0.01 − 0.04. Moreover to avoid large gradients we will switch on the coupling linearly, with the switch-on time T switch = 20.
Initial Green function is found by numerically solving Lorentz-time equilibrium Dyson-Schwinger equation as described in Appendix C. The bath Green function is obtained in a similar fashion. In order to avoid large speads in the bath's Green function, β B will be in the range 70 − 100.
To compute the integral in KB equations we will use trapezoid method and for the time propagation we use predictor-corrector scheme. The same techniques have been used in [15,16]. For KB equations one has to be careful with propagating the Green function along the diagonal. Fortunately, for Majorana fermions there is a simple relation: However, for the Green function obtained by numerically solving the DS equation the diagonal value is not exactly −i/2, so on a discrete lattice we will just propagate this value: For the mixed G LR we do not have a simple relation like that. So we will use again the predictorcorrector scheme. The value on the diagonal can be found by either propagating along t 1 or t 2 . We will take the avarge of these results. Let us define the "corner slice" Green functions as G > T,ab (x): G > T,ab (x) = θ(x)G > ab (T − x, T ) + θ(−x)G > ab (T, T + x) Then the energy at time t = T can be computed analogously to eq. (12): Computing the time derivatives using the KB equtions one arrives at the following integral form: We expect that they should match for frequencies much less than the inverse discretization timestep: n/L 1/dt. In on this graph this UV cutoff is at βn max /L ∼ 400. The region used for the fit, β|n|/L ≤ 10, is within the dashed lines. Now let us return to our problem with the µ term. So now we have Neumann fermions plus the interaction term µψ L ψ R coupling the two boundaries. The partition function can be easily found: where the matrixμ is given by eq. (10). For large µ we have To conclude that Z µ 1 ∝ Z D we need the matrix relation In fact this relation coincides with the conformal(i.e. neglecting the time derivative) Dyson-Schwinger equation (8).