From Black Hole to Qubits: Evidence of Fast Scrambling in BMN theory

BMN Matrix theory admits vacua in the shape of large spherical membranes. Per- turbing around such vacua, the setup provides for a controlled computational frame- work for testing information evolution in Matrix black holes. The theory realizes excitations in the supergravity multiplet as qubits. These qubits are coupled to matrix degrees of freedom that describe deformations of the spherical shape of the membrane. Arranging the ripples on the membrane into a heat bath, we use the qubit system as a probe and compute the associated Feynman-Vernon density matrix at one loop order. This allows us to trace the evolution of entanglement in the system and extract the characteristic scrambling timescale. We find that our numerical analysis is consistent with this time scaling logarithmically with the entropy of the qubit system, in tune with suggestions by Sekino and Susskind.


Introduction and Highlights
Tracing the evolution of information falling into a black hole has recently become the focus of increased research activity. The premise naturally tests the limits of quantum information theory and string theory, and necessitates a peek into black hole horizon mechanics. It is now believed that, in one way or another, information falling past the black hole horizon gets deeply and quickly entangled with the degrees of freedom of the black hole [1,2,3]. This process of information scrambling is not well understood. Complex physical systems scramble information over times scaling as a power law in the entropy. This phenomenon is generic, dictated by ergodic, random walk-like dissipation of information. Motivated by various paradoxes and thought experiments, the authors of [4,5] have proposed that a black hole is the fastest scrambling machine, diffusing information on timescales given by where S is entropy and T is temperature. Indeed, any faster scrambling would lead to a violation of the no-cloning principle of quantum mechanics. To date, there is no known physical system proven to achieve such a fast rate of scrambling. Contrived non-physical models have been developed in the literature [6,7], but none that are related to black holes or gravitational dynamics, let alone being realized within string theory. If the conjecture of [4] is correct, string theory -or any candidate theory of quantum gravity -should be able to demonstrate this fast scrambling phenomenon. In [5,8], it was suggested that Matrix theory, as a framework for M theory, may provide for a fast scrambling mechanism.
In a previous work [9], a controlled setting for addressing the fast scrambling conjecture was developed within Berenstein-Maldacena-Nastase (BMN) Matrix theory [10]. The latter is related to the Banks-Fischler-Shenker-Susskind (BFSS) Matrix model [11] and is also conjectured to describe dynamics of M-theory -in particular light-cone M-theory in a planewave background. The BMN setup comes with a infra-red cutoff (related to the background plane wave) which allows for vacua in the form of Bogomol'nyi-Prasad-Sommerfield (BPS) spherical membranes (M2 branes of light-cone M-theory), or fuzzy spheres (D2 branes of IIA theory) [12,13]. In the proper regime of the parameters, gravitational geometries have been identified that are dual to these Matrix vacua [14]. On the Matrix theory side, one can perturb around these vacua and effectively study the dynamics of these objects.
These perturbations of the fuzzy spheres come in two types: geometric deformations of the spherical membrane, appearing as bosonic degrees of freedom in the Matrix theory; and harmonics of the eleven dimensional supergravity multiplet -arising from fermionic degrees of freedom in the Matrix theory. In [9], it was shown that, for the vacuum of BMN theory corresponding to a single spherical M2 brane with N units of light-cone momentum, qubits correspond to the spherical harmonic modes of the eleven dimensional supergravity multiplet on a sphere. The system has both an infrared (IR) cutoff (the size of the sphere) and an ultraviolet (UV) cutoff (related to N , the regularization of the fuzzy sphere). The remaning 3 N 2 bosonic degrees of freedom describe the deformations of the shape of the M2 sphere in a three dimensional subset of the full space. (b) Our setup makes a heat bath out of the membrane ripples, in the background of which the qubits evolve. Integrating out the bosonic heat bath generates a dense network of interactions between the qubits. We look for entanglement evolution in this qubit system to determine whether BMN Matrix theory has the right ingredients to be a fast scrambler.
the perturbation dynamics can be mapped onto the dynamics of 8 N 2 qubits coupled to 3 N 2 bosonic variables that describe the shape of the M2 brane (see Figure 1). It was then demonstrated that the coupling between the qubits and the bosonic variables is necessary if the BMN theory is to achieve fast scrambling in the qubit sector.
In this work, we focus on an explicit computation in BMN Matrix theory with aim to ascertain whether the theory is indeed a fast scrambler. We start with the system of [9] and arrange the bosonic degrees of freedom in a heat bath at a fixed temperature T . This represents for us a background Matrix black hole model [15,16,17] through which we can trace the evolution of the qubits. We then compute the density matrix of the qubit system and measure entanglement as a function of time at one loop order. Since the qubits represent the graviton, 3-form, and gravitino modes of eleven dimensional supergravity, we are then tracking the evolution of information stored in the supergravity multiplet near the surface of a spherical rippling M2 brane. We are able to obtain a closed form expression for the entanglement entropy. We then extract the timescale of information scrambling from this result. Remarkably, we confirm that the scrambling time scales as the logarithm of the entropy of the qubits and is given by (1)! We believe this is the first known physical system exhibiting fast scrambling. It also is fully embedded in string theory and closely related to Matrix black holes.
In section 2, we develop the perturbation theory about the BPS membrane vacua of BMN theory. In Section 3, we present the core technology, the Feynman-Vernon technique of integrating out a heat bath; and we apply it to our system leading up to identifying the fast scrambling timescale. In Section 4, we present our conclusions and a discussion of future directions. The Appendices collect some of the more detailed results: Appendix A shows samples of combinatorial terms arising from the Feynman diagrams used in the computation; Appendix B shows a justification of the chosen coupling constant for a perturbative expansion; and Appendix C shows the main result -part of the computed entanglement entropy as a function time -used in identifying fast scrambling.

Setup
The BMN system is a 0 + 1 dimensional supersymmetric Matrix theory describing the dynamics of N D0 branes. The Lagrangian is given by [10,12] The X i s and Y a s are N × N Hermitian matrices with i = 1, 2, 3 and a = 4, . . . , 9. The Ψ is also a matrix in the adjoint of U (N ) with entries that are Majorana-Weyl spinors in ten dimensions. Two variables parameterize the system, R and µ. The limit where µ → 0 leads to the BFSS Matrix model, that is D0 branes in a flat background or M theory in the light cone frame with light cone momentum p + = N/R. For µ = 0, the theory describes D0 branes in a plane wave background, or light-cone M-theory in a plane wave background with again p + = N/R. We fix the static gauge A t = 0, with the covariant derivative D t becoming simply the time derivative -at the cost of the constraint where the Πs are the canonical momenta of the Xs. In our conventions, the length units of the various variables are given by The BPS vacua are described by the configuration The τ i s are the Pauli matrices, with [τ i , τ j ] = 2 i ε ijk τ k . Considering N dimensional representation of SU (2), these vacua are fuzzy (sometimes called non-commutative) spheres of radii 3 This setup is our starting point. In a proper regime of the paramaters and at finite temperature, it can serve as a model for a black hole within a fully quantized theory of gravity. We proceed by perturbing the vacuum as follows We freeze out the perturbations in the six transverse directions, Y a = 0, essentially focusing on a compactification to three target space dimensions. We also rescale time 3 In the language of [13], this corresponds to N 2 = 1 and N = N 5 1 for a single spherical membrane. The theory can also describe spherical M5 branes in other regimes of the parameters [13].
Henceforth, all variables are then dimensionless, and time or energy are measured in units of µ.
We then expand all matrix structures in matrix spherical harmonics Y j m with j = 1, . . . , N − 1 and m = −j, . . . , j [12,18]: note that we have 4 as expected. And the Y j m 's are N × N matrices satisfying the algebra [19] where written in terms of 3j and 6j symbols. We also have the normalization condition For the bosonic variables, the decomposition that diagonalizes the quadratic part of the Hamiltonian looks like [12] x The bosonic degrees of freedom are hence written in terms of physical complex variables α jm (t) and β jm (t). This decomposition diagonalizes the quadratic part of the Lagrangian. The indices j and m are summed over such that the range of the indices on the α jm s is j = 1, . . . , N −2 and m = −j, . . . , j, while the range on the indices of the β jm s is j = 1, . . . , N and m = −j, . . . , j.
Note that due to the conjugation property the α jm s and β jm s satisfy the conditions We can then easily check as required. Hence, the physical degrees of freedom are mapped onto all of the α jm s and β jm s. Similarly we find 8 N 2 fermionic degrees of freedom amongst the η jm s (4 N (N + 1) in total) and χ jm s (4 N (N − 1) in total) . Substituting these decompositions into the Lagrangian, we obtain where g is the effective coupling constant for large N g ≡ R µN 3/2 (21) and the notation is such that L BBB is cubic in the bosonic variables, L BF F is linear in the bosonic variables and quadratic in the fermionic ones, etc. L BBB and L BF F are linear in the structure constants f jm j m j m , while L BBBB is quadratic. In Appendix B, we demonstrate that these f jm j m j m s are N independent for large N . That is, all explicit N dependence in the Lagrangian have been factored into our definition of the effective coupling constant g.
For small g and large but finite N , we have a well-defined perturbative expansion in g. The radius of the membrane scales as g −2/3 as seen from (6). This assures that the membrane is large in Planck units. However, for g N ∼ 1, we can instead view the setup as a spherical D2 brane with a non-commutative worldvolume theory [13]. From this perspective, the theory has an IR cutoff set by µ. The scale of non-commutativity is given by And N tunes the ratio of the IR and UV cutoffs. The dimensionless coupling on the worldvolume is given by at the IR scale. The temperature in the system, 1/β, should be between the IR and UV cutoffs 1 where temperature is measured in units of µ. And the effective dimensionless coupling at the temperature scale of interest is then g 2 ef f (1/β) = g 2 N 2 β 2 . Thus, for we can capture the strong coupling regime g 2 ef f (1/β) 1 of this non-commutative D2 brane worldvolume dynamics through a perturbative expansion in g in BMN Matrix theory 5 . To probe this interesting regime, we can choose for example g ∼ 1/10, β 1, and 10 < N < 100.
The reason we are interested in the regime delineated by (25) is that strongly coupled non-commutative dynamics of a spherical D2 brane is rather reminiscent of the Matrix black hole proposals of [15,16,17]. Hence, through perturbation in g we may expect to capture the essence of black hole horizon physics -keeping in non-local dynamics in the picture yet scaling out superfluous complications from gravity, as is typical in non-commutative field theory settings.

Pure bosonic terms
The bosonic sector of the dynamics involves properly diagonalized quadratic terms [12] identifying two sets of masses in units of µ. It will be useful to write these terms in a more compact notation where we package the N (N − 2) components of the α jm s as with j = 1 . . . , N − 2 and K = 1, 2, 3; and we package the N (N + 2) components of the β jm s as where j = 1 . . . , N and K = 1, 2, 3. Equation (26) then takes the form where the dot represents sum over the m-tagged components, and with an implicit sum over K = 1, 2, 3. The remaining terms of the Lagrangian, L BBB and L BBBB , involve cubic and quartic interactions between these modes -scaled respectively by g and g 2 . In the limit of small g we will be focusing on, we will see that we do not need to consider these interactions for the purpose of demonstrating fast scrambling. Hence we do not show the explicit and rather complicated forms of L BBB and L BBBB .

Fermionic terms
The fermionic sector involves diagonalized quadratic terms [12] where the dot implies sum over the SU (4) indices of the spinors. Quantizing this diagonalized form leads to the Hamiltonian with the commutation relations Hence, we have a system of 8 N 2 disjoint qubits. However, there is also coupling to the bosonic sector through the terms where we have used the compact notation introduced in the previous section; and the ωs are defined by The coefficients shown as W (i)K jmj m ;j , with i = 1, 2, 3, 4, are rather complicated numerical expressions. Samples are shown in Appendix A. They satisfy the identities which assure that the Lagrangian is real. They consist of linear combinations of the structure constants f jm j m j m of equation (11), and hence involve the combinatorics of combining angular momentum modes through 3j and 6j symbols. Appendix B shows numerical profiles of these structure constants, demonstrating that the f jm j m j m s scale as N 0 for large N . All explicit N dependence of the action is hence factored into our definition of the coupling constant g through equation (21).

Perturbations
For weak coupling g 1, the fermionic sector consists of 8 N 2 qubit coupled to the bosonic sector through the L BF F term. The maximum entropy of the qubits scales as N 2 . We want to see whether the coupling between the fermionic and bosonic sectors can lead to the fast scrambling expected of a quantum theory of gravity as proposed in [4], with scrambling timescale τ ∼ ln N/T .
Our approach is to compute the density matrix for the qubits as a function of time -using the bosonic sector as a heat bath at a fixed temperature. We then construct the reduced density matrix by integrating over half of the qubits and obtain the entanglement entropy as a function of time. Analyzing this, we extract the characteristic scrambling timescale as a function of N . To manage the computation, we consider the limit of small coupling constant g 1 and large N 1. In this small g regime, our leading order computation amounts to considering the schematic Feynman diagrams shown in Figure 2. We would then be integrating out the bosonic dynamics to order g 2 leading to a four qubit Fermi interaction. These quartic interactions are bound to generate a dense network of connections between the qubits, and the question we are to address is whether this high density of connections is enough to lead to fast scrambling of qubit information. In this regard, note that entanglement can only be generated by a boson-fermion-fermion vertex since this is the only dynamical mechanism for flipping and entangling qubits. In particular, there is no direct scrambling at a work in the last two diagrams of Figure 2(a).

Density matrix 3.1 Setup
Our first goal is to compute the density matrix of the qubit system in the background of a bosonic thermal bath. We do this by using the Feynman-Vernon path integration technique [21]. In this language, the η jm s and χ jm s become Grassmanians. To simplify the Integrating out the bosonic heat bath leaves a Fermi four-qubit interaction that creates a large number of interaction links between the 8 N 2 qubits of the system. The scrambling mechanism arises entirely from diagrams I, II, and III, all involving one or more L BF F couplings. This high density of effective qubit interactions originates from the fact that many f s of (11) are non-zero -i.e. an angular momentum mode couples to almost all other angular momentum modes. Interestingly, the strongest couplings arise between the largest and smallest angular momentum modes.
notation, we represent them as F with I being the SU (4) index. The F s then satisfy The bosonic variables are similarly written as In subsequent path integral expressions, F and B represent all of the fermionic and bosonic degrees of freedom.

Integrating out the thermal bath
The novelty that arises in the Feynman-Vernon technique -as opposed to the traditional path integration approach to scattering amplitude computations -is that the in-out qubits (row-column entry labels) of the density matrix specify the boundary conditions in the path integration. In addition, we need to adapt the original Feynman-Vernon technique to a system involving fermions. As depicted in Figure 1(a), the computation is expected to be at the one-loop level to order g 2 . Our starting point is the Feynman-Vernon density matrix, properly adapted to Grassmanians and using a measure consistent with the Grassmanian coherent state formalism [22] ρ ρ(F i , F i , 0) is the density matrix at t = 0 which we take for simplicity to be the vacuum of the unperturbed qubit system Hence, we evolve ρ(F i , F i , 0) to ρ(F f , F f , t) using the Feynman-Vernon propagator which is given by 6 6 see [23] for the fully bosonic scenario.
Putting the bosonic degrees of freedom in a heat bath at temperature 1/β, the Feynman-Vernon action S F V takes the form where As can be seen from diagram V of Figure 2(a), the S BBBB terms do not arise in a context that flips qubits and hence cannot contribute to the boson-qubit coupling at leading order in g in (43). They can however arise in the computation of the heat bath terms Z R and H R as we shall see shortly. The S BBB terms arise in diagrams III and IV. But once again diagram IV is not involved in the boson-qubit coupling dynamics. This leaves the S BBB terms arising in diagram III only, while the remaining two diagrams, I and II, involve the S BF F terms. Due to the complexity of the computation, we will attempt to simplify things further: we drop diagram III with the following reasoning: considering this diagram in addition to diagram II can only add interactions between the qubits that may lead to faster but not slower scrambling; and if we can show that diagram II is enough to achieve fast scrambling without considering diagram III, our goal is achieved. On the other hand, if we were to find that the system does not fast scramble, we would need to come back and add diagram III and check once again. Phrased differently, we are trying to find an upper bound on the scrambling timescale; and if considering either diagram II or III demonstrates a timescale bound scaling logarithmically with entropy, we have proved that the system is a fast scrambler. This additional simplification then allows us to drop the S BBB terms as well from (43). We then have instead of equation (43). The thermal bath at temperature 1/β comes in through the expression [23] where S E is the usual Euclidean action And Z R is the partition function. For a collection of oscillators of the forṁ we have (48) We can think of the effect of a small g as shifting the oscillator frequencies λ k where λ (0) k are the frequencies determined at zeroth order from (26). The higher order corrections λ (i) k with i ≥ 2 involve contributions from diagrams IV and V of Figure 2(a)that is they come from the S BBB and S BBBB terms of the action. When computing the entanglement entropy in the fermionic sector, we will see that these higher order corrections to λ (0) k drop out since they correspond to disconnected diagrams. Incorporating all the bosonic and fermionic variables through this Feynman-Vernon process, we finally obtain the effective action with arising from the masses of the bosonic oscillators of equation (26); and we define The dot product signifies sum over the m components of the ωs. To leading order in g 2 , the λ i j s appearing in F i j can be replaced by j/3 and (j + 1)/3 from (51). The higher order corrections to the λ i j s do however contribute to the bosonic part B(t, N, β) in (50). However, the entire B(t, N, β) factor corresponds to disconnected diagrams and needs to be divided out when writing the density matrix. Hence, the relevant part of the Feynman-Vernon action becomes with λ 1 j = j/3 and λ 2 j = (j + 1)/3. Thus throughout we need not consider contributions from the L BBB and L BBBB terms of the original Lagrangian -to leading order in g. We then arrive at an effective action for the qubits, equations (52) and (53), that is quartic in the fermions and non-local in time, as expected.
To proceed further we note that S F V is proportional to g 2 , and hence we may write since g 1. This allows us to compute the Feynman-Vernon propagator using the standard technique of a generating functional where we define by introducing the appropriate Grassmanian source terms J and J . We can then proceed by computing G F V , and we get (see for example [23]) where we have defined and In these expressions, a is the mass of the corresponding fermionic mode from equation (30). Putting everything together, we arrive at the density matrix for the qubit system given by where the C (i) (t) are rather complicated numerical coefficient that we have explicitly computed. The total number of terms in ρ is around 20, 000. However, only C (1) and C (6) will end up contributing to the entanglement entropy. The first involves diagram I of Figure 2(a) and includes about 2, 000 combinatorial terms; the last involves diagram II and includes about 1, 000 terms.

Entanglement entropy
Once the density matrix is obtained, computing the Von Neumann entropy is rather straightforward. We first trace over about half of the qubits by integrating out the η variables to generate a reduced density matrix ρ = Tr η ρ .
This operation takes the integral form over Grassmanians This gives an expression of the form K (1) (t) arises from diagram I of Figure 2(a) or C (1) in equation (62); whereas K (2) j 1 m 1 j 2 m 2 (t) arises from diagram II or C (6) in equation (62). We can then obtain the Von Neumann entropy through S(t) = −Trρ lnρ .
However, there is a subtlety one needs to be careful about. Our computation at one loop renormalizes the qubit wavefunctions. Hence, we need to make sure that the proper counter terms are added so that Tr ρ = 1 .
We write where ρ is the 'bare' density matrix computed as a series expansion in g ρ → ρ = ρ 0 + g 2 ρ 2 + g 4 ρ 4 + · · · (69) and Z is the (wavefunction) renormalization factor and we can see that one needs This then leads to the following expression for the entropy in terms of the bare density matrix Note that the leading contribution of the entanglement entropy arises at order g 4 and there is no contribution to this from any part of the density matrix beyond order g 2 . Hence, to obtain the entanglement entropy at leading order in g, we must compute Tr ρ and Tr ρ 2 -both of which we have done. However, for the purposes of identifying the scrambling timescale, we have found that looking at the (Tr ρ 2 ) 2 term is enough to demonstrate fast scrambling. Hence, we do not write the very lengthy full result and instead only quote the result for Tr ρ in Appendix C.

Analysis
In this section, we collect the results of analyzing the evolution of the Von Neumann entropy. We have computed the entire expression given by (72) in closed form. The complexity of this computation arises from the combinatorics of the diagrams, i.e. the many sums over angular momentum modes. The first term of (72), (Tr ρ 2 ) 2 , involves for example eight such sums and the number of terms scales as N 8 . The second term, Tr ρ 2 2 , contains sixteen sums or ∼ N 16 terms. We have tried to use various asymptotic methods to tame these expressions without much success. We are hence unable to present analytical results. But since we have a closed form expression for S(t) -albeit a long and complicated one -we can evaluate it numerically. As mentioned earlier, we need to consider large N . This makes even the numerical evaluation challenging. For the Tr ρ 2 2 term, we have not been able to numerically evaluate the result for N > 4 in a reasonable amount of time. However, we are able to tackle the (Tr ρ 2 ) 2 term for N up to around 40 using parallelized algorithms. As we shall see, looking at only the (Tr ρ 2 ) 2 term for such large values of N is enough to demonstrate that the scrambling time indeed scales logarithmically with entropy. Henceforth, we focus then on a numerical analysis of the first term of (72), the (Tr ρ 2 ) 2 term, that arises from diagram I of Figure 2(a). Figure 3 depicts the time evolution of the Von Neumann entropy due to the first term in (72). On the left is a plot of the logarithm of the entropy as a function of time. Noting that the initial state was the vacuum of the unperturbed theory, and that the scheme of obtaining the reduced density matrix involves a rather symmetric partitioning of the full system, we can see that the evolution of the entanglement entropy is not chaotic. We also may be seeing the phenomenon of Poincaré recurrence since the evolution demonstrates a characteristic cyclic period. The right figure shows a plot of entropy as a function of time zoomed onto a range of timescales for which the entropy evolves exponentially with time. Typically, the process of entropy evolution is characterized by several timescales. The timescales associated with the exponential growth can be identified as the scrambling timescale [24,25,9], which we now focus on. We find it using τ ∼ Ṡ S .
(73) Figure 4 shows a plot of scrambling time as a function of N . Removing the first few small N values, we fit the data to a logarithmic profile and find excellent agreement. To contrast this with alternative entropy evolution models, we show two more fits to the data in Figure 5. Both are power law fits τ ∼ N c , with the plot on the left having no restrictions on the fit parameters while the plot on the right restricting c to be positive. While the model on the left gives as good of a fit as the logarithmic model of Figure 4, it involves a pathological conclusion: τ ∼ constant − 1/ √ N . Intuitively, we must have scrambling time scaling at worst as a positive power of N . The plot on the right shows a fit with τ ∼ N 0.2but, as is apparent, the fit is a very poor one. We are then lead to conclude that the BMN matrix model, in the regime of the parameters we have explored, scrambles information on timescales τ ∼ ln N ∼ ln S .
That is, the BMN matrix model is indeed a fast scrambler, saturating the quantum nocloning bound. Stated more carefully, our analysis truncated away diagram III of Figure 2.
In this sense, it can only be used to find an upper bound on scrambling timescale. Hence, we are able to conclude that τ ln S. The additional diagram that we have ignored can only make the scrambling faster. Yet, we do not expect to violate the quantum no-cloning bound since the framework of computation is quantum mechanical. Instead, we must conclude that this system saturates the scrambling bound as in (74).
Finally, we look at the temperature dependence of the scrambling time. Figure 6 shows the scrambling time as a function of inverse temperature β. As discussed earlier, the system has a IR cutoff of order one in units of µ. Hence, at low temperatures β 1, we expect finite size effects to kick in and a freeze out of the heat bath temperature. Correspondingly, we see from the Figure that the scrambling timescale becomes independent of temperature -that is, there is no dynamics in larger thermal wavelengths than the size of the box in which the system is placed. For low enough temperatures β 1 above the IR cutoff scale, we should however expect τ ∼ β. The inset in the figure shows a zoomed view of this temperature regime, demonstrating a linear dependence on β as expected.

Conclusions and Outlook
Given the connections fast scrambling has to black hole dynamics [7,4,5,26], our result is further circumstantial evidence that BMN Matrix theory is a theory of quantum gravity. However, the mechanism of fast scrambling we identified is also present in the BFSS Lagrangian: a dense network of connections between qubits arising from the commutator of matrices, essentially from the density of the structure constants of equation (11). Hence, the more general implication of this work is a confirmation of the suggestion by [5] that M-theory-related Matrix theories are indeed fast scramblers. Furthermore, to achieve this fast scrambling, the problem relies on the combinatorial complexities involved in coupling representations of SU (2) -that is, a dense web of interactions between all spherical harmonic modes -which in turn makes an asymptotic analytical treatment most challenging. It is interesting to note that the structure constants of equation (11) which determine the density of the effective qubit-qubit network are peaked for couplings that mix the highest and lowest angular momentum modes -as if a stringy UV-IR mixing is central to the highly efficient scrambling dynamics.
The model we employed is very reminiscent of the Matrix black hole of [15,16,17]. In a sense, we have made this original Matrix black hole narrative more concrete: a heat bath associated with the shape of a spherical membrane -a structure akin to the membrane paradigm of the black hole horizon -to which one tethers qubits that represent the supergravity degrees of freedom. The equilibrium state would correspond to the two sectors of the model -the membrane shape degrees of freedom and the qubits -being in thermal equilibrium. In this work, we artificially arranged the qubits in their ground state and tracked the subsequent evolution towards thermalization. This then naturally allowed us to identify the scrambling timescale. The energy regime we focused on corresponds to non-local fuzzy D2 brane dynamics, strongly coupled yet without gravitational dressing. The implication is that gravitational physics is a red herring in the information scrambling problem: the relevant ingredient is non-local interactions near the black hole horizon 7 . This setup is only one step of exploration in an otherwise very rich system that admits many more interesting black hole-related applications.
For future explorations, there are many interesting venues to pursue. One can complete the full calculation by including all Feynman diagrams to construct the full form of the entanglement entropy as a function of time. While we do not expect that this will add to the punchline of fast scrambling, analyzing the full form of the evolution of the entropy may help identify other interesting timescales that are typically at work in thermalization processes [25]. And while we were forced to fall back on a numerical analysis of the final entropy expression, there is the potential for a successful analytical analysis: the asymptotics of the structure constants for N → ∞ can be worked out using Edmonds' relations and their extensions [27]. The analogue of the thermodynamic limit in our system is indeed the regime N → ∞; and hence techniques used in traditional statistical systems can perhaps be used to tame the complicated entropy expression.
Black hole evaporation should correspond to qubits freeing themselves from the thermal bath, the membrane at the horizon. The entanglement between Hawking radiation and the black hole then translates to the entanglement between evaporated qubits and membraneattached qubits. Toy models of qubit-based evaporation mechanisms for black holes have already been explored in the literature [29,31,30], albeit with too much freedom in choosing arbitrary mechanisms of qubit-qubit interactions. However, this Matrix theory model provides for a concrete string theory-embedded mechanism for determining the details of Matrix black hole evaporation. Along the same line of thought, one can attempt to address the in-fall problem by initially arranging for a small subspace of the matrices away from the spherical configuration and then tracking the merging of the probe into the membrane structure. Qubits can then once again be used as probes of the evolution of information into the Matrix black hole. Perhaps all this can lead to a better understanding of a horizon firewall [32,33,34]. BMN Matrix theory also admits other interesting regimes in its parameter space. In particular, spherical five brane vacua can be realized in the theory, along with a similar perturbation problem yet to be explored.
Finally, vacua of BMN Matrix theory have dual holographic spacetimes that are valid geometries in complementary regimes to the one we have focused on [14,13]. Our setup corresponds instead to a setting without gravity but with non-local dynamics -the realm of non-commutative field theories. It would be interesting to develop a more detailed map between qubits or perturbations in the Matrix theory language and corresponding supergravity excitation modes. This would help to directly connect the scrambling dynamics we analyzed to the supergravity black hole scrambling phenomenon.

Appendix A: W coefficients
Here we give a couple of samples for the type of expressions that arise in computing the boson-qubit coupling terms of the action, packaged into what we label as the W In writing W (i)K j 1 m 1 j 2 m 2 j 3 , we have i = 1, . . . , 4 and K = 1, · · · , 6. For K = 3 and K = 6, the W (i)K j 1 m 1 j 2 m 2 j 3 is a scalar, but for other Ks it is a vector in the m 3 index. As shown in Appendix B, the f s are N independent. The interactions packaged in the W s are essentially a packaging of the combinatorics of combining SU (2) representations.

Appendix B: Structure constants
In this appendix, we show that the structure constants f j 1 m 1 j 2 m 2 for all values of its indices. The horizontal axis is scaled from 0 to 1 to map onto the various possibilities of (j 1 , m 1 , j 2 , m 2 , j 3 , m 3 ) -sorted according to the value of the corresponding f j 1 m 1 j 2 m 2 j 3 m 3 . The different curves correspond to different values of N , from N = 5 to N = 21. We see a clear convergence pattern, with the maximum and minimum values of f capped off about ±18, independent of N . Indeed, the asymptotic form of this profile for large N can be obtained in closed form through Edmonds formula [27].
The peaks of f arise when two of the js are maximized and the third is minimized. For example, for N = 11, maxima arise around (j 1 , m 1 , j 2 , m 2 , j 3 , m 3 ) = (1, 0, 10, −10, 10, 10) and permutations of these (j 1 , m 1 ), (j 2 , m 2 ), and (j 3 , m 3 ). This implies that the predominant contribution in our computation comes from dynamics that mixes the highest and lowest angular momentum modes. The N independence of the f s inspires the definition of the coupling constant given in (21) for large N .

Appendix C: Entropy evolution
In this appendix, we write the full expression for Tr ρ , the square of which gives us part of the entropy as a function of time.