Matrix Thermalization

Matrix quantum mechanics offers an attractive environment for discussing gravitational holography, in which both sides of the holographic duality are well-defined. Similarly to higher-dimensional implementations of holography, collapsing shell solutions in the gravitational bulk correspond in this setting to thermalization processes in the dual quantum mechanical theory. We construct an explicit, fully nonlinear supergravity solution describing a generic collapsing dilaton shell, specify the holographic renormalization prescriptions necessary for computing the relevant boundary observables, and apply them to evaluating thermalizing two-point correlation functions in the dual matrix theory.


Introduction
In recent years, gauge/gravity duality, also known as "holography," has emerged as a rare tool for the study of strongly coupled systems far from equilibrium. Originally motivated by the creation of a quark gluon plasma in ultrarelativistic heavy ion collisions, many authors have used the AdS/CFT correspondence to study what happens when energy is suddenly injected in a strongly coupled quantum field theory. Interesting results include thermalization times short enough to be compatible with experiments [1][2][3][4][5][6][7][8][9][10], a thermalization pattern in which short-wavelength modes thermalize first [7,11], and new insights in the spreading of entanglement entropy after a 1+1d quantum quench [5,7,[12][13][14].
It is interesting to ask whether holography can be used to make predictions for the thermalization of systems that can also be studied using conventional techniques. If so, this would provide a framework in which holography can be quantitatively tested in a far-fromequilibrium regime. With this question in mind, we will study holographic thermalization in BFSS matrix theory [15,16], a quantum-mechanical model of N × N matrices which, in the N → ∞ limit, has been proposed as a nonperturbative definition of M-theory in asymptotically flat backgrounds [16]. Our considerations will revolve around the relation of this model [17] to a non-conformal version of the AdS/CFT correspondence [18][19][20][21]. It has also appeared in recent discussions of "fast scrambling" [22], the fast spreading of information that is added "locally" (e.g., in one matrix component). While our focus will be on far-from-equilibrium processes driven by energy injection, a simpler holographic setup involving the same matrix theory has been previously studied, in a way involving extensive numerical simulations, in a sequence of papers including [23][24][25][26][27][28]. In those considerations, a stationary black hole was introduced in the gravitational bulk, corresponding to thermodynamic equilibrium, rather than thermalization, in the dual matrix theory. Further analytic considerations of matrix theory thermodynamics can be found in [29,30]. In [31,32], dynamics of moduli fields has been explored as a tool to probe thermal properties of higherdimensional super-Yang-Mills theories, though applications to the case of matrix theory are less straightforward.
Despite the apparent simplicity of the BFSS matrix theory, which involves quantum mechanics rather than quantum field theory, real-time evolution of this model in the appropriate strong coupling regime presently appears to be out of reach of conventional techniques, even for small values of N . The nine scalar matrices and their fermionic partners contain too many degrees of freedom to allow direct diagonalization, and the interactions between the various matrix elements appear to be too nonlocal for variational techniques such as tensor network methods to be directly applicable. It would be really nice if these or other techniques could be developed up to the point where they can capture matrix theory, first for small N and later for larger values of N , in order to allow detailed comparison with holography. (We would like to mention an intriguing attempt to tackle the quantum dynamics of a simpler bosonic matrix theory undertaken in [33].) In the meantime, numerical simulations have been carried out in another regime, where the matrix theory can be treated classically [34][35][36]. (A similar study of the related BMN matrix model can be found in [37].) This is a simplification which would not arise in higher dimensions; see [36] and references therein. One motivation is that, according to numerical simulations, there is no phase transition between the different regimes [23][24][25][26][27][28], so some qualitative features can be expected to be similar [36]. Further studies of semiclassical processes in the matrix theory revolving around the idea of continuity from weak to strong coupling can be found in [38,39].
In this paper, we use D0-brane holography [18][19][20][21] to study far-from-equilibrium evolution of matrix theory after a sudden injection of energy. In higher-dimensional AdS/CFT, a simple way to inject energy in a holographic field theory is by briefly turning on and off a homogeneous source, for instance for an anisotropic component of the stress tensor [3], for an electric current [40] or for a scalar operator [41]. In the bulk, this corresponds to turning on nontrivial boundary conditions for the corresponding bulk field. For a small-amplitude scalar perturbation, an approximate AdS-Vaidya spacetime was found [41], and this has become an often-used toy model for homogeneous, isotropic energy injection. Interestingly, the electric field perturbation of [40] yields an exact AdS-Vaidya spacetime.
For D0-brane holography, if one restricts to an ansatz that is spherically symmetric in the "internal" directions (transverse to the D0-brane worldvolume), the supergravity field equations simplify to those of a dilaton-gravity model [20,21] coupled to an additional scalar (the "breathing mode" of the internal sphere) [20]. We will explicitly solve the dilaton-gravity equations (with the breathing mode set to zero) with an arbitrary boundary profile for the dilaton (corresponding to an arbitrary source for the dual operator). As a consequence of the lack of dynamical degrees of freedom in 2d dilaton-gravity, if one turns a source on and off, the late-time bulk metric agrees with the early-time bulk metric, and no net energy was injected. We will find, however, that one can inject energy in the system by considering a boundary condition for the dilaton that is constant at early times and evolves to a different constant value at late times. In field theory language, this corresponds to starting with a thermal state and ending with a thermal state at a different temperature (and with a different value of the coupling constant).
Concretely, in Appendix B we derive the following exact analytic solution of IIA supergravity, expressed in a dual frame in which the 2d metric is asymptotically AdS (see Section 2, where more details will be provided): where M 0 is a mass parameter and φ (0) (v) is a function that one is free to choose as Dirichlet boundary condition for the dilaton field φ and which we also call dilaton source. The metric (1) describes a black hole with mass M = M 0 e 4 3 φ (0) and is asymptotically AdS 2 × S 8 for x → 0, where the timelike boundary of AdS 2 is located. Provided that we have M 0 = 0, a non-constant dilaton boundary value φ (0) (v) will effectively result in a non-constant mass term in the metric.
Even though the above collapsing solution allows arbitrary energy injection patterns, in this paper, we will mostly consider the thin-shell limit in which a black hole spacetime of some initial temperature is glued to a black hole spacetime of higher temperature at a null surface v = v 0 . This can be achieved by assuming the following profile for the dilaton source: with φ 0 a negative constant. For computational simplicity, we will often consider the φ 0 → −∞ limit in which the initial temperature vanishes and the early-time geometry is vacuum v > v 0 : At least within an energy range to be discussed in the next section, this solution is holographically dual to matrix theory excited ("quenched") away from equilibrium through energy injection. However, the solution (1)-(2) does not describe propagating degrees of freedom and will not be sufficient for computing non-trivial correlation functions. As a dynamical probe of this background, we then consider fluctuations in the size of the compact S 8 , i.e., the breathing mode. This mode has already been considered in previous holographic works [19,20] and has been identified in [19] to be dual to a matrix theory operator T −− , to be defined in (35), by matching of generalized conformal scaling dimensions [42][43][44]. Our setup will allow us to holographically compute its retarded two-point correlation function in the quenched dual state, thereby providing a first non-trivial observable which, in the future, one may hope to compare with direct matrix theory computations. Predictably, the late-time behavior of this correlation function is dominated by the lowest quasinormal mode of the final state black hole.
AdS 2 holography has a reputation for being very subtle and relatively poorly understood (see [45][46][47][48][49][50][51][52] for a sampling of the literature, with an emphasis on recent discussions), so one might wonder why we did not run into problems when considering AdS 2 backgrounds and excitations thereof. To see the difference between our D0-brane holography and what is usually referred to as AdS 2 holography, note that our AdS 2 solution arises in the dual frame, in which the effective dilaton-gravity action takes the form (43) with constant b = 25/4. This action has a nontrivial dilaton kinetic term, the removal of which would require a dilaton-dependent rescaling of the metric. After such a rescaling, the metric would no longer be asymptotically AdS 2 . This should be contrasted with conventional AdS 2 holography, which considers asymptotically AdS 2 solutions in the frame without a dilaton kinetic term, which turns out to be more subtle. More specifically, subtleties such as the absence of finiteenergy excitations for fixed asymptotics arise for AdS 2 solutions with constant dilaton. 1 In the theory we are considering here, the dilaton field depends at least on the radial coordinate for all solutions of the equations of motion, and holography works in a way similar to the running dilaton solutions considered in [52].
The paper is organized as follows. In Section 2 we review the duality between matrix theory and IIA supergravity originally conjectured in [18], setting thereby the conventions that are going to be used throughout this work. In Section 3, we study the bosonic part of type IIA supergravity with asymptotically AdS 2 × S 8 geometry in the dual frame. In order to simplify the problem as much as possible, we only consider spherically symmetric solutions, leading to a two-dimensional effective theory describing the metric, the dilaton and the breathing mode accounting for the S 8 size dynamics. This mode is the only propagating physical degree of freedom in that system, and will be our probe for computing non-trivial correlations functions in the quenched dual state. It is important to note that the breathing mode cannot be considered nonperturbatively as noted in [20], because it deforms the boundary away from AdS 2 (see Appendix C). In terms of matrix theory, its dual operator T −− is irrelevant and cannot be sourced nonperturbatively. Nonetheless, a proper perturbative treatment of the breathing mode is expected to correctly reproduce the correlation functions of the dual matrix theory operator [20,53]. In Section 4 we perform the holographic renormalization procedure [54]. Knowing the fields' asymptotics near the AdS 2 boundary as well as the on-shell effective action, local boundary counterterms are added in order to cancel boundary divergences. This is part of the precise holographic description of matrix theory. In earlier works holographic renormalization has been performed for various cases, including non-linear gravity-dilaton solutions [21] and breathing mode perturbations around pure AdS 2 × S 8 [20]. Related work also includes [55]. A general discussion of holographic renormalization in the presence of irrelevant operators deforming the AdS boundary can be found in [53]. Here we consider breathing mode perturbations around non-linear gravity-dilaton background solutions, allowing in particular for time-dependent backgrounds of the form (1)- (2).
In Section 5 we compute the retarded boundary-to-bulk propagator of the breathing mode in the case of pure AdS 2 and in the more interesting case of the thin-shell solution (1)- (6), which is dual to a quenched state in matrix theory. For the retarded propagator in the latter case, we use numerical evolution and show that its asymptotic value near the boundary is rapidly dominated by a single decaying and oscillating mode after crossing of the shell located at v = v 0 . We also show that the associated single complex frequency dominating the retarded two-point function in this quenched state with final temperature T coincides with the lowest quasinormal mode frequency of breathing mode fluctuations around a static black hole at the same temperature T . The first and second quasinormal mode frequencies are therefore computed in Appendix G using a numerical shooting method. Using these results, we holographically derive in Section 6 the retarded non-equal-time two-point function of T −− . Earlier computations of holographic two-point functions in equilibrium states (as opposed to our far-from-equilibrium setting) can be found in [19,56,57].

Review of IIA Supergravity -Matrix Theory Duality
We start by reviewing the duality originally presented in [18], looking only at terms relevant for the present work and setting our conventions. Useful references include [19,58,59]. The bosonic part of the 10d type IIA supergravity action in string frame is with g s and √ α = l s being the string coupling and the string length, respectively. This action involves the metric, a scalar dilaton φ and a gauge potential C M with field strength F M N = ∂ M C N − ∂ N C M and density F 2 ≡ F M N F M N . This system admits a solution representing N coincident electric D-particles at the origin [60]: where H is a single-centered harmonic function on the Euclidean space labeled by Cartesian coordinates x i , given by It has been conjectured that the near-horizon limit or decoupling limit of the above D-particle background is dual to matrix theory [18]. Explicitly, this decoupling limit is where the energy is kept fixed while taking the limit, and the Yang-Mills coupling of matrix theory is identified with Performing a Weyl transformation on the string frame metric (8) in which the action reads For further simplification, we apply the following fields redefinitions, where we define β 1 ≡ 5×5 4/5 4×2 1/10 (3π) 3/10 , bringing the action to the form Finally, by performing the coordinate redefinition the D-particle background in the decoupling limit (15) becomes The geometry of this dual frame metric is therefore manifestly that of AdS 2 × S 8 , with asymptotic boundary located at z = 0. Moreover, it has been argued in [61] that the coordinate z of the dual frame is proportional to the inverse energy scale of the boundary theory, making the dual frame a natural holographic frame. It is also important to note that all relevant expressions appearing in (19)-(23) depend on quantities held fixed in the decoupling limit (12). The validity of the supergravity description requires the string frame curvature and dilaton to be small [18]: These conditions can be conveniently summarized by We see that near the AdS boundary (in the UV regime of the dual field theory), supergravity loses its validity because of strong curvature; in fact, very near the boundary perturbative matrix theory becomes valid. On the other hand, string loop corrections invalidate classical supergravity for large values of z (in the IR regime of the dual field theory). The regime of validity of supergravity can be made parametrically large, however, by working in the large N limit and at large values of the 't Hooft coupling g 2 Y M N . There, the gravitational description can be trusted for matrix theory observables that are neither dominated by the extreme UV or IR parts of its spectrum. Indeed, in [62] Monte Carlo computations of various two-point functions in the matrix theory ground state have been shown to agree well with earlier holographic predictions, at least when neither of the two extreme regimes mentioned above are considered.
On the other side of the duality, matrix theory is the non-abelian U (N ) gauge theory describing the decoupled low-energy dynamics of N probe D-particles on a flat background. Its action can be written as [63] where X i with i = 1, ..., 9 are N ×N hermitian matrices, θ is a 16-dimensional SO(9) spinor whose components are N × N Grassmann matrices, γ i are the associated SO(9) gamma Although in the present discussion we considered a supergravity background with pure AdS 2 × S 8 geometry, the duality applies in principle to any background with asymptotically AdS 2 × S 8 geometry (in the dual frame). The black D0-brane solution giving rise in the decoupling limit to an AdS 2 black hole has for example been considered in the original paper of Itzhaki et al [18], and associated holographic computations have been done in [56].

Dimensional Reduction of IIA Supergravity
We are now going to study IIA supergravity for asymptotically AdS 2 × S 8 geometry in the dual frame, and we will restrict our attention to solutions that are spherically symmetric (l = 0 Kaluza-Klein modes). For this, we use the dual frame action (19), dropping for simplicity the tilde from all symbols together with the fixed overall factor, which can be easily restored: where γ is the induced metric at the asymptotic boundary whose extrinsic curvature is K, S GH being the standard Gibbons-Hawking term necessary to have a well-posed variational principle with Dirichlet boundary conditions on variations of the metric. The equations of motion obtained from varying this action together with the Bianchi identity associated to the gauge field read We now postulate that the 10-dimensional spacetime is the warped product of a 2dimensional Lorentzian spacetime and a Euclidean S 8 , with the following block diagonal metric form: We use Greek indices for labeling the two-dimensional spacetime and Latin indices for the S 8 part. This metric ansatz almost decouples the two subspace metrics, except for the breathing mode b(x ρ ). This mode has been identified in [19] as being dual to the matrix theory operator [64] T where we define In (35), STr stands for the symmetrized trace operation, consisting of full symmetrization over all possible orderings of the F -operator factors, followed by taking the ordinary trace. This identification has been performed by matching of generalized conformal scaling dimensions [42][43][44]. It is important to keep in mind that, in the conventions we are using, g Y M in (35) is given by the constant value (13), which does not include a dilaton dependence, and hence does not depend on time.
We also choose an ansatz for the gauge field energy-momentum tensor F , often referred to as a Freund-Rubin ansatz [65]: for which A(x R ) is a priori a function of all coordinates, and µν is the Levi-Civita symbol of the two-dimensional spacetime associated with coordinates x ρ . This ansatz amounts to assuming a purely radial electric field. Furthermore we assume that the dilaton is independent of the S 8 coordinates, The function A appearing in (39) is actually expressed by equations (32)-(33) through b and φ. In this situation, it is advantageous to develop an effective theory involving the two-dimensional metric g µν , b and φ. While we give the details of this dimensional reduction in Appendix A, we point out here a subtlety that arises when performing this dimensional reduction at the level of the action.
In an asymptotically AdS setting, g µν , b and φ satisfy Dirichlet boundary conditions at the conformal boundary of AdS. If (39) is imposed (with A expressed through b and φ), this will imply that the gauge field strength F (rather than the gauge potential C) satisfies Dirichlet boundary conditions. If we want to produce an effective action for g µν , b and φ we should thus start with a ten-dimensional action which allows a consistent variational principle based on variations of F satisfying Dirichlet boundary conditions. This requires adding the following boundary term to (28): with n M being the outward-pointing normal unit vector of the boundary. Similar boundary terms have been considered in [66,67]. Equation (42) highlights the fact that the boundary term effectively corrects the coefficient of the Maxwell term when the action is evaluated on field configurations satisfying the gauge field equations of motion. This feature will correspondingly be reflected in the action of the dimensionally reduced theory, in which the gauge field is expressed through g µν , b and φ. (One can straightforwardly verify that naively substituting the two-dimensional ansatz in the ten-dimensional action without taking into account the boundary term subtlety would have resulted in an action that no longer reproduces the correct equations of motion of the dimensionally reduced theory.) Putting everything together (see Appendix A for the details) one finds that the action (28) is effectively reduced to two dimensions, where all geometric quantities are constructed from the two-dimensional metric, and C = 1 2 ( 14 5 ) 2 ( 25 4 ) 8 . Furthermore one can check that this effective action leads to the correct equations of motion. We are effectively considering two scalars φ and b living on a twodimensional spacetime (with non-minimal couplings).

Holographic Renormalization
In this section, we give the result of the holographic renormalization procedure in the usual Fefferman-Graham gauge [54]. This procedure consists in adding local, gauge invariant boundary counterterms to the on-shell action such that divergences cancel. The renormalized on-shell action is then interpreted as the generating functional of correlation functions of the boundary matrix theory. (The discussion in [54] is for Euclidean correlation functions. In Lorentzian signature, there are various correlation functions and extra care is needed; see for instance [68][69][70][71]. ) We recall that we are now effectively working with a two-dimensional metric, whose Fefferman-Graham form is Consider field fluctuations about a background, where we can recognize the background value of the breathing mode from (21). As mentioned in the Introduction, one cannot consider a dynamical breathing mode at the nonlinear level. Indeed, its presence strongly deforms the boundary away from AdS 2 , as is manifest in the near-boundary asymptotic expansion of Appendix C. On the other hand, treating this mode perturbatively produces the holographic correlation functions and will be sufficient for our purposes. In the path integral formalism applied to matrix theory this amounts to considering sources in the generating functional which allows for the computation of n-point correlators by functional differentiation, and setting these to zero at the end of the calculation.
In order to perform holographic renormalization, we only need to know about the nearboundary asymptotics of these fields, which can be easily deduced by studying the equations of motions in the limit z → 0. By expanding the fields in (possibly fractional) powers of z, one usually finds up to two undetermined coefficients per bulk field. The first one is typically the leading coefficient in that expansion and is usually referred to as the field source. It needs to be specified as Dirichlet boundary value. The second undetermined coefficient is a subleading one and is directly related to the expectation value of the operator dual to the bulk field under consideration, in the corresponding matrix theory state. A necessary ingredient for computing correlation functions in a specific state is to know the solution deeper in the bulk (not only near the boundary), which fixes uniquely the subleading coefficient. We postpone this to Section 5, as it is not required for carrying out holographic renormalization.
In Appendix C one can find details of the near-boundary asymptotics of the fields. The outcome is that there are sources for each of the physical fields considered here, but only the breathing mode has a subleading coefficient left undetermined by this asymptotic analysis. This comes from the fact that the metric and dilaton are non-propagating fields in this two-dimensional system, whereas the breathing mode does propagate.
As already said, holographic counterterms are boundary terms added to the action such that, when evaluated on-shell, the small z cutoff (corresponding to a UV cutoff in the dual field theory) can be removed without generating boundary divergences. This means that the divergences coming from the on-shell action are to be cancelled against those coming from the counterterms. In Appendix D we give the gauge invariant form of the on-shell action up to quadratic order in the expansion parameter . In this work we are mainly interested in non-equal-time two-point correlation functions, and we will therefore only give the relevant boundary counterterms. These correlators encode non-trivial dynamical information about the dual matrix theory state. For the interested reader, we comment on one-point functions in Appendix E and explain why the information they encode is somewhat trivial. Other counterterms are needed to ensure finiteness of the renormalized on-shell action, but these will not contribute to the values of one-point functions and nonequal-time two-point correlators. 2 Neither should we care about terms local with respect to sources (those are often referred to as contact terms), since they contribute to equaltime correlators only. Here we therefore give the gauge invariant expression of the relevant boundary counterterms: One can check that leading divergences are thereby canceled against those present in the on-shell action (131) and that the finite renormalized on-shell action is given by where in the first line we implied a limit removing a small z cutoff. The final expression explicitly depends on the Fefferman-Graham expansion coefficients left undetermined by the asymptotic analysis, as well as on the mass parameter M 0 characterizing the most general solution to the gravity-dilaton system (1)-(2); see Appendix C. The leading terms in the asymptotical expansions, which are also called sources, are f is related to the expectation value of the matrix theory operator T −− as we will see in Section 6. In the following sections, we will use the notation b source , hoping that this will make it easier for the reader to follow the logic of the computations. The renormalized action is to be understood as the generating functional of matrix theory, and will be used in Section 6 to compute non-equaltime correlators.
2 For example, the full counterterm at order zero in is given by where differential operators are boundary ones. Note that the first term is identical to (48) at order zero in .

Retarded Boundary-to-Bulk Propagators
In this section we compute the retarded boundary-to-bulk propagators of the linearized breathing mode fluctuations, for two different backgrounds. The first one is simply pure AdS 2 . The second one is the thin shell limit (5)-(6) of the exact analytic solution (1)- (2) presented in the Introduction, whose dual representation can be thought of as matrix theory quenched from its ground state to a non-zero temperature state by sudden energy injection. For convenience we recall here the background solution. The metric in ingoing Eddington-Finkelstein gauge adapted to the description of collapsing spacetimes is In this gauge, the most general background solution is given in closed form by We provide a derivation of this solution in Appendix B. The free function φ (0) 0 needs to be specified as Dirichlet condition and is interpreted as source for the scalar operator dual to the dilaton.
We then consider fluctuations about this general background solution, allowing for the breathing mode which, as mentioned earlier, is the only propagating degree of freedom in the reduced action (43). Using the linearized equations of motion, one can derive a decoupled equation that this mode must satisfy, A boundary-to-bulk propagator is a solution to this equation which reduces to a delta function profile at the boundary. We now compute the retarded boundary-to-bulk propagator, which vanishes outside the future lightcone of the boundary point where the delta function has support, for both backgrounds under consideration.
Vacuum AdS The pure AdS 2 background solution is obtained by setting M 0 = 0 and φ (0) 0 = 0 for all times in (52)- (53). We derive the retarded boundary-to-bulk propagator (with delta function support at t = 0 at the boundary) 3 in Appendix F; its expression in Eddington-Finkelstein gauge is Note that the ingoing lightlike coordinate v reduces to the boundary time t in the near boundary limit z → 0. In that limit one can show that G AdS R (t, z) z −14/5 δ(t) which is the appropriate asymptotic behavior for a delta source at the boundary; see Appendix C. Therefore a generic source profile b source Thin Shell Collapse Assuming a fixed non-zero mass parameter M 0 , the choice of a timedependent dilaton source profile φ (0) 0 as Dirichlet boundary condition offers the possibility of having collapsing solutions. Here we consider a thin shell limit in which vacuum AdS is glued to a black hole background at a null surface v = v 0 . This is achieved by choosing the following profile for the dilaton source: Note that we can also use any dilaton source profile after the gluing surface v = v 0 , resulting in a generic time-dependent black hole background. Note also that the limit (58) is perfectly acceptable since one expects supergravity to be valid for small value of e φ (except if one were to consider observables with strong support in the deep infrared regime, as discussed in Section 2). Consider a delta function source b source 1 (t) = δ(t − t s ) on the boundary before the shell collapse, at t s < v 0 . The retarded boundary-to-bulk propagator on the pure AdS segment is already known, and one simply has to numerically continue it across the shell located at v = v 0 . Therefore, by using G Shell is obtained, one can extract the asymptotic coefficient that was left unfixed by asymptotic analysis for Dirichlet boundary condition. This coefficient is the one appearing at order z 28/5 near the boundary; see Appendix C for more details. In particular, in the limit z → 0 and for non-equal times one has and the leading coefficient G is directly related to the retarded two-point function of the matrix theory dual operator, as we will soon see. This coefficient as well as subleading ones can be extracted from the numerical solution by fitting a power series of the form (125) appropriately expressed in coordinates (v, x), close to the boundary x → 0. Note that the first three towers of terms in the expansion (125) are effectively missing whenever the breathing mode source is turned off, which is presently the case for t = t s .
The main result coming out of the numerical evolution from the gluing surface at v = v 0 is that this leading coefficient is very rapidly dominated by a single oscillating and decaying mode of the form G The more massive the black hole, the more rapidly this mode dominates. The timescale for perfect matching (within numerical errors) is of the order of the inverse horizon radius r h = M 5/14 0 of the associated final black hole. We show an example of this behavior in Figure 1. In Figure 2, we give the extracted values of ω R and ω I for various black hole masses. This leads to the following frequency-to-temperature ratios: where the Hawking-Unruh black hole temperature is given by [56] T H = 7 10π r h .
Importantly, the single complex frequency ω = ω R − iω I for a choice of final black hole temperature T H can be matched to the lowest quasinormal frequency of breathing mode fluctuations around a static black hole at the same temperature. Indeed, using a numerical shooting method we compute in Appendix G the first and second quasinormal frequencies whose values are given by Second QNM: Similar results in a related context have been found in [72,73]. In Appendix G we also derive the universality of the frequency-to-temperature ratio, implicitly indicated in (63)-(69).

Linear Response in Matrix Theory
We are now ready to compute non-equal-time retarded two-point functions of the matrix theory operator T −− , as probe of thermalization. We assume a δ-function source profile b source 1 (t) = δ(t − t s ) on the boundary and recall the relation (56) between this source and the induced retarded bulk solution. We note that the field dual to T −− is naturally given by e − 6 7 φ b rather than simply b. Indeed, couplings of the boundary operators to their dual bulk fields can be recovered from considering the effective D-brane action in a given bulk field background. For example, at the level of a single D-particle (matrix theory being the low-energy limit of the U (N ) generalization of this action), the metric-dependent part of the action in the dual frame defined by (14) is given by Linearizing around the Minkowski background, one gets where we have defined v i ≡ẋ i . One can therefore conclude that fields coupling to operators like v i v j or v 4 (the latter being proportional to T −− in this U (1) case) are e − 6 7 φ 0 h rather than the simple metric perturbations h. (Of course, if the background dilaton φ 0 is timeindependent, this distinction becomes unimportant.) According to the usual holographic dictionary [54] and focusing on the order 2 piece of the renormalized on-shell action (50) that will lead to a non-equal-time correlator, we then have where . . . stands for terms analytic in b source 1 (t) which only contribute to equal-time twopoint functions. Since the one-point function in absence of source is trivial, O(t) s=0 = 0, the response function is simply and we identify the non-equal-time retarded two-point function of the matrix theory as Ground state Two-Point Function As a first special case we consider a pure AdS 2 background, which is dual to the matrix theory ground state. The retarded boundaryto-bulk propagator is given in (55) while the dilaton source is φ One can check that the power decay agrees with previous work [19].
Thermalizing Two-Point Function The collapsing thin shell setup that we studied in Section 5 is dual to matrix theory initially in its ground state and suddenly excited by energy injection. The retarded thermalizing two-point function in this quenched state with insertion times t s < t quench ≡ v 0 < t is rapidly dominated by a single complex frequency for t > t quench , as analyzed in Section 5. Its expression is given by where G In Figure 3, we compare thermalizing and ground state two-point functions. One can see that the former relaxes more rapidly, as can be generally expected from the exponential decay of quasinormal modes. (See also equation (62) for the generic behavior.) Of course, this is consistent with the general picture of how thermalizing observables should evolve.

A Details of the Dimensional Reduction
In this appendix we work out the details leading to the effective two-dimensional action (43). We also show how one gets the associated equations of motions. With the assumptions (34), (39), (40) of Section 3, equations (32)- (33) impose with C a constant. A simple way to identify the value of C that supports asymptotically AdS 2 × S 8 solutions of (30)-(33) is by computing Any asymptotically AdS solution must agree in the asymptotic region with (21)- (23), which implies that that the second one gives the unique equation while the third one is trivially satisfied. Instead of working with the fields {g M N , φ}, one can work with {g µν , b, φ}. First one notices that all quantities in the previous equations can be written in terms of {g µν , b, φ}.
In particular, we have where (2) indicates that all geometric quantities are constructed from the two-dimensional metric g µν only, and where n µ is the outward-pointing normal vector to the timelike boundary. These relations can be inserted in (85)-(87) in order to get the effective two-dimensional equations of motions. Using (83),(88),(89) together with √ −g = − g (2) √ g S 8 b 4 , one can show that the action (28) becomes Note that the bulk term −Cb −4 in the effective action is obtained not only from the term − 1 4 e 6 7 φ F 2 in the bulk part of the action (28), but also from the term 1 2 e 6 7 φ F 2 in the boundary part of the action S extra (41). One can further check that this gives the correct effective EOM's by extremization of S ef f . Although S extra is explicitly written as a Maxwell bulk term in (42) for fields configurations satisfying the gauge field EOM (32), one would usually think that boundary terms should not contribute to the EOM's. In this case, there is an underlying reason why it does contribute. Indeed, fixing the gauge we deduce from (39)-(82) that A possible gauge potential is therefore given by which is a non-local functional of the fields g (2) µν , b, φ. In particular bulk variations of those fields induce boundary variations of C 0 . Therefore, in order to get the EOM's for those fields by action extremization, it is not sufficient to consider variations of g (2) µν , b, φ and C 0 localized in the bulk. An alternative way of seeing this is that for the subclass of fields configurations (97)-(98), one can show that the boundary term S extra effectively produces a bulk term through the non-locality of C 0 : where n M = −zδ M z is the outward-pointing unit normal vector to the boundary, and γ its induced metric. The same result is of course directly obtained from the expression of S extra given in bulk form in (42).
In the main body of this paper, we will exclusively use the two-dimensional effective system. We therefore omit the (2) superscript there, since all geometrical quantities will refer to the two-dimensional metric.

B Derivation of General Background Solution
We provide a derivation of the general background solution (52)-(53) in ingoing Eddington-Finkelstein gauge In these coordinates, the equations of motion for the background fields are given by 16 21 Integrating (105) readily gives where φ (0) (v) and A(v) are arbitrary functions of v. Knowing φ(v, x), equation (102) is just an ordinary differential equation for g with variable x. The solution is The two functions C 1 (v) and C 2 (v) are then constrained by (103)-(104), with M 0 a free mass parameter. One can check that (106) is also satisfied. The function A(v) can be gauged away through the replacement bringing the solution to the form (52)-(53), φ(v,x) = φ (0) (v) + 21 10 logx.

C Bulk Field Asymptotics
In this appendix we derive the near-boundary asymptotics of the bulk fields (45)- (47), by analyzing the equations of motions in Fefferman-Graham gauge (44) in the limit z → 0. We start with the background fields. For illustrative purposes, we display their equations of motion, which can be derived by variation of (94) or by use of (85)-(92) together with b(t, z) = 25 4 : 4

8
Let us assume the following leading behavior, Plugging this ansatz into (113)-(117) and solving at leading order in the limit z → 0, we get Moreover f (t) and φ ( 14 5 ) 0 (t) coming from such subleading terms in the equations of motion. The resulting asymptotic expansion is 0 (t) . . .
Let us comment on the structure of these expansions. Except for the logarithmic term in the dilaton expansion that can easily be recognized from (22), there are five towers of terms, each of which has an internal spacing of two in powers of z. Such a spacing is very familiar from AdS-CFT and is referred to as Fefferman-Graham expansion [54]. It follows from the equations of motion (113)-(117) being second order differential in z. The first tower starts at order z 0 with some undetermined coefficients f 0 called sources that need to be specified as Dirichlet conditions. Then the second tower starts at order z 14/5 with coefficients that are related to the expectation values of the operators dual to the considered bulk fields. In contrast to familiar examples of AdS-CFT dualities, their dependence on sources is local: with M 0 an undetermined constant that can be interpreted as a mass parameter characterizing the most general solution of the gravity-dilaton system (1)-(2) derived in Appendix B.
The reason for this complete locality is that the dilaton and metric are non-propagating fields in this two-dimensional system. Finally the remaining towers are generated by the non-linearities of the equations of motion. 5 A similar structure has been displayed for a somewhat simpler system in [53]. One should note that all coefficients depend locally on the sources f 0 , although we will refrain from giving the precise relations. These can be easily recovered by inserting (120)-(121) in the equations of motion and solving order by order in z.
The asymptotics of fluctuations are similarly obtained. From the linearized equations of motion, a decoupled one can be derived for the breathing mode, Performing a similar analysis as for the background fields, we find the following asymptotic expansion: Again there are five towers of equally spaced terms in this expansion. The first one starts at order z −14/5 with a source that needs to be specified as Dirichlet boundary value. The fourth tower starts at order z 28/5 with a coefficient left undetermined by the asymptotic analysis. (In Euclidean signature, it is also a functional of the sources {f }, but the dependence is nonlocal and cannot be determined from a boundary asymptotic analysis only; specific solutions need to be known deeper in the bulk in order to fix its dependence in sources, which is in general not local or analytic. In Lorentzian signature, it also depends on the choice of state, or initial conditions, and encodes information on a propagating degree of freedom.) The remaining towers are generated by the non-linearities of the equation of motion (124) with respect to the background fields. A particular point to note is that the coefficient b (0) 1 is also undetermined but needs to be set to zero, as otherwise it would change the breathing mode background value in (47).
Still at linear order in fluctuations, the dilaton and metric being non-propagating fields are constrained by the remaining equations of motions. In particular they enjoy the same power expansion, with running index given by j ∈ {− 14 5 , . . . , 2, . . . , 14 5 , . . . , 28 5 , . . . , 42 5 , . . .}. One can again check that all coefficients in (125)-(127) can be expressed as local functionals of the sources } and the remaining undetermined coefficient b . 6 A striking feature of expansions (125)-(127) is that fluctuations seem to dominate over background fields in the limit z → 0, and deform the geometry from being asymptotically AdS 2 . This is the reason why we cannot consider a dynamical breathing mode nonperturbatively, and is mirrored by the fact that its dual matrix theory operator T −− is an irrelevant one.

D On-Shell Evaluation of the Action
Considering fluctuations of the form with constant b 0 = 25 4 , we give the gauge invariant expression of the action (43) when evaluated on-shell, up to quadratic order in the expansion parameter : 6 Local dependences on b ( 28 5 ) 1 arise at order j = 28 5 and further in the above series expansions. with theory state is somewhat trivial in the sense that it is non-dynamical. As mentioned in the Introduction, sources for the irrelevant operator T −− can be used in the generating functional of the matrix theory for computing n-point correlators, but have to be set to zero afterwards. From the holographic point of view, the breathing mode source b source plays this role; see Appendix C. Using the expression for the renormalized on-shell action (49) we derive the various one-point functions, One can check that the following trace and diffeomorphism Ward identities are satisfied [21]: This actually is a simple consequence of (122)-(123). Let us briefly comment on (134)-(136). As can be seen, all these one-point functions simply follow the dilaton source profile φ (0) 0 (t). Looking at a quenched solution of the type (5)-(6) that is dual to sudden energy injection in the dual matrix theory groundstate, we see that one-point functions simply jump between zero and a constant set by the mass parameter M 0 (related the to final temperature through equation (65)). Expectation values instantaneously take their thermal value and no dynamical thermalization can be seen from these observables. In order to see dynamical thermalization in matrix theory, we therefore need to look at higher-point correlators like the one we study in Section 6.

F Boundary-to-Bulk Propagators in Vacuum AdS 2
We review the derivation of breathing mode boundary-to-bulk propagators in pure AdS 2 . For convenience, we start in Poincaré coordinates and then give the expression in Eddington-Finkelstein gauge.
In the Euclidean version of vacuum AdS 2 described by the metric the unique breathing mode boundary-to-bulk propagator is given by Indeed, one can show that it satisfies the equation of motion (54) and is normalized to a delta function source at the boundary, It is well-known that this unique Euclidean boundary-to-bulk propagator corresponds to the Feynman time-ordered one after Wick rotation. For this we perform the analytic continuation τ = e iθ t from θ = 0 to θ = π 2 (see [72] for a similar treatment) and obtain iG AdS where the overall minus sign is needed in order to satisfy We then write the retarded propagator as which is correctly normalized at the boundary. Note that the retarded propagator vanishes outside the future lightcone as expected. By performing the variable change v ≡ t − z, we find the expression of the Feynman and retarded boundary-to-bulk propagators in ingoing Eddington-Finkelstein coordinates,

G Lowest Quasinormal Modes
In this appendix we numerically compute the first and second quasinormal frequencies of breathing mode fluctuations around static black holes at temperature T H and show that the lowest one matches the single complex frequency dominating the late-time behavior of the leading coefficient (62) of the retarded boundary-to-bulk propagator when considered on the collapsing thin shell background solution (58)- (59). We also derive the universality of frequency-to-temperature ratios that the results (63)-(64) were indicating, using an argument given in [75]. For reviews of quasinormal modes, one can consult [76][77][78].
We therefore need to determine the lowest QNM frequencies of linear breathing mode perturbations around a static black hole solution. The discussion will closely follow that of [75] since the asymptotic geometry in both cases is that of AdS.
The static black hole background considered here is the near-horizon limit (or decoupling limit) of the black D-particle background found in [60]. Its expression in the dual frame is given by 7 together with the associated background dilaton φ(r) = − 21 10 log r.
In these coordinates, the AdS boundary is located at r → ∞. The horizon radius of such a black hole is related to its mass through r h = M 5/14 , while its Hawking-Unruh temperature is given by (65). In this coordinate system the decoupled equation of motion satisfied by the linearized breathing mode perturbation is
(155) 7 Note that the coordinate r here is different from that defined in equation (11).
It is obvious that (154) has the exact same form as (153) up to the replacement r h → r h ≡ ar h . This has the important consequence that the quasinormal frequency spectrum scales with the horizon radius r h , or equivalently with the black hole temperature, which is the only scale in the problem. To see this, it suffices for example to consider the timedependent part e −iω[r h ]t = e −iω[r h /a] at of such a quasinormal mode, which has to be equal to e −iω[r h ]t from the form invariance of the equation of motion, and therefore of its solutions. We conclude that the frequency must be linear in r h , or equivalently where the constant is different for each of the quasinormal modes forming a discrete spectrum.

G.2 Frobenius Expansion Near the Horizon
Using the following ansatz for the breathing mode perturbation: b 1 (t, r) = r −9/10 χ(r)e −iωt , and defining the tortoise coordinate r * through dr dr * = N (r), one arrives at a Schrödinger-like equation In order to study numerically the behavior of ingoing modes at the horizon, we factorize the associated singular part in the solution: and end up with the equation of motion where we are now using Fefferman-Graham coordinates with z = 1/r and where we define the inverse horizon radius h ≡ 1/r h . One can check that the functions s(z), t(z) and u(z) are regular between the horizon (z = h) and spatial infinity (z = 0), so that ζ admits a Frobenius series expansion close to the horizon: where α is a complex parameter which should have two possible values, corresponding to ingoing and outgoing solutions. In order to determine α and the coefficients {a n (ω, h)} n∈N , we also expand the known functions s(z), t(z) and u(z) around z = h: Then, by plugging (165)-(168) in (161) and looking at the lowest order in (z − h), one gets the indicial equation s 0 α(α − 1) + t 0 α + u 0 = 0.
One can indeed identify those solutions as ingoing and outgoing ones by noting that close to the horizon [75] such that ingoing and outgoing modes indeed match the two solutions we have found: ingoing solution : b 1 (t, r * ) ∼ e −iω(t+r * ) , outgoing solution : b 1 (t, r * ) ∼ e −iω(t−r * ) = e −iω(t+r * ) e 2iωr * e −iω(t+r * ) (z − h) iω/κ Since we are interested in quasinormal modes, we impose that nothing comes out of the horizon by choosing α = 0. With this choice, we can return to (161) and solve iteratively for the a n (ω, h), leading to the recursive formula [75] a n = − 1 n(n − 1)s 0 + nt 0 ∞ k=0 [k(k − 1)s n−k + kt n−k + u n−k ] a k . (175)

G.3 Numerical Shooting Method
We will determine numerically the lowest quasinormal complex frequency with ω R , ω I ∈ R. From the positivity of V (r) outside the horizon radius, it can be proven that ω I must be a strictly positive real number, such that only decaying solutions are allowed [75]. A simple method is to evaluate the Frobenius expansion (165) of ζ(z) with α = 0, at a point close to the horizon where it is valid, 8 and then integrate (161) numerically up to spatial infinity, for various values of ω. Then quasinormal frequencies are those for which the solution is normalizable (does not diverge) at spatial infinity. Of course one has to truncate the Frobenius expansion at some order N and check for the convergence of the results for increasing N . For the purpose of this work, it was sufficient to truncate the series at order N = 40.
Using this numerical shooting method, we give results for the first and second lowest quasinormal frequencies: Second QNM: As one can see, the value of the first quasinormal frequency is indeed in good agreement with the single complex frequency (63)-(64) dominating the leading coefficient of the retarded boundary-to-bulk propagator, whose time evolution is computed in Section 5. 8 The radius of convergence of the Frobenius series is at least equal to |h − p|, where p is the singular point of (161) which is the closest one to h in the complex z-plane. In this case the closest singular point is e iπ/7 h, such that the radius of convergence is at least 0.445 h. We have chosen z ini = 0.75 h as initial point for our numerical integration, which therefore lies in the convergence disk centered around z = h.