Replica wormhole and information retrieval in the SYK model coupled to Majorana chains

Motivated by recent studies of the information paradox in (1+1)-D anti-de Sitter spacetime with a bath described by a (1+1)-D conformal field theory, we study the dynamics of second R\'{e}nyi entropy of the Sachdev-Ye-Kitaev (SYK) model ($\chi$) coupled to a Majorana chain bath ($\psi$). The system is prepared in the thermofield double (TFD) state and then evolved by $H_L+H_R$. For small system-bath coupling, we find that the second R\'{e}nyi entropy $S^{(2)}_{\chi_L, \chi_R}$ of the SYK model undergoes a first order transition during the evolution. In the sense of holographic duality, the long-time solution corresponds to a"replica wormhole". The transition time corresponds to the Page time of a black hole coupled to a thermal bath. We further study the information scrambling and retrieval by introducing a classical control bit, which controls whether or not we add a perturbation in the SYK system. The mutual information between the bath and the control bit shows a positive jump at the Page time, indicating that the entanglement wedge of the bath includes an island in the holographic bulk.


Introduction
The black hole information paradoxes refer to various kinds of obstruction in combining black hole gravitational physics and quantum mechanics. As a well known example, Hawking's calculation [1] for the entropy of the radiation from a pure state evaporating black hole leads to a monotonically growing result, which is inconsistent with an unitary evolution where the late time entropy is expected to follow a Page curve [2].
To compute entropy in holographic systems, a powerful tool is provided by the Ryu-Takayanagi (RT) formula [3][4][5]. The formula was first proposed for stationary asymptotic anti-de-Sitter (AdS) spacetime, and has been subsequently generalized to time-dependent cases, known as the Hubeny-Rangamani-Ryu-Takayanagi (HRRT) formula [6]. After taking the contributions from the bulk quantum fields into account [7,8], the general proposal states that the von Neumann entropy of a boundary region A is determined by finding all extremums of the generalized entropy: among all possible bulk surfaces γ A that are homologous to A and then look for the minimal one. In the formula, S bulk is the entropy of the bulk quantum fields in the region bounded by the quantum extremal surface γ A and the boundary.
Recently, new insights on the information paradox have been brought by discovering a new quantum extremal surface for an evaporating AdS black hole [9,10]. As emphasized in [11,12], when applying eq. (1) to calculate the entropy of the radiation, one must consider possible solutions involving entanglement "islands" in the bulk. In most of these papers, the set-up is to allow the black hole evaporating into a bath by gluing the AdS boundary with an auxiliary spacetime with no gravitational degree of freedom. Before the Page time, the quantum extremal surface of the radiation in the bath is trivial, and the entropy of the radiation agrees with Hawking's field theory calculation. In contrast, after the Page time the quantum extremal surface becomes nontrivial and bounds an isolated island in the bulk. The first order transition between these two solutions gives the Page curve. In [13,14], these new solutions are explained as coming from the replica wormhole solutions in the gravitational path integral derivation of the entropy, and the transition becomes a cross-over if one sums over all geometries in the model of [14]. See also recent discussions on the information paradox, islands and replica wormholes in [15][16][17][18][19][20][21][22][23]. Since the quantum extremal surface of the radiation bounds an island in the bulk, the entanglement wedge reconstruction proposal [24] implies that the information inside the island, which covers part of the black hole interior, should become accessible to the bath after the Page time [25]. Concrete ways to recover the information inside the island have been proposed via the Petz map [14] or the modular flow [22].
In this paper, we study the black hole evaporation problem by considering a Sachdev-Ye-Kitaev (SYK) model coupled with a (1 + 1)-dimensional free fermion bath. The SYK model [26,27] is a (0 + 1)-dimensional strongly correlated fermion model with emergent nearly conformal dynamics at low temperature. The low energy dynamics of the SYK model has a holographic dual theory which is the Jackiw-Teitelboim gravity in AdS 2 [26,27]. Previously, the physics of a SYK model coupled to a large SYK bath has been studied in [28][29][30], while the entropy dynamics of the SYK model have been studied in [14,31] using coupled SYK model with equal number of modes. A new saddle point solution after the Page time, which corresponds to the replica wormhole, has been discovered in the micro-canonical ensemble in [14]. In this work, we instead model the bath by free Majorana chains. Having the simple free dynamics in the bath is helpful in simplifying the problem, and is also closer to the setups with AdS black holes coupled with non-gravitational flat space bath.
For simplicity we consider a thermofield double state of the SYK model coupled to free fermion bath, which then contains two SYK models (left and right) and two baths (also left and right). Denoting the Hamiltonian of each side as H L and H R respectively, the thermofield double state is invariant under the time evolution of H L − H R . We consider the time evolution by H L + H R which changes the state and leads to increase of the entanglement between the SYK sites and the baths. This corresponds to a setup with a two-sided eternal black hole in equilibrium with two flat baths, as has been discussed in Ref. [12]. The set-up is discussed in detail in section 2. Using the Schwinger-Dyson equations in the large N limit, we study the second Rényi entropy of the union of the two baths in the time-evolved thermofield double state. Equivalently, this can be expressed in terms of the correlation function of two twist operators. As explained in section 3, for small system-bath coupling, numerically we find that the second Rényi entropy S (2) χ L ,χ R shows a firstorder transition. The short-time saddle can be studied using perturbation theory (section 4) and the long-time solution can be explained by twist operator factorization (section 5). By introducing a classical control bit, we could ask whether the information thrown into the SYK system can be extracted from the bath by looking at the mutual information between the bath and the control bit. This is discussed in section 6. We find that the mutual information has a jump at the Page time, which signals that there is an "island" outside the horizon in the gravity picture. Finally, in section 7, we discuss the holographic picture for our model, which is different from the previous gravity models where matter fields are conformal in the joint spacetime.

Set-up
The SYK model [32,33] describes N Majorana fermion modes χ i labeled by i = 1, 2, ..., N with random interaction. We consider coupling each mode χ i to an individual (1+1)-d free Majorana chain ψ i (x) with constant hopping Λ/2 and periodic boundary condition. These Majorana chains serve as a thermal bath for the SYK system. In most of the following discussion, we will refer to χ as the system (or the black hole in the gravity analogy), while ψ as the bath. The Hamiltonian of the coupled system is the following: where x ∈ Z labels different sites in the Majorana chain, and the interaction term H int will be specified later. We choose the convention that {χ i , χ j } = δ ij and {ψ i (x), ψ j (y)} = δ ij δ xy . J ijkl are Gaussian random variables with the mean and variance: One can perform a Fourier transform on the Majorana chain: where N L is the total number of sites. This gives the summation is over a half of the first Brillouin zone. Here a is the lattice spacing of the Majorana chain. In the low energy limit, the bath contains a left-moving Majorana mode with k ∼ π/a and a right-moving Majorana mode with k ∼ 0. For simplicity, we choose a = 1/Λ. As a result, in the continuum limit Λ → ∞, we have k ∼ ±k near the gapless points. The central charge c of the conformal field theory for the bath in the continuum limit, which contains N copies of the Majorana chains, is then N/2. In this work, we choose the interaction term H int as a hopping of the SYK fermion to the center site x = 0 of the Majorana chain: Here we have introduced a factor of corresponds to the continuous fermion operator in the limit Λ → ∞, with the Dirac δ-function anti-commutator {η i (x), η j (y)} = δ ij δ(x − y).
The set up of the problem is as follows. We first introduce two copies of the coupled system -Left: χ L , ψ L and Right: χ R , ψ R , and prepare them in a thermofield doubled (TFD) state [34] with inverse temperature β. When we only look at the left or right system, it is in a thermal density matrix with Hamiltonian given by (2), while the whole system is in a pure state. The definition of the thermofield double state in this model is not unique. Without losing generality, we make the following explicit choice. We begin by constructing the state |I χ L ,χ R , which satisfies |I χ L ,χ R is a maximally entangled state between the χ L system and χ R system. Similarly, we construct a maximally entangled state |I ψ L ,ψ R between the ψ L system and ψ R system with spatial locality, which satisfies The thermofield double state is then given by where H L and H R are the Hamiltonian (2) defined on the left and right system. After we have the state |T F D , we evolve the system in time using H L + H R . One important property of the thermofield double state is that it is annihlated by Figure 1: The graphical representations of (a) the |T F D(t) state and (b) the reduced density operator ρ χ L ,χ R (t).
H L − H R , and thus we can pull all the evolution on the right system onto the left system, and write the time-evolved TFD state as: The |T F D(t) state can be represented graphically as in fig. 1 (a). The inner/outer line represents the χ/ψ system. We have suppressed the extra spatial dimension for ψ. The half circle with length β/2 corresponds to the Euclidean preparation for the |T F D state, followed by a real time evolution of 2t represented by the horizontal lines. The dotted lines between χ and ψ denote the interaction in the system. In this paper, we will focus on calculating the second Rényi entropy of subsystem χ L ∪ χ R . This is in analogy of calculating the entropy for the black holes in [12]. Due to the inherent unitarity here, this is the same as calculating the second Rényi entropy of the baths, since we started from a pure state. This is different from the gravity story, where a priori one cannot assume unitarity when trying to address the information paradox.
The reduced density matrix of subsystem χ L ∪ χ R is given by: We draw the graphical representation of the density matrix in fig. 1(b). We take two copies of the state in fig. 1(a), one for the ket and one for the bra, and then trace out the ψ L ∪ ψ R system (denoted by the dashed lines in Fig. 1(b)). The n-th Rényi entropy of a density matrix ρ is defined as Specifically, we are interested in the second Rényi entropy of the density matrix ρ χ L ,χ R (t), given by Equivalently, we could express the right hand side as the expectation of twist operators on two copies of the coupled system on state |T F D(t) ⊗ |T F D(t) : Here the twist operator T L/R operates on the two copies of the χ L/R system by swapping their states: We can formulate the calculation of the second Rényi entropy in terms of a pathintegral over a replicated contour C with twisted boundary conditions. The contour C is shown in fig. 2(a), where we take two copies of the density matrix in fig. 1(b), and join the open ends of the χ systems in a twisted way (denoted by the dashed lines). In this figure, we've also marked how we parametrize the contour using a real parameter s: s ∈ [0, β + 4t) covers the upper part of the contour in the clockwise direction, while s ∈ (β + 4t, 2β + 8t] covers the lower part of the contour, also in the clockwise direction. The parameterization will be needed later in the presentation of our numerical results. Another equivalent way to picture the contour C is in fig.  2(b), where it makes clear that the replica contour has the topology of four circles.
The path integral has the following form (17) Because we are using a single real parameter s to label the contour, we need to introduce an extra factor f (s) to account for whether we are doing imaginary time evolution (f (s) = 1), forward real time evolution (f (s) = i) or backward real time evolution (f (s) = −i) 1 . The expression (16) applies for a single realization of the SYK Hamiltonian. However, in order to apply the standard large N technique of the SYK model, one has to average over the disorder coupling, and approximate S (2) χ L ,χ R (t) by the disorder-averaged value: 2 The second approximation comes from the assumption that the dominant saddle point remains replica diagonal. It should be noted that there are two different kinds of replica discussed here. In computing (tr(ρ 2 )) k , there are 2k replica labeled by s = 1, 2, α = 1, 2, ..., k. We assume the dominant saddle point is diagonal in α, so that (tr(ρ 2 )) k (tr(ρ 2 )) k . In general, the solution is off-diagonal in s which labels the two replica we discussed above in computing tr(ρ 2 ).
After the standard procedure of introducing the bilocalG,Σ fields, and integrating out the fermion fields χ and ψ(x), one arrives at where we've introduced the factor F (s, if both s, s lie on the same contour in fig. 2(b) and otherwise zero.
is the Green's function for the bath fermion ψ, without coupling to the SYK system. Since χ only couples to bath ψ at x = 0, only g(s, s ) ≡ G 0,ψ (0, s, s )Λ appears in the second line of (21). For completeness, the explicit expression for g(s, s ) is given in Appendix A. We calculate the large N leading order result of S (2) χ L ,χ R by doing saddle point approximations to the above path integral. The saddle point equations are This set of equations can be solved using iteration numerically. Generally, there could be several different solutions to the saddle point equation and the dominating solution in the large N limit is determined by comparing the action. The on-shell action I C can be written as 1 More explicitly, by the parametrization in fig. 2(a), f (s) is defined as 2 There have been many recent discussions in gravity about the role of disorder/ensemble average and its relation with replica wormholes, see [14,23,35] for examples. The "replica wormhole" solution that we will discuss below does not rely on the disorder average, since different replicas are directly coupled together as in fig. 2(b).
Here the log det term for ψ does not depend on the saddle point solution G χ and cancels with the normalization Z(β). Here we take t/β = 6 as an example. Here and in latter figures we have removed tiny matrix elements |G(s, s )| < 10 −3 in the plot to make the plot clearer.
1. For the large coupling case V 2 /J = 0.25, as shown in Fig. 3(a), the entropy is a smooth function of time.
2. For small coupling V 2 /J = 0.05 case in (b), there is a first-order transition between two different saddle point solutions. The entropy initially grows almost linearly in time, and then switches to be almost time-independent, governed by a different saddle point. The two saddle points coexist for a finite time interval, which can be reached by choosing different initial conditions for the iteration. The transition time is the analog of Page time in evaporating black holes [12].
We will study the solutions in more detail in latter sections. Here we just mention several properties of the solutions: Firstly, the short-time solution for the small V 2 /J case is almost replica diagonal: the correlation between two SYK Majorana operators on different solid contours in Fig. 2(b) is small. On the other hand, the long-time saddle is highly nondiagonal. The anti-diagonal peak for the long-time saddle is from the approximate cancellation between the forward and backward evolution for two halves of the χ contour interacting with the same ψ system. There is a change of the paring between the forward and the backward evolution, which has been found in [36] for a coupled SYK system. Note that the forward evolution and the backward evolution on the same solid contour of Fig. 2(b) could not cancel exactly since they interact with different ψ systems.
Secondly, for large V 2 /J, the off-diagonal terms are of the same order as the diagonal terms. The short-time solution described above is smoothly connected to the long-time solution. This is similar to the equilibrium problem of coupled SYK model [36], where for a small coupling there is a first order transition, while for a large coupling the free energy is smooth.
Thirdly, the entropy for the long-time solution is almost constant and close to 2S th is the second Rényi entropy between χ and ψ in a thermal density matrix with the Hamiltonian in (2).

The short-time solution
In this subsection, we discuss the analytic calculation of the Rényi entropy. Without turning on interaction between the system and the bath in the real-time evolution, the entropy S (2) is time-independent. As a result, focusing on the time dependence of S (2) , we perform a perturbative calculation in V 2 starting from the V = 0 replica diagonal solution for short-time. The calculation is similar to [31], where one first calculates the action in Euclidean time and then continue the result to Lorentzian time.
The conformal limit solution for the χ system is Similarly, at the low-temperature limit Λβ → ∞, we expect G 0,ψ (0, τ 1 , τ 2 ) to be dominated by low energy modes. We can use the conformal two-point function with scaling dimension ∆ = 1/2 for the bath: To calculate the change of entropy, we rewrite the n-replica action as where is a small UV regulator. We take τ 0 = β 2 − 2it in the end. The details are provided in Appendix B. For time t β, the entropy shows a linear growing behavior: Here in the slope of the linear growth, there is an additional factor √ β compared to the results in [14,31]. This is due to the fact that the coupling V is relevant with scaling dimension 1/4, whose effect becomes larger when we lower the temperature. In Fig. 4, we compare the analytic formula (53) with the numerics for V 2 /J = 0.005, βJ = 4 and Λ = 5J. There we have chosen the constant piece in (28) to match the t = 0 numerical result.
The linear entropy growth in (28) can not persist to time much long than t ∼ √ βJ V 2 , since the entropy cannot exceed the thermal entropy at the same energy density. At this time scale, a first order transition occurs for small V 2 /J. Note that here the Page time is finite at N → ∞, because the central charge of the bath is proportional to N and there is an extensive entropy flow between the bath and the system. This is also consistent with the gravity analysis in [12] with c ∝ N .

The long-time solution
In this section, we discuss the structure of the long-time solution and give arguments for the long-time saturation value.
To gain some insight of the problem, in fig. 5(a), we show how the long-time solution of G looks like (it is similar to the one shown in fig. 3(d), but the time is longer, and we set the small elements to zero). In fig. 5(b), we show the solution with the same parameters, but on the contour without inserting any twist operators and thus being two disconnected circles. The important observation is that the solution of (a) is very similar to (b), with the only significant differences locating at the places near the twist operators.
What this tells us is that for the long-time solution, the backreaction of the twist operators is only local. Thus we expect the correlation function of the twist operators to be approximately factorized, i.e.
The solutions with only one insertion of twist operator are shown in fig. 6 (a)(b). A more careful argument for the relation between local back-reaction of twist operators on the on-shell solution and the factorization in (29) is given in Appendix C.
For the computation of each one-point function T L/R , the other side of the system (R/L) is traced out within replicas. As a result, the one-point function of the twist operator computes the second Rényi entropy between one copy of the system and the bath in a thermal state: one prepares the system χ and bath ψ into a thermal ensemble with density matrix ρ th = exp(−βH)/Z, and then compute the second Rényi entropy of the system χ, denoted by S (2) th . Note that for the one-point function in a thermal ensemble, the time evolution is trivial since the density matrix commutes with the unitary evolution. Consequently, by the above argument, the saturation value of the second Rényi entropy should be S th . This is also consistent with the gravity calculation in [12]. In fig. 7, we see that this estimation agrees quite well with the numerical results.

Information retrieval from the bath
By the gravity picture in [9][10][11][12][13][14], after the Page time, there is an island in the bulk which belongs to the entanglement wedge of the bath. This means that if one throws in a particle after the Page time and wait for a scrambling time, one should be able to tell the difference between the state with and without the particle by only accessing the state of the bath. In other words, the information carried by an infalling particle can be retrieved from the bath after a scrambling time. In comparison, a particle thrown in before Page time can only be retrieved after the Page time. From the quantum mechanical point of view, this suggests that an information initially scrambled in the SYK system should emerge in the bath after the Page time. In this section, we will study this phenomena in our model.
In the quantum mechanical system, throwing a particle into the black hole cor- χ L ,χ R after the Page time (the dots) and twice the second Rényi entropy S (2) th between χ and ψ in the thermal density matrix of a single system (the dashed lines). responds to adding a perturbation δU to the evolution of the TFD. For simplicity, we choose the perturbation to be δU = √ 2χ 1 , applied at time t 0 . We then study its effect on the state of the bath at a later time t ≥ t 0 . The set up is illustrated in fig.  8.
We denote the bath reduced density matrix at time t in the perturbed case as ρ B,1 (t), and the unperturbed one as ρ B,2 (t). Obviously, ρ B,1 (t 0 ) = ρ B,2 (t 0 ) since tracing over the system cancels δU with δU † . For t > t 0 the effect is generically nontrivial. The informational retrieval depends on the distinguishibility of ρ B,1 and ρ B,2 .
We would like to define an appropriate quantum information measure for the distinguishability of the two density operators. For this purpose, we introduce an ancilla qubit A with two internal states |0 and |1 . We initialize A in a maximally mixed state ρ A = 1 2 I 2 , where I 2 is the identity matrix with dimension 2. We then perform a classically control operation: if the A system is in the state |0 (|1 ), we prepare a system without (with) perturbation δU , respectively. Tracing out the system χ, the reduced density matrix for the ancilla qubit A and the bath is given by We can then compute the second Rényi mutual information I A,B (t) between the bath and the control bit A, which quantifies whether we are able to reconstruct the Figure 8: We prepare two density matrices ρ B,1 and ρ B,2 , where for ρ B,1 we inserted a fermion operator χ i at time t 0 (represented by the red crosses). In the figure, the denominators are the proper normalization factors, while the factor 1/2 comes from perturbation by measuring the bath. 3 The relevant entropy quantities are given by: exp −S As a result, we have A,B (t) = − log Right after the insertion of the Majorana operator, we have ρ B,1 (t 0 ) = ρ B,2 (t 0 ) and I (2) A,B = 0. As a result, the bath is not affected by the perturbation. On the contrary, from the gravity picture, we expect that in the long time limit, density matrix with and without the perturbation becomes orthogonal tr (ρ B,1 (t)ρ B,2 (t)) = 0, which leads to I A,B = log(2).
To study (33), we could further express both tr (ρ B,1 (t)ρ B,2 (t))/tr (ρ B,2 (t) 2 ) and tr (ρ B,1 (t) 2 )/tr (ρ B,2 (t) 2 ) in terms of the Green's functions of Majorana fermions on the contour in fig. 2. As an example (see fig. 9), we have: Note that, although the mutual information is more physical, this overlap, which is equivalent to the Green's function, also measures the difference between the density matrices.  Similarly, tr (ρ B,1 (t) 2 ) /tr (ρ B,2 (t) 2 ) can be expressed as a four-point function, the leading order contribution of which is given by the factorized result: From the discussion in the previous sections, we learned that across the phase transition at the Page time, the Green's function changes discontinuously. In fig.  10, we show an example of how the the mutual information and the overlap looks like if we add the pertutbation at t 0 = β. It should be noted that the Rényi mutual information is nonzero even before the Page time, which is because the particle created by χ 1 (t 0 ) has a finite probability to directly hop to the bath. In other words, χ 1 (t 0 ) creates a particle that is not entirely infalling, but has a finite probability of going out into the bath. The contribution of this direct coupling to the Rényi mutual information is proportional to V 2 .
At the Page time, both the mutual information and the overlap are discontinuous. The jumps correspond to the fact that after the Page time, the information carried by the infalling particle is encoded in the bath in a nonlocal way. If we throw in a particle before the Page time, most information will only be retrievable after Page time. In comparison, if we throw in the information after the transition, we only need to wait for a scrambling time t sc ∼ β, as shown in fig. 11. Note that because the bath has a central charge N , the information retrieval time is order 1 rather than order log N . In the bulk picture [12], this corresponds to an island that is finite distance outside the horizon. In Fig. 10, we have also plotted the perturbative result for the Green's function in (34). Focusing on a single replica, the change of the Green's function comes from the effective action (27). We copy it here for convenience: To the V 2 /J order, we have the change of the Green's function being: Note that the connected part of four-point function is 1/N , which cancels the factor N in front. Instead of using the imaginary time and applying an analytical continuation at the end, here we directly include the real-time contour in the calculation. We set τ 1 = β/2 + /2 + it 0 , τ 2 = β/2 − /2 + it 0 , τ 3 = β/2 + /2 + it 3 and τ 4 = β/2 − /2 + it 4 . The main time dependence is from: (38) Figure 12: Bulk dual interpretation of the boundary perturbation of adding a χ i fermion. We sketch the Penrose diagram of an AdS 2 eternal black hole coupled with the flat space bath, discussed in Ref. [12]. The blue solid curve indicates the worldline of the boundary of the island. (a) When backreaction is neglected, this perturbation creates a superposition of infalling fermion mode (blue arrow) and outgoing fermion mode (red arrow). The blue diamond is contained in the entanglement wedge of the bath at the Page time. A particle created earlier than the Page time can be retrieved from the bath soon after the Page time, while a particle created after Page time can be retrieved after a finite scrambling time. (b) Considering the backreaction, the infalling fermion induces a change of the boundary location, described by the Schwarzian action in low energy. This leads to a decreasing overlap between the perturbed and unperturbed states of the bath even before the Page time (see text).
Here we use the four-point function F c from the Schwarizan theory (for β = 2π) [27,33]: We evaluate the integral in eq. (38) numerically since there is no closed analytic form. We choose the cutoff = c 0 2πJ and adjust c 0 to match with the numerical result. For c 0 = 4, the result is shown in fig. 10 (a) as the black dashed line, which suggests that the estimation works reasonably well. In the long time limit t β, the first order perturbative result is linear in time: with θ = 2π β = c 0 βJ . The linear decrease comes from the fact that for δG(τ 1 , τ 2 ) with τ 1 = β/2 + /2 + it 0 , τ 2 = β/2 − /2 + it 0 , the four-point function is finite even if t 3 ∼ t 4 t 0 . Physically, the reason is that the infalling particle carries an SL(2, R) momentum, and the boundary will gain an opposite momentum, required by the overall SL(2, R) symmetry. Consequently, the backreaction induced by the infalling particle does not decay with time.
For comparison, we also study the result if we approximate the four point function by only the disconnected part: In the bulk interpretation, this corresponds to neglecting the backreaction and considering a free fermion problem. In this case, it is easy to see that F c → 0 in the limit t 3 ∼ t 4 t 0 . The free fermion result is also plotted in Fig. 10 (b) by the red dashed line. Instead of linear t dependence, the overlap saturates to a finite value (until the Page time), corresponding to a finite probability of the initial particle moving outwards. Comparing the two approximations, we see that the change of the overlap before Page time is mainly due to the backreaction.
Comparing the Rényi mutual information and the overlap, we see that the mutual information changes much slower for short time. This is because in the short time limit, V 2 /J contribution for the mutual information vanishes: The decrease of I A,B then comes from higher order corrections including correlation between different contours.

Conclusion and Discussion
In conclusion, we have studied the SYK model coupled to a free Majorana fermion bath, as a toy model to investigate the physics of black hole evaporation. We studied the time evolution of the two-point function and the second Rényi entropy of the thermofield double state of this coupled system. For low coupling with the bath, we found a first order transition in the second Rényi entropy, which corresponds to the formation of a "replica wormhole", similar to the results in the section 5 of Ref. [14]. We also studied the information retrieval from the black hole by creating a single fermion on the boundary. By comparing the perturbed and unperturbed reduced density operators of the thermal bath, we see a combination of two kinds of effects. Before the Page time, the bath already knows partially about the perturbation to the black hole system, because the boundary fermion has a finite chance to directly leak to the bath, and also because the backreaction of the infalling fermion. The latter effect makes dominant contribution. At the Page time, the information available to the bath about the perturbation has a finite jump, which is consistent with the expectation that the black hole almost saturates to its maximal entropy state after the Page time and therefore almost cannot preserve any information about the perturbation.
With this concrete model, there are many open questions. Although we've shown that the information is in principle retrievable through quantum information argument, we've not provided an explicit construction of the form of the bulk fermion operator. It will be interesting to study the bulk fermion operators more explicitly, especially the fermions behind the horizon. If the bulk fermion operators can be identified, it may become possible to more explicitly investigate the black hole information paradox such as the firewall paradox [37]. It would also be nice if one can see explicitly how the proposals for recovering operators in the island [14,22] work in this set-up. Another question about the information retrieval is whether the backreaction effect we observed for two-replica calculation should vanish if we take the von Neumann limit. For example if we compute the relative entropy S(ρ B1 |ρ B2 ), will the backreaction effect still be significant?
ψ 0 ψ 0 We are mainly interested in the low-energy modes. Consequently, we further make the approximation by using a linear dispersion k ≈ k with a cutoff of the order of Λ. The integral can then be carried out explicitly which gives Here B z (a, b) is the incomplete beta function defined as and the factor of 2 comes from the summation over the left-moving and right-moving modes. At low-temperature limit Λβ → ∞, using and (49), we recover the conformal Green's function as expected.
Appendix B Perturbative analysis in the short time limit In this appendix, we give details for the calculation of the short-time action (27). Defining θ i = 2π β τ i and using the explicit formula (24) (25) for the Green's function, we find with θ = 2π /β. The integral can be carried out explicitly. After continuation to real time by τ 0 = β 2 − 2it and taking → 0, we get (53) Here E (φ| m) is the elliptic integral of the second kind and p F q (a 1 , ...a p ; b 1 , ...b q ; z) is the generalized hypergeometric function. Since we are mainly interested in the time dependence, we do not give the explicit formula of the constant term. Taking the leading order contribution with t β, we arrive the result quoted in the main text:

Appendix C Factorization of the two-point function of twist operators
In this appendix we give more detailed argument for the factorization of the correlation function of the twist operators. In numerics, we observed that the backreactions of the twist operators are local. This means that we can approximate the Green's function of χ by where G I is the solution with no twist operator inserted ( fig. 5(b)), where I stands for identity, while δG T L and δG T R corresponds to the backreaction of the twist operator T L and T R . The support of δG T L and δG T R are separated by a real time evolution t much greater than β on the contour. When we evaluate the action, this leads to Here we have separated out crossing terms of δG(T L ) and δG(T R ) into δI C . There are two kinds of diagrams in δI C . Terms from G 4 is of the form: Here we have a > 0, b > 0 and a + b ≤ 4. Terms from the log det term is also suppressed by e −t/β . As an example: Since G −1 (s, s ) decays exponentially for large real-time separation. Combining these results, we have Realizing that exp(−I C (I)) = Z 2 is just the partition function, we find the saturation value of the entropy is just twice of the thermal Rényi entropy for subsystem χ, which means the factorization of twist operators.

Appendix D Entropy dynamics with an SYK bath
In this appendix, we present results for a related model by replacing the Majorana chain bath by an large SYK bath. This is an extension of the results in [14], where the authors studied the case with equal number of modes. We find similar results as the chain bath case. We now consider two SYK systems χ and ψ with different number of modes N χ and N ψ . The Hamiltonian of the system is written as: with the variances for J χ ijkl and J ψ ijkl being: We add an interaction term H int that couples χ and ψ system. In this model, we consider two types of interaction. One is of the "χ 2 ψ 2 " form: and another is of the "χψ 3 " form: We define the ratio of fermion number to be r ≡ N ψ /N χ . We again compute the second Rényi entopy (13) after the evolution of a TFD state. For the χ 2 ψ 2 case, the path-integral formalism now reads e −S (2) χ L ,χ R = 1 Z 2 DΣ χ DΣ ψ DG χ DG ψ exp(−S C [Σ,G]) (64) which gives the saddle point equation Similar effective action and saddle point equation can be worked out straightforwardly for the χψ 3 case. The r = 1 case with χ 2 ψ 2 interaction has been studied in [14,31], where no transition is found in the canonical ensemble. The transition appears if we instead consider the micro-canonical ensemble [14].
Here we instead focus on large r limit, whose equilibrium physics and quench dynamics have been studied in [28][29][30]. Since the qualitative features (the short time linear growth and tirst order transition regardless of the strength V 2 /J 2 . In contrary, for the χψ 3 interaction, similar to the chain case, there is no consistent exact replica diagonal solution and the transition only appears for small V 2 /J 2 .