Target space entanglement in Matrix Models

We study target space entanglement in gauged multi-matrix models as models of entanglement between groups of D-branes separated by a planar entangling surface, paying close attention to the implementation of gauge invariance. We open with a review of target space entanglement between identical particles, which shares some important features (specifically a gauged permutation symmetry) with our main problem. For our matrix models, we implement a gauge fixing well-adapted to the entangling surface. In this gauge, we map the matrix model problem to that of entanglement of a $U(1)$ gauge theory on a complete or all-to-all lattice. Matrix elements corresponding to open strings stretching across the entangling surface in the target space lead to interesting contributions to the entanglement entropy.


Introduction
The study of quantum entanglement between spatial regions in strongly coupled, large-N theories, starting with the classic work [1][2][3], has led to crucial insights regarding the emergence of the bulk spacetime in holographic theories. However, we expect that other means of separating degrees of freedom and computing their entanglement are crucial for a variety of reasons: • The existence of entanglement shadows (see for example [4][5][6]), regions of the bulk that are not in the entanglement wedge of any spatial region.
• The fact that the entangling surfaces for spatial entanglement in holographic conformal field theories do not resolve the large "internal" manifolds transverse to the anti-de Sitter factors of the gravitational duals. Possible duals of extremal surfaces in this direction have been discussed in [7,8]. Note that these surfaces are not minimal [6]).
• The existence of large-N gauged matrix quantum mechanics models such as [9,10] with holographic duals, and no spatial regions to entangle. For the one-matrix model a direct calculation with a clear dual interpretation was provided in [11,12], and there has been some discussion of bulk surfaces and possible dual boundary entanglement quantities in the BFSS model [13][14][15], but the question remains open. 1 A natural avenue to investigate is the entanglement between matrix degrees of freedom in large-N gauged matrix quantum theories. It has been argued for many years that these degrees of freedom are crucial in exploring local physics at scales below the bulk radius of curvature. An early hint comes from [16], from which it is clear that in global anti-de Sitter space, a patch of size R AdS has N 2 degrees of freedom, when the dual holographic theory is a large-N gauge theory. Recent work on resolving the (AdS-scale) directions transverse to anti-de Sitter space has indeed focused on entanglement between sectors charged under different unbroken gauge symmetries in the Coulomb branch phase of these theories [7,8].
The first significant challenge is isolating degrees of freedom and defining their entanglement with the rest in a gauge-invariant fashion. For the quantum mechanics of a single N × N Hermitian matrix transforming in the adjoint of the U (N ) gauge symmetry, a physical of entanglement and its dual interpretation in 2d string theory was described elegantly in [11,12,17]. However, in that case, the gauge symmetry leaves only O(N ) gauge-invariant degrees of freedom, namely the eigenvalues of that matrix up to permutations. Here we will consider matrix quantum mechanical examples with two or more matrices transforming in the adjoint representation, for which there are O(N 2 ) physical, gauge-invariant degrees of freedom. Through a straightforward gauge fixing we can readily define a physical notion of "target space entanglement". Crudely speaking the diagonal elements of the matrices correspond to the positions of D-branes, and the off-diagonal modes to open strings stretching between them. For planar entangling surfaces, we can choose the matrix representing the direction transverse to this surface, use the gauge invariance to diagonalize it, and define an appropriate notion of target space entanglement. This leaves two issues: • This procedure leaves a residual U (1) N S N gauge symmetry unfixed. The latter is the Weyl group of U (N ) and permutes the D-branes.
• We must make a physical decision as to the inclusion and treatment of off-diagonal matrix elements corresponding to open strings stretching across the entangling surface.
In this note, we attack the above questions by studying bosonic models with two and three matrices. These can be thought of as toy versions of the BFSS model [9]. The residual permutation symmetry is taken care of by considering "target space entanglement" [17], in which we use spacetime position as a quantum number to distinguish particles, much as one would in an EPR experiment involving correlated photons or spin- 1 2 particles. Once this is done, we can think of matrix excitations as living on a complete (all-to-all) lattice, with each vertex corresponding to the location of the diagonal matrix elements, and each link to the strings stretching between these locations. The U (1) N gauge invariance ensures that the number of strings and anti-strings entering each node have total charge equal to zero. The nodes can be separated according to which side of the entangling surface they reside on. We must then decide how to treat the open strings stretching across the entangling surface. Since we cannot observe one end of these strings, we treat them as unobserved.
The resulting reduced density matrix breaks up into superselection sectors according to the observed U (1) charge carried by these strings. This is closely analogous to the discussion in [18][19][20] of lattice gauge theories in the electric center definition. Furthermore, for lowenergy configurations, the dependence of the state of these strings on the position of the hidden and visible branes induces additional entanglement. This note is organized as follows: In section 2 we review EPR-like entanglement for identical particles. This serves as a primer/review of target space entanglement in the presence of a gauged permutation symmetry, in a simple set of examples. Next, in section 3 we study toy models with two and three matrices at weak coupling, to understand entanglement between matrix degrees of freedom. We draw the connection to entanglement in lattice gauge theories, and point out distinct features involving the off-diagonal degrees of freedom in the two and three-matrix cases which we think captures important qualitative features of entanglement in the full BFSS model. Finally, in section 4 describe some open problems and directions of future work. Appendix A contains details of the construction of the wavefunction used in the 2-matrix case while appendix B describes leading order corrections to this wavefunction in the Born-Oppenheimer approximation.

Relation to recent work
While this paper was in preparation, two excellent papers appeared with substantial overlap [14,15]. As in our work, the authors study the connection between target space entanglement in Dp-brane field theory and bulk entanglement entropy across surfaces separating Dp-branes: first using the same gauge as our work, [14], and then in a gauge-invariant fashion [15]. These works are complementary and we develop different points, both in the overall explication and in the details of our calculations. To begin with, the focus of that work is on states near to origin of the "Coulomb branch" in which the eigenvalues of the matrices (aka D0-brane locations) are close together. In this case, given an entangling surface, the wavefunction breaks up into components with different numbers of D0-branes on each side: the reduced density matrix for a spatial region then breaks up into terms with different numbers of D0-branes in that region (as in [17]). The coefficients of each of these terms contribute a "classical" or "disorder" component to the entanglement entropy. In this work we are interested in well-separated clumps, so we project onto a fixed number.
A second difference is in our treatment of the off-diagonal matrix elements connecting degrees of freedom on each side of the entangling surface. In [14,15] the authors suggest two different prescriptions to compute entanglement: in the subalgebra of observables.
Our prescription for computing entanglement entropy agrees with 2), and we argue that this is a physical choice. We show that these degrees of freedom lead to the reduced density matrix breaking up into superselection sectors given by the gauge charge carried by strings stretching across the entangling surface, providing a "classical" component to the entanglement very much in analogy with the "classical" components that arise in some treatments of lattice gauge theory [18][19][20]. Furthermore, in studying the low-energy states for well-separated branes in the Born-Oppenheimer approximation, we show that the dependence of the state of the off-diagonal modes on the relative separation between the eigenvalues induces further quantum entanglement.

Review: target space entanglement
In this section we will prepare the way for discussing matrix models by revisiting the question of entanglement for identical particles. The essential reason is that permutation symmetry also will arise in our treatment of matrix models. As a simple example, well discussed in the literature, consider the quantum mechanics of a single N × N Hermitian matrix which transforms in the adjoint of a U (N ) gauge group. We can partially fix the gauge by diagonalizing this matrix and writing a Lagrangian for the eigenvalues. This leaves a residual permutation symmetry which permutes the eigenvalues. This model can then be mapped to a system of non-interacting fermions (see for example [12,21]), in which the residual permutation symmetry of the matrix theory is the permutation symmetry of the identical fermions.
The essential point both for our discussion of matrix models and for other identical particle systems is that the permutation symmetry is a gauge symmetry: if one begins with a Hilbert space that is the product of identical factors, one projects the Hilbert space onto states that are invariant under permutations of these factors up to a sign. The problem is then to find a physically meaningful gauge-invariant notion of entanglement. To this end, the target space entanglement discussed in [17] and hinted at in a specific example in [22] introduces a natural type of entanglement which extends clearly to the matrix models we will study here. We will this focus on this definition.
To review, the essential confusion that arises with identical particles is this: Given N particles each described by a one-particle Hilbert space H, the physical Hilbert space is H ⊗N asym , where "asym" denotes a projection onto states antisymmetric with respect to the permutation group S N exchanging factors of H ⊗N . This projection amounts to imposing a gauge symmetry. One might have been tempted to break up the collection into two groups with k and N −k particles and compute a reduced density matrix. The antisymmetrization would seem to give a state which is entangled due to the antisymmetrization. But this split is not gauge invariant as S N permutes particles between these clumps.
There are of course clear notions of spatial entanglement for the quantum field theory of bosonic and fermionic fields, whose excitations correspond to identical particles: for example, one can defines entanglement between regions of the space on which the particles propagate. In essence that is what we will do here. However, the matrix model naturally lends itself to a first quantized picture, so we wish to understand entanglement in this language. Indeed, in discussions of identical particles, there are circumstances for which a natural first-quantized discussion is desirable. For example, entanglement or lack thereof between two identical spin-1 2 particles can be described in terms of realizable experiments, of the type attributed to Einstein, Podolsky, and Rosen (more precisely, the version developed by Bohm). This should be describable without quantizing the fermionic field [23].
To this end, several frameworks have been discussed in the literature to define some notion of entanglement, based on identifying a subclass of observables. Our work will follow from constructions of entanglement based on a subalgebra of observables. We will follow the specific discussion in [17] (see also the references in that work) which leads to a notion of entanglement that naturally include EPR-type experiments. A formal definition of entanglement based on identifying operator subalgebras is also given in [22,24].
The essential lesson is this: given a subalgebra A of the algebra of observables, one can prove that the Hilbert space decomposes into a direct sum of tensor factors, for which the observables A act as the identity operator on the second factor. Given a (pure or mixed) state in H, we can construct a reduced density matrix that encodes all expectation values of observables in A, such that it is a sum over reduced density matrices ρ i in each subfactor, weighted by a a probability p i . One can the define an entanglement entropy with contributions from a 'classical' piece − i p i ln p i , and a 'quantum' piece S q = − i p i Trρ i ln ρ i which comes from a weighted sum over von Neumann entropies in each subspace. In [17] in particular, the authors describe observables which are constructed as combinations of |x x | where x, x lie in a specified subregion of the configuration space of the fermions. In the first-quantized language, the index i then corresponds to the number of fermions that sit in this region. If we further project the state down onto a specific i before computing an entanglement entropy, we clearly describe EPR-type experiments in which each observer is observing a definite number of particles. This is the approach we will take. Since our goal is to describe matrix models which encode D-brane dynamics, it would encode entanglement between clumps of D-branes separated in space, with fixed number of D-brane in each clump. In this setup one can discuss entanglement in an intuitive way, if we understand that we are looking at the entanglement of measurements in two spatial regions. In essence we are using the position basis to "identify" the identical particles. Note that this approach, with a discrete "position" label, was also discussed in an example in [22].
In this section we will review and discuss the target space entanglement of fermions in this way, in the case that space is a finite lattice. We show how we factorize the Hilbert space for a fixed number of fermions at each position. This is basically working out simple examples of target space entanglement, an exercise we have found useful in laying the groundwork for entanglement in gauged multi-matrix models.
An alternative treatment of identical particles [25] makes use of an appropriate subspace spanned by a set of observables, rather than a full subalgebra. This has some antecedents in the quantum chemistry literature. The provides a basis for "entwinement" in quantum field theories with a gauged permutation symmetry, a quantity which has an appealing holographic dual [25][26][27]. We leave for the future any application of this approach to gauged multimatrix quantum mechanics.

Examples
The general results discussed above state that the Hilbert space of particles with a permutation symmetry breaks up into a direct sum of factorized Hilbert spaces; thus, by projecting onto a summand, we can discuss entanglement in the canonical way. On the way to applying this to the study of gauged matrix models, we have found it useful to work very explicitly through some simple examples involving identical particles on a finite lattice.

Two spin-1 2 particles
We begin with the simple case of two indistinguishable spin-1 2 particles. Our goal is to provide a concrete example in language that can be generalized to systems with a larger number of particles. Thus, we will be somewhat pedantic in the hopes that this will pay off for the reader when we discuss more general examples.
In the usual EPR-type thought experiment one considers two spin-1 2 particles which are prepared locally in some state, and then physically separated. The observation that the particles are or are not entangled is then based on first making local measurements of each spin, and then comparing the results. This gives an operational meaning to entanglement between indistinguishable particles: by comparing successful experiments made in each region, we have projected the state of the two particles onto a subspace in which two separate regions each have a definite particle numbers (that is, 1). In essence we have distinguished the particles by tagging them with their position, and we can construct reduced density matrices in the standard way.
We next restate this in more precise mathematical language. First, we can get at the essence of physical separation of particles by considering two particles each with spin-1 2 and with two possible position eigenstates (this model also appears in [22]). The position quantum numbers are written as X = ±1 ≡ ±. A basis of the four-dimensional Hilbert space of single-particle states is |αa , where α =↑, ↓ denotes the eigenvalue of S z / (with ↑−→ S z = 2 , ↓−→ S z = − 2 ), and a = ± the eigenvalue of the position operator X.
Before imposing the permutation symmetry, the Hilbert space of the "ungauged" twoparticle system is 16-dimensional. Projecting onto an antisymmetric wavefunctions leads to a smaller six-dimensional subspace. We can write as a basis |αa 1 |βb 2 − |βb 1 |αa 2 , although this includes zero vectors, and so is not the most compact notation. On the other hand, we can separate the spin and position degrees of freedom, so that beginning with basis vectors |ab |αβ ≡ |αa 1 |βb 2 , the antisymmetric states are: Let us consider further the subspace of states for which the two particles are physically separated, which is the image of the projection P = P pos ⊗ 1 spin where P pos = | + − + − | + | − + − + |. This leaves a four-dimensional Hilbert space spanned by ψ 3,4,5,6 . At this point we could simply fix the gauged permutation symmetry by placing particle 1 at position − and particle 2 at position +. Nonetheless, we find it illuminating to provide a more manifestly S 2 -invariant setup. Working with the basis ψ 3,4,5,6 , we seek a natural factorization C 4 = C 2 ⊗ C 2 corresponding to the usual EPR notion of entanglement in which we measure spins that are associated to specific position quantum numbers. To this end, we introduce the operators where we have written these in the same basis as (2.2), in which the first two factors label the position and the second two factors label the spin. These operators measure the spins associated to specific position states. They are diagonal in the basis spanned by: in which they can be written as O +z = diag(1, −1, 1, −1) and diag(−1, 1, 1, −1) respectively. One can check quickly that the operators do not have support outside this subspace: they annihilate ψ 1,2 . Within this subspace, O ±,z are each doubly degenerate but together form a complete set of commuting observables: their eigenstates thus naturally describe It is straightforward to check that the states |k are the four states in which each position can be associated to a definite value of the spin. From these we can take linear combinations to construct nontrivial EPR pairs. We can simplify the notation further by labeling the states as |O +,z , O −,z according to the eigenvalues of O ±,z . Here the first label is the spin at the position + and the second the spin at the position −. In this language, and the basis is that of two distinguishable spin-1 2 particles in C 2 ⊗ C 2 , for which standard notions of entanglement can be defined. In particular linear combinations such as describe the entanglement of separated particle in a way that is gauge-invariant and intuitive. It is clear, for example, that observations with O +,z alone are described by a mixed state with density matrix equal to the identity operator on C 2 ; this operator can be computed via a standard partial trace.
2.2.2 3 spin-1 2 particles As a second example, let us work out the example 3 spin-1 2 particles on a three-site lattice. The position basis is |X = 1, 2, 3 ; for spin we again take the basis |m = ± 1 2 of eigenstates of σ z . The full Hilbert space is Following the discussion above of "target space entanglement", we wish to consider the entanglement of a particle at X = 1 with two other particles at X ∈ {2, 3}. The corresponding subalgebra of observables is that generated by where µ ∈ {0, 1, 2, 3}, σ 0 is the identity operator on C 2 , σ 1,2,3 are the Pauli matrices acting on C 2 , and P 1 = |1 1| is the projector onto X = 1 in the Hilbert space of positions. The operators O 1µ generate 1− and 2− particle observables for fermions located at X = 1. The Hilbert space H can be repackaged as where H (k) X=i 1 ...ip corresponds to the k-particle Hilbert spaces for particles restricted to the sites i 1 . . . i p . This factorization is manifest in the occupation number basis, in which we label the states by the number of particles for each value of X and σ, with the total number of particles set to N = 3 and the Pauli exclusion principle imposed. In this basis, it is convenient to work in second quantized language, with fermionic annihilation and creation operators a iα , a † iα , i 1 . . . 3, α = ± 1 2 . A complete basis of states consists of all non-vanishing combinations where |0 is the fermionic vacuum. Note that in this language, The basis states |i 1 α 1 ; i 2 α 2 ; i 3 α 3 will have zero entanglement between particles localized at X = 1 and particles located at X = 2, 3: entangled states are linear combinations of these basis states.
The unentangled basis |i 1 α 1 ; i 2 α 2 ; i 3 α 3 can also be written in position space as: This latter presentation is closer to the presentation that naturally emerges from considering gauged matrix models. Note that if we were to embed this state into a state of three distinguishable particles in the Hilbert space , and compute the entanglement of the particle in the first factor with the others, we would find that the first particle is entangled with the other two. This is not target space entanglement: for each term in the sum the particle in the first factor of C 3 is at a distinct location. More importantly, it is not gauge invariant (where we take the action of S 3 as a gauge symmetry). A gauge-invariant quantity is the target space entanglement of a single particle at X = 1 with two particles at X ∈ {2, 3}. We could compute this within the space C 3 position ⊗ C 2 spin ⊗3 = C 15 by embedding the σ = 1 term in (2.13), considering linear combinations of such states, and computing the entanglement between the first particle and the other two. This amounts to a gauge fixing with respect to S k=3 , to a gauge in which the unentangled nature of the basis vectors (2.13) is clear. In practice this is what we will do below, but we emphasize that we could phrase our entanglement in such a way that the S k gauge symmetry is manifest.
In physical systems energy eigenstates of N particles may, depending on the Hamiltonian, be more naturally constructed from states with definite properties under S N for the position and spin degrees of freedom separately. The symmetry types are coupled by the demand that the total wavefunction be antisymmetric. To see how this works out, we follow the textbook construction of these states via Young diagrams. (See for example [28] for a clear discussion). We will see that such wavefunctions are entangled in the target space.
Recall that an irrep R of S N is determined by a Young diagram with N boxes. From each of these diagrams one may construct a Young operator Y acting as a projection operator Given a basis |i = 1, . . . n of H, the elements of H R correspond to N -rank tensors of fixed symmetry type: To construct a basis in H R , we choose k numbers from {1, . . . n} and fill in the Young diagram to create a standard tableaux so that the numbers are nondecreasing along rows and strictly increasing along columns. Reading these numbers from left to right along the top row, then the second row, etc, we assign them in turn to i 1 . . . i n ; we then act on |i 1 . . . |i N with the associated Young operator.
In our examples with two quantum numbers, the position and spin wavefunctions must be constructed from conjugate Young diagrams, with the rows and columns transposed, for the total wavefunction to be antisymmetric. The specific formula for combining the wavefunctions can be found, for example, in [28] (section 7-14).
Let us consider the case that the positions of the particles are distinct, and the individual particles have spin- Here the counting is understood by the number of ways of filling up the boxes with α = 1, 2 (corresponding to spins m = − 1 2 , 1 2 ) for spin or n = 1, 2, 3 (corresponding to x 1 , x 2 , x 3 ) for positions, so that the numbers are non-decreasing along rows and increasing along columns. For example, in the first line of 2.16, the boxes corresponding to completely symmetric spin wave function can be filled up in four different ways (1, 1, 1), (1, 1, 2), (1, 2, 2), (2, 2, 2) where 1 refers to spin up and 2 refers to spin down.
Let us again focus on states for which each particle is at a distinct position. The completely asymmetric state in line 2 of 2.16 survives and gives rise to four wavefunctions. Among the eight states of mixed symmetry, only two -(1,2,3) and (1,3,2)-survive and contribute a total of four wave functions so that the dimension of the resulting Hilbert space is 8 = 2 3 . We can use the basis described in the previous section to write these wavefunctions; this is a nontrivial change of basis, so that the states constructed as above from Young tableaux are naturally entangled from the point of view of the target space. We first define the completely antisymmetric position wavefunction and a partially antisymmetrized wavefunction respectively as follows: |r 1 , r 2 , r 3 as = 1 √ 6 (|r 1 , r 2 , r 3 − |r 2 , r 1 , r 3 + |r 3 , r 1 , r 2 − |r 3 , r 2 , r 1 −|r 1 , r 3 , r 2 + |r 2 , r 1 , r 3 ) where we the position labels take the values r 1 , r 2 and r 3 . Similarly, we can also define the position eigenfunctions symmetric and antisymmetric respectively w.r.t exchange of particles 2 and 3, Using these basis states for the position and spin wavefunctions we can write a complete basis for the full Hilbert space of spin-1 2 particles in distinct positions: Let us consider the state ψ 3 . In the occupation number basis discussed above this corresponds to the state The corresponding reduced density matrix for the particle at x 1 is: That is, there is nontrivial target space entanglement between a particle at any fixed position, and the two particles at other positions. We could continue in this vein with more particles, and more generally partitioning the "target space" between two groups, with fixed numbers of particles in each group; these particles may or may not be coincident in position within either group. The upshot will be the same: the occupation number basis will naturally yield states that are unentangled from the standpoint of target space entanglement, while the states constructed from irreps of the permutation group acting on each of the quantum numbers will naturally yield entangled states.

Entanglement in gauged matrix models
In this section, we set up a framework to compute target space entanglement in low-energy states of simple bosonic gauged matrix models of two and three matrices. We view this as a step to understanding entanglement in large-N theories relevant for holography, for example those which arise as the strong-coupling, low-energy limits of D-brane dynamics. The gauge theories under anything like computational control on either side of the holographic duality are typically supersymmetric, contain more bosonic matrices, and often extended spatial directions. In particular, supersymmetry allows for low-energy states of well-separated Dbranes, and as such ensures the consistency of the Born-Oppenheimer approximation we use to describe these configurations. For the sake of clarity, we will work with the reduced bosonic models, and assume the cancellation of the ground state energy occurs. At the lowest order in which we would require this effect, the extra degrees of freedom will simply contribute additively to the entanglement .

Outline of an entanglement entropy calculation
In order to motivate the specific form of our simplified models, we start by recalling the action [9] for the BFSS matrix model, which can be attained by dimensional reduction (in flat space) of the 10-d maximally supersymmetric SYM action to 0+1 -dimensions: This describes the massless sector of open strings ending on N D0-branes in 9 + 1 flat spacetime dimensions, with spatial coordinates on this spacetime as Y i=1,···9 , and with target space and worldline time identified. It is well known that at sufficiently low energies, closed strings as well as oscillator modes of open strings decouple. The theory has a U (N ) gauge symmetry, with all fields transforming in the adjoint of the gauge symmetry. The index (i) in parentheses is a "flavor" index, identified with the index for target space coordinates. Thus the bosonic fields transform as SO(9) vectors. We will use unbracketed subscripts a, b = 1 . . . N to refer to the matrix elements, on which the gauge symmetry acts. Thus in index notation we can write the bosonic fields as X (i)ab . The 0+1-dimensional gauge coupling e 2 = g s , where g s is the dimensionless string coupling. The natural dynamical length scale in the problem is l 11 p = (g s ) 1 3 l s [29,30], the elevendimensional Planck scale; this denotes the size of bound states at threshold.
We will need to pay careful attention to two features of this theory. First, the gauge symmetry must be treated properly. For the specific questions we will answer, we will do this by choosing a gauge appropriate to answering those questions, and imposing the residual gauge invariance on quantum states. Secondly, the fermionic variables in the theory are crucial to understanding D0-brane dynamics. In particular, they ensure the existence of low-energy states of well-separated branes. That said we will for the most part ignore them, and focus our studies on entanglement in reduced bosonic models. we believe the lessons we draw from the bosonic sector should generalize when fermions are included. While fermions give a negative contribution to the vacuum energy, they generally contribute positively to entanglement.
We imagine low-energy states which are well localized about configurations far along the flat directions of the potential in Eq. (3.1), corresponding to well-separated D-branes. The separation of the D0-branes amounts to a Higgsing of the gauge symmetry U (N ) → U (I) N . It is well known that in this 'classical' regime, the diagonal entries of these matrices can be interpreted as the coordinates of the D0-branes in the 9+1-dimensional bulk theory. The off-diagonal variables represent the fields for strings stretching between the D0-branes. Fermionic variables will cancel the zero-point energy of the off-diagonal strings, and ensure there is no static potential between the D0-branes.
We must be careful about gauge invariance, but this identification is close to the gauge-invariant reality if the branes are well-separated. We could consider explicitly gaugeinvariant variables by writing where U is a unitary matrix sometimes denoted the "angular components" and Λ (i) is the real, diagonal matrix of eigenvalues. In the limit of well-separated branes, the degrees of freedom denoted by U (i) become massive with masses equal to the difference between eigenvalues; the eigenvalues of Λ (i) are gauge-invariant up to discrete permutations, and can be approximately identified with the diagonal elements X ii . The discrete permutations reflect the fact that the D0 branes should be treated as identical particles.
We will work with a specific target space entanglement calculation that lends itself to a simple gauge fixing. We consider states in which n D0-branes are localized in the region Y 1 < 0 and the remaining N − n localized in the region Y 1 > 0. We can effectively use the position of the branes to label them, as in §2. We can then ask the question: 'how entangled are the branes and strings within the region X (1)ii < 0 with the strings and branes in the region X (1)11 > 0?'.
We begin by fixing to the axial gauge A t = 0. Once done, one still has to impose the condition that δL/δA t = 0; this condition is identical to the demand that the wavefunction Ψ(X (i) ) is invariant under U (N ) transformations, where U is a U (N ) matrix: to see this, consider infinitesimal transformations in Eq. (3.2).
For the question at hand, it makes sense to use this invariance to diagonalize X (1) : the resulting components of the matrix are precisely the eigenvalues of that matrix in the full theory. Note that this does not completely fix the remaining gauge invariance, as this choice for X (1) is unchanged by a residual U (1) N S N when the eigenvalues are distinct. For simplicity we will focus on this case. Under the residual gauge group, the matrix elements Physical states in this theory will have charge zero under this residual U (1) N . So long as our entangling surface is planar (and the target space flat), our gauge choice is natural.
As is well known, this gauge fixing of X (1) means that the measure of the path integral includes a Vandermonde determinant for the eigenvalues of X (1) . This determinant can be included in the wavefunction so that the measure on the eigenvalues is flat. The result is that under actions of the S N Weyl group, the wavefunction should transform with a sign equal to the sign of the permutation. (See appendix A as well as [14,31].) If the eigenvalues of the matrices are distinct, we can fix the permutation symmetry by ordering the eigenvalues so that X 11 < X 22 < . . . X N N .
Our desire is to compute the entanglement entropies (of von Neumann or Rényi type) between degrees of freedom localized on each side of Y 1 = 0. The negative elements of X (1) define an n × n block (after relabeling the indices), corresponding to strings stretching between the visible branes. We consider the entanglement of all of the matrix elements of X (i) within this block with the remaining degrees of freedom.
In doing so, there are off-diagonal matrix elements which couple the visible and "hidden" blocks of the matrix fields. These correspond to string fields for open strings stretching across this entangling surface. We must decide how to treat them. In string language, these strings are labeled by which D-branes they begin and end on. In the full string theory, the open strings would have internal degrees of freedom controlling their vibrational and rotational states, about which at least partial information could be gained via closed string scattering experiments done in the Y (1) < 0 region. 2 In the low-energy limit we consider, such information is invisible to us. What we can record is the total U (1) n charge carried by strings stretching across the entangling surface. This must be cancelled by additional strings stretching between the D0-branes localized in the Y (1) < 0 region, in order to satisfy the Gauss' law constraint that the total gauge charge be zero. As we will discuss in detail below, the full Hilbert space will break up into superselection sectors labeled by this U (1) n charge. It is possible that one could get additional information about these string fields and their coupling to the other hidden degrees of freedom, by measuring the force they exert on the D-branes at Y (1) < 0. We leave this possibility aside for now. The upshot, as we will describe in more detail in §3, is that the density matrix for the visible degrees of freedom will break up into pieces with different external U (1) n charge, leading to a "classical" contribution to the entanglement entropies entirely which is analogous to the classical contribution that appears in the Extended Hilbert Space formulation of entanglement in lattice gauge theory [18][19][20].
While we have discussed our problem in terms of target space entanglement of D0-branes, which gives a clear geometric picture, we can phrase the story more directly in terms of the gauge theory. Starting with a U (N ) gauge theory with adjoint matter, we consider configurations which spontaneously break the gauge symmetry according to the pattern We are interested in the entanglement of the degrees of freedom that transform as adjoints of the U (n) factor. The density matrix will break up into superselection sectors based on the total U (1) n charge carried by degrees of freedom which are "bifundamental" under the group We will now proceed to a more explicit discussion in terms of reduced models.

Entanglement in a two-matrix model
We first consider a bosonic theory with just two Hermitian N × N matrices, given by the action, where D t is the gauge covariant derivative for a U (N ) gauge group which acts on the bosonic Hermitian matrices X, Y in the adjoint representation. This is a crude model for N D0-branes in 2 spatial dimensions. Up to gauge transformations, the degrees of freedom of this theory are the eigenvalues and angular directions of these matrices. Following the discussions in §2 and §3.1, we wish to compute entanglement entropies for a given state after projecting onto the subspace for which k < N eigenvalues of X have negative values, corresponding to k D-branes localized at Y (1) < 0.
As above, we choose axial gauge and further solve for the remaining gauge invariance by gauging away the angular directions of X. We restrict our states by demanding they are invariant under the residual U (1) N gauge invariance, and sum over the images of this state under actions of the Weyl group. For a given image, the k negative diagonal components X define a block; we then construct the reduced density matrix for all degrees of freedom of X, Y that live within that block, breaking it up into superselection sectors corresponding to the U (1) k charge of the strings that are bifundamental with respect to the decomposition G k × G N −k with G n = U (1) n .
We will be interested in constructing the reduced density matrix for states that are models of low-energy states of D0-branes. To this end we proceed by writing out the Hamiltonian more explicitly. We will use lower case letters to denote the diagonal matrix elements of X and Y , and denote them by a single index. For example, (X) 11 := x 11 = x 1 .
The Hamiltonian corresponding to Eq. (3.4) can now be written as: where we have used the Hermiticity of Y to equate Y ij = Y * ji , π Y ij = π * Y ji . We have separated the Hamiltonian into "slow" and "fast" parts. The interaction terms make the off-diagonal elements of Y massive, with a mass that scales as the separation between the two branes. For the states localized around sufficiently well-separated values of x i , modeling D-branes separated by sufficiently large distances, the diagonal degrees of freedom have slow dynamics compared to that of the "heavy" off-diagonal matrix elements, which justifies the use of the Born-Oppenheimer approximation to construct low-energy states (though as we will see, we will need to include fermionic degrees of freedom to make this approximation quantum-mechanically consistent.). We also note that the off-diagonal components Y ij and Y ji are not independent degrees of freedom as they are complex conjugates of each other.
Note that this Hamiltonian has a large global symmetry group, corresponding to the phase rotation of each field Y ij . Recall that these complex fields are string fields; we can take excitations with positive charge under this summetry to be "strings" and with negative charge to be "anti-strings".

Example: N = 2
Let us begin with the simple case of N = 2, modeling 2 D0-branes. Here H f ast describes a single complex oscillator: (3.6) The interaction term shows that the potential energy is independent of y i , while for x i , only the relative separation x 1 −x 2 matters. These are consistent with the expectation that the physics of D0-branes in flat space should be invariant under target space translations. Note that the y i -independence of the action is just a feature of our gauge choice. We could have chosen a gauge where we diagonalize the matrix Y instead, and potential would be independent of x i but dependent on y 1 − y 2 . The next step is to construct interesting states with respect to which we compute entanglement entropies for appropriate subsets of the degrees of freedom. We will consider models of low-energy states of D0-branes. For such fixed x i = x j , Y 12 is a complex oscillator which has charges ±(1, −1) under the action of U (1) 1 ×U (1) 2 . When |x 1 −x 2 | is sufficiently large, the variables x i have slow dynamics compared to Y 12 , and we may use the Born-Oppenheimer approximation. For fixed x i the instantaneous energy eigenstates can be labelled as |n + , n − ; x 1 − x 2 , where n ± are the oscillator numbers with equal and opposite U (1) × U (1) charge.
After the aforementioned gauge fixing, we must still impose the residual U (1) 2 × S 2 gauge invariance. In particular, we need to write down wavefunctions that are odd under the S 2 action In this simple case, U (1) 2 invariance simply means that n + = n − = n, and we will collapse these to describe the instantaneous states as |n, x 1 − x 2 . This also guarantees that the offdiagonal sector is even under S 2 ; thus, the diagonal sector is odd. In essence the D-brane can be treated as fermions. Thus, we can write low-energy states of the total Hamiltonian Eq. (3.5) as: Again, we have absorbed both the Vandermonde determinant and the integral over the gauge group in the amplitudes ψ 0 ,δψ n (see appendix A for more details). These amplitudes can be general 2-body wavefunctions with unit norm; we write out the antisymmetrization explicitly, with the factor of 1/ √ 2 out front yielding the correct norm for the full theory. That said, for the purposes of computing target space entropy, o, we can simply fix gauge and keep only the first term in the above sum so long as x 1 = x 2 , dropping the factor of 1/ √ 2 In the limit that |x 1 − x 2 | l s , we expect ψ n>0 to be small compared to ψ 0 , of order ∼ e 3/2 s δx 3/2 for eigenvalues separated by δx ≡ |x 1 − x 2 (see Appendix B). At leading order, we set c n = 0. To find the dynamics of the slow modes, we first solve the Schrödinger equation for the fast modes with fixed values of x i , y i and compute the energy E 0 (x 1 − x 2 ) of the instantaneous ground state; this appears as an additional potential term in the effective Hamiltonian for ψ 0 , At this point we have a problem. The linear potential that appears in our purely bosonic two-matrix model means that there are no good wellseparated states with energies low compared to the oscillator modes. This problem can be cured by supersymmetry; the fermionic degrees of freedom contribute terms to the effective potential which cancel the static potential, leaving a weak velocity-dependent force. We will assume that such a cancellation exists, without including the degrees of freedom that lead to it. In a weakly coupled system, we expect fermions to contribute additively to entanglement, and that the qualitative lessons regarding the role of "fast" modes in entanglement calculations will carry over to a more precise treatment. Now that we have set up the Born-Oppenheimer framework, we will describe our prescription to compute entanglement entropy for branes localized at Y (1) < 0. We project our space onto a subspace of the Hilbert space for which one of x i is negative and the other positive. Once this is done, the upshot of the discussion in our previous section is that we can fix the permutation gauge and do our calculation with one image of the permutation group S 2 for which x 1 < 0. We next need to make a decision as to whether or not the degrees of freedom Y 12 are to be included in the observable sector. Following the discussion in §3.1, we assume they are not visible.
For the case of general N , after we project on a fixed number of X (1) eigenvalues satisfying x i < 0, the density matrix breaks up into superselection sectors corresponding to the U (1) n flux generated by excitations Y ij with x i < 0, x j > 0. We will discuss this further below. For the case of N = 2 with one negative eigenvalue, U (1) 2 invariance means this flux is always vanishing. In this sense the N = 2 example misses an important part of the entanglement calculation for large-N theories. Nonetheless, it isolates and captures another important feature, entanglement induced by Y 12 , so we will proceed.
Let us now try to study the entanglement between the "first" and the "second" brane when the fast modes are placed in the instantaneous ground state |0, x 1 − x 2 , so that we are (assuming the existence of additional fermionic modes which cancel the zero point energy) studying low-energy states of the system at leading order in the Born-Oppenheimer approximation. We will suppress the Y -coordinates throughout. We further assume that we can write low-energy states for which ψ 0 (x i , y i ) = ψ 1 (x 1 , y 1 )ψ 2 (x 2 , y 2 ). Our essential point is that the x 1 − x 2 dependence of the instantaneous ground state of the heavy modes already induces nontrivial entanglement..
Tracing over both the center of mass coordinates X 2 , y 2 of the second brane, and the string field Y 12 , the reduced density matrix for x 1 is (We ignore the y 1 -dependence as in our gauge, y 1 does not couple to Y 12 and the wavefunctions for x i , y i are unentangled). We can easily check that the trace of the reduced density matrix defined above is one for appropriately normalized wavefunctions ψ. An additional factor of 2 arising from the two terms in the antisymmetrized wavefunction give the same contribution to the density matrix is cancelled by a factor of 1/2 = 1/ √ 2 2 from the normalization of the appropriately symmetrized wavefunction.
, arises from the x 1 − x 2 dependence of the state |0, x 1 − x 2 of the "heavy" modes that we have declared unobservable. In the present case this yields non-trivial entanglement even for the case that ψ 0 factorizes. In order to see this let us compute the second Renyi entropy, also known as the purity of the state. This is simpler to compute than the von Neumann entropy but is a good a probe of the existence of entanglement. We choose for simplicity . Technically, note that ψ 1 has support on both sides of the entangling surface; however, this part of the wavefunction is exponentially small and we will ignore its effects. Also, technically the delta function for x 2 carries infinite energy; we are simply using this as a model for a well-localized wavefunction to simplify our analysis. The essential point is that the wavefunction for the diagonal components has support over a range of values of x 1 − x 2 . Since the state of Y 12 depends on this separation, Y 12 becomes entangled with x 1 .
In the end, in our approximation, the purity becomes: The integral can be evaluated using the saddle-point approximation if σ d. Performing thex 1 integral first, We have an overall factor of two as we are integrating over half the real line. Upon performing the integral over x 1 , again in the saddle point approximation, we find We have retained only the leading saddle x 1 ∼ −d in this integral. Note that this expression for the purity is something we could have guessed by dimensional analysis except for the numerical factor in the second term. We see that the the state is mixed as long as one of the branes is not held fixed at a location (if both branes are held fixed, σ → 0 and there is nothing to decohere). This entanglement is induced by the off-diagonal degrees of freedom Y 12 . From the point of view of the D-branes, entanglement between clumps of D-branes will be induced in part by the dependence of the state of the "off-diagonal" strings, stretching between these clumps, on the location of said clumps.
As we have stated, our states are only parametrically lighter than the mass of the offdiagonal modes if we are able to somehow cancel off the contribution of these modes to the effective potential for x i . In the BFSS model, this is a consequence of supersymmetry. In these cases, the dependence of the instantaneous ground state of the off-diagonal fermionic degrees of freedom will also contribute to the kernel K. We expect this will not cancel the entanglement induced by K. In simple free field examples, fermions contribute additively to entanglement, for example in the case of spatial entanglement in a 2d CFT, which scales with the central charge of the theory [33][34][35].
A more accurate calculation would include corrections at sub-leading order in the Born-Oppenheimer approximation: that is, corrections to the adiabatic approximation to the state of the off-diagonal modes. As explained in the appendix B, these corrections lead to even, pairwise excitations of the strings which leads to zero net flux on each string; a simple consequence of imposing invariance under the U (1) 2 gauge group which remains after we diagonalize X (1) .
Qualitatively, new effects occur when we extend our analysis to the U (N > 2) ,with k < N visible branes in the region Y (1) < 0. In this case, the density matrix will generally break up into "superselection sectors" as we have discussed above, corresponding to components of the quantum state that contain different configurations of unobserved off-diagonal matrix elements which carry nontrivial U (1) k charge. To understand these better, we will pass in the next section mapping our target space entanglement calculations onto a computation of entanglement in a gauge theory living on a complete lattice. Once we have explained this map, we return to the discussion of three matrices in BFSS theory which we will argue captures all the qualitatively essential features of target space entanglement in the full BFSS theory. .

N > 2: superselection sectors
For N > 2, we can have excitations of off-diagonal matrix elements such that the state of a given matrix element carries non-vanishing charge under a U (1) × U (1) subgroup of U (1) N , while the total U (1) N charge is vanishing. In particular, this means that the density matrix for k 'visible' eigenvalues breaks up into terms labeled by the total U (1) k charge carried by excitations charged under U (1) k × U (1) N −k . Let us explore this phenomenon further. Before projecting onto the S N -symmetric states -or if the eigenvalues of X 1 are distinct, after fixing the permutation symmetry by ordering the eigenvalues -the matrix quantum mechanical model can be treated as a U (1) lattice gauge theory on a 'complete' (all-to-all) N-site lattice, in the Extended Hilbert Space formalism of [19]. A natural physical treatment of this system, motivated by theories of D0-branes, leads to an 'electric center' prescription for treating the degrees of freedom that connect the 'hidden' and 'visible' degrees of freedom. We begin with N lattice points and N 2 links between them. The diagonal matrix elements X (i)kk 'live' on the kth lattice site, while the off-diagonal terms X (i)kl live on the links joining the kth and lth lattice sites. Each site k has an unbroken gauge symmetry U (1) k associated to it, with the diagonal matrix elements being neutral and the complex fields X (i)kl having charge ±(1, −1) under U (1) k × U (1) l . Each link can carry nontrivial U (1) 2 , so long as adjacent links also carry charge in such a way that the total U (1) N charge is zero. We can clearly organize the Hilbert space according to the amount of charge flowing through each link, with 'charged' links joining to form closed loops. Note that for the example discussed in the previous section, we have a lattice with two vertices representing the coordinates of the two branes; the total U (1) charge associated to the edge must have n + − n − = 0 (in string theory language, excitations will come in string-anti string pairs).
The upshot is that the Hilbert space is the U (1) N S N -invariant subspace of the "extended" Hilbert space where H V i is the local Hilbert space at each vertex and H ij is the local Hilbert space on each link. Invariance under U (1) N transformations leads to the conservation of total U (1) charge at each vertex. We must further impose the S N permutation symmetry of our matrix theory. We can impose this by writing states in a basis of S N orbits of U (1) N -invariant wavefunctions on H. This extends the "unentangled" basis in §2. Alternatively, if the eigenvalues of X 1 are distinct, we can choose a particular element of this orbit, for example the obe such that X 11 < X 22 < . . . X N N . Target space entanglement is gauge-invariant and we can use either treatment of S N to calculate it. Let us now consider k of the N − k lattice sites, and the associated links between them, to be 'visible', modeling for example D-branes that lie on one side of an entangling surface, and the string fields for open strings stretching between them. We can now trace over the degrees of freedom on the N − k 'hidden' vertices, as well as the corresponding links, as well as the string fields living on the links connecting the hidden and visible vertices. Since this a complete lattice all vertices are also boundary vertices. The resulting Hilbert space and the corresponding reduced density matrix break up into superselection sectors where the reduced density matrix ρ in is block diagonal and q = q 1 , · · · q k is the vector that denotes the total U (1) k charge due to excitations on links joining the visible and hidden vertices. The entanglement entropy associated to this state is given by The first term, which we refer to as the classical piece, is given by the entropy due to the probability to be in a particular superselection sector. The second piece, which we refer to as the quantum piece, is given by a weighted sum of the entanglement entropy of each superselection sector. No gauge invariant operator can act to change the superselection sector unless it acts on degrees of freedom both inside and outside the spatial cut (namely a Wilson loop crossing the boundary); in particular no observable supported on the "visible" degrees of freedom can change q.
We now proceed to illustrate our map in a simple example.

N = 4 example
Let us return to the case of two matrices. We might ask whether corrections to the leadingorder Born-Oppenheimer approximation would induce states with string excitations forming loops on the lattice, such that the net U (1) 2 gauge charge on each link is nonzero. However, as noted previously, the Hamiltonian 3.5 has an enhanced global symmetry corresponding to Y ij → e iφ ij Y ij . A single "string" associated to a given link carries nonzero charge under this symmetry. If the leading term in the Born-Oppenheimer approximation corresponds to the instantaneous vacuum, which carries vanishing charge, corrections to adiabaticity from x i being diagonal will not break this symmetry, and the aforementioned loop configurations will not be induced. In the language of D0-branes, open string would be created in string-antistring pairs, so that the excitations on each link carry vanishing gauge charge. As we will see below, for three or more matrices, the Hamiltonian no longer has this symmetry and even the instantaneous ground state for the heavy off-diagonal modes will contain components with loops of strings.
However, for illusrative purposes, we simply construct a state in the two-matrix thwory, for which the reduced density matrix for visible branes breaks up into superselection sectors labeled by the charge carried by the hidden degrees of freedom. This will be some highly excited state in the 2-matrix theory. Let us consider 4 "D0-branes" which are split into two groups separated by a large distance d, as illustrated in Figure (1). In order to focus on the classical contribution we will consider the situation where we have pinned down the branes to unique locations Y 1 , . . . , Y 4 , where Y k = (X kk , Y kk ) is a two-dimensional vector of diagonal matrix elements which we can think of as D0-brane position, and we are working in the gauge with diagonal X as before. We take Y 1 , Y 2 to lie on on the visible side of the entangling surface: that is, X 11 , X 22 < 0. In the figure, the sites of the lattice correspond to the D-branes; we have denoted the associated location by the number 1, . . . , 4 corresponding to the subscript of Y I , The links denote the "off-diagonal" matrix elements, labeled by position subscript at each end. The arrows between sites i and j define the meaning of the sign of the charge α ij ; an excitation with charge α ij = +1 on a given link carries positive charge at site j at the head of the arrow, and negative charge at site i at the tail of the arrow. The Figure denotes a particular configuration of off-diagonal excitations which satisfies U (1) 4 gauge invariance; that is, the total charge at each vertex is zero.
If the vertex labels correspond to target space position, the graph in Figure (1) then corresponds to a specific U (1) 4 × S 4 -invariant state. As before we can describe this in a manifestly S N -invariant fashion, or by fixing the S N gauge symmetry. In the former case, we start with the state |ψ = | Y 1 , . . . , Y 4 |α 12 = 3, α 13 = −2, α 14 = −1, α 23 = 1, α 24 = 2, α 34 = −1 (3.16) in the "extended" Hilbert space 3.13 and then sum over the S N orbit that includes this state, with appropriate signs: For example, the contributions to the sum from σ = 1 and σ = (12) are: Alternatively, we can fix the S N symmetry, for example by imposing the ordering X 11 < X 22 < X 33 < X 44 . With this ordering we can just work directing with (3.16) and compute entanglement in the standard fashion.
States of the form (3.17) have vanishing target space entanglement. On the other hand, if we begin with the extended Hilbert space state and sum over the orbit of S 4 , which acts on each term separately, or equivalently impose an ordering on X kk , we have a state which has target space entanglement. We perform the partial trace over all degrees of freedom except Y 1,2 , α 12 to find the reduced density matrix: where the (1 ↔ 2) appears if we have not fixed the S 4 symmetry. Each term in this sum corresponds to a distinct superselection sector with different amounts of charge flowing into the visible vertices. Each superselection sector appears with probability 1 2 ; the total classical entropy is thus S = ln 2.

Entanglement with three matrices
In the previous sections, we discussed the two-matrix model at length and introduced our framework for computing entanglement between two clumps of well-separated branes. Starting from the Born-Oppenheimer approximation for low-energy modes, with the offdiagonal terms in the instantaneous ground state, we found (see B for more details) that non-adiabatic corrections to the Born-Oppenheimer approximation included only stringanti string pairs associated to each off diagonal mode (a link of the complete lattice, in the language of §3.2.2.) Next, we consider the case of three bosonic matrices, for which new qualitative features arise in the entanglement structure of low-energy states. 3 Specifically, corrections due to interactions induce components of the quantum state with excitations that carry U (1) charge from the hidden to the visible degrees of freedom, inducing "classical" terms in the entanglement entropy. We will sketch the manner in which such components appear.
We start with the following action, As before, we can use the U (N ) symmetry to diagonalize the matrix X. Passing to the Hamiltonian, and assuming the eigenvalues of X are well-separated, we separate the dynamics of the "slow" diagonal elements and the "'fast", heavy off-diagonal elements: The interaction terms of H f ast within the round braces induces masses of order O(|x i −x j |) for the off-diagonal terms in the Y and Z matrices. The important new feature is the commutator term tr[Y, Z] 2 . In the limit of large |x i − x j |, it can be treated perturbatively. A simple analysis of how the perturbative corrections scale -see Appendices B, B.2shows that this term will induce corrections to the instantaneous energy eigenstates of order √ g ≡ el 3/2 s δx 3/2 and of order g 2 . These small parameters also control non-adiabatic corrections to the Born-Oppenheimer approximation, although the order √ g terms are further suppressed by a ratio of energy scales. A complete treatment will have to take care of both effects simultaneously. The new commutator term breaks the global symmetry of separate phase rotations of the fields on each link of our complete lattice. t Thus it can and does induce corrections to the instantaneous ground state with excitations of the off-diagonal modes of Y , Z around 3-and 4-link loops of our complete lattice: at leading order, these have single excitations on each link of the loops. When such loops connect vertex sites on each side of the entangling surface, the correction to the instantaneous ground state leads to a term in the reduced density matrix for the "visible" degrees of freedom in which non-trivial U (1) charge flows into visible lattice sites across the entangling surface. Figure 2 gives a visual picture of such corrections in the case of N = 4, with each numerical label on the vertex denoting a target space position.
In a generic excited state, the computation of entanglement breaks up into superselection sectors labeled by the total U (1) charge at each vertex. As explained in section 3.3, we get both a quantum piece and a non-trivial classical piece. In a computation of entanglement in the full BFSS model, we expect all the features described in this section to appear. We hope to understand these aspects at a more quantitative level in a future communication.

N scaling
We close with some preliminary comments on the scaling of the target space entanglement entropies we have been discussing, when the eigenvalues are well separated across an entangling surface. We imagine N 1 eigenvalues of X are negative, and N 2 = N − N 1 positive.
First, let us revisit entanglement induced by off-diagonal matrix elements which connect matrix eigenvalues on each side of the entangling surface. We imagine a case for which the state of these modes depends on the relative separation of the diagonal modes which give them mass, as in the N = 2 2-matrix example given above, and not on the "visible" off-diagonal modes. Let us further imagine a situation for which in the full quantum state, the eigenvalues on each side of the entangling surface have wavefunctions of the form ψ 1 (x 1 , . . . x N 1 )ψ 2 (x N 1 +1 , . . . , x N ), and the O(N 2 1 ) off-diagonal modes Y ij , Z ij for 1 ≤ i, j ≤ N 1 are similarly uncorrelated with those for N 1 + 1 ≤ i, j ≤ N . In this case the entanglement is carried entirely by the O(N 1 N 2 ) off-diagonal modes with 1 ≤ i ≤ N 1 < j ≤ N , From the standpoint of the visible branes, these modes entangle only the O(N 1 ) eigenvalues on one side of the entangling surface with the O(N 2 ) eigenvalues on the other. For a quantum system in a pure state, with M visible degrees of freedom and P hidden degrees of freedom, the von Neumann entropy is bounded from above by min(M, P ). In this case that minimum will be N 1 , even if N 1 N 2 . In other words, if the strings stretching across the entangling surface are not observed, they can provide a large amount of entanglement even if there are few eigenvalues with support in the "hidden" region.
Next, we consider the "classical" entanglement induced at weak coupling from components of low-energy states with excitations of strings that carry charge under the U (1) N 1 gauge group for the visible degrees of freedom as well as the U (1) N 2 gauge group for the hidden degrees of freedom. We work at lowest order in the Born-Oppenheimer approximation for which states of the off-diagonal matrix elements are taken to be their ground states. For the case of 3 matrices, the commutator term tr[Y, Z] 2 in Eq. (3.22) induces corrections to the ground state away from the simple harmonic oscillator vacuum. As explained in section 3.4, this term creates triangular and square loops which contribute to the first-order corrections to the wavefunction: where δx is the characteristic separation between eigenvalues, and we have used the results in Appendix B.2 to estimate the size of the corrections.  The pattern of excitations of off-diagonal matrix elements that appear in corrections to the ground state at leading nontrivial order are shown in Figure 3. Tracing over everything except the visible degrees of freedom (I, J), we obtain a density matrix of the form, where we have introduced the 't Hooft coupling λ 2 = e 2 N 2 . The reduced density matrix has O(N 2 1 ) sectors, thus, in the limit λ 2 1, δx l s , the leading "classical" correction the the von Neumann entropy will be of order

Future directions
In this work we have discussed the computation of target space entanglement for gauged multimatrix models, in particular at weak coupling when the eigenvalues are well-separated. There are many important directions for further research, of which we list a few below: N scaling A more complete and precise accounting of the N -scaling of target space entanglement entropy for two separated "clumps" of eigenvalues is needed, in particular for understanding the holographic duals of this entropy in the situations studied for example in [6][7][8].
Entanglement and bulk surfaces An obvious question is whether, in the limit of strong coupling and large N , the target space entanglement of matrix degrees of freedom has a dual holographic interpretation as the area of a bulk surface. This issue has been explored in several papers including [6][7][8][13][14][15], but there is as yet no conclusive answer. This is an area in which working with theories containing extended spatial directions, far out on the Coulomb branch, may be a useful strategy.
SUSY A more precise treatment of the fermionic degrees of freedom in a genuinely supersymmetric matrix model is an important direction for future work.
Dynamics In this paper, we have mainly focused on static entanglement measures in states with well separated D0-branes. It would be interesting to study entanglement dynamics in a scenario where D0-branes scatter off each other [37]. In the scattering region, the Born-Oppenheimer approximation breaks down and open strings stretching between D0-branes are created, introducing additional entanglement between groups of D0-branes.
Multipartite entanglement In this paper we have focused only on bipartite entanglement; it would be valuable to consider multipartite versions of target space entanglement in these models.
We have successfully used the U (N ) gauge symmetry to diagonalize one of the matrices. This leaves a residual U (1) N × S N symmetry, where the permutation symmetry S N acts on the labels of each matrix element under which the wavefunction Ψ is invariant. Imposing this implies that Ψ picks up the sign of the permutation. For instance, for a permutation σ, Ψ(λ 1 · · · λ N ; X 2 , · · · X N ) = sign(σ)Ψ(λ σ(1) · · · λ σ(N ) ; X σ(2) , · · · X σ(N ) ) (A.4) Let us now specialize to the case of N = 2 and two matrices. Expanding the state using Born-Oppenheimer approximation, where we have denoted the diagonal entries of the two matrices X and Y with a single index and in small cases letters to emphasize the fact that we think of them as coordinates. The above expression can be easily generalized to the case of N branes, It is understood that we need to take a tensor product of the both the fast and slow basis under the usual product. The harmonic oscillator states are normalized as usual n ij ; x i − x j |m kl ; x k − x l = δ ik δ jl δ m ij ,n kl and thus, we can define a normalized state.

B Corrections from nonadiabatic and finite coupling effects
In this Appendix, we estimate the size of non-adiabatic and finite-coupling corrections to the leading-order Born-Oppenheimer treatment, for low-energy states of the two-matrix model. We will see that the natural small parameters are: • ≡ E slow /E f ast , measuring the strength of the adiabatic approximation. Here E slow is the energy of low-energy scattering states of the matrix eigenvalues, while E f ast is the mass of excitations of the off-diagonal components.
• g ≡ e 2 ( s /δx) 3 , the strength of quartic interactions, where δx is the separation between eigenvalues on each side of the entangling surface.
• We will also find corrections of order √ g .
We take these parameters to be small, and show below that the corrections to the leading Born-Oppenheimer approximation are self-consistently of order g, , √ g .

B.1 Sub-leading Born-Oppenheimer corrections
We write the n = 0 and n > 0 terms in Eq. (3.8) as: |ψ = |ψ 0 + |δψ n>0 . The first term is the state constructed at lowest order in the Born-Oppenheimer approximation; the second term captures the "non-adiabatic" corrections. We will consider solutions to the timeindependent Schrödinger equation, assuming that the low-lying states lie in a continuum of scattering states. We will thus fix the low-lying energy; the goal of computing corrections to the Born-Oppenheimer approximation is to find a more accurate projection onto states with energies smaller than the masses of the "fast" modes. The time-independent Schrödinger equation Eq. (3.5) acting on the first term in Eq. (3.8) leads to: (H − E)|ψ 0 = 1 √ 2 dx 1 dx 2 dy 1 dy 2 e 2 l s 2 (π 2 x 1 + π 2 x 2 + π 2 y 1 + π 2 y 2 ) + 1 l 2 s |x 1 − x 2 | ×ψ 0 (x 1 , x 2 , y 1 , y 2 )|x 1 , x 2 , y 1 , y 2 s |0; Here, we have already taken H f ast to lie in the instantaneous ground state; the term proportional to |x 1 − x 2 | ≡ δx in brackets is the eigenstate of H f ast in Eq. (3.5). As discussed in the text, this term would spoil the assumption that we can have low-energy states for well-separated eigenvalues of X. In supersymmetric models such as [9], this vacuum energy is cancelled by fermions; while we do not include the fermions here, we will drop this term in the rest of our discussion. (The presence of fermions should not change our estimate of the size of corrections to the Born-Oppenheimer approximation, so this is the extent to which we will consider the effects of fermions here). Next, inserting resolutions of identity in the above equation, we find: dx i dy i dz i dw i ψ 0 (x 1 , x 2 , y 1 , y 2 )|z 1 , z 2 , w 1 , w 2 s × w 2 , w 1 , z 2 , z 1 | e 2 l s 2 (∂ 2 x 1 + ∂ 2 x 2 + ∂ 2 y 1 + ∂ 2 y 2 )|x 1 , x 2 , y 1 , y 2 s |0; δx f − (1 ↔ 2) −E|ψ 0 = −(H − E)|δψ .

(B.2)
We can now use the fact that s z i | − ∂ 2 x |x i s = −∂ 2 x i δ 2 (x i − z i ) and integrate by parts to find: where ∂ 2 x i := ∂ 2 x 1 + ∂ 2 x 2 , ψ 0 (x i , y i ) := ψ 0 (x 1 , x 2 , y 1 , y 2 ) and ∂ x i ψ 0 ∂ x i := ∂ x 1 ψ 0 ∂ x 1 + ∂ x 2 ψ 0 ∂ x 2 . If we take the inner product of this equation with |x i , y i antisym |0, δx := 1 √ 2 (|x 1 , x 2 , y 1 , y 2 |0, δx − (1 ↔ 2)), the middle line gives the Berry connection. When the ground state wavefunction for the heavy modes is real, it is well known that this term then vanishes. 4 The final term is of higher order in the Born-Oppenheimer approximation. To leading order, then, the wavefunction ψ 0 (x, y) satisfies the free wave equation. We can write these low-energy modes in terms of eigenstates of the target space momentum operator −i∂ x i , −i∂ y i . If the characteristic wavenumber is k, the energy is E slow = e 2 s k 2 /2, characteristic of a free particle with mass M = 1/(e 2 s ). The Born-Oppenheimer approximation should be valid when E slow E f ast ∼ δx/l 2 s Imposing this, the remaining terms are higher-order in the Born-Oppenheimer approximation (reflecting corrections to adiabaticity). Upon imposing the leading-order equations we are left with the subleading terms: We wish to find δψ n at the first subleading order in the corrections to adiabaticity. Thus, we ignore the terms in the second line of (B.4) for which the derivatives act on |n; δx . Taking the inner product of this equation with δx, n| x i , y i | antisym , and solving formally for δx, we find: We can further simplify our estimate by replacing the denominator by its dominant term. We claim that to leading order the term nδx/l 2 s ∼ E f ast dominates. It is larger than E ≡ E slow by a factor of . The derivatives will act either on ψ to bring down a factor of k ∼ 2E slow /e 2 l s or on the matrix elements of the fast modes, bringing down a factor of 1/δx. This gives terms of order g, , √ g relative to the dominant term. Thus, to leading order in these parameters, we find: δψ n = 1 E f ast e 2 l s 2 ψ 0 (x i , y i ) n, δx|(−∂ 2 x i − ∂ 2 y i )|0, δx +e 2 l s (∂ x i ψ n, δx|∂ x i |0, δx + ∂ y i ψ 0 n, δy|∂ y i |0, δx ) (B.6) Applying the same principles as above to estimating the size of derivatives, the first term will give a correction of order gψ 0 /n, and the second, a correction of order √ g ψ 0 /n.

B.2 Perturbative corrections to "fast" modes
In computing the ground state of the fast modes for 3 or more matrices, we have in addition to the harmonic oscillator terms the commutator term tr[Y, Z] 2 in Eq. (3.22). This contributes the following classes of terms: • Additional contributions to the masses of off-diagonal terms: δH mass = 1 e 2 l 5 s i<j (y i − y j ) 2 |Z ij | 2 + (z i − z j ) 2 |Y ij | 2 (B.7) • Terms cubic in off-diagonal matrix elements, such as • Terms quartic in off-diagonal matrix elements.
For simplicity, we will assume that all of the eigenvalues are well-separated by distances of order δx l s . Ignoring the cubic and quartic interactions, the off-diagonal modes have the Hamiltonian of simple harmonic oscillators with mass M = 1/e 2 l s and frequencies ω ∼ δx/l 2 s . The characteristic spread of the wavefunction for the off-diagonal modes is Recall that for a give perturbation δH to the Hamiltonian, the ground state acquires a component proportional to the unperturbed excited state |n with coefficient where |0 is the ground states, and E n −E 0 is the energy difference between the unperturbed energy eigenstates: here the energy difference is ∼ nω. We estimate matrix elements of the form n|ZY Z|0 ∼ L 2 , n|ZY Z|0 ∼ L 4 . With this, the cubic interaction give corrections of order O( √ g), and the quartic interactions give corrections of order O(g)