Black holes Entangled by Radiation

We construct three models to describe the scenario where two eternal black holes are separated by a flat space, and can eventually be entangled by exchanging radiations. In the doubly holographic setup, we compute the entanglement entropy and the mutual information among the subsystems and obtain the dynamic phase structure of the entanglement. The formation of entanglement between the two black holes is delayed by the space where the radiations must travel through. Finally, if the two black holes exchange sufficient Hawking modes, the final state is characterized by a connected entanglement wedge; otherwise, the final entanglement wedge contains two separated islands. In the former case, the entanglement wedge of the two black holes forms at the time scale of the size of the flat space between them. While in both cases, unitarity of the evolution is preserved. When the sizes of two black holes are not equal, we observe a loss of entanglement between the smaller black hole and the radiation at late times. In the field theory side, we consider two Sachdev-Ye-Kitaev (SYK) clusters coupled to a Majorana chain, which resemble two black holes connected by a radiation region. We numerically compute the same entanglement measures, and obtain similar phase structures as the bulk results. In general, a time delay of the entanglement between the two SYK clusters is found in cases with a long Majorana chain. In particular, when the two SYK clusters are different in size, similar entanglement loss between the smaller SYK cluster and the Majorana chain is observed. Finally, we investigate a chain model composed of EPR clusters with particle exchanges between neighboring clusters, and reproduce the features of entanglement observed in the other models.


Introduction
Understanding black hole physics in the context of quantum mechanics is a notoriously difficult task. It leads to the famous problem of whether a black hole evaporates in a unitary fashion, which is also called the black hole information paradox [1][2][3]. On one side, the nearly thermal spectrum of Hawking radiation requires the entropy to keep growing, according to Hawking's earlier calculation. On the other side, however, the quantum nature of black holes demands that the entropy of radiation should decrease at late times and its dynamical evolution should obey the rules described by the Page curve [4][5][6]. See [7][8][9][10] for further debates.
Recently, a breakthrough was made to understand the unitary evolution in the context of Gauge/Gravity duality [11,12]. By extremizing the generalized gravitational entropy [13,14], the entanglement contributed by Hawking radiation can be captured by the quantum extremal surface (QES), through the island formula [15] S[R] = min I ext I Area [∂I] 4G N + S[R ∪ I] . (1.1) Here R and I represent the radiation and the islands, respectively. At the beginning of evolution, the entropy increases due to the accumulation of Hawking modes. While at late times, the growth of entropy is pinched off by the emergence of the island in the gravitational region. As a result, the entropy of radiation during evolution is in line with the Page curve [16]. See  for more recent discussions on islands. The island formula (1.1) was firstly proposed in the model of two-dimensional gravity coupled to twodimensional CFT matter sectors [15]. In addition to this time-dependent case, the island also emerges in the two-dimensional static geometry, where the eternal black hole is in equilibrium with two flat baths [39]. The full-time behavior of entropy for this model can be obtained by considering the holographic dual of this two-dimensional theory, known as double holography, where the QES becomes a standard RT/HRT surface in the threedimensional geometry [15,17]. The existence of islands in higher dimensions was first demonstrated in [40], where the lower-dimensional black hole is replaced by a brane with tension where Neumann boundary conditions are applied [41][42][43]. In this context, however, the dynamical Page curve can be obtained only under the condition that there are sufficient degrees of freedom (DOF) on the brane [44][45][46][47], where the prescriptions of inputting enough DOF have been proposed as well. See [45][46][47][48][49][50][51][52][53][54][55] for more brane-world construction in higher dimensions.
From the aspect of path integral formalism, the emergence of islands corresponds to the dominance of another saddle point, known as the wormhole saddle point, where Euclidean wormholes connect the different replicas [56,57] (See also [58,59]). Inspired by this thought-provoking viewpoint, the notion of "wormhole" was extensively discussed as the bridge between disjoint universes [60][61][62][63][64], where the wormhole behind horizons is lengthened, as well as in baby universes [65][66][67][68], which connect to parent universes through a wormhole.
To better understand the black hole information paradox, measuring alone the entanglement entropy of radiation is not sufficient. On one side, other measures were introduced to investigate the entanglement properties of radiation in [69][70][71][72][73][74][75]. On the other side, the entanglement properties inside the gravity system were studied in [76], where the emergence of the island trimmed off the growth of reflected entropy.
Further, an intriguing topic is to investigate two black holes entangled by exchanging Hawking radiation, based on the setup of thermofield double (TFD) states [39,[77][78][79]. We consider two black holes separated by a finite space. Intuitively, the Hawking radiation emitted by one black hole can eventually be absorbed by the other, if they travel through this space. By this process, two black holes will develop more and more entanglement. Due to unitarity, the entanglement will eventually reach an upper bound determined by the DOF of subsystems. In this paper, we strive to model these phenomena by analyzing the entanglement properties between the separate black holes and envisaging the formation of entanglement, from both the bulk gravity and boundary perspective. (See [80][81][82][83] for other discussions on the formation of wormholes.) For these purposes, we propose the following three kinds of models to simulate the entanglement between black holes and radiation: First, we construct a doubly holographic model with two Planck branes in (3 + 1) dimensions. The tension and Dvali-Gabadadze-Porrati (DGP) terms [84] are imposed on both branes to enhance the DOF. Second, inspired by the work in [85], we propose a quantum mechanical model with two separate SYK clusters coupled to a Majorana chain connecting them. Third, we will further simplify the model to be one composed of several EPR clusters which form a chain structure, with the outermost clusters being black holes and the inner clusters being radiation. See also [86][87][88] for other models in BCFT, BTZ black hole and black string with branes.
The paper is organized as follows: in section 2, we generalize the doubly holographic setup in [44] by introducing two Planck branes with Neumann boundary conditions on them. Then we apply the Einstein-DeTurck formalism to solve this gravitational system.
In section 3, we elaborate on the description of the holographic entanglement entropy and mutual information in different phases. More importantly, we find a novel phase concerning the appearance of a wormhole.
In section 4, we illustrate the dynamical phase structures and transitions of entanglement during the evolution in the probe limit. Subsequently, we take the backreaction into account and explore the effect of DOF in black holes on the phase structure.
In section 5, we propose a quantum mechanics model that connects SYK models with Majorana chains. The eigensystems are calculated by exact diagonalization and the related measures in quantum information theory are computed.
In section 6, we construct a toy model which consists of several EPR clusters. By allowing them to exchange particles with the nearby clusters, we qualitatively simulate the process of black holes entangled by radiation.
Finally, the conclusions and discussions are given in section 7.

The doubly holographic setup
In this section, we will present a general setup for two charged eternal black holes exchanging Hawking modes. Consider two d-dimensional eternal black holes which are asymptotic AdS d , coupled to two finite-sized flat baths, where Hawking modes living in this combined system are described by the d-dimensional CFT matter sector, as shown in Fig. 1 Fig. 1(a). The second is the bulk gravity perspective, where the matter sector is dual to a (d + 1)-dimensional bulk and B 1 and B 2 are described by Planck branes pl 1 and pl 2 in the bulk, as shown in Fig. 1(c). In the first model, the setup will be constructed from the bulk gravity perspective, while the dynamical processes will be described from both the bulk gravity perspective and the brane perspective. From the bulk gravity perspective, the action of the (d + 1)-dimensional gravity theory is specified as Here K ∂ is the extrinsic curvature on the conformal boundary ∂. The electromagnetic curvature is F = dA. K i is the extrinsic curvature and the parameter α i is proportional to the tension on the brane pl i , which will be fixed later. The last term in the second line is the junction term at the intersection of the brane pl i and the conformal boundary ∂, where θ i is the angle between the brane pl i and the boundary, while Σ i is the metric on pl i ∩ ∂. The first term in the last line is the DGP term [45], where G b,i is the additional Newton constant on each brane and R h i is the intrinsic curvature on each brane. The second term in the last line is the junction term at the intersection pl i ∩ ∂ of the brane Figure 2. A simple setup of Planck branes. Here Planck branes are anchored on the conformal boundary at (z, w) = (0, 0), (z, w) = (0, w 0 ) respectively and penetrate into the bulk with angles θ 1 and θ 2 . and the conformal boundary, where k i is the extrinsic curvature on pl i ∩ ∂. Taking the variation of the action, we obtain the equations of motion as where T is the trace of the energy-stress tensor T µν .

Boundary conditions on the Planck branes
In AdS/CFT setup with infinite volume, the (d + 1)-dimensional bulk is asymptotic to AdS d+1 which in Poincaré coordinates is described by with the conformal boundary at z = 0. Let θ i be the angle between the i-th Planck brane and the conformal boundary as shown in Fig. 2. Then Planck branes pl 1 and pl 2 are described by hypersurfaces respectively near the boundary. One should cut the bulk ending on these two branes and restrict the region between them [41,84,89,90]. We will also impose these constraints (2.5) deep into the bulk and find its backreaction to the geometry. For the boundary terms in (2.1), we impose Neumann boundary conditions on each Planck brane, which are where h i is the induced metric on the brane pl i , and λ i ≡ can be regarded as the effective coupling of the DGP term on this brane. The parameter α i in action (2.1) is fixed to be a constant by solving (2.6) near the conformal boundary to concrete the tension term on the brane. In addition to (2.6), we also impose Neumann boundary conditions for the gauge field A µ on the brane pl i [41][42][43]91], which is where n µ is the normal vector to that brane and the subscript B denotes the coordinates along the brane.

The Einstein-DeTurck formalism
The line element of the standard RN-AdS d+1 (d ≥ 3) geometry is where µ is the chemical potential of the boundary theory, the outer horizon has been scaled to be z = 1 and one can recover it by transforming coordinates In order to restrict the region between the Planck branes (0 ≤ x ≤ 1), and outside the event horizon, we transform the coordinates to That is, for x = 0 we have w+cot θ 1 z = 0, while for x = 1 we have (w−w 0 )−cot θ 2 z = 0 (see Appendix A for more discussions). With the presence of Planck branes, the geometry is not precisely RN-AdS, and the deformation of the ambient geometry due to the backreaction of Planck branes can be described by introducing the Deturck method [92]. Instead of solving (2.2) directly, we solve the so-called Einstein-DeTurck equation, which is where is the DeTurck vector andḡ is the reference metric. Hereḡ is required to satisfy the same boundary conditions as g only on Dirichlet boundaries, but not on Neumann boundaries [40].
Now we introduce the metric ansatz in the four-dimensional case and derive the boundary conditions in the doubly holographic setup. Since the translational symmetry along the x direction is broken, the most general ansatz of the ambient geometry is with {Q 1 , Q 2 , Q 3 , Q 4 , Q 5 , Q 6 } being the functions of (x, y). All the boundary conditions are listed in Tab. 1, with the tension on the brane (2.17) (see Appendix B for derivation). Moreover, the boundary conditions at the horizon y = 0 also imply that Q 1 (x, 0) = Q 2 (x, 0), which fixes the temperature of the black hole as The reference metricḡ is given by Q 1 = Q 2 = Q 5 = 1, and With the metric ansatz (2.14), the ambient geometry is numerically solved by the Newton-Raphson method. Specifically, we discretize the Einstein-DeTurck equations (2.13) both on the x and y directions with the Chebyshev Pseudo-spectral method.
Distinguishing the DeTurck soliton and the solution of Einstein equations is always important. Note that the possible solution of (2.13) with ξ = 0 is called DeTurck soliton, which is not the solution of the Einstein equations (2.2). In literature, it is shown that the DeTurck soliton does not exist for static geometries with Dirichlet boundary conditions, while for Neumann boundary conditions, this conclusion has never been proved [40,93,94]. Although there is no rigorous proof, the problem we address is well-posed elliptic, where the local uniqueness of solutions is guaranteed. That is to say, a soliton must be distinguishable   from any genuine solution. Therefore, in practice, we just need to monitor the gauge vector ξ, and make sure it will converge with the increase of the grids -see Appendix C for details.
Owing to the Neumann boundary conditions on the brane, there is no analytical expression of the metric. See Fig. 3 for a concrete example, where the values of {Q 1 , Q 2 , Q 3 , Q 4 , Q 5 , Q 6 } are obviously deviated from those ofḡ.

Quantum extremal surfaces and phase structures
From the brane perspective, the entanglement entropy of the radiation R is measured by a quantum extremal surface. Further from the bulk gravity perspective, the QES is equivalently described by a (d − 1)-dimensional HRT surface and a lower-dimensional area term on the brane pl [40,45,49,50], namely where γ I∪R is the corresponding HRT surface sharing the boundary with I ∪ R. In the ordinary doubly holographic setup with one Planck brane, there are two possible phases characterized by distinct configurations of HRT surfaces - Fig. 4. The first phase γ tv occurs at early times with the absence of an island I, while another phase occurs at late times, and a non-trivial island I emerges to keep the entanglement entropy (1.1) from divergence after the Page time [39]. For simplicity, we consider θ 1 = θ 2 and λ 2 ≥ λ 1 throughout this paper, and the involvement of a second Planck brane renders four possible phases for the configuration of HRT surfaces at most 1 -which are illustrated in Fig. 5.
As proposed previously in [44,76], we will still call the phases as shown in Fig. 5(a) and 5(c) the trivial phase and the island phase respectively, since the configurations of HRT surfaces are similar to the cases with one Planck brane. In addition, it is interesting to notice that the presence of the second Planck brane brings two more possible phases during evolution. The first possible phase is called the half-island phase - Fig. 5(b), at which only one island emerges on the brane due to the fewer DOF on it. Another possible phase may occur at the saturation of the whole system, which we call it the wormhole phase - Fig. 5(d). In this phase, the entanglement wedge connecting B 1 and B 2 appears in the bulk, representing the entanglement between them. We will call it a "wormhole" connecting these two black holes, following the idea of ER = EPR conjecture [10] (see also [95] for some constraints on experiments). The dynamics of the wormhole formation will be analyzed in the next section.
Technically, the endpoints of HRT surfaces should locate near the branes. Consider these HRT surfaces anchored at x = x b and x = 1 − x b respectively on the conformal boundaries - Fig. 5. They measure the entanglement between two subsystems: one consists of B 1 , B 2 and part of baths within the region Conventionally, we still call the former the black hole subsystem, and the latter the radiation subsystem.

The entropy and mutual information in different phases
From now on, to avoid cumbersome statements we will call the trivial, half-island, island and wormhole phase simply as Ph-T, Ph-H, Ph-I, and Ph-W respectively. From (3.1), the entropy density of the radiation subsystem R can be determined by Here V is the infinite volume of the relevant spatial directions. For instance, for d = 3 we is the area of the corresponding extremal surface G. Moreover, in the second line, we have set L = 1. While for each black hole subsystem B i (i = 1, 2), the formula of entropy density is the same as that in the ordinary double holography with one Planck brane -see Fig. 4, which is Moreover, the density of mutual information between any two subsystems K and K can be expressed as For simplicity, all these quantities are dimensionless. Following the transformation (2.11), their dimensions can be recovered as Next, we will derive the expressions for densities A/V in each phase in parallel. First, for the extremal surface γ T i at t = 0, we introduce the parameterization x = x i (y) (i = 1, 2) -see Fig. 6(a) for a concrete example, which leads to the corresponding density Then, for γ I i together with ∂I i (i = 1, 2), we introduce two different ways of parameterization in different intervals. In the (x, y) plane, for the curve in y ∈ [y c , 1], we introduce Fig. 6(b). Note that the anchoring point y 1 (0) = y I 1 on the brane can be read off from the integration procedure. Thus, the surface with the extremal area can be found by running the anchoring point all over the brane. Subsequently, the corresponding density associated with that anchoring point can be expressed by . (3.8) Similarly, for γ I 2 anchored at y 2 (1) = y I 2 , we have . (3.10) Finally, for the surface γ W , we introduce the parameterization y = y(x)- Fig. 6(c), which leads to the corresponding density With these expressions at hand, we are now ready to investigate the entanglement properties of the system. In the next section, we will elaborate on the description of the entanglement entropy and mutual information in the probe limit, which renders complicated phase structures. Further, we will take the backreaction of branes into account and study its effect on both the entanglement properties and phase structures.

Page curves and entanglement phase structures
Page curve plays a significant role in understanding the black hole information paradox. In our case, there are also Page curves caused by exchanging hawking modes of separate black holes. However, usually, it is difficult to obtain the Page curve in the doubly holographic model since the backreaction of the brane to the ambient geometry is very hard to solve such that one lacks the data inside the event horizon during evolution [44].
Fortunately, this difficulty can be circumvented by considering the special case with the probe limit [44]. In such cases, Planck branes pl i apply no backreaction to the ambient geometry. This allows us to consider the entanglement properties with respect to time. The strategy to calculate the Page curves is as follows: first, we analytically calculate the growth rate of entropy densities, and then numerically find the saturation values of the Page curves by properly choosing the distance w 0 between B 1 and B 2 . As a result, the full Page curves can be obtained by this hybrid method.
The upshot is that the unitary evolution can be preserved not only by the emergence of islands but also by the formation of a wormhole. In the next subsection, we will first calculate the growth rates of quantum information measures such as entanglement entropy and mutual information.

The growth of entropy and mutual information densities
In the probe limit, i.e. |π/2 − θ i | 1 together with λ i 1, we have g →ḡ, and hence, the corresponding ambient geometry in any dimensions can be treated as the standard RN-AdS black holes as shown in (2.8) with 0 ≤ w ≤ w 0 . Now we obtain the data of the geometry inside the horizon, which allows us to keep track of the growth of γ T i with time. To see this, we firstly express the density functional of γ T i in the Eddington-Finkelstein coordinates {v, z, w, w j } as where Ξ i is the intrinsic parameter on γ T i , and the corresponding time on the boundary system is Due to the cyclic coordinate v in (4.1), the integral of motion can be obtained as Furthermore, since the action in (4.1) is invariant under reparametrizations, we can freely choose the integrand as We subsequently substitute both (4.2) and (4.3) into (4.1), and then the result is d dt Here z max denotes the turning point of the trivial surface γ T i and the derivation is presented in Appendix D. Moreover, the relation between z max and the integral of motion C is given by At late times, the extremal surface γ T i tends to surround a special slice, with z max = z M [79]. Define it is easy to show that C 2 = F (z max ) 2 will keep growing until approaching an extremum at z max = z M , where we have the following relation By substituting the solution of (4.8) into (4.4), we finally get the growth rates at late times as As a result, the growth rate of the entropy density of R at late times is Ph-I and Ph-W. (4.10) Note that the first is always twice as much as the second since the entropy on one of the branes reaches saturation in Ph-H. Specifically, after recovering the dimension, for the neutral cases where µ → 0, we have 2d . Moreover, a = 2 if the system is in Ph-T, a = 1 if the system is in Ph-H, and a = 0 if the system is in Ph-I or Ph-W. While for the extremal cases where T h → 0, we have and the linear behavior with respect to T h indicates that the near horizon geometry of the near-extremal black hole is AdS 2 × R d−1 spacetime. In this case, the evolution is nearly frozen due to the low temperature and the entropy density barely grows. Further in Ph-W, the entropy density of R (3.2) saturates, while the mutual information density (3.4) between B 1 and B 2 might not. Similarly, the late-time growth of the mutual information density can be expressed as Hence, the linear growth of mutual information density is similar to that of entanglement entropy density as shown in (4.11) and (4.12). The similar derivations for remaining s [K] and I[K : K ] will not present here. These quantum information measures are crucial for analyzing the evolution process including the formation of wormholes.

Dynamical evolution and phase structures in the probe limit
In this subsection, we will elaborate on the entanglement phase structures during evolution. First, we investigate the cases with equal-sized black holes, where B 1 and B 2 contain the same DOF, which is simple but enough to generate a wormhole. Then, we study the cases with different-sized black holes, where the smaller (larger) black hole B 1 (B 2 ) contains fewer(more) DOF.
For numerical convenience, here we continue to work with the coordinate set {t, y, x, w 1 } in 4-dimensional spacetime. In addition, the UV cut-off near the conformal boundary is fixed to be = 1 − 1/100. We also use scale-free parameters in the following discussion, such as

Black holes of equal size
Different phases are undergone by the combined system during evolution - Fig. 7, at different distancesw 0 between B 1 and B 2 as well as temperatureT h . As shown in Fig. 7(a), there are two critical surfaces in orange and blue, which separate Ph-T, I and W. The reason for no Ph-H is that the black hole subsystems with the same DOF are saturated at the same time.
For instance, at the beginning of evolution, the hot systems with subsystems B 1 and B 2 being far apart are always in Ph-T, which indicates that black hole subsystems are not entangled with each other in this phase. While the cold systems with adjacent subsystems are generally in Ph-W, which manifests that in these cases B 1 and B 2 are entangled with each other and hence, there is a wormhole connecting them.
Based on these observations, we extract three characteristic types of evolution - Fig. 7(b). Intuitively, these types are distinguished by the formation of wormholes. First, for adjacent black holes where the distancew 0 is less than a critical scalew c1 (T h ), wormholes have emerged from the beginning. Second, for middle-ranged black holes wherew c1 (T h ) < w 0 <w c2 (T h ), wormholes are generated in the intermediate stage. Finally, for distant black holes wherew 0 is greater than another critical scalew c2 (T h ), wormholes never appear. Again, to avoid cumbersome statements we will name these types of evolution simply as Type-I, Type-II, and Type-III respectively.
Next, we will elaborate on these types of evolution by tracking the entanglement entropy density and mutual information density of each subsystem, and the detail is given as follows:  • Type I: the evolution of adjacent black holes,w 0 <w c1 (T h ) For adjacent black holes, the whole system is always in Ph-W during evolution - Fig. 8(a).
At the beginning, a wormhole has already been generated between the adjacent black holes B 1 and B 2 , so they share nonzero mutual information densityĨ[B 1 : B 2 ] at t = 0. Whereafter, the entanglement between B 1 and B 2 becomes stronger owing to the accumulation of exchanged Hawking modes between them. By the entanglement wedge reconstruction, the information on both Planck branes can only be encoded in B 1 ∪ B 2 , but not in the radiation subsystem -see Fig. 5(d).
Untilt =t 13 , black hole subsystem B 1 (B 2 ) is fully correlated with R ∪ B 2 (B 1 ), and hence, the mutual information densityĨ[B 1 : B 2 ] approaches saturation. In this circumstance, in addition to the former reconstruction, we can also extract information on brane pl 1 (pl 2 ) from the combined system R ∪ B 2 (B 1 ).
These phenomena are different from those in the standard double holography supporting one Planck brane [15,39,40,44,57,96], and are also direct consequences of the formation of the wormhole.
• Type II: the evolution of middle-ranged black holes,w c1 (T h ) <w 0 <w c2 (T h ) For Middle-ranged black holes - Fig. 8( As a result, we obtain a dynamical curve ofs[R] - Fig. 8(b). Although it exhibits similar behaviors as in the standard double holography, in this circumstance, the unitary evolution is preserved by the formation of wormholes rather than the emergence of islands.
• Type III: the evolution of distant black holes,w 0 >w c2 (T h ) For black hole subsystems which are far apart - Fig. 8(c), the evolution also tends to begin with Ph-T. That is, B 1 ∪ B 2 starts to be entangled with R as the process goes on.
Aftert =t 13 , the radiation subsystem R is fully entangled with B 1 ∪ B 2 , wheres[R] saturates, as well ass[B i ], but the growth ofĨ[B 1 : B 2 ] is suppressed by the long distance.
Note that we acquire a Page curve here without the formation of wormholes. Therefore, the unitarity during evolution is preserved by the emergence of islands. In this sense, the situation is much similar to those in literature.

Black holes of different sizes
In the following, we will generalize the discussion to the cases with different-sized black holes, where the smaller(larger) black hole B 1 (B 2 ) contains fewer(more) DOF. The phase diagram is demonstrated in Fig. 9, by varying both the distancew 0 and Hawking temperaturẽ T h . We remark that one more possible phase (Ph-H) appears in this case, due to the different instants of the saturation ofs B i . Nevertheless, we point out that the addition of Ph-H will not change the classification of evolution. Thus we still envisage the evolution with the following types: • Type I: the evolution of adjacent black holes For adjacent black holes, the full system still tends to stay in Ph-W - Fig. 10(a). Nonzero mutual information density att = 0 indicates the existence of the wormhole connecting B 1 and B 2 . As time passes by, the accumulation of Hawking modes directly enhances the entanglement between them.
Att =t 12 , the entropy densitys[B 1 ] of the smaller black hole subsystem B 1 saturates first, while that of the larger one B 2 does not. As a result, the entanglement between B 1 and R will continue to "pass" into B 2 (whereĨ[B 1 : R] decreases whileĨ[B 2 : R] increases). From the brane perspective, in this phase, we are allowed to reconstruct the information inside the black hole on brane pl 1 from the boundary region R. Finally, att =t 23 , all measures under consideration reach equilibrium.
• Type II: the evolution of Middle-ranged black holes For middle-ranged black holes, the system undergoes either Ph-T,W or Ph-T,H,W, and then reaches the saturation. The desired wormhole will be generated during the evolution in both routes.
Let us firstly discuss the evolution undergoing Ph-T,W phases- Fig. 10(b). At early times, B 1 and B 2 start to exchange Hawking modes with R. Att =t 14 , R is fully entangled with the surroundings, and hence, B 1 begins to be correlated with B 2 . Untilt =t 12 , the entropy density of B 1 saturates, while that of B 1 does not. Thus the entanglement between B 1 and R passes into B 2 all the way until the equilibrium is reached att =t 23 .
As for the second case - Fig. 10(c), the main difference from the former is that the entropy density of B 1 saturates prior to that of R. After the saturation of B 1 at t =t 12 , B 2 continues to interact with R. Until the saturation of R att =t 24 , the entanglement between B 1 and R is reassigned into B 2 . We remark that in this case, the wormhole occurs later and the maximal entanglement is much lower, reflecting the intuitive fact that the farther the distancew 0 , the harder it is for B 1 and B 2 to be correlated.
We also point out that in both cases unitarity is preserved by the formation of a wormhole rather than an island.
• Type III: the evolution of distant black holes For distant black holes, the system is in Ph-T at the beginning and finally saturates in Ph-I.
In this type of evolution, the entropy density of B 1 saturates first att =t 12 , while those of B 2 and R saturate later att =t 23 . As a result, the entanglement between B 1 and B 2 continues to grow untilt =t 23 .
Similar to the Type-III of black holes with equal size, the wormhole is never generated if only considering the leading order of mutual information and unitarity is preserved by the emergence of islands as described in literature.
In this subsection, we have demonstrated that when the amount of DOF on the brane is small such that the backreaction of the brane can be ignored, the evolution of the system can be classified into three universal types according to the distance between the black holes. Then a natural issue arises and remains unanswered. That is, whether and how these types of evolution will be affected by the amount of DOF on the branes. Therefore, inspired by the work in [44], we tend to take the backreaction of branes into account and investigate its effect on the existence of the wormhole in the next subsection.

The wormhole phase in backreacted spacetime
To illustrate the universality of the wormhole phase, or more precisely, the universality of the evolution characterized by three types, in this subsection, we take the backreaction into account, which is equivalent to joining DOF on the black hole subsystems B 1 and B 2 . Moreover, in the probe limit, we need an artificial nonzero scale w b to produce a dynamical Page curve, but in backreacted cases, we can safely eliminate it. Therefore, we change the definition of symbol "tilde", which results in new scale-free parameters, such as For simplicity, here we explore neutral cases with the same DOF on branes. That is to say, we specifyμ = 0 and θ 1 = θ 2 = 3π/16, while vary λ := λ 1 = λ 2 andw 0 . The corresponding HRT surfaces and types of evolution are shown in Fig. 11(a) and 11(b), respectively. It can be seen that there are also two critical scalesw c1 andw c2 which classify the evolution into three types. For λ > 0.7,w c2 grows with λ, whilew c1 is rarely changed with λ. This result indicates that with more DOF, B 1 and B 2 are more likely to be correlated, and thus the wormhole is more likely to be generated during evolution.

SYK clusters coupled to Majorana chains
In this section, a similar setup in the context of SYK models is constructed by identifying B 1 and B 2 as distinct SYK systems, and radiation R as Majorana chains connecting them. This combined system is studied by exact diagonalization (ED). The entanglement entropy and mutual information of subsystems are investigated during the evolution, and the results are compared with those obtained in the previous section.

The setup
To be specific, the SYK systems labeled by χ 1 and χ 2 are composed of N χ 1 and N χ 2 Majorana fermions with random couplings among them, respectively, while the Majorana chains are described by ψ(x) on N ψ separate sites with hopping Λ/2, and the fermions on the first and last sites are coupled with nearby SYK systems respectively. The combined system is sketched in Fig. 12. The corresponding Hamiltonian of the whole system is given by where J i (i = χ 1L , χ 2L ) are in Gaussian distribution, which satisfy Starting with a thermofield double (TFD) state, with six subsystems χ 1L , χ 2L , ψ L and χ 1R , χ 2R , ψ R , the time evolution of the TFD state can be constructed as where |I K L ,K R is a maximally entangled state between K L and K R subsystems, where K ∈ {χ 1 , χ 2 , ψ}. On the SYK systems, the state is defined as On the Majorana chains, the state is defined as Next, the density matrix of the subsystem K can be obtained by tracing out the DOF of its complementK, which is given by For instance, when K = χ 1 , we haveK = χ 2 ∪ ψ. Then the entanglement entropy of a system K is defined by S[K] = −TrK(ρ K log ρ k ), (5.8) and the mutual information between K and K is defined to be

The entanglement properties of subsystems
In this subsection, we investigate different measures in quantum information to keep track of the entanglement properties during evolution. Precisely, χ 1 and χ 2 serve as the black hole subsystems B 1 and B 2 respectively, while ψ is regarded as the radiation subsystem R. The amount of DOF in B 1 and B 2 are referred to as the number of fermions in χ 1 and χ 2 , namely N χ 1 and N χ 2 , and the distance between B 1 and B 2 is captured by the number of sites N ψ in Majorana chains.
In principle, we hope the combined SYK system can not only resemble a classical gravity system, but also reach local thermalization quickly for each subsystem, which specifies the hierarchies between the parameters in the above Hamiltonian. The former requires 1 J i β N χ i and the latter requires J i β V i √ Λβ. In practice, although it is difficult to realize all the large hierarchies in ED, we still maintain the proper gaps between the parameters. The numerical computations of the entanglement entropy and mutual information are illustrated in Figs. 13 and 14. Their evolution still shares many qualitative features with the entanglement evolution in the gravity system, as discussed below, which suggests the universal evolution of the entanglement of subsystems. At the end of this subsection, we will also discuss some distinctions between the present setup and the gravitational one. Fig. 13 • For adjacent SYK subsystems with small N ψ , the mutual information I[χ 1 : χ 2 ] is positive and grows immediately at the early stage - Fig. 13(a) and 13(d). This indicates that the SYK systems have already been correlated at t/β = 0, which agrees with Type-I evolution - Fig. 8(a), where the wormholes have been generated for adjacent black holes.

Equal-sized cases with
• With the growth of the number of sites N ψ in Majorana chains, the growth of I[χ 1 : χ 2 ] is delayed, and I[χ 1 : χ 2 ] saturates at a lower level. The first agrees with the result that the appearance of Ph-W is delayed with the increment of the distance of black holes in Type-II evolution - Fig. 8(b), where the unitarity is preserved via this phase transition. While the second implies that the evolution approaches Type-III - Fig. 8(c), with the growth of distance.
• For all of N ψ , S[χ 1 ] are always close to S[χ 2 ], since χ 1 and χ 2 have the same number of Majorana fermions and their couplings obey the same Gaussian distribution. In the gravity system, we always finds B 1 =s B 2 , when B 1 and B 2 contain the same amount of DOF.  χ 1 and ψ starts to "pass" into χ 2 . This phenomenon is similar to that in the gravity system, where the entanglement transmits into B 2 at late times, as shown in Fig. 10(b).

Different-sized cases with
• The growth rate of I[χ 1 : χ 2 ] has fallen nearly by half at the scale of the bifurcation - Fig. 14(b), which also agrees with the result in Fig. 10(b).
• S[χ 1 ] always saturates earlier and at a level lower than S[χ 2 ], due to fewer DOF in χ 1 compared with those in χ 2 . In the gravity system, this result is also manifest.
At the end of this section, we remark that aside from these similarities, there are also some distinctions in this quantum mechanical setup compared to the gravity system, which mainly Cstems from the deviation from both the large N limit and strongly coupled limit, especially in the Majorana chain. First, in the combined SYK model, χ 1 and χ 2 will always be entangled within the range of parameters, and it is reasonable to expect the existence of entanglement in any choice of parameters. While in the gravity system, since we only observe the results in the leading order, the entanglement between B 1 and B 2 always vanishes as long as the distance is far enough. Second, the saturation of entropies is always smooth in the combined SYK model, which is distinct from the first-order phase transition of entropies in the gravity system. Third, some decaying oscillations occur in the growth of the entanglement entropy and mutual information, due to the lack of thermalization in the Majorana chain.

A toy model for the evolution of entanglement: EPR clusters
We propose a toy model to mimic the evolution of two-body entanglement in the previous systems. Considering several clusters of EPR pairs, each pair consists of one left and one right qubit, which together build a maximally entangled state. Collecting all the left and right qubits respectively, the resulting huge maximally entangled state resembles a TFD state, which is the model we are mainly interested in in this section.
Specifically, we denote M clusters of EPR pairs as {G i |i = 1, 2, .., M }, each of them containing N i EPR pairs, and denote the total number of EPR pairs as N = M i=1 N i . Then, consider the time evolution only on the left qubits under a simple interactionexchanging one left qubit in G i with one left qubit in G j at the frequency of J ij , which plays the role of a coupling constant. Note that the right qubits remain static during evolution. We define N ij as the number of left qubits in cluster G j entangled with the right qubits in cluster G i , which satisfies the constraints N i = M j=1 N ij = M j=1 N ji . The evolution equation for N ij in mean field theory is In the unit of 2 ln 2, the entanglement entropy of G i is S i = N i − N ii , and the mutual information between G i and G j is I i;j = N ij + N ji . Now we divide all clusters into three subsystems, namely B 1 = G 1 , B 2 = G M , and R = ∪ M −1 i=2 G i , which respectively serve as two black hole subsystems and radiation subsystems in the previous holographic model. Then, we specify the coupling constants J ij as J 12 .., M −2}, with others being 0, and N α = N α+1 = n R , with α ∈ {2, · · · , M − 2}. As a result, these clusters together construct a chain model, which resembles the previous holographic setup - Fig. 15. In this combined system, we have the following identifications: • N 1 , N M , n R are regarded as the DOF in the black hole and radiation subsystems.
• M is referred to as the distance of distinct black hole subsystems.
• V characterizes the strength of interactions between B 1,2 and R.
• Λ characterizes the speed of propagation of qubits in R.
In this context, the entanglement entropies and mutual informations can be expressed as The evolution of entanglement for typical parameters of M 1 and n R < N 1 < N M is illustrated in Fig. 16. The results are also similar to those obtained in the previous holographic model: • Overall, entanglement entropy and mutual information grow at early times, because of the exchanging process, and finally reach equilibrium.
• The growth of I[B 1 , B 2 ] exhibits a time delay related to M , because the qubits that are exchanged between B 1 and B 2 have to traverse R.

Conclusions and discussions
In this paper, we have constructed three different kinds of models to demonstrate the dynamic process that two black holes can be entangled by emitting radiation, which exhibits a generic phenomenon in a quantum system with black holes. First of all, from the bulk gravity perspective, the doubly holographic model has been constructed numerically, where the holographic CFT matter sectors are dual to a higher dimensional bulk and the black holes B 1 and B 2 are replaced by two Planck branes in the bulk. The bulk geometry is solved by the Einstein-Deturck equations with Neumann boundary conditions on the branes. The DOF on branes pl i are specified both by the angle θ i and the DGP coupling λ i . Restricting B 1 to be a smaller black hole with fewer DOF, we have found four possible phases of QES during evolution, including a novel wormhole phase, which corresponds to the formation of wormholes owing to the exchange of Hawking modes between black holes. In the probe limit, we have elaborated on the entanglement phase diagrams and dynamical evolution for black holes of both the same and different sizes. In general, the entanglement entropy of each black hole subsystem increases at early times due to the exchange of Hawking modes with the environment, while saturates at late times. Furthermore, three universal types of evolution can be distilled: first, for adjacent black hole subsystems, a wormhole has been generated from the beginning, and no dynamical Page curve of radiation is found. Second, for middle-ranged black holes, the instant of wormhole formation is delayed by the distance, and unitarity of the evolution is preserved. Third, for distant black hole subsystems, the entanglement between black holes vanishes up to the leading order of quantum corrections in the bulk, because of the long distance. Consequently, unitarity in this type will only be preserved by the emergence of islands. Moreover, for the cases with B 1 being a smaller black hole, the entanglement between the smaller black hole and radiation will leak into the larger one at late times.
The second model we have proposed is a quantum mechanical system that is composed of two separate SYK clusters χ 1 and χ 2 connected by a Majorana chain ψ. In this setup, SYK clusters are regarded as black hole subsystems and Majorana chains serve as radiation. Concretely, the DOF both in SYK clusters and the Majorana chains are specified by the number of fermions they contain, namely, N χ 1 , N χ 2 and N ψ . To resemble the former gravity setup and consider the limitation of the numerical calculation, we further require the weak hierarchies N χ i > βJ i > 1 and βJ i > V i √ Λβ, where βJ i and V i √ Λβ are the dimensionless couplings inside the SYK system χ i , and between the SYK system χ i and the Majorana chain ψ. With the finite number of fermions, the system is solved by exact diagonalization, and we have numerically computed the same entropy measures in this coupled SYK model. We have investigated both the equal-sized SYK clusters with N χ 1 = N χ 2 and the differentsized SYK clusters with N χ 1 < N χ 2 . In summary, with short Majorana chains, the SYK clusters will be entangled from the beginning. With the increment of the chain length, the instant of the entanglement between the SYK clusters is delayed and the maximum of the entanglement is suppressed. In addition, in the cases of the SYK clusters with different sizes, the entanglement between the "smaller" SYK cluster χ 1 and the Majorana chain ψ will pass into the "larger" SYK cluster χ 2 .
The third model consists of several EPR clusters {G i |i = 1, 2, · · · , M } forming a chain structure including two outermost clusters G 1 , G M referred to as black holes B 1 , B 2 and the inner clusters serving as radiation R. Concretely, the DOF in each subsystem is specified by the number of EPR pairs it contains, and the distance of B 1 and B 2 is referred to as the number of clusters M . By exchanging particles with nearby clusters, the entanglement properties of the subsystems can be explored. Focusing on the case with N 1 < N M , and in the large distance limit M 1, we find the growth of entanglement between B 1 and B M exhibits a time delay. Also, B 1 will lose some entanglement with R at late times due to the monogamy of entanglement.
Nevertheless, the deviation from the large N and strong coupling limit in our quantum mechanical setup will generally lead to some discrepancy in the entropy and mutual information in comparison with the results in a gravity system. So we are looking forward to studying the entanglement by solving the Schwinger-Dyson equations of the coupled SYK model in the future.
It is also interesting to compare our model to [97], where an eternal traversable wormhole was constructed in both nearly-AdS 2 gravity and the SYK model by coupling two entangled subsystems. These two subsystems share the same Hamiltonian and are coupled by a double-trace interaction which matches their entanglement structure of the ground state, a nearly-TFD state. The ground state is eternally traversable, where a signal injected from one subsystem can travel into the holographic bulk and appears in the other subsystem [98]. Turning on this interaction does not affect the entanglement entropy between these two sides [99]. While in our construction, the subsystems B 1,2 or SYK 1,2 are not directly coupled, and their entanglement increases gradually as Hawking modes travel through the bath. In other words, the clusters SYK 1,2 in our system get entangled and reach the wormhole phase in a different manner.
In addition to the above discrepancies, our models also share some similarities with [97]. On one side, our phase structures at late times are similar to those in [97] at finite temperature. From the bulk gravity perspective, the entanglement wedge of B 1 ∪ B 2 is connected (disconnected) in the wormhole phase (island phase). This feature is similar to the connected (disconnected) bulk in the eternal-traversable-wormhole phase (two-blackhole phase) in [97]. On the other side, the transitions between these two phases are also similar. To see this, we trace out the bath and consider the purity of the remaining two subsystems 2 , which is negatively correlated to the second Renyi entropy of the bath. When the bath is small (huge), the purity is high (low), then the two subsystems tend to end in the wormhole phase (island phase). In this sense, this transition is similar to the Hawking-Page transition in [97], where at low (high) temperature or high (low) purity, the eternal-traversable-wormhole phase (two-black-hole phase) dominates. Another interesting direction is to consider the dynamics in the full backreacted spacetime, which cannot be covered by our setup in this paper, due to the lack of metric data inside the event horizon. Nonetheless, by applying the ingoing Eddington coordinates in the Einstein-DeTurck formulation [100], this problem may be conquered. It must be pointed out that, the Einstein equations become a mixed elliptic-hyperbolic system, and hence the local uniqueness of solutions is not guaranteed in this setup [92]. Beyond the stationary cases, it is also interesting to further generalize the numeric setup into the nonequilibrium cases [17].
Furthermore, observing various entropy computed in the coupled SYK model indicates that the result of each individual realization seems to have some noise around the averaged result, although our computation is limited by the number of fermions. This could be related to the recent "half-wormhole" interpretation of the noise of the spectral form factors [101][102][103][104][105]. It is interesting to see in our case that the half-wormhole also contributes to the entanglement entropy. This is in fact a natural result from the holographic point of view; holographically the entropy is computed by the bulk gravitational path integrals that involve bulk geometries compatible with the replica trick in the boundary field theory. Then it is very likely that (half-)wormhole-like off-shell bulk geometries could also contribute to this computation, which could be the origin of the fluctuations of the entropy for each realization of the random coupling. It is interesting to show how various "halfwormhole" contributions enter the holographic computation of various entropies explicitly. In particular, the current case involves two asymptotic boundaries connected by a pair of branes, which could be translated to a boundary with two pairs of matter trajectories and the half-wormhole analysis in [106], see also the follow-up discussions in [107] could be applied directly. It is interesting to further optimize our numerical code, probably borrowing some numerical techniques from [108], and verify this behavior in the future.

Note added
While this work was being completed Ref. [88] appeared, where the authors considered a similar doubly holographic setup but in AdS 3 /BCF T 2 duality and found a similar dynamic phase structure of entanglement. Our discussion applies to general dimensions, and also to the case where the black holes are different in size, realized via different tensions and DGP terms. Furthermore, we also construct two microscopic many-body models to study the evolution of entanglement. Finally, we discuss the time delay of entanglement due to the propagation of radiation and the loss of entanglement at late times due to the different sizes of the black holes.

B Derivation of the free parameter α i
Near the boundary, the geometry is asymptotic to AdS 4 which in Poincaré coordinates is described by For the Planck brane at w + cot θ 1 z = 0, the outward normal co-vector is n = n µ dx µ = L z −1 1 + cot 2 θ 1 (cot θ 1 dz + dw) . In this appendix, we examine the numerical convergence of ξ 2 to distinguish the standard solution of Einstein equations (2.2) and the corresponding DeTurck equations (2.13). By adjusting the number of Chebyshev grid points both along x and y directions, which we denote as N x and N y respectively, we monitor the norm of DeTurck vector ξ 2 := ξ µ ξ µ -see Fig. 17(a). Specifically, for fixed grid points {N x , N y } and other parameters such as {L, θ 1 , θ 2 , µ, w 0 }, we take the maximum of ξ 2 within domain {(x, y)|0 ≤ x ≤ 1∩0 < y < 1}, which we denote as ξ 2 max . Then varying {N x , N y }, we finally obtain ξ 2 max as a function of grid points, as illustrated in Fig. 17(b). The key point is that even though there is no proof of the non-existence of the DeTurck soliton with Neumann boundaries, with sufficient Chebyshev grid points, we are still capable to obtain a solution that is arbitrarily close to that of Einstein equations.

D Derivation of the growth rate of entropy density
Recall that the entropy density of Planck branes penetrating the horizon is written as where i = 1, 2 denote the different Planck branes. In addition, the conserved momentum conjugated to v and the invariance of the reparametrization can be expressed respectively as where we have set L = 1.