Quantum Black Hole Formation in the BFSS Matrix Model

We study the various head-on collisions of two bunches of D0-branes and their real-time evolution in the BFSS matrix model in classical limit. For a various matrix size N respecting the 't Hooft scaling, we find quantitative evidence for the formation of a single bound state of D0-branes at late time, which is matrix model thermalization and dual to the formation of a larger black hole.


Introduction
Black holes have always been at the centre of theoretical interests in quantum gravity and string theory. Through the non-perturbative formulation of string theory and quantum gravity via gauge/string correspondence [1], significant progresses have been made in our understanding of black hole microstates, especially its entropy counting. However, we still need to understand better various aspects of quantum gravitational system, and one of the important aspects is understanding their real-time evolution.
One model, which is more tractable to analyse the real-time dynamics compared with many other theoretical models in the gauge/string correspondence, is a BFSS (Banks-Fischler-Shenker-Susskind) matrix model [2]. This model is a 0 + 1 dimensional large N matrix quantum mechanics, and conjectured [3,4] to be dual to either ten-dimensional string or eleven dimensional M theory depending on the various parameter limit. Since BFSS matrix model is dual to the higher dimensional spacetime, it would be very interesting to ask how the higher dimensional black hole formation can be understood from the dual matrix model. In this set-up, we consider the various head-on collisions of two bunches of D0-branes and their real-time evolution, which is expected to become a single bound state of D0-branes at late time. This is a matrix model thermalization and dual to the formation of one large black hole in either ten-or eleven dimension. Our main interest in this paper

JHEP07(2015)029
is, by solving the matrix model classical equation of motion, to study quantitatively how such a bound state is formed in the real-time evolution of the large N matrix model.
In order to study the complicated time-evolution of large N matrix model, we make several simplifications; First, we consider only the classical limit in this paper, namely we solve the classical equation of motion numerically and find its time-dependent solution. At first thought, this simplification might sound very boring limit. However, there are several reasons why classical limit of late time evolution is interesting and capture aspects of non-perturbative physics. Even though classical limit is justified in the weak coupling limit, since what we conduct is not the perturbative analysis near the trivial vacuum but rather the analysis seeking the time-dependent soliton configuration, this can capture non-perturbative physics. Remember that thermalization is a very complicated but also universal phenomenon which one cannot see in various solvable models, nor in perturbation analysis from the trivial vacuum [5]. It has a property that all memories of the initial data are lost in late time, and this late time thermalization may occur only in the large N limit. Therefore even though we simply solve the classical equations of motion, we are looking some non-perturbative physics through the late time evolution.
Another simplification is, that we neglect effects of fermions for simplicity. Since supersymmetry is expected to be more crucial at low energy but not in high energy deconfined phase, we expect that, for late time thermalization, the analysis with/without supersymmetry does not change much. In our analysis, we have conduct N up to N = 24, and regard it as large enough to perform an extrapolation to the large N limit.
Before we close the introduction, we comment on several related references. Thermal equilibrium (i.e., static) properties of the BFSS matrix models for black holes are very well studied. Thermodynamic simulations of BFSS matrix model was initiated by Kabat et al. [6,7] by the mean-field method [8], and recently by Monte Carlo method [9][10][11]. These results confirmed that the duality is valid not just at supergravity level but also at stringy level; at finite-coupling [12] and finite-N [13]. For non-equilibrium properties, a probe limit thermalization was studied in [14][15][16]. For so-called BMN matrix model, 1 Berenstein et al. studied the thermalization numerically in the series of papers [18][19][20] in detail. BFSS matrix model with initial conditions different from ours have also been studied in [19,20]. Finally the black hole formulation at the correspondence point for the BFSS matrix model was analytically studied recently in [21,22]. In this paper we consider the similar setting to these works, though the parameter range is different, which is possible since we solve the model numerically.
2 The model and simulation method BFSS matrix model [2] is the maximally supersymmetric matrix quantum mechanics, which is the dimensional reduction of 4d N = 4 supersymmetric Yang-Mills theory to one dimension (time). In this paper, we concentrate on the classical dynamics, and set fermion to zero for simplicity. Henceforth, we consider only the bosonic part of the Lagrangian of JHEP07(2015)029 BFSS matrix model and it is given by Here X M (M = 1, 2, · · · , 9) are N × N Hermitian matrices and (D t X M ) is the covariant derivative given by where A t is the U(N ) gauge field. This matrix model is a gravity-decoupled theory on the N D0-branes [23]. Large N bound state of D0-branes with finite temperature T with large 't Hooft coupling is dual to a black hole (black 0-brane) with Hawking temperature T [3]. Such a bound state is described by highly non-commutative matrices. On the other hand, if matrices are close to block-diagonal, it describes multi-D0-branes and each of which is highly quantum 2 black hole.
In the A t = 0 gauge, the classical dynamics is described by the equation of motion with the Gauss's law constraint In the rest of this paper we solve this classical equation of motion. From the action (2.1) it is clear that 't Hooft coupling λ ≡ g 2 Y M N appears only as an overall factor. Therefore the classical limit, i.e., → 0 limit, is equivalent to the effective coupling constant λ eff ≡ λ/U 3 goes to zero limit in the 't Hooft scaling limit, where N → ∞ and U is the typical VEV scale in the matrix model associated with the temperature, or quantum fluctuations of the matrix degrees of freedom. Therefore, our analysis of solving classical equation of motion is justified in hight temperature/energy limit. 3

Discretization
In order to study the time evolution, we discretize the equation of motion (2.2) while preserving the constraint (2.3) exactly. For that purpose, we write the EOM as and (2.5)

JHEP07(2015)029
The constraint becomes Under the infinitesimal time evolution t → t + δt, X M and V M change as as one can easily check by using the standard Taylor expansion. We terminate this expansion at the order of (δt) 2 , (2.10) Then, a simple but tedious calculation shows that the constraint (2.6) is satisfied at t + δt if it holds at t, without any additional higher-order terms: Therefore, we use (2.9) and (2.10), and choose the initial configuration to satisfy (2.6) at t = 0,

Initial condition -collision of two bunches of D0-branes
We consider a collision of two bunches of D0-branes, which is analogous to two black zero-branes, or two black holes (the left of figure 1). The initial condition is 4 and JHEP07(2015)029 Figure 1: If a single bound state of eigenvalues is formed, it is analogous to a formation of one large black hole, and this is the dual to the thermalization process of the whole matrix system.
for µ = 2, 3, · · · , 9, i = 1, 2, · · · , N/2 and j = N/2 + 1, · · · , N , where a µ,ij and b µ,ij are Gaussian random numbers with the weights e −a 2 /2 and e −b 2 /2 . Here σ is the source for off-diagonal elements and therefore non-commutativity, which is crucial ingredients for the bound state formation. If σ = 0, then commutator square potential vanishes and two bunches just pass by. For nonzero σ, the quantitative details of the initial condition can depend on the choice of random numbers, especially when N is small. Note that we consider only even values of N . This initial condition satisfies the Gauss-law constraint We consider the 't Hooft large-N limit, where λ ≡ g 2 Y M N is fixed, L and V are fixed, and then, the energy scales as N 2 . When two bunches are separated, the off-diagonal element σ has larger mass as σ and so quantum mechanically path-integrated out perturbatively in λ/σ 3 and suppressed. This describes the quantum fluctuation of open strings stretched between two bunches. Because the potential energy at t = 0 is We will see that the energy actually scales as E ∼ N 2 then.

Formation of a single bound state
Let us first ask 'when a single bound state of eigenvalues is formed'. This is the process where two bunches of D0-branes, or equivalently two quantum black holes merge to form one large quantum black hole (figure 1).
The easiest way to see this process is to plot TrX 2 M , especially TrX 2 1 , as a function of t. In figure 2, TrX 2 1 /N for N = 8, for L = 5.0, V = 0 and √ N σ = 0.12, with several different values of time step dt, is shown. Here we set the 't Hooft scale λ = 1. At late time, errors associated with discrete time steps become non-negligible. We can see that behaviours at t ≤ 70 can reliably be studied with dt = 0.00010. TrX 2 1 /N bounces a few times and then stays small, which suggests the formation of a single bound state. Note that, if instead two bunches oscillated as shown in figure 3, TrX 2 1 /N would oscillate without decaying. At late time, we expect that the system goes to the the "typical" configurations and that such typical states are rotationally symmetric after taking the average over time and/or different initial configurations. In figure 4, the average T rX 2 M /N (M = 1, 2, · · · , 9), by using 0 ≤ t ≤ 10, 10 ≤ t ≤ 20, · · · , 60 ≤ t ≤ 70, for 100 different initial configurations, JHEP07(2015)029 are plotted. We can see that the average values converge to the same value at late time, suggesting the restoration of the rotational symmetry.

Insensitivity to initial conditions for thermalization
In order to add further evidence for thermalization, we studied two completely different initial conditions with the same value of the energy. In order to tune the total energy, we JHEP07(2015)029 set V M = 0 and rescale X M to αX M . Then the energy scales as E → α 4 E. By using this, we can tune the energy to any value. Here we consider (i) N = 4, L = 10.0, V = 0, √ N σ = 0.12, and (ii) N = 4, all matrix components are generated Gaussian weight. Then we rescale X M so that E/N 2 = 0.5. Note that (i) is much more anisotropic initial condition than (ii). We employed dt = 0.00001 and dt = 0.000001 in this calculation, which give consistent results at t < 90. From figure 5, we can see that 9 M =1 TrX 2 M /9N with two different initial conditions behave the same manner in late time when the energies are the same, and it is hard to tell the initial condition unless we explicitly solve the equation of motion to go back to the past. These suggest typical thermalization, i.e., whatever initial conditions it starts with, a system ends up with similar typical states.

Large-N limit and thermalization
In order to study the statistical nature of the BFSS matrix model, it is not essential to take the step size dt very small, as long as the energy is conserved and we take enough late time evolution. In the following, we take dt = 0.0005.
Let us first consider a large-N limit with fixed L, V and √ N σ. As a concrete example, we consider L = 5.0, V = 0 and √ N σ = 0.12. Similar results were obtained for other values too.
In figure 6, we plot TrX 2 1 /N for N = 8, 12, 16. We can see qualitatively similar behaviours, and the fluctuation at late-time becomes smaller as N becomes larger. When a bound state is formed, the Virial theorem relates the kinetic energy K and the potential energy V as K = 2 V = 2 3 E, where · stands for the time average and E = K + V is the total energy, which is conserved. Therefore, (K − 2 3 E)/N 2 is a good measure for the fluctuation. As we can see from figure 7, this quantity fluctuates around zero and suppressed more as N becomes larger.
In order to see the statistical property at late-time, we take the time average at 50 ≤ t ≤ 100 and then take an average over random initial configurations (50 samples for N = 4, 8,

Large-N limit with fixed E/N 2
Let us consider another, perhaps more natural large-N limit, in which E/N 2 is fixed exactly. 5 We first generate initial configurations with L = 5.0, V = 0 and √ N σ = 0.12, and then set E/N 2 to 1.5, by rescaling scalars as we did in the end of section 3.  samples for N = 24. We use 50 ≤ t ≤ 100 again. In principle, since the energy is fixed, the time average for all samples with different initial configurations should give the same average value within errors for a sufficiently long time interval. (Still, we terminate the simulation at t = 100 so that the conservation of the energy is not violated due to the discretization error.) Therefore, statistical fluctuations of the average 9 M =1 TrX 2 M /9N become much smaller compared to figure 8. The results are plotted in figure 9. We can see that the correction to 9 M =1 TrX 2 M /9N also starts with 1/N 2 .

Large-N limit and thermalization
As we have seen, the system shows typical thermalization process and the fluctuation of macroscopic quantities disappears at large N . Therefore, the macroscopic quantities converge to the same value at late time, irrespectively of the initial condition, as long as the energy is the same. The information of the initial condition is hidden in the 1/N correction. This is exactly what we expect. Since the large N limit corresponds to the thermodynamic limit there is no distinguish between micro canonical ensemble (this is what we did) and canonical ensemble where temperature is fixed. The difference only appears in 1/N expansion. Note that an exact thermalization occurs only in the large N limit.   1 /N has a well-defined large-N limit. Remember that large N limit is a thermodynamic limit, where canonical ensemble and microcanonical ensemble are distinguishable only by looking at the 1/N suppressed order. Therefore, we call the system get thermalized at the late-JHEP07(2015)029 time if the system becomes indistinguishable "typical states", and the deviations for the macroscopic quantities from the typical states are only of order 1/N . Then, to read-off the thermalization time, i.e., the time-scale system get thermalized, we have to read-off the decay of the deviation at late-time in great detail to the order of 1/N . However this is hard generically since the error bar is too large and we have initial-condition-dependence for the late-time average.

JHEP07(2015)029
To avoid this issue, we define the thermalization time by another way. This is defined by the time scale for a small perturbation added on top of the bound state to become invisible. Note that if N is not sufficiently large, then the thermalization time defined in this way would depend on the detail of the perturbation. Let us consider two-point function along time direction which does not require the perturbation, and from that, we will read off the thermalization time. Let us consider a gauge-invariant operator where W s,s+t is the Wilson line,

JHEP07(2015)029
In our setup we took the A = 0 gauge, so Let us consider where the time average · is taken at late-time. We have already seen that f (0) is of order one. When t becomes large, f (t) decays because X M (s + t) "forgets" the information of X M (s). Therefore, it should be possible to determine the thermalization time from the time scale for the decay of f (t). The numerical result is shown in figure 11. It turned out that f (t) can be fitted well as where the overall constant c = M 1 N TrX M (t) 2 , the decay width Γ and the frequency ω are of order N 0 . The width Γ naturally gives the time scale for the thermalization, We can also consider connected two-point functions of gauge invariant operators. We consider one of the operators studied by Berenstein and collaborators, Tr(X M (t)X N (t)) (M = N ). Although we did not find a simple fitting function for this, the decay is consistent with e −Γ t (figure 12). Then the thermalization time defined by this operator is (Thermalization time) ∼ 1/Γ ∼ N 0 .
(3.8) Figure 12 suggests that the thermalization time, defined from the exponential decay, converges to finite value in the large N limit. This is quite interesting since it is consistent with finiteness of quasinormal modes in the black hole background, although this convergence is more subtle at late time due to finite N effect. As is discussed in [24], although the correlation function keep decaying to zero at any finite order in the 1/N -expansion, the decay should disappear once full finite-N corrections are taken into account. The fat tail at long time scale in figure 12 would demonstrate this property.

Under what initial conditions, does thermalization occur?
Let us consider at what values of V and σ the thermalization, or equivalently the formation of the bound state, can be realized. For simplicity we fix L to be 5.0. Firstly let us note that the off-diagonal elements are the source for the non-linearity, which is responsible for the formation of the bound state. If we set the off-diagonal elements to be zero (σ = 0) at t = 0, then the off-diagonal elements are not generated by the classical equations of motion and no interaction appears at all. In such case, two bunches just pass through each other. Even if σ is nonzero, if it is too small, the non-linearity cannot grow large enough before the collision. Therefore, it is clear that σ must be sufficiently large for the bound state formation. However, as we will see shortly, it turns out that σ should not be too small but also must not be too large in order to form the bound state. This indicates that if the σ accelerates the two bunches before the first collision too much, they tend to pass by without forming the bound state.

JHEP07(2015)029
In order to illustrate these, we show the evolution of TrX 2 1 /N for N = 16, for L = 5.0, V = 1.0 and various values of √ N σ in figure 13. For √ N σ = 0.16, 0.12 and 0.08, the latetime behaviours at t 20 look qualitatively the same. A fact that a bounce at √ N σ = 0.12 is larger than those at √ N σ = 0.08 and √ N σ = 0.16 suggests very complicated dynamics. At √ N σ = 0.04, on the other hand, two bunches just pass through with each other. It might be possible that the bunches comes back later after very long time and form a bound state; in fact sometimes we observe that TrX 2 1 /N becomes as large as O(100) and then come back and form a single bound state. As we can see from the early-time behaviour in figure 13, as σ increases the bunches are accelerated more before the collision (i.e. TrX 2 1 /N decreases more rapidly). So whether a single bound state corresponding to thermalised state is formed or not depends on a very subtle competition of the acceleration before the collision and deceleration after the collision, both of which are enhanced by a larger σ. The green "×" symbols in figure 14 show (V, √ N σ) in which two bunches passed through with each other and TrX 2 1 /N became at least twice bigger than the initial value. A single-bunch state would be formed after long time, but not immediately. (See e.g. √ N σ = 0.22 in figure 13.) The red "+" symbols show (V, √ N σ) where a single bound state is formed without bouncing too much. We can see that a single bound state is not formed when √ N σ is too large, which means the acceleration before the collision beats the deceleration after the collision. We cannot get any conclusive pattern from these, probably the result depends on the choice of the random numbers. We left this issue as open questions.

Summary and discussions
In this paper, we considered the head-on collisions of two bunches of D0-branes and study its real-time evolution in the BFSS matrix model. Especially our interests in this study are how time evolution differs and under what conditions the big bound states corresponding thermalized state are formed and what is the its time-scale.
We have seen that all of TrX 2 M /N where M = 1, · · · , 9 converges at late time to the same values and rotational symmetry is restored irrespective to the anisotropic initial conditions in section 3.1. Especially the insensitivity to the initial conditions strongly JHEP07(2015)029 suggest the thermalization and formation of a large bound state dual to quantum black holes in late time. In section 3.2, for various scaling limit, we have checked that at the thermalized stage, Virial theorem, which is expected to hold at the equilibrium state, is in fact satisfied and deviation from that occurs only at the subleasing order in 1/N expansion. In section 3.3, we estimated thermalization time from the time-direction two point functions and its exponential decay at late time, and shows that the timescale is N 0 order, which seems to converge in the large N limit. This result is very interesting since it is consistent with quasi-normal mode in the dual black hole. Finally in section 3.4, we studied the initial condition dependence for the thermalization. Clearly for small σ, the off-diagonal sources for non-commutativity, two D0-brane bunches pass by while we notice that for too large σ, two bunches also pass by in such cases. We however have no clear interpretation for this results. Our analysis is conducted under several assumptions. The big assumptions we take in our analysis is that we can completely neglect the quantum effects and fermion effects. Quantum effects are important as the effective coupling becomes larger, and fermions are also very important since they change the interactions between two bunches of D0-branes. It is quite possible that once we take into account these effects, then the time-scale we estimated at section 3.3 can be modified significantly. Obviously more detail study is necessary to get conclusive results.
To better understand under what initial conditions does the system thermalized, let us discuss what is the typical states in the BFSS matrix model. If the final states is made JHEP07(2015)029 √ N σ) in which two bunches passed through with each other and TrX 2 1 /N became at least twice bigger than the initial value. A single-bunch state would be formed after long time, but not immediately. The red "+" symbols show (V, √ N σ) where a single bound state is formed without bouncing too much. We used the same set of random numbers a ij µ and b ij µ for all V and √ N σ.
up by sets of almost commuting matrices, then the entropy of the gas is made up by N D0-branes, which has only O(N ) entropy. Off-diagonal elements are highly suppressed in such case and does not contribute to the entropy. On the other hand, a single large black hole, which is a single bound state of the D0-branes, has the entropy of order N 2 , because all of the off-diagonal elements are excited. We can also consider multi-black hole states. Clearly a single black hole is entropically and therefore free-energy-wise dominant, since all of the off-diagonal elements are equally excited and contributes to the entropy.
However there are subtle issues. In the BFSS matrix model, because moduli space is not bounded, when the two bunches pass through each other, it is possible, at least classically, that they just pass by and run to infinity. This is all due to the fact that all of the adjoint scalars have zero mass. 6 However by taking into account quantum fluctuations, 6 There is a subtlety if the moduli space is infinity, how the ergodicity can hold. Probably the appropriate limit is that, first we introduce a proper IR cutoff for matrix eigenvalues (for examples, by putting D0branes in a box or by introducing small but nonzero mass for the adjoint matrices) so that the phase space becomes bounded, the ergodicity can hold for generic initial conditions. And then we evaluate spectrum and correlators and in the end, we remove the IR cut-off. For the BMN matrix model [17], this IR problem is avoided.

JHEP07(2015)029
this is unlikely, because in the real system involving quantum fluctuation, there are always small but nonzero quantum fluctuation and when the non-linearity from the interaction is large enough at late time, two bunches are likely to be pulled each other and merge to a single bound state. If the entropy of the black hole is large enough (i.e. N is large enough), then it is unlikely to see a large deviation from the black hole within a finite simulation time.
Note that since our calculation is classical, it is crucial to take small effective coupling limit i.e., λ eff → 0 to justify our results. Even though effecting coupling is small, the equation of motion can accumulate large non-linearity at late time, which can capture non-perturbative soliton formation physics. A large N effect (which plays the same role as thermodynamic limit) and the late time universality effects (which capture enough non-linearity of the interaction and scramble) are crucial ingredients for the thermalization. There is an argument that the classical time evolution of the BFSS matrix model is stochastic [25]. Such natures are probably crucial for the system get thermalized and show the universal behaviour at late time.