Entanglement and Confinement in Coupled Quantum Systems

We study some general properties of coupled quantum systems. We consider simple interactions between two copies of identical Hamiltonians such as the SYK model, Pauli spin chains with random magnetic field and harmonic oscillators. Such couplings make the ground states close to the thermofield double states of the uncoupled Hamiltonians. For the coupled SYK model, we push the numerical computation further towards the thermodynamic limit so that an extrapolation in the size of the system is possible. We find good agreement between the extrapolated numerical result and the analytic result in the large-$q$ limit. We also consider the coupled gauged matrix model and vector model, and argue that the deconfinement is associated with the loss of the entanglement, similarly to the previous observation for the coupled SYK model. The understanding of the microscopic mechanism of the confinement/deconfinement transition enables us to estimate the quantum entanglement precisely, and backs up the dual gravity interpretation which relates the deconfinement to the disappearance of the wormhole. Our results demonstrate the importance of the entanglement between the color degrees of freedom in the emergence of the bulk geometry from quantum field theory via holography.


Introduction
Quantum entanglement provides us with various curious phenomena, such as the quantum teleportation [1]. From the point of view of holography, there is an interesting connection between quantum entanglement and wormhole [2,3]. Recently it has been proposed that, by coupling two quantum systems appropriately, it is possible to engineer a traversable wormhole [4]. This provides a dual gravitational description of the quantum teleportation [4,5] and also helps deepen our understanding of the puzzle of information loss during the evaporation of black holes [6,7,8,9].
A simple and explicit model that realizes the setting in this context is the "coupled SYK model" constructed in [10]. As we will review in Sec. 2, this model has several features of interest both from quantum mechanical and gravitational points of view. We would like to know whether results obtained within this model are generic, and if yes, try to give prescriptions useful for theoretical considerations and experimental studies. Therefore, we investigate analogous prescriptions for spin and bosonic systems, and study their properties. Furthermore we study coupled gauged systems with gauge singlet constraints in order to understand the connection with a wider class of models of quantum gravity on a firmer footing.
A key player in this context is the thermofield double state (TFD). In order to define the TFD, we introduce two copies of the identical Hilbert space -the "left" and "right" Hilbert spaces, denoted by H L and H R -on which two copies of identical Hamiltonian, denoted byĤ L andĤ R , act. The thermofield double at temperature T = β −1 is given by Here the sum runs over all energy eigenstates. The symbol * denotes the complex conjugate. The TFD state is a purification of the thermal state, in the sense that the reduced density matrixρ L = Tr R |TFD; β TFD; β| is the same as the thermal density matrix of the left system. When the quantum theory admits a dual gravity description, the TFD state is dual to the eternal black hole, which is a maximally entangled black hole connected by the Einstein-Rosen bridge [2,3]. In general, the TFD state is not the ground state of H L ⊗ 1 + 1 ⊗Ĥ R . The Einstein-Rosen bridge in the gravity dual, namely the eternal black hole, is not traversable. Alternatively, one can add a certain coupling between the two copies, such that the ground state of the coupled system mimics the TFD state of the uncoupled theory. Such systems provide us with useful setups for studying the TFD state both numerically and experimentally. For instance, the ground state can be studied more easily than the TFD using numerical techniques such as the Markov Chain Monte Carlo. Also, it appears possible to realize such system experimentally by engineering the coupled Hamiltonian, with for instance suggested protocols using quantum wires or graphene flakes bilayers [11].
When the system is in the regime that allows a gravitational dual description, the coupled model admits an interpretation as an eternal traversable wormhole [4,5,10].
In this work, we focus on this coupled system protocol [10], but we note that other ways for approximating the TFD state have also been proposed, e.g. [12,13]. Our first goal is therefore to understand which kind of couplingĤ int can achieve this objective for a variety of theories. This will be the topic of Sec. 2, Sec. 3 and Sec. 4. Our second goal, pursued in Sec. 5, Sec. 6 and Sec. 7, is to further improve the understanding about the relation between the quantum entanglement and spacetime. We will study how the loss of the entanglement and disappearance of the traversable wormhole are related. In the past, this problem has been studied for the SYK model [10]. We will consider gauge theories, for which precise calculations and simple gravity picture are available based on the knowledge about the confinement/deconfinement transition. 1 The detailed organization of this paper is as follows. In Sec. 2, we review the coupled SYK model, setting the strategy for finding the prescriptions for other systems. We additionally present numerical results for larger systems than previously available (up to N = 52 fermions for q = 4), as well as for a different value q = 8. In Sec. 3, we propose a prescription for spin systems, and check the validity by numerically studying a model of coupled spin chains in various parameter regimes. In Sec. 4, we study a similar prescription for the harmonic oscillators. We derive analytic expressions for the ground state, which are useful for later sections. We also perform a numerical calculation to demonstrate the loss of entanglement in the excited states. Technical details on the numerical computations performed in the above sections are given in the Appendices. In Sec. 5, the prescription considered in Sec. 4 is applied to matrix models. The gauged Gaussian matrix model is studied in Sec. 5.1, and it is explained how deconfinement and loss of entanglement are related. In Sec. 5.2, we briefly comment on the interacting matrix models. In Sec. 6, we perform similar analyses on the O(N ) vector model. The dual gravity interpretation is provided in Sec. 7. The notion of partial deconfinement [15,16,17,18,19] enables us to estimate the quantum entanglement precisely, and backs up the dual gravity interpretation which relates deconfinement to the disappearance of a traversable wormhole.

Coupled SYK Model
As a concrete example of a coupled Hamiltonian (2), we consider the coupled SYK models [10] by takingĤ Hereχ i L andχ i R are Majorana fermion operators acting on the left and right copies, respectively, which satisfy the anticommutation relation {χ i α ,χ j β } = δ ij δ αβ . The random couplings J j 1 ...jq are taken uniformly from a normal distribution of zero mean and variance We further set J = 1 as an energy scale. We used similar notations and normalizations as in [10] and [20] (the later uses the notation k for the coupling instead of µ). There are N 2 fermions on each copy; N is the total number of fermions in the coupled model.
From two Majorana fermions a Dirac fermion can be constructed. Let us useĉ j = It is straightforward to see that [Ĥ,P ] = 0. Note also thatQ mod 4 commutes withĤ; see e.g. Ref. [20]. This choice ofĉ j for the Dirac basis is useful for numerical calculations, as explained in Ref. [11], as the Hamiltonian is real for q = 4 and q = 8 and block diagonal with respect toQ mod 4.
We reproduce previously published overlap results [20,11] for q = 4 and extend them to larger systems, and also present new results for q = 8. First, we compute the ground state of the coupled Hamiltonian (2) and then compare it with the TFD state (1), which is obtained either directly from eqn. (1) or computed by |TFD, β ∝ e − β 4 (Ĥ L +Ĥ R ) |I . Thanks to improved sparse linear algebra methods and using the practical representation introduced in [11], we are able to push the calculations up to N = 52 Majorana fermions for q = 4 and N = 40 for q = 8. Details on the numerical methods are presented in Appendix A.  The results for the average overlap O (average over at least 50 realizations of disorder for each value of N, µ) are presented in Fig. 1. As expected, the overlap reaches unity in the two limiting cases of µ → 0 and µ → ∞. For intermediate values of µ, the overlap differs from unity, and this deviation increases with a system size. It would be natural to expect the overlap to be identically vanishing in the large-N limit at any non-zero, non-infinite µ, even though the numerical results above are not enough to conclude this. However it is quite remarkable that for systems with very large Hilbert space (i.e. with a relatively large value of N ), the TFD at the optimal β differs only by a very small amount from the ground state of the coupled model: at most 5.5% for q = 4 (up to dimH = 2 26 6.7 × 10 7 ) and 8% for q = 8.
Another interesting feature in Fig. 1 is that the overlap displays a minimum as a function of coupling strength µ, and the value µ * at which this occurs moves considerably with system size, in particular for q = 4. While for N = 28 and 32, we have µ * 0.15, 0.1 for q = 4 (as observed in previous studies [20,11]), we obtain µ * 0.05 for N = 48, 52. The largest sizes for q = 8 develop a minimum around a similar value µ * 0.04. However, it is hard based on our data to extrapolate to the thermodynamic limit behavior, in particular for q = 8 where the minimum µ * does not vary monotonously with system size. Based solely on the q = 4 data, it is possible that for larger N , µ * continues to decreases with a decreasing minimum overlap.
As we adopted the same normalization convention as in [10], we can directly compare our results (taking J = 1) as presented in Fig. 2. The agreement is very good when µ is sufficiently large, both for q = 4 and q = 8. We find a larger discrepancy in the small-µ region, which is below roughly the same scale where the overlap starts to decrease below unity, even when N is large. For q = 4 and for a range 0.2 µ, we attribute this to finite-N effects. Indeed in this range, we have been able to extrapolate our results to the thermodynamic limit using a 1/N fit for large-enough N (typically ranging from N = 32 to N = 40, see inset of Fig. 2): the extrapolated values are in good agreement with the analytical expression. For lower values of µ, we cannot extrapolate correctly the data due to the larger error bars caused by the greater fluctuations of β from sample to sample. This is also the region where the overlap starts to deviate much more significantly from unity, and therefore we expect the agreement to be less good. In fact, it is also possible to identify the deviation from unit overlap with the possible 1/q corrections that one has to take into account in the large-q analysis [10]. It is argued in [10] that such a deviation is due to the possible excitations due to the stronger left-right coupling to the coupled model so that the ground state of the coupled system is not accurate to mimic the TFD state at some temperature. 2 We finally note that the finite-N corrections are more complex to interpret in the q = 8 case. We observe that in the same range of values of µ where we could extrapolate the q = 4 data, the dependence of β(µ) on N is not always monotonous: for the very small set of data that we have access currently, β(µ) tends to first increase with N (in the opposite direction than the one expected from the analytical prediction), and then it tends to decrease again for larger N . Note that this effect is hardly visible on the scale of Fig. 2. We thus think that, despite the fact that the analytical expressions should be closer to the numerical data for large q, the range of N that we can reach with q = 8 is too small to be in an asymptotic regime where finite-N effects are easily interpreted.

Coupled Spin System
Next let us consider Hamiltonians consisting of Pauli spin operators. The Hilbert space consists of up and down spins at each lattice site. We do not specify for the moment any detail of the interaction; it can be local or nonlocal, and any lattice structure in any dimensions is allowed.
Our strategy is simply to mimic the coupled SYK model. We introduce the coupling such that, when the coupling is large, the TFD state at infinite temperature where BecauseΣ (| ↑ | ↑ * + | ↓ | ↓ * ) =Σ † (| ↑ | ↑ * + | ↓ | ↓ * ) = 0, the ground state at µ = ∞, which is obtained from GS|Ĥ int |GS = 0, is |GS = 1 2 L (| ↑ | ↑ * + | ↓ | ↓ * ) ⊗L/2 indeed. We turn to concrete numerical experiments to check whether this simple type of coupling plays a similar role as in the coupled SYK models. We study one-dimensional spin chains with random magnetic field along the z-direction: where the random magnetic field is chosen uniformly random in [−W, +W ] and α = L, R denotes the left / right chain. This model offers the opportunity to probe the adequateness of the TFD to mimic the ground state of different phases of coupled matter. For W = 0, the model (Heisenberg model) is integrable. For 0 ≤ W < W c , the system is in an chaotic/ergodic phase, while for W > W c , the system is in a Many-Body Localized (MBL) phase. Current estimate of the critical disorder (in the middle of the spectrum) is [22] W c 3.7. It is important to remark that the nature of the ground state of the coupled system does not necessarily reflect the underlying phase of the uncoupled system. This is indeed the case for the disorder-free model (W = 0) where the coupled system is known to be a gapped paramagnet for any µ > 0 [23] (whereas the uncoupled system is a critical integrable liquid). While we are not aware of specific predictions for coupled systems for W > 0, we expect a similar behavior (gapped paramagnetic ground state) for all systems in the limit µ 1. We study the coupled model Eq. (2) with the leftĤ L and rightĤ R taken with the same disorder realization of (13). As in the SYK case, we measure the average overlap O between the ground state of the coupled system and the TFD. We additionally compute the Kullback-Leibler (KL) divergence KL = Trρ L logρ L −ρ L logσ(T ), wherê ρ L = Tr R (|GS GS|) is the reduced density matrix andσ(T ) is the thermal density matrix of the uncoupled theory. The KL divergence is another measure of how different the TFD and the coupled ground state are: it vanishes in case they are equal. We compute both the inverse temperature β(O), which maximizes the overlap, and β(KL), which minimizes the KL-divergence. We average our results over more than 100 realizations of disorder for each data set (W = 0,L).   Our numerical results on the overlap are summarized in Fig. 3. From the results, we observe the same qualitative tendency, irrespective of whether the uncoupled system is located in the ergodic phases, the MBL phases (small and large W , respectively) or integrable (W = 0). The data are very similar to what is observed for the coupled SYK model: the overlap is close to unity in both limits µ → 0 and µ 1, and displays a deviation from unity which increases with the system size (modulo a non-trivial effect for small disorder strenghts depending on the even-odd parity of L/2). This is corroborated by considering the small values of the KL divergence both in the ergodic and MBL cases of the uncoupled system: see bottom panel of Fig.4, where data for disorder W = 1, 8 are presented. The KL is maximal in the range where the overlap is minimal, as naturally expected, but its scale is overall quite small, considering the sizes of the vectors involved. The left top panel of Fig.4 presents the inverse temperature β as a function of coupling strenght for these two values of disorder: again the data is very similar to the SYK case, even though there is no precise analytical prediction to compare with in the spin case (we also expect that β scales as 1/µ in the large µ limit, based on the dimensional counting).
We close these numerical observations by showing in the right top panel of Fig.4 that the inverse temperatures as determined either by minimizing the KL divergence or the maximization of the overlap are very similar: the relative difference is at most a few percent, and non-negligible difference appears only in the very low µ regime where the error bars on β are larger.
The conclusion of this numerical study of the coupled spin chains system is that the results appear qualitatively very similar to the ones obtained for the SYK model. In particular, one can find the inverse temperature β(µ) such that the ground state of the coupled model is very close to the TFD, at any value of the coupling strength µ. Moreover, this appears to be the case irrespective of the underlying nature of the uncoupled system (even if the spin chain system is not as chaotic/ergodic as the SYK model). This may not be so surprising after all since the inclusion of the coupling can clearly change the nature of the ground state, as we mentioned previously.

Coupled Harmonic Oscillators
In the previous sections we considered models whose fundamental degrees of freedom are fermionic (the 1d spin chain can be recasted in terms of a fermionic chain). In the following we consider models with bosonic degrees of freedom. In order to construct coupled bosonic systems with interesting properties, we start with the simplest but very instructive example, i.e. the coupled harmonic oscillator [24], with the following uncoupled Hamiltonian Following the strategy adopted for the SYK model and spin chain, we introduce a coupling between two identical copies of the hamonic oscillators such that the ground state of the coupled Hamiltonian mimics the TFD state at infinite temperature. The coupling we introduced is parameterized by two independent coupling constants C 1 and C 2 : The physical intuition that suggests a direct connection between the ground state of this model and a thermofield double state is the following. We first set C 2 to zero and   vary C 1 , the coupled model is equivalent to a pair of uncoupled oscillators with frequency √ ω 2 − 2C 1 . As C 1 approaches ω 2 2 , the wave function spreads out, and the ground state approaches to Next we turn on C 2 to a large value. Then x L and x R are forced to become close, and hence the normalized ground state of the coupled system becomes close to ∞ −∞ dx|x |x . This is nothing but the TFD state at infinite temperature. In the rest of this section we provide a quantitative analysis that verifies such intuition explicitly.
For the convenience of later analysis, we reparametrize the deformation terms such that Notice that in principle we could separate out the C + −C − 2 (x 2 L +x 2 R ) piece in the left-right coupling from the genuine "interaction" termx LxR . We choose not to do so because, as we will show below, the ground state of the coupled system is identical to a TFD state if the coupling satisfies 1 + 2C to be the 'original' frequency, the ground state would then be the TFD state of the 'shifted' frequency ω. 4

Ground state as the TFD
The Hamiltonian of the uncoupled theory is given by Eq. (14), with the canonical commutation relation [x,p] = i. The creation and annihilation operators are defined byâ † =x −iωp √ 2ω andâ =x +iωp √ 2ω . The Hamiltonian can be written in terms of the number operatorn =â †â asĤ = ω n + 1 2 . The vacuum |0 is defined byâ|0 = 0, and the normalized excited states are constructed as |n =â †n √ n! |0 . The coupled model is defined by (16), which can be rewritten aŝ where Note that both ω + and ω − are different from ω for any nonvanishing C ± . The creation operators areâ † where r ± = ω ± ω . The ground state that satisfiesâ + |0 coupled =â − |0 coupled = 0 is where the normalization factor N is given by The ground state (20) can be rewritten as where This is a TFD state when A 1 = 0, or equivalently With this condition, A 2 simplifies to the following expression: We can rewrite this ground state into a form that resembles the TFD state in a manifest manner where and From this expression we can read off the effective temperature for any given C ± that satisfy (24) By construction, A 2 can take the values between −1 and +1. At A 2 = ±1, ω + or ω − becomes zero, and the theory becomes ill-defined with a continuous spectrum; the effective temperature T eff becomes infinite there. For all the other values of A 2 , we have shown that the ground state of the coupled model (16) satisfying (24) is identical to a TFD state at temperature (29). Such a TFD state is analogous to the one originally discussed in [2]. Further notice that all the above discussion assumes ω 2 ± > 0, and together with the condition (24) they are the incarnation in this model of the fact that there is only one good choice of the sign of the coupling in order to make the wormhole traversable in the general discussion [4].

Decay of entanglement at high temperature
In this section, we consider the excited states of the coupled model, instead of the ground state, where the coupling constants still satisfy 1 + 2C + /ω 2 1 − 2C − /ω 2 = 1. It is expected that the quantum entanglement between the two sides is washed away when sufficiently large amount of energy is added to the system. We check this expectation quantitatively in this section.
We consider the excited state where the energy of the system is We would like to study how does the entanglement between the left and right sectors of the system, which is in the state (30), decay at high temperature/high energy. In the following, we carry out some of our computations numerically imposing a cutoff Λ to the Hilbert space. For details of the numerical methods, see Appendix A.2.

Overlap with TFD
We can calculate the overlap with the TFD state, For each (n + , n − ), the value of β is chosen so that the overlap is maximized. The results for ω = 1, C + = 0.1 and 1.0 are shown in Fig. 5. We can see that the overlap with TFD becomes smaller as the energy (or equivalently the temperature) is increased. Note that the temperature of the TFD state which maximizes the overlap depends on the given excited state, and except for the ground state this temperature is rather high. Note also that (n + , n − ) = (a, b) and (b, a) give the same overlap.

Mutual Information
With the general discussion [4] and the particular example of the SYK model [10] in mind, we would like to relate the correlation of the two sides with the existence of possible wormhole phases. As known previously [25,26,27], the left-right propagator is not alway a good diagnose of this connection. In the rest of this section, we consider a refined mutual information, namely S EE,L + S EE,R − S diag , to probe the quantum correlation of the left and right sides and demonstrate that the connection between the coupling of the two sides and the quantum entanglement between them. The mutual information (MI) is defined by Here S EE,L and S EE,R are the entanglement entropy obtained from the reduced density matricesρ L,R = Tr R,Lρ , while S therm is the thermal entropy of the entire system. As we can see from Fig. 6, the MI does not decay significantly, even when C + is as small as 0.1. It is because the MI picks up both the classical and quantum correlations. At high temperature, the MI should be dominated by the classical correlation. On the other hand, the decrease of MI at low temperature reflects the decay of the quantum entanglement. In order to study the quantum correlation more clearly, we need to subtract the classical correlation. To do so, recall that the thermal density matrix of the coupled system can be written in the following form: When the coupling parameters C + and C − are small, tiny off-diagonal elements are generated and contribute to the entanglement. If we keep only the diagonal part and define a separable stateρ diag = n L ,n R ρ n L ,n R ;n L ,n R (|n L n L |) (|n R n R |) , it should not capture the entanglement, rather it should capture the classical thermal correlation between left and right sectors. If S EE,L + S EE,R − S thermal is not different from S diag − S thermal , where S diag = −Tr (ρ diag logρ diag ), it is natural to expect that the MI is not picking up the quantum entanglement that is not due to thermal effects. Therefore we use the quantity to characterize the quantum correlation between the left and right sides. As we can see from Fig. 7, the difference S EE,L + S EE,R − S diag decays significantly when C ± are small. We interpret it as the evidence that the quantum entanglement decays as temperature goes up. Note that S diag manifestly depends on the choice of the basis of the Hilbert space. It is possible that a more elaborate choice of the basis could lead to a better estimate of the classical correlation. This is indeed aligned with the fact that there does not seem to be a canonical measure of the classical contribution to general entanglement entropy. We consider the quantity (36) because the entropy so defined monotonically decreases as temperature increases, as shown in Fig. 6, which is what we expect since raising temperature generally destroys quantum entanglement.

4.3
The case of C + = C − ρ diag , S diag and the overlap with the TFD state explicitly depend on the choice of 'uncoupled' and 'interaction' parts, unlike the entanglement entropy and mutual information. Therefore, let us consider another natural example: (16) In this case, the ground state of the coupled system is not the TFD state in the L-R basis, since we can check that (24) is not satisfied in this case for all C > 0. The system is well-defined when ω 2 ± = ω 2 ± 2C ≥ 0, namely − ω 2 2 ≤ C ≤ ω 2 2 . As we can see from Fig. 8, the ground state is close to the TFD as long as the coupling does not become too large, and S EE,L + S EE,R − S diag is small.

Coupled matrix models
In this section, we apply the results obtained in Sec. 4 to gauged matrix models, which is closely related to physics of black holes via gauge/gravity duality. The models exhibit a deconfinement transition, which corresponds to the formation of black hole in the gravity side via holography.

Gauged Gaussian matrix model
Let us consider the simplest, analytically calculable example: the gauged Gaussian matrix model. The Euclidean action is given by The covariant derivative D t X I is defined by where A t is the gauge field. When the gauge field is integrated out, the gauge-singlet constraint emerges. Other than the singlet constraint, this system is nothing but a bunch of non-interacting DN 2 harmonic oscillators with m = ω = 1. The coupled version is We use the same gauge field for the left and right copies, so that X and Y transform as the adjoints under the same SU(N ) gauge group, because otherwise the coupling term is not gauge invariant. Having the results in Sec. 4 in mind, let us take C − such that is satisfied. Then the ground state can be interpreted as a product of TFD's discussed in Sec. 4. The gauged Gaussian matrix model exhibits the confinement/deconfinement transition because of the gauge-singlet constraint (see e.g. [19,17]). The free energy of the original, single-copy uncoupled theory is 6 where x = e −β and u n = 1 N TrP n , where P = diag(e iθ 1 , · · · , e iθ N ) is the Polyakov line. Strictly speaking, this is the effective action in terms of θ 1 , · · · , θ N ; the free energy is obtained by minimizing it with respect to θ's. The term N 2 D 2 log (det (−D 2 0 + 1)) is the contribution from D scalars, while N 2 2 log (det (−D 2 0 )) is associated with the gauge fixing.
There is a first order phase transition at T c = 1 log D , where |u 1 | jumps from 0 to 1 2 ; see Fig. 9 (It cannot go beyond 1 2 because the density distribution ρ(θ) = 1 2π (1 + 2u 1 θ) must not be less than zero). This leads to the first order phase transition without hysteresis. At T > T c , u 2 , u 3 , · · · become nonzero as well, while u 1 becomes larger, so that the free energy is minimized while ρ(θ) remains non-negative.
If we consider the microcanonical ensemble (i.e. use the energy as a parameter, rather than the temperature), a rich structure can be found at the phase transition. The key concept is the partial deconfinement [15,16,17,18,19]: an SU(M ) subgroup of SU(N ) deconfines, and M N increases from zero to one as energy grows, as E ∝ M 2 . The value of M is M = 2N |u 1 | in this case. When the SU(M ) subgroup is deconfined, N 2 − M 2 degrees of freedom remain confined, namely they remain as the ground state.  (38). Blue, orange and red lines are identified with the confined, partially deconfined and completely deconfined phases, respectively. These figures are taken from Ref. [17].
For the coupled model, we obtain where ω = √ 1 + 2C + , x = e −βω and x = e −β/ω . As temperature is raised, the deconfinement phase transition takes place at 1 − D(x + x ) = 0. The critical temperature is T c = 1 log 2D for C + = 0 and T c = 0 for C + = − 1 2 , ∞. The theory is ill-defined at C + < − 1 2 , because the energy is not bounded from below. At T = T c , the coefficient in front of |u 1 | 2 becomes zero, and hence, |u 1 | can take any value between 0 and 1 2 . In the confining phase (T < T c ), the system is indistinguishable from the ground state up to the 1/N corrections. Therefore, the confining phase should be the TFD, up to the 1/N corrections. Associated with the deconfinement, the harmonic oscillators are excited, and hence, the quantum entanglement decreases. Hence we expect the phase diagram shown in Fig. 10. There is one subtlety though; when T c is small (i.e. C + → − 1 2 or ∞), the jump of the energy is also small, while the entanglement in the ground state is large. Therefore, the entanglement cannot be washed away immediately. The same holds also when D is large. It may be an artifact of the free nature of the theory. Just like the uncoupled theory, we obtain an interesting phase diagram with SU(M )deconfinement when we consider the micorocanonical ensemble. When the SU(M ) subgroup is deconfined, N 2 − M 2 degrees of freedom remain confined, namely they remain as the ground state. Therefore even when the coupling parameter is small, the quantum entanglement survives until all the degrees of freedom deconfine, see Fig.11. The amount of the entanglement can easily be estimated by counting the number of confined degrees of freedom: Here we have assumed that the coupling is small and the entanglement in the deconfined sector is washed away by the thermal excitation. Because E ∝ M 2 , we can also express the amount of the entanglement as where E deconf is the energy needed for the complete deconfinement. Note that the estimate above is different from the naive entanglement entropy at finite energy; we have omitted the contamination coming from the nonzero entanglement entropy in the deconfined sector which does not correctly measure the quantum entanglement. In Sec. 7, we will consider the dual gravity description based on this estimate. Figure 11: Partial deconfinement phase in the matrix model, at small but finite deformation parameter C + 1. When the energy is not large enough such that all degrees of freedom are deconfined, only a part of degrees of freedom deconfine. The deconfined sectors in the left and right Hilbert spaces are not entangled, while the confined sectors are entangled.

Gauge-invariance in terms of Hilbert space
The states in Hilbert space have to satisfy the gauge-singlet constraint. This is automatic if we construct the states by acting gauge-invariant operators to the gauge-invariant vacuum.
The gauge transformation is generated bŷ where f αβγ is the structure constant of the SU(N ) algebra. By usinĝ we can rewrite it asĜ Hence the ground state of the uncoupled theory |0 L |0 R , or more explicitly ⊗ I,α (|0 LIα |0 RIα ), is gauge-invariant:Ĝ we obtainĜ Hence the ground state of the coupled Hamiltonian |0 + |0 − is also gauge-invariant.

Yang-Mills Matrix Model
As a concrete example with interaction, let us consider the Yang-Mills matrix model. The action of the uncoupled model we consider is where I and J run from 1 to D. The simulation data [31] is consistent with the partial deconfinement. The coupled version is We can take only C I+ ≥ 0 and C I− ≤ 0, because otherwise the potential is not bounded from below (if we introduce the mass term, for example by considering the plane-wave deformation [32], C I+ < 0 and C I− > 0 are allowed). We expect that the X-Y coupling introduces entanglement between X and Y sectors, although it is not easy to see analytically if the TFD naturally appears or not. Below the deconfinement transition there is no temperature dependence up to the 1/N -corrections, and hence, the entanglement should survive up to the deconfinement temperature. When the X-Y coupling is not too large, the jump of the energy at deconfinement transition should be large enough to eliminate the entanglement in the ground state. Note that, as in the case of the Gaussian matrix model, the entanglement can survive until all the degrees of freedom deconfine.

Coupled vector model
Let us consider the coupled free vector model on a three-dimensional sphere, with the gauge-singlet constraint. We introduce N -component vectors φ f , where f = 1, · · · , N f are the flavor index. In the path integral formalism, the singlet constraint can be imposed by introducing the Chern-Simons gauge field [30]. Non-constant modes along the sphere have bigger 'effective mass' ∼ √ m 2 + ∆ 2 , where m is the mass and ∆ 2 is the eigenvalue of the Laplacian. Due to the lack of the interaction, each mode behaves as harmonic oscillator, whose frequency is the same as the 'effective mass'. Therefore we can use the findings of Sec. 4. The left-right coupling C f makes each pair the TFD with certain shifted frequency and effective temperature. Note that the values of the shifted frequency and effective temperature depend on the pair. The ground state is the tensor product of such TFDs.
When m and coupling parameters are much smaller than the curvature of the sphere, the phase structure is close to the massless theory with 2N f flavors. There is a Gross-Witten-Wadia (GWW) type phase transition [30] at T = 3N 2π 2 N f [30], which separates the partially deconfined phase and completely deconfined phase [17].
As in the case of the matrix models discussed in Sec. 5, the entanglement survives all the way up to the GWW transition. It is rather surprising that the entanglement survives up to such high temperature, which is of order √ N . When the temperature is T = (54)

Geometric interpretation
In this section, we propose connections between the mechanism of deconfinement, loss of the entanglement and the dual gravity interpretation (the disappearance of the traversable wormhole). As mentioned in Sec. 5 and Sec. 6, it has been proposed [15,16,17,19] that the deconfinement transition takes place gradually, namely there is a partially deconfined phase in which a part of the SU(N ) gauge group, SU(M ), is deconfined and the rest is still confined. The size of the partially deconfined sector M increases gradually toward N . Via AdS/CFT duality, the natural candidate of the gravity dual of the partially deconfined phase is the black hole with negative specific heat, which is smaller than the AdS scale [15,17]. While the proposal has been made for the usual, uncoupled gauge theories, it is natural to expect the same mechanism for the coupled theories, at least when the coupling term is small. Then what would be the dual gravity interpretation?
The coupled SYK model has a first order transition similar to the deconfinement, and an intermediate state, which resembles the partially deconfined phase, exists [10]. The dual gravity interpretation proposed in [10] is that the low-temperature phase (∼ confinement phase, strong entanglement) is dual to the traversable wormhole, the high-temperature phase (∼ completely deconfined phase, (almost) no entanglement) is dual to two separate black holes without wormhole, and the intermediate phase (∼ partially deconfined phase) describes gradual formation of black hole and disappearance of the traversable wormhole.
If this picture is indeed true for the gauged model, it is confirming the following gravity interpretation [15,16,17]: in the partially deconfined phase, the deconfined SU(M )-sector describes a small black hole, while the confined sector describes the rest of the geometry in the dual gravity theory. When M N 1, the black hole is so small that the wormhole geometry is not much affected. As M increases, the geometry is gradually filled by the black hole and the wormhole becomes thinner. In the completely deconfined phase, all degrees of freedom are excited, and the two copies of deconfined phases are not entangled. In the language of gravity, the dual black hole is so large that it is almost completely filling AdS and the traversable wormhole is destroyed.

Conclusion
In the introduction, we have set two goals of this paper: (i) introduce appropriate couplings such that the the ground state of the coupled Hamiltonian mimics TFD, and (ii) improve the understanding about the relation between the quantum entanglement and spacetime. As for the first goal, we have confirmed the validity of the proposal for the SYK model in the reference, and generalized it to the Pauli spin systems and bosonic systems including matrix model and vector model. As for the second, we pointed out the partial deconfinement enables us to estimate the entanglement at finite temperature and leads to a natural geometric interpretation. The wormhole in the gravity side arises from the entanglement between the color degrees of freedom in the QFT side. We regard this mechanism as a good demonstration of the importance of the entanglement between the color degrees of freedom in the emergence of the bulk geometry from quantum field theory via holography. Institute for Theoretical Science at the University of Chinese Academy of Science, Tsinghua University, the Songshan Lake Materials Laboratory, and The Weizmann Institute of Science for warm hospitality during the various stages of the project. The work of  [39,40,41] and SLEPc [42] libraries.

A.1 Coupled Spins and SYK models
For the spin chain and SYK models, we use similar exact diagonalization techniques. The size of the Hilbert space for the coupled system is 2 N/2 for SYK in the standard complex Dirac fermion basis. It was shown in [11] that one could use an adapted Dirac representation where the Q mod 4 symmetry of the coupled system is explicit: the Hamiltonian is block-diagonal with 4 blocks, reducing the dimension of the Hilbert space by about 4. The ground state of the coupled system is always located in the Q mod 4 = 0 sector for any µ = 0. On top of this, the Hamiltonian becomes real in this basis which furthermore eases computations by limiting memory requirements. For the coupled spin chains, the coupling preserves the U(1) symmetry of each chain (the magnetization σ z α = L/2 i=1 σ z i,α is conserved for α = L, R), and furthermore, we find that the coupled ground state is always in the sector where these two magnetizations are of opposed sign. The total Hilbert space size in this sector is then To reach systems larger than those accessible with full diagonalization (limited to L ≤ 16 for the spin chain, N ≤ 32, 36 for SYK), we take advantage of the sparse nature of the matrix representation of the Hamiltonians. Indeed, each matrix line contains only on average ∝ 3L/4 (∝ (N/2) 4 /192 in the complex basis) non-zero elements for the coupled spin chain (SYK) model. Note the favorable prefactor for the coupled SYK system compared to a single SYK model which has ∝ (N/2) 4 /24 non-zero matrix elements per line. This overall small number of non-zero matrix elements permits to use sparse linear algebra iterative algorithms such as the Lanczos algorithm to obtain the ground state of the coupled system. We can reach systems of size up to L = 32 spins for the spin chain (even though we only present data up to L = 28 for computational resources reasons), and up to N = 52 for the SYK model. To reach these large systems, we need to use large-scale parallel computations, which is possible with sparse iterative methods. For the computation of the KL divergence in the case of the spin system, we furthermore perform a singular value decomposition (SVD) of the coupled ground state once it has been obtained. For this, we find it beneficial to construct the basis of the coupled system using Lin tables [33] for the right and left subsystems respectively.
The thermofield computation can become also very costly for large L, N . Indeed, in the method detailed in Ref. [20], one first need to compute all eigenstates |E of a single system with N/2 particles (which is easy for the considered sizes using full diagonalization) and then perform the outer product between all these eigenstates in Eq. 1, which is the most demanding part for large systems, before computing the overlap with the coupled ground state. We do not keep in memory all matrix elements of outer-product but only compute the overlap line-by-line (as explained in Ref. [20]). The outer-product is performed using optimized shared-memory parallel routines of the BLAS library. The U (1) symmetry of the spin chain eases this computation as only corresponding sectors of magnetization of the single spin chain (with L/2 spins) are matched to form the TFD.
For the coupled SYK model, we did not find an efficient way for computing the TFD double using the full diagonalization of a single SYK model using the real representation of Ref. [11] (which is needed to reach large N for the ground state of the coupled model). There we use a different method, also taking advantage of sparsity. As emphasized in Ref. [10,11], the TFD at β = 0 has a simple representation (a single basis state |I = |00000 in the coupled basis of Ref. [11]). The TFD at finite β can be obtained (up to normalization) by applying exp(−βH(µ = 0)/4) to this state. While this requires to use the full coupled Hamiltonian (even though systems are not coupled), we can use Krylov expansion techniques [34] to perform this application, without the need of forming explicitly the matrix exp(−βH(µ = 0)/4). In practice, all the computations are done in the same run: we start by computing the ground state of H(µ), store it in memory, and modify the diagonal of the matrix to set µ = 0. We then apply recursively exp(−δβH(µ = 0)/4) to |TFD(β = 0) in small steps (typically chosen as δβ = 1/(100µ)). Starting from β = 0 (where the overlap with the coupled ground state is minimal), the overlap grows as β increases until it reaches a maximum, which is the value we seek. We stop the iterative application of exp(−δβH(µ = 0)/4) as soon as we observe a drop in the overlap.
We compute the maximal overlap for each disorder realization (of h i or J j 1 ···jq ), then average over disorder to obtain the results presented in the main text. We use from 30 to 1000 realizations of disorder (depending on the system size) for each value of µ.

A.2 Coupled harmonic oscillators
We writeρ thermal explicitly in the L-R basis, with the explicit cutoff, n L , n R , n ± < Λ. Then the dimension of the truncated Hilbert space is Λ 2 .
By using them, we can constructâ ± andâ † ± explicitly. By utilizing the sparseness, we can restrict the cost of the multiplication ofâ ± andâ † ± on a generic state to be O(Λ 2 ). Note that, if we simply define |n L , n R =â †n L Lâ †n R R √ n L !n R ! |0 coupled , the norm can deviate from 1 due to the cutoff effect. We multiplied a positive number so that the norm becomes 1. They are not exactly orthogonal, due to the regularization effects.
The thermal entropy of the total system can easily be calculated analytically, so we did not use the numerical method in the results shown in this paper. 7 In order to obtain the entanglement entropy, we constructed the reduced density matrix numerically, calculated the eigenvalues of the reduced density matrix p 1 , p 2 , · · · , p Λ , and then determined the entropy as − i p i log p i . In order to avoid numerical singularities, we have omitted the eigenvalues smaller than 10 −15 .  7 We used it to check the debugging, namely we checked that the entropy of the total system is correctly obtained.

A.2.1 Comparison with analytic results
Thermal entropy of the coupled harmonic oscillators is where ω + = ω 2 + 2C + , ω − = ω 2 − 2C − . To make the ground state to the thermofield double state, we take r + = 1 r − , where r ± = ω ± ω . The entanglement entropy of the ground state is where T eff is defined by T eff is zero when C + = 0, and diverges as C + → ∞ or C + → − 1 2 . We can use these relations to check the validity of the numerical calculations. In Fig. 13, we have plotted the thermal entropy calculated numerically and the analytically. We can observe a good agreement, when the cutoff Λ is sufficiently large.