Matrix Entanglement

In gauge/gravity duality, matrix degrees of freedom on the gauge theory side play important roles for the emergent geometry. In this paper, we discuss how the entanglement on the gravity side can be described as the entanglement between matrix degrees of freedom. Our approach, which we call 'matrix entanglement', is different from 'target-space entanglement' proposed and discussed recently by several groups. We consider several classes of quantum states to which our approach can play important roles. When applied to fuzzy sphere, matrix entanglement can be used to define the usual spatial entanglement in two-brane or five-brane world-volume theory nonperturbatively in a regularized setup. Another application is to a small black hole in AdS5*S5 that can evaporate without being attached to a heat bath, for which our approach suggests a gauge theory origin of the Page curve. The confined degrees of freedom in the partially-deconfined states play the important roles.


Introduction
The importance of quantum entanglement is getting more and more widely appreciated. In the context of quantum gravity, a proper understanding of quantum entanglement is important for the resolution of the black hole information puzzle. (See e.g., Refs. [1,2] for recent reviews.) Furthermore, quantum entanglement might be a key to understanding the emergent space in holography [3]. In this paper, we would like to introduce a new tool having such a big picture in mind.
In the Ryu-Takayanagi prescription [4] for holographic entanglement entropy, quantum field theories with nonzero spatial dimensions are assumed. Then, the space of dual quantum field theory is split into two parts. Let us call entanglement defined in this way as 'spatial entanglement'. Spatial entanglement is not the only way to define quantum entanglement. between these objects.
For such a definition of entanglement to make sense, we have to find a natural way of splitting matrix degrees of freedom into two or more pieces. Such a partition is obtained by understanding the meaning of 'eigenvalues' in gauge theory properly [23]. An important example is the partially-deconfined states dual to small black hole, in which an U(M )-subgroup of U(N ) (0 < M < N ) is deconfined [24][25][26][27][28]. The confined and deconfined sectors in the partially-deconfined states provide us with a natural partition, and gauge fixing can be regarded as a superselection [27]. One of our motivations is to define the entanglement between black hole and Hawking radiations solely in terms of dual QFT or matrix model without referring to gravity, as a tool to study quantum gravity based on the first principles.
This paper is organized as follows. In Sec. 2, we explain how the geometry on the gravity side can be encoded into matrix degrees of freedom in the matrix model or QFT. It is important to note the relationship between the gauge-invariant Hilbert space H inv and the extended Hilbert space H ext that contains gauge-non-singlet states, discussed in Sec. 2.1. The latter is suitable for geometric interpretations. In Sec. 2.2, we use wave packets in the extended Hilbert space to encode geometry into matrix degrees of freedom. Based on such a geometric interpretation, the matrix entanglement is introduced in Sec. 3. For a class of quantum states we consider, geometric interpretation becomes clear after gauge fixing, and we can define the matrix entanglement in a straightforward manner. In Sec. 4, we consider the matrix entanglement associated with the fuzzy-sphere background of the BMN matrix model. In appropriate large-N limits, 3d SYM or 6d N = (2, 0) SCFT are realized, and matrix entanglement can be used to regularize the spatial entanglement in these QFTs. In Sec. 5, we consider the matrix entanglement for the evaporating small black hole. We show that the confined degrees of freedom play an important role. Sec. 6 is devoted for conclusions and discussion of open problems.

Geometry encoded into matrix degrees of freedom
In this section, we explain how the bulk geometry can be encoded into matrices, following Ref. [23]. We consider the matrix model as a starter and then discuss the generalization to QFT.
We consider the BFSS matrix model [29] and simpler matrix models. One of the simplest is the bosonic part of the BFSS matrix model, I, J run from 1 to 9. The details of models do not matter unless we consider quantitative correspondence with gravity. Later, we will also consider the Gaussian matrix model as an even simpler model which can be solved exactly. Although the details of the Hamiltonian is not important, for concreteness let us consider the Hamiltonian interpolating the Gaussian matrix model and the bosonic part of the BFSS matrix model: The Gaussian matrix model is obtained at g 2 = 0, and in the strong coupling limit g 2 → ∞ the mass term is negligible and the bosonic part of the BFSS matrix model is obtained.
The mass term is analogous to the conformal mass term in 4d maximal super Yang-Mills compactified on three-sphere. By using the generators τ α of U(N ) normalized as Tr(τ α τ β ) = δ αβ , the matrix-valued operators are expressed aŝ with the canonical commutation relation

Gauge-invariant Hilbert space and extended Hilbert space
To understand the matrix model in terms of quantum states, it is important to appreciate the difference of gauge-invariant Hilbert space H inv consisting only of gauge-singlet states and extended Hilbert space H ext that contains non-singlet states as well. We can use the coordinate eigenstates |X that satisfyX I,ij |X = X I,ij |X for all 9N 2 combinations of (I, i, j), as the basis of H ext . Equivalently, we can use the momentum eigenstates |P as the basis. Note that operatorsX I,ij andP I,ij are defined on this extended Hilbert space. Later in this paper, we will use wave-packet states as an over-complete basis of H ext . The use of the extended Hilbert space is naturally justified if we use the projector from H ext to H inv that corresponds to the integration over Polyakov loop [23,28]. The partition function at finite temperature is written in two different but equivalent ways. The first expression is Based on this expression, people often say 'physical states are gauge-singlets'. The second expression uses the extended Hilbert space that contains non-singlet modes: Here G = U(N ) is the gauge group, g is a group element, andĝ is the representation of g acting on the extended Hilbert space H ext . The Haar measure is used for the integral. As shown in Appendix A, the second expression is directly related to the path integral with temporal gauge field A t , and g ∈ G corresponds to the Polyakov loop. In this second expression, we can say that 'states connected by a gauge transformation are identified', i.e., |φ andĝ|φ are indistinguishable, just as we do in the classical theory or in the path-integral formulation. Note thatP is a projection operator from H ext to H inv that satisfiesP 2 =P. By usingP, (5) can also be written as From a state |φ ∈ H ext , which can be non-singlet, we can obtain a singletP|φ ∈ H inv . By using this correspondence, we can describe the same physics by using the singlet states or non-singlet states. BecauseP andĤ commute with each other, this singlet-projection and Hamiltonian time evolution commute.

Geometry from wave packet in R 9N 2
Low-energy states can be written as superpositions of wave packets. We consider a wave packet around Y I in coordinate space and Q I in momentum space. Such a wave packet naturally describes geometry formed by D-branes and open strings [23,30], along the lines of Ref. [31]. A natural way to obtain a low-energy wave packet is to find |Y, Q that satisfies constraints The state minimizing the energy is of particular interest. Note that such a |Y, Q is an element of the extended Hilbert space H ext , but not necessarily in the gauge-invariant Hilbert space H inv . It transforms as via the gauge transformation; see Fig. 1. States that transform to each other via gauge transformation should be identified. Note that the shape of the wave packet does not change via gauge transformation. The only 'matrices' that can be diagonalized are Y I and Q I . Each wave packet smoothly spreads in the coordinate space R 9N 2 , and also in the momentum space R 9N 2 . Along each of 9N 2 dimensions, the size is at most of order N 0 . Two wave packets centered around Y and Y in R 9N 2 have vanishingly small overlap if the centers of wave packets are well separated, i.e., Tr(Y − Y ) 2 1. See Fig. 2.
The simplest example: Gaussian matrix model and coherent state As an instructive example, let us consider the Gaussian matrix model that is obtained by setting the coupling in the toy model Hamiltonian (2) to zero: Namely, we consider 9N 2 harmonic oscillators subject to the U(N ) gauge symmetry. The ground state is simply the tensor product of 9N 2 Fock vacua, This ground state is U(N )-invariant. The wave packet is defined above is the coherent state, Figure 1: Wave packets in the coordinate space R 9N 2 . Black points are the centers of wave packets (Y , Y (U ) and Y (U ) ), and gray disks are the wave packets (more precisely, the regions where the wave functions are not vanishingly small). Under the gauge transformation, the center of the wave packet moves as Y → Y (U ) = U −1 Y U , but the shape of the wave packet does not change. The center of the wave packet gives 'slow modes' that describe the geometry consisting of D-branes and strings. Figure 2: Each wave packet smoothly spreads in R 9N 2 . Along each of 9N 2 dimensions, the size of the wave packet is at most of order N 0 . Two wave packets centered around Y and Y in R 9N 2 have vanishingly small overlap if the centers of wave packets are well separated, i.e., The ground state is a wave packet about Y = 0 and Q = 0, Under the gauge transformationX I ,P I → UX I U −1 , UP I U −1 , this wave packet transforms as The wave function is well localized around Y I in the coordinate space and Q I in the momentum space. For example, when Q I is zero, The overlap between two wave packets is They have vanishingly small overlap if the centers of wave packets are well separated, i.e., Tr(Y − Y ) 2 1. In typical situations considered in holography, Therefore, typical separation is much larger than 1, and overlap between typical wave packets is practically negligible.

Geometric interpretation
As explained in Refs. [23,30], matrices Y I are related to the geometry consisting of D-branes and open strings via Witten's interpretation [31]. a,b Let us list a few features justifying this: • Ref. [31] considered low-energy states, which should be wave packets. The only 'matrices' associated with a wave packet that can be diagonalized are Y I and Q I .
• Suppose all 9 matrices can be diagonalized simultaneously: Around this background, fluctuation of (i, j) component has mass proportional to | y i − y j |. c This allows the interpretation that y i is the location of the i-th D0-brane and off-diagonal entries describe the open strings stretched between D-branes.
a The original proposal [31] concerns low-energy effective description of D-branes and open strings, and did not claim that such a interpretation is valid in gauge/gravity duality; actually Ref. [31] appeared before the duality was proposed. Matrix Theory proposal [29] incorporated this viewpoint to the description of gravity. In AdS/CFT, the location of probe D-brane is described in this way [32].
b In the past, people tried to apply Witten's interpretation to a generic coordinate eigenstate in a wave packet, and encountered pathological issues. One has to consider a wave packet as a whole to extract the geometry.
c More precisely, the transverse modes have mass, while the longitudinal modes corresponding to the motion of D-branes do not. It is straightforward to check it analytically in the weak-coupling region of a toy model Hamiltonian (2).
• If the matrices cannot be simultaneously diagonalized but can be taken to be simultaneously block diagonal, then each block describes an extended bound state of D-branes and open strings. We can construct such states by choosing Y and Q appropriately.

Superselection at large N
In the standard model of particle physics, SU(2) L ×U(1) Y is 'spontaneously broken' to U (1) EM by the vacuum expectation value of the Higgs field. More precisely, the global part of gauge symmetry is broken [33] by a specific choice of boundary condition, although the local part is not. Note that quantum states corresponding to different boundary condition are separated by superselection, because no gauge-invariant local operator can connect them. Formally, we can construct a gauge-invariant vacuum by taking a superposition of all possible boundary conditions. However, due to the superselection, we cannot distinguish such a state from a gauge-symmetry-breaking vacuum with a specific boundary condition. Essentially the same mechanism for superselection applies to many states in large-N gauge theories [27]. The wave packet |Y, Q defined above is not U(N )-invariant (unless Y and Q are proportional to the identity matrix). The gauge-invariant counterpart is obtained by symmetrizing over all possible U(N ) transformations, i.e., where N is a normalization factor. In the large-N limit, the size of the wave function along each direction (I, α) is of order N 0 or less, while typical eigenvalues of Y I and Q I we consider are of order √ N . Even for a small U(N ) transformation (specifically, g = e i α τα , where α ( α ) 2 is small but of order N 0 ), the overlap between |Y, Q and |U −1 Y U, U −1 QU is parametrically suppressed. In the free limit, the wave packets are coherent states, and we can explicitly see that the overlap behaves as e − 1 8 I Tr(Y I −U −1 Y I U ) 2 . In typical situations we are interested, TrY 2 is sufficiently larger than O(N 0 ). Therefore, the overlap is essentially the delta function. As observables, we consider gauge-invariant operators which does not change the energy of the system significantly, e.g., Tr(X I 1X I 2 · · ·X I k ), k = O(N 0 ). Such operators cannot connect |Y, Whether we take the expectation value of such an operator by using |Y, Q or its gauge-invariant counterpart, we obtain the same value. In this sense, different gauge choice corresponds to different superselection sector. |Y, Q = 0 and its gauge-symmetrization give the same expectation values, and hence, we cannot distinguish them.

Generalization to QFT
Generalization to QFT (i.e., theories with nonzero spatial dimension) is straightforward. The partition function is the same as (6), except that we need to consider all possible local gauge transformation instead of G that can be schematically written as as G = ⊗ x G x where G x is the gauge group G acting on a spatial point x. Thermal partition function can be written as As a concrete example, let us consider (p + 1)-dimensional maximal super Yang-Mills theory. If we use a lattice regularization for the spatial dimensions, wave functions are defined on R (9−p)×N 2 ×n site , where 9 − p is the number of scalar fields and n site is the number of spatial lattice sites. (Strictly speaking, we have to take into account the spatial directions of the gauge field, whose discretized version live on the links.) The center of a wave packet is specified by Y I ( x) and Q I ( x).
In the case of p = 3, D-brane geometry on R 6 = R >0 × S 5 can be generated from six scalar fields. R >0 is interpreted as the radial coordinate of AdS 5 . In the ground state, Y I = 0 and Q I = 0, and hence, all D-branes are sitting at the origin of AdS 5 , and no open string is excited. This is exactly the expected property of the BPS limit of black three-brane.

Examples
Let us list a few examples of the encoding of geometry into matrix degrees of freedom. We consider matrix model for simplicity. Generalization to QFT is straightforward.

Probe D-branes in BPS black hole geometry
The ground state is given by Y = 0, Q = 0. d All D-branes are sitting on top of each other at the origin, no string is excited, and the area of the horizon is zero.
We can separate a few D-branes from others and probe the geometry. Let us take the BFSS matrix model as an example, and take the (N, N ) components y ≡ (Y 1,N N , · · · , Y 9,N N ), q ≡ (Q 1,N N , · · · , Q 9,N N ) nonzero. Such wave packet describes the BPS black zero-brane and a probe D-brane. Typical distance scale considered for the duality to type IIA black zerobrane is small but of order √ N (i.e., √ N , where is small but of order N 0 ). When | y| is of this order, (1 + 9)-dimensional geometry is probed. As | y| becomes smaller, string coupling becomes larger, and eventually (1 + 10)-dimensional geometry sets in [34]. If we also take the (N − 1, N − 1) components nonzero, the scattering of two D-branes in the black zero-brane geometry can be described.

Two-black-hole geometry and thermofield double state
Let us consider states corresponding to two stacks of D-branes separated by a distance L. We write that the center of wave packet Y I and Q I to be block diagonal, d In the case of BFSS matrix model [29], YI = cI · 1, where 1 is the N × N unit vector and cI is any real number, is equally fine. If we impose the traceless condition, cI has to be zero.
where Y (1) and Q (1) are N 1 × N 1 , and Y (2) and Q (2) are N 2 × N 2 . Furthermore, we assume where the eigenvalues of Y and Q are negligibly small compared to L. Let us use the following notation: If L is sufficiently large,Ŵ s acquire mass of order L and decouple. More precisely, these off-diagonal blocks, together withR I=2,··· ,9 , can be treated as harmonic oscillators with large frequency of order L, and they get frozen to the ground state.Ẑ does not acquire mass, but it enters the Hamiltonian only via the interaction withŴ I=2,··· ,9 , and hence, it also decouples. Then we have copies of two matrix models with SU(N 1 ) or SU(N 2 ) gauge group. (Note that, when off-diagonal blocks decouple, L is so large that dual gravity description in terms of black zero brane geometry is not valid.) This setup is analogous to the ones studied in Refs. [35,36]. If diagonal blocks have sufficiently large energy, each of them can be regarded as a black hole.

Thermofield double
If the two blocks are infinitely separated, off-diagonal entries completely decouple from the dynamics. Then we have two decoupled Hilbert spaces. The Hamiltonian turns to the sum of the SU(N 1 ) Hamiltonian and SU(N 2 ) Hamiltonian. When N 1 = N 2 , this is exactly the setup of two copies of same Hilbert spaces with the same Hamiltonian, considered in Ref. [37]. By taking certain linear combinations of such two-stack wave-packet states, the thermofield double states can be constructed.

Small black hole (partially-deconfined phase)
A small black hole can be described by using partially-deconfined states [24][25][26][27][28], in which a subgroup of U(N ) is deconfined. Roughly speaking, the SU(N 1 )-deconfined states can be constructed as linear combinations of wave packets of the form (23), with Y (2) = Q (2) = 0. Namely, we consider We assume Y , while other blocks are set to zero. Both N 1 D-branes and open strings between them are excited and form an extended object. Other N 2 = N − N 1 D-branes are sitting at the origin and no open strings are excited between them, or between those N 2 D-branes and the other N 1 D-branes. Note that all D-branes contribute to the exterior geometry. Such states can be interpreted as a small black hole, with no radiation emitted yet.
To describe the Hawking radiation, we can take Y (2) and Q (2) to be nonzero, and add small blocks whose sizes are of order one, and interpret them as gravitons. The idea behind it is that the large and small blocks are long and short strings, which are black hole and graviton, respectively. To understand it in terms of operators, let us considerÔ IJ ≡ Tr(X IXJ ) as an example. By definition,Ô IJ = i,jX ij IX ji J . Terms with i = j is negligible at large N , and hence,Ô IJ =X 12 IX 21 J + (permutations). Therefore, with an appropriate gauge fixing, such an operator excites a 2 × 2 block which can be seen as a closed string made of two open strings. In the same manner, from a trace of product of L operators, we can get a closed string made of L open strings. In this way, the block-diagonal structure naturally appears after gauge fixing. Note that N 1 decreases as the black hole shrinks, and the number of small blocks (number of gravitons emitted from black hole) grows.

Fuzzy sphere
Another important class of bound states of D-branes is fuzzy sphere. Fuzzy sphere appears through the Myers effect [38] when D-branes are coupled to background flux. As a concrete setup, here we consider fuzzy sphere in the BMN matrix model [16]. In this case, fuzzy spheres can describe D2-branes, M2-branes, NS5-branes or M5-branes, depending on the way the large-N limit is taken [39][40][41]. Intuitively, fuzzy spheres are analogous to a lattice regularization. Spatial geometry of such higher-dimensional objects are encoded into matrix degrees of freedom, and hence, the matrix entanglement in the BMN matrix model can be reinterpreted as the spatial entanglement in the world-volume theory on the higher-dimensional objects; see Sec. 4.
The bosonic part of the BMN matrix model is written aŝ where I, J = 1, 2, · · · , 9, i, j, k = 1, 2, 3, a = 4, 5, · · · , 9. ijk is the structure constant of SU(2) normalized as 123 = 1. The fuzzy-sphere background is a wave packet localized about the SU(2) algebra, where There are various representations of SU (2), and hence there are many fuzzy sphere backgrounds that are BPS ground states. Let us take J i to be N 2 copies of N 5 -dimensional irreducible representations, where N = N 2 N 5 . When N 2 is fixed and N 5 is sent to infinity, the system of N 2 -coincident D2-or M2-branes are described. (Coupling constant g and flux parameter µ should be scaled appropriately with N .) On the other hand, when N 5 is fixed and N 2 is sent to infinity, the system of N 5 -coincident NS5-or M5-branes are described.

Entanglement between matrix degrees of freedom (Matrix Entanglement)
In this section, we define the matrix entanglement. We utilize the geometric interpretation reviewed in Sec. 2. Although we consider only 'matrix' entanglement here, the generalization to theories with generic matter content, e.g., theories with a field in the fundamental representation such as QCD, is straightforward.

Matrix model
Let us start with a rather formal argument. Suppose N 2 generators of U(N ) are split into two groups A = {α} andĀ = {α }. The extended Hilbert space H ext factorizes as where Therefore, by tracing out H A or HĀ, we can define the entanglement entropy. A nontrivial issue is if there is a natural separation of colors into A andĀ. The answer is positive, as we have seen in Sec. 2. Namely, there are states that admit geometric picture relating the matrix eigenvalues and bulk coordinate, based on the wave-packet states. One might be worried about an apparent need for the gauge fixing. However, such a 'gauge fixing' is automatic at large N thanks to the 'spontaneous breaking' of U(N ), or more precisely the emergence of superselection sectors.
When the superselection takes place, the situation is analogous to the system of identical particles which are well separated. Let us consider the system of two indistinguishable bosons with wave function where the one-particle wave functions φ and χ have vanishingly small overlap due to the separation in space or energy barrier. In such a case, symmetrization of the wave function does not have physical effects, and apparent 'entanglement' coming from the symmetrization of a product state should not have physical effect either. (See e.g., [42] regarding this point.) By replacing the permutation group with U(N ) and repeating the same argument, we should not allow apparent 'entanglement' due to the symmetrization over U(N ). A natural procedure is to take a superselection sector and define the entanglement in that sector.
As an example, let us consider two stacks of D-branes described by a wave packet in which Y I 's and Q I 's are block-diagonal as in (23). Off-diagonal blocks mediate interaction between diagonal blocks. Whether the off-diagonal blocks are treated separately or not is a nontrivial issue. One option is to combine them with one of the diagonal blocks. Then the factorization of the Hilbert space is schematically expressed as For this, we can define the entanglement entropy as whereρ is the density matrix for the entire system. Equivalently, this is the von Neumann entropy of region A. If the density matrix describing the entire system is pure, then S A = SĀ.
Another option is to treat the off-diagonal blocks separately, Then, we can first trace out H C to obtain the reduced density matrix for diagonal blocks aŝ The reduced density matrixρ A∪B is not necessarily a pure state, i.e., can be nonzero. A natural measure for the entanglement between region A and region B would be the mutual information When the theory under consideration has weakly-coupled gravity dual and the partitioning of matrix degrees of freedom corresponds to geometric partitioning on the gravity side, it would be natural to expect that the matrix entanglement is calculated by using holographic fine-grained entropy formulas [4][5][6][7][8]. Depending on the geometry encoded into matrices, there can be multiple ways that matrix degrees of freedom are split into pieces. Then, different partitioning of matrix degrees of freedom should correspond to different extremal surfaces on the gravity side. In Sec. 5, we will consider black hole and Hawking radiations as such an example.
Note that we do not find (or at least have not found yet) a natural way of splitting colors into two sectors when a given state is gauge-invariant even in the extended Hilbert space, such as the confining vacuum that corresponds to Y = Q = 0. This is very different from the target-space entanglement whose main motivation includes the application to the ground state [19][20][21]; see Sec. 3.3.

Generalization to QFT
In the past, the definition of the spatial entanglement in lattice gauge theory has been considered. Refs. [43][44][45] introduced the extended Hilbert space H ext because otherwise the factorization of the Hilbert space does not hold, that is, while where A is a spatial region andĀ is its complement. We use the same extended Hilbert space H ext to define the matrix entanglement. The use of H ext is crucial not just for the factorization of the Hilbert space, but also for the geometric interpretation, as explained in Sec. 2.
With a physically-well-motivated gauge fixing, we can define the entanglement entropy by splitting the matrix degrees of freedom (the internal space) into two pieces. In general, different partitionings of matrix entries are allowed at different spatial regions. The factorization (41) holds for such generic partitions as well. In this paper, we consider only the partitioning of color space into two pieces.

Target-space entanglement vs matrix entanglement
To understand the difference between matrix entanglement and target-space entanglement, and also some similarities in special situations, let us briefly review target-space entanglement.
Suppose there are N particles in R 3 . We fix N and work in the first-quantized setup. There are 3N coordinate operatorsx i ,ŷ i andẑ i (i = 1, 2, · · · , N ), where a set of three operators (x i ,ŷ i ,ẑ i ) defines the location of i-th particle in 'target space' R 3 . When N particles are indistinguishable, the states transforming to each other by S N permutation are identified. Namely, the S N -permutation is a gauge symmetry. Up to the gauge redundancy, the N -particle Hilbert space can be written as We can split the target space R 3 into two pieces A andĀ, and consider a subspace H where and Note that the gauge symmetry is fixed down to S m ×S N −m in each H m . The full Hilbert space H (N ) can be written as where H A ≡ ⊕ ∞ m=0 H A,m and HĀ ≡ ⊕ ∞ m=0 HĀ ,m . Then we can obtainρ A by tracing out HĀ. The von Neumann entropy defined by using thisρ A is the 'target-space entanglement entropy' [18]. It is possible to define such quantity in a gauge-invariant manner.
Such a definition is motivated by the expectation that partitioning of target space captures important physics. This is often true, and hence the target space entanglement can provide us with useful insight. However, the partitioning of target space is not always the right way to go. For example, superfluid and normal fluid are not spatially separated. Still, if there is a natural way to separate the system into two pieces such as some sort of superselection, then entanglement between these two sectors can be defined. In the 'particle entanglement' approach [46], superselection associated with particle number is used.
Next, let us consider U(N )-gauged matrix models. For one-matrix model [47][48][49], a natural approach is to fix U(N ) by diagonalizing the matrix. Then the system is interpreted as a system of N particles with the residual S N symmetry, and hence, we can define targetspace entanglement just as above. Things get complicated when there are two or more matrices, say X 1 , · · · , X d where d ≥ 2. In the past, it has been proposed that one should diagonalize one of the matrices, say X 1 , and use its eigenvalues to specify two regions A and A on R [19][20][21]. Although such a definition is well-defined, it is not clear if the entanglement entropy defined in this way correctly captures emergent geometry in the context of holography. In fact, at least for low-energy states including the ground state, the approach based on wave packets captures the geometry of string theory motivated by the probe-D-brane picture more precisely, as we have seen in Sec. 2.

Application 1: Noncommutative geometry, M2-branes and M5-branes
As the first example of application of matrix entanglement, we consider fuzzy sphere in the BMN matrix mnodel, which was introduced in Sec. 2.5.4. We consider small excitations above the fuzzy-sphere background Then, in the appropriate large-N limit, D2/M2brane theory (3d maximal SYM) or NS5/M5-brane theory (6d N = (2, 0) theory) can be described. The BMN matrix model can be used to regularize these theories, similarly to the lattice gauge theory. As we will see shortly, geometric information is encoded into matrix degrees of freedom. Therefore, by using matrix entanglement, we can define the regularized version of spatial entanglement.
Firstly, let us consider the large-N limit that describes N 2 -coincident D2-or M2-branes, i.e., fixed N 2 and N 5 → ∞. In this limit, 3d maximal SYM on two-sphere is obtained. More precisely, the radius of sphere is 3 µ , the coupling constant is g 2 3d = g 2 1d µ 2 N 5 , and spatial coordinate has noncommutativity characterized by the noncommutativity parameter θ = 1 µ 2 N 5 . If the noncommutativity parameter is sent to zero, the ordinary 3d SYM on commutative sphere is obtained. (The commutative limit is not singular because of maximal supersymmetry [50,51].) Geometric information on two-sphere is encoded into matrix degrees of freedom via the fuzzy spherical harmonics; see e.g., Ref. [52]. Therefore, geometric partitioning on two-sphere can be mapped into a partitioning of matrix degrees of freedom. Refs. [53][54][55] applied this idea to scalar theories on fuzzy sphere. They considered one-matrix model whose kinetic term is proportional to Tr(Φ 2 + µ 2 Note that J i 's are introduced by hand. Our matrix-entanglement approach generalize it to gauge theory, in which J i emerges dynamically. Specifically, let us fix the gauge such that J 3 is diagonalized in such a way that J (N 5 ) 3 = diag(s, s − 1, · · · , −s + 1, −s) and where 2s + 1 = N 5 . Then, if we split fuzzy sphere into two spatial regions as shown in the left panel in Fig. 3, the corresponding partitioning of matrix degrees of freedom becomes like the right panel of Fig. 3 [53][54][55]. Note that there is no ambiguity in the partitioning of the matrix side, including the off-diagonal entries. Next, let us consider another large-N limit that describes N 5 -coincident NS5-or M5branes, i.e., fixed N 5 and N 2 → ∞. In this case, it is expected that the eigenvalue distribution of Y 4 , · · · , Y 9 becomes five-sphere that corresponds to spherical NS5-or M5-branes [40,41], modulo the assumption that Y 4 , · · · , Y 9 are simultaneously diagonalizable. If this expectation is true, then the geometric partitioning on five-branes can be related to the partitioning of matrix degrees of freedom. Unlike the partitioning of fuzzy sphere (Fig. 3), in this case it is not immediately clear how to treat the off-diagonal entries. Both (35) and (37) seem to be natural.

Application 2: Small black hole and Hawking radiations
In this section, we consider small black hole and Hawking radiations as an application of the matrix entanglement.
As shown by Horowitz [17], a small black hole in AdS 5 ×S 5 described by dual 4d SYM on S 3 can evaporate into gas of gravitons if the coupling constant on the QFT side is sufficiently large (λ tHooft N 8/17 ) e,f and energy is sufficiently small (E N 20/17 ). The same might be true for 11d black hole described by Matrix Theory [12][13][14][15]. Note that there is no need for external heat bath, and hence such systems provide us with a theoretically clean setup to study the black hole information problem.
As discussed in Sec. 2.5.3, small black hole is described by partially-deconfined states. A simple mechanism of the emergence of Page curve naturally follows from partial deconfinement. (See Ref. [56] for an early suggestion that partial deconfinement should be important to resolve the information puzzle.) Below, we first consider a small but not too small black hole that does not evaporate (Sec. 5.1), and then move on to the case of sufficiently small black hole that evaporates completely (Sec. 5.2).

Bekenstein-Hawking entropy from entanglement
Let us consider a small black hole without radiations g that is represented by a superposition of wave packets of the form (27). We can define the reduced density matrix and von Neumann entropy S A by keeping only the deconfined sector (region A), as in (36). We would like to know whether S A obtained in this way can be the Bekenstein-Hawking entropy.
Let N 1 be the size of the deconfined sector, N 2 = N − N 1 N 1 , and the coarse-grained entropy be S coarse grained = log K.
Then, there are K wave packets |Y [a] , Q [a] (a = 1, 2, · · · , K) that give the same values for the macroscopic observables (energy, coarse-grained eigenvalue distribution, etc) up to 1/Ncorrections and the coarse-grained reduced density matrix is written aŝ e For the stringy corrections to type IIB supergravity to be small, g 2 YM ∼ gs 1 and α R 2 AdS ∼ (g 2 YM N ) −1/2 1 are required. In the 't Hooft large-N limit ('t Hooft coupling λ = g 2 YM N fixed), the former condition is automatically satisfied, and the latter can be satisfied as well if λ is taken to be large. To realize λ tHooft N 8/17 without having large stringy corrections, we should take g 2 YM ∼ N −α , 0 ≤ α ≤ 9 17 . f If the energy of black hole is not sufficiently small, graviton gas becomes dense before black hole completely evaporates and the equilibrium state consisting of black hole and graviton gas is realized. If λ tHooft N 8/17 , the stringy effect is not negligible for such a small black hole. This is the reason that λ tHooft N 8/17 is required.
g More precisely, we consider a situation that black hole does not evaporate completely and the density of radiations is negligibly small.
Only the upper-left block (region A) of (Y [a] , Q [a] ), denoted by (Y (1) [a] , Q [a] ), is nonzero, as in (27). At weak coupling, explicit computations show K ∼ e N 2 1 [27]. On the gravity side, this coarse-grained entropy should be the Bekenstein-Hawking entropy. Therefore, our task is to see whether this coarse-grained entropy coincides with the von Neumann entropy S A .
We consider a superposition of K wave packets as a typical state corresponding to the small black hole: We would like to claim that, by tracing out HĀ, we obtain [a] Y [a] , Q [a] |.
Then the von Neumann entropy is also S A ∼ log K. Intuitively, we would like to claim that regionĀ behaves like a heat bath. For (54) to hold, has to be negligible. One might think this is trivial because the overlap between low-energy wave packets are exponentially suppressed at large N . However we have to be a little bit careful here: we have to confirm that the overlap is small in HĀ, where Y and Q are zero. We argue that this is indeed the case, and [a] , Q [a] Y [b] , Q with c e −N 1 N 2 .
To motivate this scaling, let us consider a simplified situation where Y [a] and Y [b] are diagonal, Q [a] = Q (1) [b] = 0, and eigenvalues are sufficiently far separated such that off-diagonal entries can be treated as harmonic oscillators with frequency ω proportional to the distance between eigenvalues. Let y i [a] and y i let's take this to be large but of order N 0 . For i = 1, · · · , N 1 and j = N 1 + 1, · · · , N , the wave functions are the normal distribution with the width (g · | y i [a] |) −1/2 and (g · | y i [b] |) −1/2 , respectively. Then, unless | (which is not the case in general), the wave functions do not overlap perfectly, and the inner product of these components in two wave packets are suppressed. There are N 1 N 2 matrix degrees of freedom contributing to such a supression, and hence the exponentially small suppression factor e −N 1 N 2 is expected. The lower-right block (i, j = N 1 + 1, · · · , N ) may give further suppressions when they are traced out. We expect similar suppression also when Y [a] and Y [b] are not simultaneously diagonalizable, because there is no reason to expect perfect overlap of wave function along any of the dimensions. Therefore, (54) should hold, and S A = log K.
So far we called region A in this picture as 'black hole', but in fact, D-branes described by regionĀ are also sitting inside the 'black hole'; they are sitting at the origin of the bulk geometry. (See Sec. 2.5.3.) It seems to be natural to interpret region A as the degrees of freedom near the black hole horizon as depicted in Fig. 4, because the deconfined sector describes an extended bound state of D-branes and strings, and hence, D-branes in region A cannot be sitting at the origin. A A Figure 4: Geometry of a small black hole encoded into a partially-deconfined state, in terms of D-brane configuration. Region A describes the degrees of freedom near the black hole horizon. D-branes described by regionĀ is sitting at the origin of R 9−p . Only region A contributes to the coarse-grained entropy.
In terms of matrix degrees of freedom, the formation process of black hole would look like Fig. 5, from left to right. In the figure, first there were multiple small objects (∼ multiple small deconfined sectors). Eventually, they merge and form a single deconfined sector. When a black hole is newly formed, we do not expect that regions A andĀ are fully entangled. As an extreme example, we can imagine a pure state consisting of just one specific wave packet |Y, Q . After some time (presumably the scrambling time), the system evolves to a generic superposition of many wave packets and then A andĀ are strongly entangled. The argument above assumes such a situation.

Two stacks of D-branes
As an analogous but slightly different situation, let us consider two stacks of D-branes considered in Sec. 2.5.2. We assume that Y I 's and Q I 's are block-diagonal as in (23). To be concrete, we assume each block is N 1 = N 2 = N 2 . Furthermore, we assume that Y 1 is close to diag L 2 , · · · , L 2 , − L 2 , · · · , − L 2 , where + L 2 and − L 2 appear N 2 times for each, and other Y I 's and Q I 's are close to zero. Still, we can consider a linear combination of wave packets of this form: In the limit of L → ∞, off-diagonal blocks decouple. Thermofield double states are special case of such two-stack states. It is convenient to trace out such decoupled degrees of freedom (region C in (37)) and keep diagonal blocks (region A and region B in (37)). In this case, matrix entries in the off-diagonal blocks are harmonic oscillators with frequency ω = gL up to negligible corrections, and hence they are indistinguishable regardless of the details of the diagonal block, unlike the case of the small black hole considered above. (Strictly speaking, we are considering L → ∞ for large but fixed N .) We obtain [a] |Y (2) [a] , Q [a] Y (1) [a] |Y (2) [a] , Q The decoupled degrees of freedom do not play nontrivial role in entanglement.

Evaporation and Page curve
A cartoon picture of the evaporation of small black hole in terms of matrix degrees of freedom is shown in Fig. 5, from right to left. Different from the black hole formation (from left to right), small blocks should correspond to gravitons; see Sec. 2.5.3.
Suppose that the small BH corresponds to M × M deconfined block (region A) and radiations are described by M × M block (region B). In addition, there is a confined sector (region C); see Fig. 6. We assume that M, M N , because otherwise the black hole cannot evaporate completely. We trace out the Hilbert space corresponding to the confined sector C, and define the reduced density matrixρ A∪B : Reduced density matrixρ A∪B obtained this way is a mixed state, even if the original density matrixρ is pure. Hence the von Neumann entropy S A∪B (entropy of the system of black hole and radiations) is nonzero. We can also define the von Neumann entropy S A (entropy of black hole) and S B (entropy of radiations) as Based on the argument along the lines of Sec. 5.1, these von Neumann entropies should be the same as the coarse-grained entropy. As the black hole evaporates, S A decreases, while S B increases. The crucial point is that M and M change with time; see Fig. 7. Note that S A∪B = S A at early time and S A∪B = S B at late time. There is no reason to expect that region A and region B are directly entangled, and hence we expect S A∪B = S A + S B . Apparently, the time dependence of S A∪B does not exhibit the Page curve. As we will see, however, there is no contradiction. One has to take into account the degrees of freedom in the confined sector to obtain the Page curve. This is different from the assumption in Page's original argument [9] that black hole and radiations form a pure state.

Page curve
Let us introduce regions A, B and C as shown in Fig. 6. The region C corresponds to Y = Q = 0 and hence it does not contribute to the coarse-grained entropy. If we trace out region C, we obtain the reduced density matrix discussed above. Note that C contains most of the matrix degrees of freedom; otherwise the small black hole cannot evaporate completely. In Fig. 6, A and B are drawn much bigger than their actual sizes. Because the entire system A ∪ B ∪ C is in a pure state, we immediately obtain and Therefore, if we interpret that min(S B , S B∪C ) is the entropy of radiations, we obtain the Page curve; see Fig. 8. Clearly, there is a close similarity to the island formula [10,11]: it appears that two ways of partitioning the matrix degrees of freedom (A and B ∪ C, or A ∪ C and B) correspond to two ways of choosing the extremal surfaces in the island formula. In Sec. 5.2.4, we will discuss whether this mechanism is related to the island proposal or not. At this stage, we will study the matrix entanglement as it is, without referring to the island proposal.

Meaning of the confined sector on the gravity side
Geometry obtained by identifying matrix eigenvalues and D-brane locations [23,30,31] is shown in Fig. 9. (Here we assumed that the radial coordinates determined from small blocks are large such that they can be regarded as gravitons outside the black hole. Explicit demonstration of this picture is desirable to strengthen our claim below.) Region A describes the degrees of freedom near the black hole horizon. D-branes in region C are sitting behind the horizon. Region C is strongly entangled with radiations (region B) at late time, and when this region and radiations are combined, the von-Neumann entropy becomes much smaller, as shown in Fig. 8.
In Page's argument, only radiations and black hole are concerned. Therefore, if we want to put our finding to his framework, we have to include the region C as a part of radiations or black hole. Time dependence of von Neumann entropy shown in Fig. 8 suggests that, before and after the Page time (the time when S A = S B holds), B and B ∪ C should be regarded as "radiations", respectively. Can we justify this interpretation?
It is easy to see that the alternative choice leads to a pathology. We assume that S A and S B agree with the coarse-grained entropy of black hole and radiations, respectively. (This is the reason that we believe that black hole can completely evaporate, based on the calculation on the gravity side [17].) Intuitively, and somewhat naively, this is the number of independent wave packets, i.e., possible choice of Y and Q, in each sector. Seen this way, region C does not contribute to the coarse-grained entropy, because Y and Q vanish there. However if we use precise definition of coarse-grained entropy (i.e., consider all possible density matrices that give the same expectation values for macroscopic observables under consideration, calculate their von Neumann entropies, and then choose the largest value as 'coarse-grained entropy'; see e.g., Ref. [1]), region C does have a visible contribution. Suppose we interpreted that B∪C describes radiations before Page time. Then the von Neumann entropy S B∪C = S A is larger than S B . Then the coarse-grained entropy cannot be S B , because the von Neumann entropy gives the lower bound for the coarse-grained entropy. We would encounter the same pathology if we interpreted that A ∪ C describes black hole (equivalently, B describes radiations) after the Page time, because the von Neumann entropy S A∪C = S B is larger than S A which we want to have as the coarse-grained entropy. Therefore, as long as we split the system to two pieces ("black hole" and "radiations") in a way consistent with semiclassical picture on the gravity side, we are forced to interpret that "Hilbert space of radiations" is H B before the Page time and H B∪C after the Page time. It would mean that, although both orange and red lines in Fig. 8 are well-defined on the gauge theory side, the dual gravity description in semiclassical language can be valid only for the gray line (Page curve).
To have a better argument, we have to specify the rule of the game precisely: which part of the Hilbert space can Bob, the exterior observer, manipulate? Clearly, we should assume that Bob can access to the region B, but not to region A. How about the confined sector (region C)? Although D-branes in the confined sector is sitting at the origin and hence they are behind the horizon of the small black hole at the center of the bulk, they interact with region B via off-diagonal entries (open strings), and due to such interactions, radiations in region B know they are in asymptotically AdS 5 ×S 5 spacetime. Furthermore, when the AdS radius is very large and black hole is very small, black hole can sit away from the origin at least for some time, if we fine-tune the initial condition. h It is as if region C prepares the background spacetime which black hole and radiations live in. See Sec. 6 for further discussion. Therefore, we find it more natural to assume that Bob has access to region C as well. In the following, we assume that Bob can access B ∪ C. Specifically, we assume that Bob knows the reduced density matrixρ B∪C .
As usual, Bob's task is to read Alice's diary which was thrown into the black hole. For that purpose, he wants to know as much information on radiations as possible. He has three options: (1) keepρ B∪C , (2) discard C and getρ B , and (3) discard B and getρ C . He wants to specify the state as precisely as possible, and hence, he wants to minimize the entropy. The von Neumann entropies associated with (1), (2) and (3) are S B∪C = S A , S B and S C = S A + S B , respectively. Therefore, before and after the Page time, he should choose option (2) and option (1), respectively. The Page curve represents the best of Bob's knowledge.

Time scales associated with the evaporation
Before the Page time, the smallest von Neumann entropy Bob can achieve is S B , that is the coarse-grained entropy of radiations. He does not know the information on the confined sector. i,j On the other hand, after the Page time he can make the entropy as small as h There are two effects which pull black hole toward the center. Firstly, there is a conformal mass term for scalar fields. Secondly, D-branes attract with each other when supersymmetry is broken by thermal excitations. The latter works in the BFSS matrix model as well.
i He can decrease the entropy further by throwing away a part of region B, but that does not add any information on C.
j At weak coupling, entanglement between region A and region C (resp., B and region C) is encoded into the off-diagonal blocks corresponding to open strings between D-branes in region A and region C (resp., B and region C). By specifically looking at those off-diagonal blocks, we can extract the information regarding black hole or radiations. However, at strong coupling, it is highly unlikely that such a simple partitioning S B∪C = S A < S B . The confined sector gave him some information. In this sense, Bob can "see" the region C only after the Page time.
It is important to note that the estimate of von Neumann entropy above is valid only after certain time needed for regions A, B and C to fully entangle with each other. When a black hole is newly formed, we do not expect that region C is fully entangled with region A. Region B does not exist or is very small at this point. As an extreme example, we can imagine a pure state consisting of just one specific wave packet |Y, Q . After some time, the system evolves into a generic superposition and then regions A, B and region C are strongly entangled. Due to such entanglement, S A and S B should coincide with the coarse-grained entropies of black hole and Hawking radiation, respectively. The time dependence of the entropies shown in Fig. 8 applies after this time scale. Such a time dependence closely resemble the time scale for the formation of the island on the gravity side.
To see when Bob can steal Alice's secret, let us slightly modify the Hayden-Preskill information retrieval protocol [57]. The original Hayden-Preskill protocol assumes maximal entanglement between black hole and a subsystem of radiations after the Page time. This assumption is satisfied if H B∪C is regarded as "Hilbert space of radiations" after the Page time. Suppose Alice throw her diary (a diagonal block in region B) into black hole (region A) after the Page time. Then region A becomes slightly bigger, but region C is not affected instantaneously. A natural time scale when the information thrown into region A is transferred to region C is the scrambling time. Although this information transfer can take place even before the Page time, Bob can have the information on the region C only after the Page time. Therefore, he can read Alice's diary only after the Page time. Figure 8: Regions A and B are interpreted as black hole and radiations, respectively. Region C describes D-branes sitting at the origin, which is inside the horizon. After the Page time, it is natural to interpret that C is a part of radiations.

Entanglement Island?
In the discussions regarding the Page curve and confined sector, there were several features of the confined sector (region C) that closely resemble the entanglement island: • D-branes in the confined sector are sitting behind the horizon (Sec. 5.2.2 and Fig. 9).
exists. Note that the gauge invariance of the confined sector restricts the possible options for partitioning rather strictly.  • Time scale for the information retrieval appears to be the same (Sec. 5.2.3).
Therefore, if the gravity side admits the entanglement island, it would be natural to expect that the confined sector corresponds to the island. Still, there are several points that appear to be different from the entanglement island proposal. Firstly, it may not be appropriate to interpret that the confined sector describes the interior of black hole. Rather, it may describe the AdS 5 ×S 5 spacetime that the black hole and radiations live in; see a discussion in Sec. 6. Secondly, the notion of the entanglement wedge, which is important in the island proposal, is not clear in terms of matrix degrees of freedom. Note also that there is a debate concerning the island proposal when the non-gravitational heat bath is not introduced. Because our setup does not involve non-gravitational heat bath, we need to introduce the interface separating black hole and radiations at such a location where gravity is dynamical. There are arguments against [2,58,59] and supporting [60,61] the use of the holographic entanglement formula and the existence of the entanglement island in such a case. Although our argument on the gauge theory side is not affected by this debate, the dual gravitational interpretation would be different depending on the conclusion of the debate. There is also an argument against existence of islands in a theory of massless gravity [62]. Perhaps matrix entanglement can provide us with a concrete setup to better understand this issue. k It would be fair to say that we do not know the precise relationship between the confined sector in gauge theory and the entanglement island in gravity. Whether they are the same or not, we expect that matrix entanglement is one of the keys to the understanding of the black hole information problem in gauge/gravity duality.

5.3
Type IIA region of D0-brane matrix model (Toy model for black hole evaporation) The small black hole discussed above was conceptually very clean on the gravity side. However, analytic calculations on the QFT side was difficult. Below, let us show a toy model [63,64] that has qualitatively similar behavior from QFT point of view, although the weaklycoupled dual gravity description is lacking. Let us consider the parameter region of D0-brane matrix model that is dual to type IIA black zero-brane. The black zero-brane is a bound state of many D0-branes and open strings. As the initial condition, let us assume that all N D0-branes described by the U(N ) theory are in the black zero-brane. At large but finite N , the black zero-brane can decay by emitting D0-branes toward flat direction at a very long time scale. (The time scale for the emission of first few D0-branes is of order e N .) A cartoon picture of such 'evaporation' is shown in Fig. 10. Emitted D0-branes will run away to infinity. They are analogous to the Hawking radiations as far as the eigenvalues are concerned, although they are not massless. The black zero-brane is described by the M × M -block, where M decreases from N to 0. As D0-branes (diagonal entries) escape to infinity, open strings (off-diagonal entries) become infinitely heavy and decouple from the dynamics. In terms of wave packet, off-diagonal entries of Y becomes zero. We can regard this process as a gradual Higgsing, When D0-branes emitted from black hole go very far, weakly-coupled gravity is not a good description. Still, some important features of black hole evaporation are captured, as we will see below.
Let regions A, B and C be black zero-brane, emitted D0-branes and off-diagonal entries connecting them (Fig. 10, top row). Off-diagonal entries can be treated perturbatively, as we did in a sample calculation in Sec. 5.1. Qualitatively, the von Neumann entropy behaves as shown in Fig. 10.
Off-diagonal elements connecting the black zero-brane and emitted D0-branes (region C) are completely frozen to the ground state due to Higgsing. In the case of small black hole described by 4d SYM, a large number of degrees of freedom are frozen due to confinement. Either way, those 'frozen' degrees of freedom (region C) are strongly entangled with 'radiations' (region B) at late time.
In this example, S A = S B∪C starts from zero, increases at first and then decreases to zero, resembling the Page curve. On the other hand, S B starts from zero, monotonically increases and saturates at finite value when the 'evaporation' ends.
Note also that the decoupling of off-diagonal degrees of freedom leads to a negative specific heat [63,64]. Because the energy per dynamical degrees of freedom increases, temperature goes up.

Conclusion, discussion and outlook
In this paper, we introduced the notion of matrix entanglement for large-N gauge theories, including matrix model and QFT. Based on previous references, we argued that matrix degrees of freedom naturally split into multiple sectors that have natural geometric interpretations in terms of string theory. When applied to fuzzy-sphere states in the BMN matrix model, matrix entanglement can be used to define the spatial entanglement in emergent QFTs nonperturbatively. We also considered the evaporation of small black hole and pointed out that the confined degrees of freedom play important role in the derivation of the Page curve. The confined sector behaves similarly to the entanglement island in several ways. Currently, we do not know whether the confined sector is the gauge-theory counterpart of the entanglement island, or they are actually different and the confined sector gives the Page curve without involving the island.
To conclude the paper, let us discuss several other open problems.
Interior or exterior?
In this paper, we discussed a possibility that the confined degrees of freedom in the partiallydeconfined states correspond to the entanglement island sitting behind the horizon. On the other hand, in the past, three of the authors (MH, AJ, CP) and collaborators speculated that the confined sector describes the exterior; see e.g., Fig. 11, which is taken from Ref. [27]. What would be the precise interpretation then?
The speculation in Ref. [27] was partly based on the assumption that the 'eigenvalues' of matrices are spread out even at low energy [65,66] and fill the significant portion of bulk geometry (for the case of AdS 5 /CFT 4 , the distance of order of AdS radius from the center of the bulk). This assumption is not necessarily valid for low-energy states, because 'eigenvalues' used there describe the location of D-branes only near the boundary of the bulk [23,30]. The eigenvalues of Y I (center of wave packet) can provide us with a better description which is valid even near the center of the bulk, and by using Y I we can show that most D-branes are well localized in the low-energy states. In particular, D-branes in the confined sector sit at the origin of the bulk. When a small black hole (deconfined sector) is formed at the center of the bulk, those D-branes in the confined sector sit behind the horizon. In this sense the confined sector describes the interior. Still, the confined sector plays a crucial role for the exterior: the bulk geometry asymptotically coincides with black p-brane (AdS 5 ×S 5 in the case of p = 3) because of the existence of the confined sector, and the curvature is determined by the total number of degrees of freedom that includes both confined and deconfined sectors. A larger fraction of the bulk is covered by black hole when the confined sector is smaller. Furthermore, it is possible to excite some degrees of freedom from the confined sector and probe exterior geometry. Therefore, we are tempted to think that the confined sector describes both interior and exterior. Note also that, when the AdS radius is very large and black hole is very small, black hole can sit away from the origin, at least for some time, then the D-branes in the confined sector are not necessarily behind the horizon of this black hole. (They can be still regarded as sitting behind the horizon of supersymmetric black hole whose horizon area is zero.) It appears that the entire black hole geometry -both interior and exterior -is encoded in the way matrix degrees of freedom interact with each other. Figure 11: Fig. 4 in Ref. [27]. Is it consistent with the interpretation provided in this paper?

More on target-space entanglement
In the target-space entanglement approach proposed in the past [19][20][21], X 1 was diagonalized. Instead we can diagonalize Y 1 and define the separation of matrix degrees of freedom. If the eigenvalues of Y 1 are sufficiently separated, these two approaches give the same result. However, in the strongly coupled region where the separation is not large, matrix entanglement captures emergent geometry more precisely. More generally, it would be better to think that Y I 's describe matrix-valued target space, because the block-diagonal structure is of crucial importance.
In Ref. [67], extremal surfaces in the black zero-brane geometry dual to the D0-brane matrix model have been investigated. It would be interesting to study if these extremal surfaces correspond to some partition of matrix degrees of freedom.
What are the geometry made of ?
We conjecture that the area of various extremal surfaces correspond to entanglement entropy associated with physically meaningful partition of matrix degrees of freedom. If this is correct, then the information regarding the metric in the gravitational geometry can be obtained from quantum entanglement, by using the duality to define quantum gravity. This approach fits into the direction of study called 'It from qubit'.
Another natural way to determine the gravitational geometry is to see the interactions. For example, if one of the matrix degrees of freedom is used as a probe D-brane, its motion should be described by the Dirac-Born-Infeld action on a geometry sourced by other matrix degrees of freedom [32]. l This approach seems to be practically useful at least outside the horizon.
Although it is natural to expect that these two approaches define the same geometry, such an agreement is far from trivial. By confirming the consistency between these two approaches, we might be able to provide strong tests of the validity of the emergent geometry in holography. For explicit confirmation, machine-learning-motivated variational Monte Carlo method [68] and quantum simulation (including quantum-classical hybrid algorithm which can run on NISQ devices) [69] may be useful.

Numerical simulations
To study small black hole quantitatively from QFT side, numerical methods will be needed. 4d SYM may be difficult at this moment, but the M-theory parameter region of D0-brane matrix model is within reach [15]. A tricky issue is that the small black hole corresponds to the maximum of free energy in the canonical ensemble, which is not efficiently sampled by standard lattice simulation methods. This problem can be circumvented rather easily by constraining the value of the Polyakov loop [70]. A gauge fixing in which confined and deconfined sectors are clearly separated can be achieved by combining the static-diagonal gauge and appropriate permutations of Polyakov line phases [71]. Once the gauge is fixed properly, the Rényi entropy can be calculated by using the Replica trick. It would be nice if we could reproduce black hole entropy quantitatively.

Other theories, other matter contents
In this paper, we considered only the matrix-valued fields. It is straightforward to generalize the construction to other representations, i.e., theories with fundamental fields such as the O(N ) vector model or QCD. It might be interesting to study the entanglement between different component of the vectors in the O(N ) vector model with singlet constraint, because this model has dual higher spin theory on AdS [73], and the mapping between the two sides of the duality is understood [74,75] better than in the duality between super Yang-Mills theory and superstring theory. Note that partial deconfinement takes place in the O(N ) vector model as well [27], which provides us with a natural way of splitting vector components into two pieces. See also [76] for the construction of the thermofield double states.
The SYK model does not have gauge symmetry and hence the definition of entanglement might be conceptually simpler. For the calculation of entanglement associated with the partition of fermions to two groups, see e.g., [77,78]. The coupled-SYK model [79] has a phase structure analogous to 4d maximal super Yang-Mills on three-sphere, in that a phase with negative specific heat resembling a small black hole exists. Perhaps there is a natural splitting of degrees of freedom similar to partial deconfinement: in the 'confined' sector fermions in the left and right copies would be entangled, while in the 'deconfined' sector the entanglement would be lost. Such a splitting, if exist, can explain the dual gravitational interpretation proposed in [79]: the traversable wormhole (confined sector) is gradually destroyed as black hole (deconfined sector) grows [80]. Such a splitting can be confirmed in toy models such as the coupled-matrix model [80].

Definition of matrix entanglement without relying on gauge fixing
In this paper, we defined the matrix entanglement after gauge fixing. Such a definition is natural at large N and for the class of states we considered, because of the emergence of superselection sectors. To consider other kinds of states or finite-N effects, a definition which does not rely on gauge fixing may be needed. See Ref. [20] for the gauge-invariant construction of the target-space entanglement entropy. rather trivial manner, the thermal partition function (6) can be rewritten as