Quantum chaos, thermalization and entanglement generation in real-time simulations of the BFSS matrix model

We study numerically the onset of chaos and thermalization in the Banks-Fischler-Shenker-Susskind (BFSS) matrix model with and without fermions, considering Lyapunov exponents, entanglement generation, and quasinormal ringing. We approximate the real-time dynamics in terms of the most general Gaussian density matrices with parameters which obey self-consistent equations of motion, thus extending the applicability of real-time simulations beyond the classical limit. Initial values of these Gaussian density matrices are optimized to be as close as possible to the thermal equilibrium state of the system. Thus attempting to bridge between the low-energy regime with a calculable holographic description and the classical regime at high energies, we find that quantum corrections to classical dynamics tend to decrease the Lyapunov exponents, which is essential for consistency with the Maldacena-Shenker-Stanford (MSS) bound at low temperatures. The entanglement entropy is found to exhibit an expected"scrambling"behavior - rapid initial growth followed by saturation. At least at high temperatures the entanglement saturation time appears to be governed by classical Lyapunov exponents. Decay of quasinormal modes is found to be characterized by the shortest time scale of all. We also find that while the bosonic matrix model becomes non-chaotic in the low-temperature regime, for the full BFSS model with fermions the leading Lyapunov exponent, entanglement saturation time, and decay rate of quasinormal modes all remain finite and non-zero down to the lowest temperatures.


I. INTRODUCTION
Our understanding of quantum chaos has significantly advanced in recent years due to numerous correspondences between chaotic systems and black holes. In particular, it was argued that physical systems which are holographically dual to black holes are maximally chaotic, with the Sachdev-Ye-Kitaev (SYK) model [1,2] and the Banks-Fischler-Shenker-Susskind (BFSS) model (supersymmetric matrix model) [3][4][5] being notable examples on the quantum field theory (QFT) side. More generally, matrix quantum mechanics provides a rather generic system for studying quantum chaos [6][7][8]. Despite this progress, many questions remain open. First of all, a direct demonstration of maximal chaos from the QFT side remains an open problem, with the important exception of the SYK model. A better understanding of the chaotic dynamics of non-Abelian gauge fields is also important for the description of early stages of heavy-ion collisions [9][10][11]. Obviously, for real QCD which should describe this process, holographic duality is not directly applicable. These problems motivate the development of numerical methods for studying quantum real-time dynamics of gauge theories [12][13][14].
Quantum chaos can be described quantitatively in terms of the exponential growth of the out-of-time-order correlators (OTOCs) CðtÞ ¼ h½ŴðtÞ;Vð0Þ 2 i ∼ expð2λ L tÞ ð 1Þ of suitable operatorsŴ,V [15][16][17]. In the semiclassical regime, the growth of OTOCs (1) is governed by the leading classical Lyapunov exponent λ 0 L of a system, CðtÞ ∼ expð2λ 0 L tÞ, at sufficiently large t. Exponential growth of OTOCs has to be contrasted with the time dependence of the conventional time-ordered correlators GðtÞ ¼ hTrðŴðtÞVð0ÞÞi, which are related to dissipative transport responses. Such time-ordered correlators typically exhibit exponentially decaying oscillations characterized by complex-valued quasinormal frequencies [18], the so-called quasinormal ringing [19,20].
While in classical systems Lyapunov exponents can be arbitrarily large, a universal Maldacena-Stanford-Shenker (MSS) bound 1 on the coefficient of exponential growth of out-of-timeorder correlators (1) can be derived in quantum theory under some mild assumptions based on analyticity properties of the OTOCs [17]. This bound is expected to be saturated by physical systems which admit a holographic dual description in terms of black holes in weakly coupled gravity (i.e., large-N, strong coupling limit of holographic QFT). This could be explicitly demonstrated in the SYK model [2], which is expected to be holographically dual to a nearly extremal black hole near zero temperature [1,2]. The BFSS model [3][4][5], obtained by reducing the (9 þ 1)-dimensional supersymmetric Yang-Mills theory down to 0 þ 1 dimensions, has a significantly richer dynamics than the SYK model and also admits a welldefined dual holographic description [21] in terms of black zero brane in the type IIA superstring theory. The BFSS model is also expected to saturate the MSS bound (2) in the strong-coupling regime at a sufficiently low temperature; actually this is the first model in the literature which has been conjectured to be a "fast scrambler" [8]. While the BFSS model is known to be classically chaotic and various aspects near the classical limit have been studied [22][23][24][25][26][27][28][29][30][31], so far not much is known about its real-time dynamics in the quantum regime because of the absence of suitable firstprinciple methods for real-time evolution of many-body quantum systems. Note that, for exactly solvable OðNÞ vector models at large N, the quantum Lyapunov exponents are parametrically suppressed as λ L ∼ T=N [32,33].
Quantum entanglement between different degrees of freedom (d.o.f.) offers a complementary language for a quantitative description of quantum chaos. It is expected that for strongly interacting chaotic systems all d.o.f. become highly entangled under quantum evolution [8,34], even if the initial state is a direct product jΨi ¼ jΨ A i ⊗ jΨ B i of states jΨ A i and jΨ B i of subsystems A and B. The entanglement entropy is expected to exhibit a rapid growth at early times followed by saturation at late times, when the system has already "scrambled" the information contained in subsystem states jΨ A i and jΨ B i [8,35]. In the semiclassical approximation, the growth rate of entanglement entropy at early times is determined by the classical Lyapunov exponents [36,37]. However, beyond the semiclassical approximation the relation between Lyapunov exponents and growth of entanglement entropy could only be demonstrated for quadratic (or approximately quadratic) Hamiltonians [36][37][38] and for models with discrete time evolution [35].
In this paper we report on numerical studies of quantum corrections to the real-time dynamics of the thermal states of the BFSS model and its bosonic sector (bosonic matrix model), addressing in particular quantum corrections to Lyapunov exponents, the relation between Lyapunov exponents and entanglement entropy generation, and quasinormal ringing. We find that quantum corrections from the bosonic sector of the model tend to make the system less chaotic and less dissipative, whereas the contribution of Majorana fermions works in the opposite direction. The characteristic Lyapunov time, entanglement saturation time, and decay time of quasinormal ringing become very long for the bosonic matrix model at sufficiently low temperatures, which roughly correspond to the confinement regime [39,40]. In contrast, for the full BFSS model with fermions these characteristic timescales remain finite even at the lowest energy accessible in our simulations. While at low temperatures at which the MSS bound is expected to be saturated, our approximation is most likely too crude to capture the full dynamics of the model, and our results suggest that quantum corrections from bosonic and fermionic sectors work in a way which is consistent with the MSS bound at lower temperatures and which evades the naive violation of MSS bound by classical Lyapunov exponents λ 0 L ∼ T 1=4 [25] at sufficiently small T. We further demonstrate that the characteristic saturation time for the entanglement entropy is in general shorter than the Lyapunov time τ L ≡ λ −1 L defined by the leading Lyapunov exponent λ L . It appears to be governed by the classical, rather than quantum, leading Lyapunov exponent. The characteristic decay time of quasinormal ringing is found to be the shortest timescale of all.
In order to simulate the real-time dynamics of the BFSS model, we approximate the density matrix of the system by the most general Gaussian function with time-dependent parameters which obey self-consistent equations of motion. Such an approach, which we will refer to as the Gaussian state approximation, is closely related to the semiclassical approximation [41,42] and is extensively used in the context of quantum chemistry [41,43]. For interacting many-body systems which admit a second-quantized QFT description, such as the tight-binding description of electron gas in solids, the time-dependent Gaussian state approximation for QFT is equivalent to the time-dependent Hartree-Fock approximation (see e.g., Chapter 12 of [44]) for the first-quantized many-body Hamiltonian. For fermionic fields interacting with classical gauge fields, this approximation is equivalent to the classical-statistical field theory (CSFT) approximation which is by now a standard tool to study real-time dynamics of fermions interacting 1 We have set ℏ ¼ 1 and k B ¼ 1.
with highly occupied soft modes of gauge fields [14]. An important property of the Gaussian state approximation is that it evolves pure states into pure states (see Appendix C for the proof), which allows one to study quantum entanglement in a consistent way.
As discussed in [45,46], for classically chaotic systems the Gaussian state approximation, surprisingly, works even better than for systems which exhibit regular classical motion and rather accurately describes the quantum evolution at the timescales of order of the classical Lyapunov time. Only some subtle late-time phenomena such as the wave-packet revival are not captured [47]. In [48] we have also compared the Gaussian state approximation with the numerical solution of the Schrödinger equation for a simple classically chaotic Hamiltonian with two bosonic d.o.f. [49] which closely resembles the bosonic matrix model, and we found a good agreement for evolution times t ≤ 2=λ 0 L less than approximately two Lyapunov times in both quantum and classical regimes. These observations suggest that the Gaussian state approximation should be at least qualitatively accurate for the description of real-time thermalization at timescales comparable with the classical Lyapunov time.
In the context of the BFSS model, one of the limitations of the Gaussian state approximation is that the gauge symmetry constraints cannot be fully respected. As a consequence, our simulations correspond to the ungauged version of the BFSS or bosonic matrix models, where no gauge constraints are imposed on the state vectors. Fortunately, the differences between the gauged and ungauged models appear to be minor at least at low temperatures, as conjectured recently in [50] and demonstrated numerically in [39]. Yet another argument in favor of accuracy of the Gaussian state approximation is that, as we will demonstrate, it reproduces the numerical results for the equation of state of the ungauged bosonic matrix model [39] within a few percent accuracy all the way from low to high temperatures.
We start our discussion in Sec. II by briefly reviewing the BFSS model and setting up the notations to be used in the rest of the paper. In Sec. III we explain the Gaussian state approximation for the real-time dynamics of the BFSS model. This approximation is rather general and can easily be extended to other models which admit Hamiltonian formulation. In Sec. IV we discuss the initial state used in our simulations, which we require to resemble the thermal equilibrium state as closely as possible. In Sec. V we present our numerical results. In Sec. VA we demonstrate that quantum corrections make Lyapunov exponents smaller than in the classical system, thus being in agreement with the MSS bound (2). We also clarify the relation of our results to out-of-time order correlators of the form (1). In Sec. V B we study real-time evolution of entanglement entropy and discuss the relation between entanglement generation and Lyapunov exponents. In Sec. V C we consider quasinormal ringing and the temperature dependence of complex-valued quasinormal frequencies.
In the concluding Sec. VI we summarize our findings and outline some directions for further work. Technical details of our simulations are described in Appendixes A-E.

II. A BRIEF REVIEW OF THE BFSS MODEL
In this paper we use the following representation of the Hamiltonian of the BFSS matrix model [3]: In this expression and throughout the paper we use the following notations and conventions: (i)X a i andP a i are canonically conjugate bosonic coordinate and momentum operators with commutation relations ½X a i ;P b j ¼ iδ ab δ ij , which have dimensions of (Mass) and ðMassÞ −1 , respectively. This is a natural convention because X a i correspond in fact to the components of the gauge field vector in (9 þ 1)-dimensional super-Yang-Mills theory.
ments of the suðNÞ Lie algebra-that is, the algebra of traceless Hermitian N × N matrices. (iv) C abc ¼ −iTrðT a ½T b ; T c Þ are the structure constants of the suðNÞ Lie algebra, with the generators T a normalized as TrðT a T b Þ ¼ δ ab . (v) λ is the t'Hooft coupling constant which is kept fixed when taking the large-N limit. The standard 't Hooft limit is realized by scaling the energy to be of order N 2 as N → ∞. λ has a dimension of ðMassÞ 3 , and without loss of generality we can set it to unity by expressing all dimensionful quantities in units of are the dimensionless Majorana fermionic operators with anticommutation relations fψ a α ;ψ b β g ¼ δ ab δ αβ . The indices α; β; … run from 1 to 16. They correspond to the 16 elements of Weyl- are the d ¼ nine-dimensional analogues of the Pauli matrices, which are traceless, real, symmetric 16 × 16 matrices with anticommutation relations σ i σ j þ σ j σ i ¼ 2δ ij (see Appendix E for explicit construction and useful identities). A nice summary of formulas for the BFSS model can be also found e.g., in [51].
Getting rid of explicit Lie algebra indices and treatingX i andψ α as N × N Hermitian traceless matrices, we can also write the Hamiltonian (3) as where the commutators and traces are understood as operations on N × N matrices, rather than quantummechanical traces, and the t'Hooft coupling λ is already set to unity. The representation (4) makes it obvious that the Hamiltonian (3) is invariant under the simultaneous unitary similarity transformations of all the matrices X i , P i , and ψ a α , which are generated by the operator acting as ½Ĵ a ;Ô b ¼ iC abcÔc on any operatorÔ a which transforms under the adjoint representation of SUðNÞ, e.g., X a i ,P a i ,ψ a α . This symmetry is a remnant of the gauge symmetry of the (9 þ 1)-dimensional super-Yang-Mills theory, from which the BFSS Hamiltonian (3) can be obtained by dimensional reduction. Correspondingly, the physical Hilbert space is defined by imposing the constraint J a jΨi ¼ 0 on physical states.

III. GAUSSIAN STATE APPROXIMATION FOR THE REAL-TIME DYNAMICS OF THE BFSS MODEL
In this section, we explain how the Gaussian state approximation is obtained by truncating the full equations of motion (Heisenberg equations) ∂ tÔ ¼ i½Ĥ;Ô for the canonical coordinate operatorsX a i ,P a i , andψ a α : Averaging these equations of motion over some density matrix, we can express the time derivatives of the expectation values hX a i i and hP a i i in terms of equal-time correlators of up to three operatorsX,P, and/orψ. The equations of motion for these correlators would include correlators with an even larger number of operators, and we would obtain an infinite hierarchy of equations similar to the Schwinger-Dyson equations which cannot be treated either numerically or analytically without further approximations.
In order to obtain a treatable approximation to the full Heisenberg equations (7) which involves only a finite number of variables, let us restrict the time-dependent density matrix hX; ψjρjX 0 ; ψ 0 i to be the most general Gaussian functional of X, X 0 , ψ, and ψ 0 with timedependent parameters [41,43] (where jX; ψi are the eigenstates of the bosonic and fermionic operatorsX a i andψ a α ). In other words, the matrix elements hX; ψjρjX 0 ; ψ 0 i can be represented as N expð−FðX; ψ; X 0 ; ψ 0 ÞÞ, where FðX; ψ; X 0 ; ψ 0 Þ is the most general quadratic polynomial of X, X 0 , ψ, and ψ 0 with a time-dependent normalization factor N . Such Gaussian density matrices can be unambiguously parametrized in terms of one-and two-point correlators of canonical variables due to Wick's theorem. Therefore, by using the Gaussian density matrix, all the equal-time correlators of the canonical variablesX a i ,P a i , andψ a α are expressed in terms of one-point and two-point correlators. Let us introduce the following concise notation: Symmetrization of the product ofX a i andP b j operators in the last definition ensures the real valuedness of the equaltime correlator ⟪X a i P b j ⟫, and it also allows us to work with Wigner functions in a more straightforward way (see below). While it is possible to introduce the mixed bosonicfermionic correlators of the forms hψXi, hψPi, and fermionic one-point functions hψi as well, one can straightforwardly demonstrate that if they vanish in the initial state, they remain zero during all the subsequent evolutions. Since states with nonzero expectation values hψXi, hψPi, and hψi are rather exotic excited states, we restrict our analysis to initial states where only the correlators (8) are nonzero. We note that these correlators can also be put in one-to-one correspondence with the Green functions G þþ ∼ hXðt þ ÞXðt þ Þi, G þ− ∼ hXðt þ ÞXðt − Þi, and G −− ∼ hXðt − ÞXðt − Þi on the Keldysh contour parametrized by time variables t þ and t − on the forward and backward branches, respectively. Throughout the paper we will often refer to X a i and P a i as the classical coordinates and momenta. This interpretation is justified when expð−FðX; ψ; X 0 ; ψ 0 ÞÞ is sufficiently localized. Averaging Eqs. (7) over the Gaussian density matrix characterized by the correlators (8) and applying Wick's theorem, we obtain the following equations for the time evolution of X a i and P a i : To make our approximation self-consistent, we also need to describe the time evolution of the two-point correlators which enter Eqs. (9). To this end let us write down the Heisenberg equations governing the time evolution of the composite operatorsX a iX b We can again average these equations over our Gaussian density matrix and apply Wick's theorem. This is straightforward for all equations except (10c), where one has to express the expectation values of the form hX b jX d iX e jP f k i in terms of two-point functions (8). SinceX andP do not commute, one cannot treat them as ordinary commuting numbers, and the application of Wick's theorem is not straightforward. Indeed, when averaging all other equations we have implicitly used a representation of the density matrix ρ in terms of the eigenstates of eitherX orP operators, which is obviously a Gaussian functional in both cases. Such a representation cannot be used for correlators which contain bothX andP operators. A simple solution to this problem is to use the Wigner transform of the density matrixρ. If hXjρjX 0 i is a Gaussian functional of X and X 0 , the Wigner transform (11) is also a Gaussian functional of X and P. Using the definition (11), one can show that the "classical" phase space integrals of the form R dX dP ρðX; PÞOðXÞP a i are related to the vacuum expectation value of the symmetrized operator product Z dX dP ρðX; PÞOðXÞP a where OðXÞ can be any operator which commutes with all operatorsX a i . Since the left-hand side of Eq. (12) is a Gaussian integral over ordinary commuting variables, we can apply Wick's theorem to the symmetrized operator products such as the ones on the right-hand side of (12).
Since we have assumed that the only nonzero correlator with fermions is hψ a αψ b β i and hψi ¼ 0, hψXi ¼ 0, hψPi ¼ 0, fermionic terms in our Gaussian density matrix completely decouple from the bosonic ones and can safely be disregarded in the above considerations. In fact, we do not even need Wick's theorem for fermions, since correlators with more than two fermionic operators never appear in our equations of motion.
Symmetrizing the operator products in the expectation values hX b jX d iX e jP f k i in (10c) and convoluting all Eqs. (10) with a Gaussian Wigner transform of the density matrixρ, we obtain the following equations for the time evolution of the two-point correlators in (8): Note that equations of motion for the bosonic two-point correlators do not contain fermionic correlators. Fermions only affect the dynamics due to the coupling to the expectation values X a i ≡ hX a i i, which enter Eqs. (13a)-(13c) via the disconnected correlators hX a iX b (9) and (13) provide a full and consistent system of equations for the time evolution of the correlators (8). In particular, one can show that these equations conserve the expectation values of the Hamiltonian (3) and the angular momentum (A1), provided these are also expressed in terms of the correlators (8) using Wick's theorem. Explicit expressions for these conserved quantities are given in Appendix A. On the other hand, supersymmetry generators (6) are not conserved; see Appendix A for a detailed discussion.
The conservation of the generators of the gauge transformations (5) is important for what follows and requires a special discussion. One can show that, similar to the energy and the angular momentum, Eqs. (9) and (13) conserve the expectation values hĴ a i ¼ hΨjĴ a jΨi of the gauge constraint, which we require to vanish in the initial state of our system. The gauge constraintĴ a jΨi ¼ 0 in the full quantum treatment is, however, much stronger and is equivalent to the vanishing of hΨ 0 jĴ a jΨi for an arbitrary state vector jΨ 0 i. It is straightforward to check that there is no normalizable Gaussian wave function jΨi which satisfies the equationĴ a jΨi ¼ 0. The Gaussian state approximation is thus only able to describe the ungauged versions of the BFSS model and the bosonic matrix model. Since the BFSS model is only supersymmetric on the space of gaugeinvariant states, we also conclude that supersymmetry cannot be preserved within the Gaussian state approximation (see Appendix A for a more detailed discussion).
Fortunately, as discussed recently in [39,50], the physics of the bosonic matrix model and the BFSS model does not strongly depend on gauging. More precisely, gauged and ungauged theories are expected to be the same up to the e −C=T correction at low temperature, where C is an order one constant. Their behavior in the high-temperature region is also qualitatively the same; in particular, the real-time aspects in the gauge singlet sector are exactly the same in the high-temperature limit, if the energies are taken to be the same. The description of the ungauged bosonic matrix model within the Gaussian state approximation appears to be rather good, as suggested by the comparison of the thermodynamic equation of state with numerical data of [39] in Sec. IV.
Another important property of Eqs. (9) and (13) is that they evolve pure states into pure states and, more generally, conserve the von Neumann entropy of the density matrixρ (see Appendix C). This allows one to study quantum entanglement between different d.o.f. in a meaningful way; see Sec. V B. Still, one has to keep in mind that the Gaussian state approximation does not describe a unitary evolution. In particular, scalar products between different Gaussian states and expectation values such as hĤ 2 i should be conserved for unitary evolution described by the operator e iĤt , but are not conserved within the Gaussian state approximation. This is because the energy eigenstates are not necessarily Gaussian.
Similar to the full Schrödinger equation ∂ t jΨi ¼ iĤjΨi, which can be obtained by extremizing the "quantum" action S q ¼ R dt hΨjð∂ t − iĤÞjΨi over all possible time histories of a unit state vector jΨi, Eqs. (9) and (13) can be obtained by restricting this extremization to the space of all possible time-dependent Gaussian states [43]. One can also interpret (9) and (13) as classical equations of motion which follow from a certain extension of the classical Hamiltonian [41,43]. This property allows one to identify a symplectic structure of these equations (see Appendix C) and devise stable leapfrog-type numerical integrators.
In contrast to the full Schrödinger equation for the Hamiltonian (3), Eqs. (9) and (13) contain a finite number of variables which scales only polynomially with the number of d.o.f., which allows for an efficient numerical solution even for large physical systems. In particular, this mild scaling is a motivation for using the Gaussian state approximation to study quantum real-time dynamics in quantum chemistry [41,43]. In our case, the most computationally intensive part of the simulations is the solution of Eqs. i i⟪X e jP f k ⟫ on the right-hand sides of (13b) and (13c). The structure of Wick contractions in these terms allows one neither to use the functions which calculate a single commutator nor to save time by contracting some of the spatial indices prior to contracting the matrix indices. Naively, index contractions in these terms require Oðd 3 N 6 Þ floating-point operations. By explicitly taking into account the structure of C abc tensors we have achieved an Oðd 3 N 5 Þ scaling, which is still significantly more dramatic than the Oðd 2 N 3 Þ scaling for the simulations of the classical dynamics, and thus significantly limits the range of accessible N values.
In this work we consider three different approximations to the full real-time dynamics of the BFSS model (3), all of which can be obtained from Eqs. (9) and (13): (1) Classical dynamics of the Hamiltonian (3). The corresponding equations of motion are obtained from Eqs. (9) by setting all two-point correlators to zero. Since there is no classical limit for fermionic dynamics, fermions are completely neglected in this approximation. The classical approximation becomes quantitatively exact for both the bosonic matrix model and the full BFSS model at asymptotically high energies/temperatures. (2) Real-time dynamics of the ungauged bosonic matrix model in the Gaussian state approximation, which corresponds to the Hamiltonian (3) and Eqs. (9) and (13) without fermionic terms. In contrast to the full BFSS model, at low temperatures the bosonic matrix model is expected to be in the confinement regime [52] with finite ground-state energy. While in the ungauged bosonic matrix model there is no strict notion of confinement and the high-and lowtemperature regimes appear to be smoothly connected, at sufficiently low temperatures physical observables in the gauged and in the ungauged models become exponentially close [39,50].

IV. THERMAL INITIAL CONDITIONS AND EQUATION OF STATE
Equations (9) and (13) which approximately describe the real-time dynamics of the BFSS model should still be supplemented with suitable initial conditions. In this paper, we are mostly interested in the real-time responses of thermal and nearly thermal states. Hence we take the initial conditions to reproduce the properties of thermal equilibrium states of the ungauged bosonic matrix model and the ungauged BFSS model as close as possible. In this section we explicitly construct such initial conditions within the Gaussian state approximation.

A. Bosonic matrix model
Within the Gaussian state approximation the thermal density matrix by definition should also be Gaussian (i.e., correspond to a Gaussian Wigner function). If the system is in contact with a thermostat which does not perform work (e.g., collisions with a hard wall), upon thermalization the von Neumann entropy of a state should reach its maximal possible value for a given energy. Based on this very general physical principle, we will approximate the thermal equilibrium states by those Gaussian density matrices which have the largest possible von Neumann entropy at a given energy.
To this end we need to know the von Neumann entropy of an arbitrary Gaussian density matrix, which can be expressed in terms of the correlators (8). This relation has been addressed in detail in [53], and more recently in [37,[54][55][56], where it was demonstrated that the von Neumann entropy of a Gaussian density matrix can be expressed in terms of the so-called symplectic eigenvalues of the block matrix f. in our system. Symplectic eigenvalues of the matrix (14) are related to the eigenvalues of the matrix ΔΩ, where Ω is the symplectic form for the canonical coordinates X a i , P b j : For any positive-definite correlator matrix of the form (14) the eigenvalues λ of the matrix ΔΩ come in complex conjugate pairs of the form λ 2k−1 ¼ þif k , λ 2k ¼ −if k . The real and positive numbers f k , k ¼ 1 Á Á Á N tot are called symplectic eigenvalues of Δ. Quantum uncertainty relations imply that f k ≥ 1=2. A necessary and sufficient condition for the correlator matrix (14) to describe a pure Gaussian state is that f k ¼ 1=2 for all k. It is easy to check that for a single bosonic coordinatex this identity implies ⟪x 2 ⟫⟪p 2 ⟫ − ⟪xp ⟫ 2 ¼ 1=4. In other words, the Heisenberg uncertainty relation should be saturated. In Appendix C we demonstrate that Eqs. (13) conserve symplectic eigenvalues f k and thus map pure states to pure states. The von Neumann entropy S ¼ −Trðρ lnρÞ of a Gaussian state characterized by the correlator matrix (14) can be expressed in terms of symplectic eigenvalues f k as [53][54][55][56] S ¼ As it should be, the von Neumann entropy is equal to zero for pure states and positive for mixed states. In the classical limit, when the f k are large, the von Neumann entropy approaches the classical entropy and can be expanded as Thermal equilibrium states should be invariant under spatial and internal SUðNÞ rotations, as well as under discrete time-reversal and parity transformations. These symmetries imply X a i ¼ 0, P a i ¼ 0 and the following form of the correlators (8): The von Neumann entropy (16) for the Gaussian density matrix characterized by such correlators is given by where f ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi σ xx σ xp p is the dðN 2 − 1Þ-fold degenerate symplectic eigenvalue of the block matrix (14) constructed from correlators (18).
Substituting the correlators (18) into the expression (A2), we also obtain the corresponding energy QUANTUM CHAOS, THERMALIZATION, AND … PHYS. REV. D 99, 046011 (2019) In order to maximize the entropy (19) at fixed energy E, we use Eq. (20) to express σ pp in terms of σ xx and E, which yields Since the entropy (19) is a monotonically increasing function of f, it has a maximum with respect to σ xx when the equation It is now convenient to express physical observables in terms of the equilibrium value of f given by (22), One can immediately check that the correlators (18) with σ xx and σ pp given by (23) provide a time-independent solution of Eqs. (9) and (13), as it should be for thermal equilibrium states. Now the only missing ingredient in our analysis of the equation of state is the temperature, which can be introduced using the standard thermodynamic relation Expressing this derivative in terms of partial derivatives with respect to f, we obtain This equation, together with (23), provides a parametric form of the equation of state and allows one to express energy and other physical quantities such as 1 N hTrX 2 i i in terms of the temperature T. Taking the high-temperature limit which is equivalent to the large-f limit, we reduce Eq. (24) to the form This immediately leads to the high-temperature asymptotics of the equation of state, which is exactly the classical equation of state for the ungauged bosonic matrix model [39,50]. It can also be obtained by replacing the quantum von Neumann entropy (19) with the corresponding classical expression S ¼ dðN 2 − 1Þ lnðfÞ þ const. The relations (23) and (25) also allow one to express the coordinate dispersion at asymptotically high temperatures as In Fig. 1 (plots on the left) we compare our equation of state given by Eqs. (23) and (24) with the numerical results of [39] for the temperature dependence of the energy and the coordinate dispersion 1 i i is such that it remains finite in the 't Hooft large-N limit.
We indeed observe a rather good agreement within a few percent accuracy for the temperature dependence of both the energy and the coordinate dispersion for all simulation parameters used in [39]. It is also interesting to note that the Gaussian state approximation also reproduces very precisely the prediction of [50] for the lowtemperature behavior of the energy of the ungauged model. Expanding Eqs. (24) and (23) to the leading order in f − 1=2, it is easy to obtain the low-temperature asymptotics of the equation of state The coefficients dðd − 1Þ 1=3 ¼ 18 and ðd − 1Þ 1=3 ¼ 2 in (28) match within statistical errors the numerical fit of the form EðTÞ − EðT ¼ 0Þ ¼ Be −C=T in [39], which yields B ¼ 20.0ð2.9Þ and C ¼ 2.043ð76Þ. Such a good agreement can probably be explained by the fact that in the dual holographic picture the difference EðTÞ − EðT ¼ 0Þ at T ≪ 1 is saturated by rather heavy open string excitations, for which the mean-field-like approximation should work rather well. It also suggests that the Gaussian state approximation is not invalidated in the large-N limit. In particular, the ground state energy P. V. BUIVIDOVICH, M. HANADA, and A. SCHÄFER PHYS. REV. D 99, 046011 (2019) in the Gaussian state approximation deviates from the large-N extrapolation of the Monte Carlo results of [39] by 1% only. For comparison, applying the Gaussian state approximation to the one-dimensional anharmonic oscillator with the HamiltonianĤ ¼p 2 þx 4 , one obtains the ground state energy with 2% precision [57]. The fact that the D ¼ ð9 þ 1Þ-dimensional bosonic matrix model is very well described by the Gaussian approximation (which is equivalent to the mean-field approximation) has been previously noticed in [58] and explained in terms of the leading order of a 1=D expansion, which very accurately describes the case of D ¼ 9 þ 1.
There are two possible ways to interpret the Gaussian state characterized by the correlators (18) with σ xp ¼ 0 and σ xx and σ pp given by (23) as initial conditions for the realtime dynamics described by Eqs. (9) and (13).
(i) A trivial way is to directly substitute the correlators (18) into Eqs. (9) and (13) which govern the realtime evolution. It is straightforward to check that for all values of f this simply yields a time-independent solution with no particularly interesting properties.
A conventional stability analysis based on linearization of Eqs. (9) and (13) also shows that this time-independent solution is stable under small perturbations. Up to corrections proportional to 1=N 2 , small oscillations of X a i and P a i around In particular, with such an interpretation we cannot extract any nontrivial Lyapunov exponents and also cannot reproduce the known chaotic behavior of the system in the classical limit. Also since X a i ¼ 0 for this solution, the fermions completely decouple and we cannot capture their influence on real-time dynamics.
(ii) In what follows we use another, physically better motivated option of interpreting mixed Gaussian states with Namely, we represent the coordinate and momentum dispersions as sums of the quantum contributions σ 0 xx and σ 0 pp which saturate the uncertainty relation σ 0 xx σ 0 pp ¼ 1=4, and the classical contributions σ c xx , σ c pp which describe classical thermal fluctuations. For the purely quantum dispersions σ 0 xx and σ 0 pp we use the values (23) with f ¼ 1=2, which correspond to the Gaussian state jΨ 0 i with the lowest possible energy (29). As long as only the total dispersions σ xx ¼ σ 0 xx þ σ c xx and σ pp ¼ σ 0 pp þ σ c pp enter the variational analysis of the equation of state, the choice of σ 0 xx is ambiguous. While for the sake of simplicity we choose the value of σ 0 xx which corresponds to the lowest-energy Gaussian state, in principle one can also make σ 0 xx a temperaturedependent quantity. We then represent the finite-temperature Gaussian density matrix characterized by correlators (18) as a mixture of pure Gaussian states, with random coordinate and momentum displacements X a i and P a i which have Gaussian distributions with dispersions where hi c denotes averaging over the classical probability distribution. We then use Eqs. (9) and (13) to individually evolve each of the randomly shifted pure states jX; Pi in time. Expectation values of physical observables are finally averaged over random initial values of X a i and P a i . As one can see from the upper right plot in Fig. 1, this representation of the initial thermal state of the system yields the correct temperature dependence of the energy with rather small statistical errors. In the lower right plot in Fig. 1 we demonstrate that in this way we reproduce also the correct temperature dependence of the coordinate dispersion 1 N hTrX 2 i i, which, unlike energy, is not conserved and has no reason to stay constant in time. Nevertheless, we find that both the early-time expectation value as well as the time-averaged late-time expectation values of this observable agree very well with the thermal equation of state. These observations justify the interpretation of a mixture of nontrivial time-dependent pure states as a dynamical equilibrium state. In Fig. 2 in Sec. V we also show the full time dependence of 1 N hTrX 2 i i. This interpretation of the classical component of the dispersions of X a i and P a i also allows one to make contact with classically chaotic behavior at high temperatures. Indeed, at high temperatures the classical dispersions will strongly dominate over the quantum ones, and the dynamics described by Eqs. (9) and (13) becomes very close to the classical one. Due to its chaoticity and ergodicity, classical matrix mechanics exhibits real-time thermalization toward a dynamical equilibrium state in which long-time averages of physical observables approach their thermal equilibrium values [25]. This thermalization process can also be interpreted as quasinormal ringing characterized by nontrivial complex-valued quasinormal frequencies [20,59,60], with real parts being quite close to our estimates (30) and (31); see Sec. V C. At the same time, classical matrix mechanics has finite Lyapunov exponents. Thus interpreting thermal states as dynamical equilibrium states we can capture quantum corrections to Lyapunov exponents and imaginary parts of quasinormal frequencies as well as the time evolution of quantum entanglement.
Of course, the two interpretations discussed above would be equivalent for unitary evolution, but yield drastically different results for the nonunitary evolution within the Gaussian state approximation.

B. Full BFSS model
To obtain the equation of state of the full BFSS model we will use the same approach as for the bosonic matrix model and find mixed Gaussian states of fixed energy which maximize the von Neumann entropy. To this end we again split the coordinate and momentum dispersions into the classical and quantum contributions, as in (32) and (34). Since according to Eq. (13d) fermionic d.o.f. only interact with the classical expectation value X a i , we assume that the fermions are initially in the ground state with fixed classical coordinates X a i . This assumption is in line with our construction of the thermal initial conditions where the thermal state was represented by a mixture of ground-state wave functions averaged over random wave packet shifts. Introducing some finite initial temperature for fermions will only increase the energy of our mixed states, which for our approximation is anyway higher than the exact value for the full BFSS model (see topmost left plot in Fig. 1), and thus will not improve our approximation.
Correspondingly, the fermionic contribution E F to the energy only depends on the dispersion σ c xx of the classical wave-function shifts X a i in (34), and it is obtained by averaging the ground-state energy of the fermionic Hamiltonian over a Gaussian ensemble of classical coordinates X a i : A detailed discussion of the spectrum and the ground state of this Hamiltonian is given in Appendix B. Since the only energy scale for the fermionic Hamiltonian is set by the classical X coordinates, dimensional analysis implies that the mean energy of fermions in (36) should scale as . General scaling arguments from random matrix theory fix the scaling of E F with N, which allows one to estimate E F up to an overall universal coefficient as While the N-and σ c xx -independent coefficient A f can be calculated exactly using the methods of random matrix theory, for the purposes of this work we obtain the value of A f numerically by averaging the fermionic energy over a sufficiently large ensemble of randomly generated X a i coordinates and fitting the dependence on σ c xx and N to Eq. (36). These fits work perfectly within statistical errors and yield A f ¼ 15.2661ð34Þ.
Since fermions are assumed to be in the ground state at fixed X a i , by virtue of Nernst's theorem their contribution to von Neumann entropy is zero. Thus in order to obtain the equation of state of the full BFSS model within the Gaussian state approximation, we have to maximize the von Neumann entropy (19) at fixed energy  [39]. The gratings show the period 2π=w XX of lowest-frequency quasinormal oscillations of the coordinate dispersion around the thermal state, with w XX given by (31).
Since the fermionic contribution depends only on the dispersion σ c xx of classical coordinates, the entropy should be maximized with respect to both σ c xx and σ 0 xx , which are now considered as independent parameters. In particular, the separation (32) of the coordinate dispersion into the quantum and classical contributions is no longer ambiguous.
In order to maximize the von Neumann entropy (19) in the space of Gaussian states with fixed energy, we now use Eq. (37) to express σ pp in terms of σ c xx , σ 0 xx , and E, which allows one to express f 2 ¼ σ xx σ pp in (19) as This function does not have a local minimum in the space of σ 0 xx and σ c xx ; that is, the equations ∂f 2 =∂σ 0 xx ¼ 0 and ∂f 2 =∂σ c xx ¼ 0 have no physical solutions. We have to remember, however, that the classical and quantum dispersions have to be all non-negative and satisfy σ 0 xx σ 0 pp ¼ 1=4. This turns the maximization of (38) into a constrained optimization problem, for which the extremum might lie on the boundary of a region allowed by constraints. To describe this region we express σ 0 pp and σ c pp in terms of f 2 , σ 0 xx , and σ c xx as Inserting the explicit expression (38) for f 2 into the above formula for σ c pp , after some algebra we can rewrite the constraint σ c pp > 0 solely in terms of σ 0 xx and σ c xx as The minimal value of the function on the left-hand side of this inequality sets the lowest value of energy at which the constraint can still be satisfied. Numerical minimization yields E > E 0 ¼ 3.9692ðN 2 − 1Þ. Thus E 0 is the ground state energy of the full BFSS model within the Gaussian state approximation. While it is noticeably lower than the value (29) for the bosonic matrix model, supersymmetry of the full BFSS model implies that the true ground state energy should vanish (see e.g., [39,52]). Again we see that supersymmetry cannot be preserved within the Gaussian state approximation (see also Appendix A). To get a lower ground state energy, one needs to include at least the three-point connected correlators of the form ⟪X a iψ b αψ c β ⟫ in the numerical analysis.
We now obtain the equation of state for the full BFSS model within the Gaussian state approximation by maximizing the von Neumann entropy (19) with respect to σ 0 xx and σ c xx within the region specified by the constraint (40) and σ 0 xx > 0 and σ c xx > 0. The value of f is now given by (38). Since f 2 given by (38) has no local maxima, its maximum lies on the boundary of the optimization region, that is, at σ c pp ¼ 0. The corresponding constrained optimization problem cannot be solved exactly, and we use numerical maximization. In this way we obtain σ 0 xx , σ c xx , and the von Neumann entropy (19) as functions of energy E. Using numerical interpolation, differentiation, and functional inversion, we then again use the relation T −1 ¼ ∂S ∂E to introduce the temperature T and express σ 0 xx , σ c xx , and E as functions of T. The resulting equation of state is illustrated in Fig. 1, where we show the temperature dependence of the energy and the coordinate dispersion h 1 N TrX 2 i i. As one can see from the plots on the left in Fig. 1, the agreement with first-principle numerical simulations of [39] is not so good as for the bosonic matrix model. Nevertheless, the Gaussian state approximation correctly captures the following features of the thermal states of the full BFSS model: (i) The ground-state energy of the BFSS model is smaller than that of the bosonic matrix model. β ⟫ are fixed by assuming that the Majorana fermions are in the ground state at fixed coordinate expectation values X a i ; see Appendix B for a more detailed discussion. Simulating the real-time evolution of these pure states, we then average the result over random shifts in the initial conditions. In the upper right plot in Fig. 1 we demonstrate that such averaging correctly reproduces the temperature dependence of the energy. On the other hand, for the coordinate dispersion 1 N hTrX 2 i i the correct temperature dependence is only reproduced by early-time averages and by late-time averages at sufficiently high temperatures. At low temperatures the latetime averages deviate significantly from their thermal values, which might be related to the conjectured real-time instability of the BFSS model with respect to spontaneous emission of D0-branes [26].
As a side remark, let us note that the temperaturedependent energies, coordinate dispersions, and entropies obtained within the Gaussian state approximation satisfy the so-called Bekenstein bound S ≤ 2πER [61,62], with R defined as R ¼ In fact, for all the models which we consider (classical matrix mechanics, the bosonic matrix model, and the full BFSS model) we have S ≪ 2πER.

V. NUMERICAL RESULTS
In this work we numerically solve Eqs. (9) and (13) with initial conditions described in Sec. IV. We use a discretization scheme described in Appendix D. Since the numerical cost of our simulations scales as N 5 , we mostly use a moderately large value N ¼ 5. We have also performed a few simulations with N ¼ 7 to make sure that our results do not change qualitatively at larger N and exhibit the proper 't Hooft scaling.
We average simulation results over several (typically, between five and seven) random initial conditions as previously discussed in Sec. IV. Where shown, error bars on our plots represent the statistical error for such an averaging. Since the number of d.o.f. in our model is sufficiently large, this statistical error is typically very small due to self-averaging, which works well even for a single instance of random initial conditions. In particular, due to high numerical cost we have used only a single instance of random initial conditions for simulations with N ¼ 7; thus the corresponding data points on our plots have no error bars.
To have a first look at the real-time dynamics described by Eqs. (9) and (13) In order to make a meaningful comparison of simulations with characteristic timescales which differ by several orders of magnitude, in Fig. 2 and in other plots in this work we express physical time in units of the classical Lyapunov time τ 0 L , which is defined as the inverse of the leading Lyapunov exponent λ 0 L for the classical dynamics of the BFSS model (3) at a given temperature. For classical matrix mechanics, the temperature dependence of the leading Lyapunov exponent is known to be [25] In Fig. 2 we also separately show the contributions of the classical dispersion hX a i X b j i c and the quantum dispersion ⟪X a iX b j ⟫. Their sum is the physical observable 1 N hTrX 2 i i. For the bosonic matrix model at all temperatures as well as for the BFSS model at high temperatures we observe that within a short time interval t < 2τ 0 L the classical contribution to the coordinate dispersion tends to decrease by a factor of roughly 4. At the same time, the quantum contribution grows in such a way that the total coordinate dispersion remains practically constant in time, up to small short-scale fluctuations with characteristic frequency close to w XX as given by (31).
The decrease of the classical contribution and the corresponding increase of the quantum contribution become particularly large at high temperatures, which indicates a rapid spread of wave functions in configuration space driven by the chaotic dynamics of their centers. In this way the system approaches a state of dynamical equilibrium. However, despite this rearrangement, it turns out that the overall coordinate and momentum dispersions which determine the von Neumann entropy (16) in our simulations remain practically constant in time and exhibit only small short-scale fluctuations. Since the energy is also conserved, we can still assign an approximate value of temperature to this dynamical equilibrium state and use the thermal equations of state derived in the previous Sec. IV.
Of course, one should keep in mind that the Gaussian state approximation probably becomes invalid at late times, more precisely, at a time of the order of several classical Lyapunov times τ 0 L . In particular, at such late times the wave function is already strongly delocalized, and one can expect that non-Gaussian terms in the Schrödinger equation become important. For example, in a simple classically chaotic model with 2 bosonic d.o.f. we have observed a good quantitative agreement with the numerical solution of the Schrödinger equation only up to t ≤ 2τ 0 L [48]. Thus while the values at which the quantum and the classical contributions saturate at late times might be not completely accurate, the enhancement of the quantum dispersion at early times should be a physical feature. For the hightemperature regime of the full BFSS model and also for larger N ¼ 7 the situation appears to be very similar, in particular, the time at which the wave function spreading starts is practically the same.
The late-time expectation value of the coordinate dispersion 1 N hTrX 2 i i deviates significantly from the thermal equation of state only for the low-temperature regime of the full BFSS model, where our approximation is most likely inaccurate due to the violation of supersymmetry. While we cannot say much about the dynamics of the BFSS model in this regime, the growth of coordinate dispersion might be related to the real-time instability of the thermal state of the BFSS model with respect to spontaneous emission of D0-branes due to repulsive fermionic forces, similar to Hawking radiation [26]. expressed in terms of the Poisson bracket fX a i ðtÞ; P b j ð0Þg ¼ ∂X a i ðtÞ ∂X b j ð0Þ , which is replaced by the commutator ½X a i ðtÞ;P b j ð0Þ in quantum theory. Since for most parity-and timereversal-invariant thermal density matrices the expectation value of this commutator vanishes, one typically considers expectation values of the form (1) which involve squared commutators [15][16][17]. In this work we consider the out-oftime-order correlator

A. Lyapunov exponents and the MSS bound
whereρ is the initial Gaussian density matrix at t ¼ 0. This definition is a direct generalization of the Lyapunov distance in X space to quantum theory. A generalization of the out-of-time-order correlator (42) can be used to define the full spectrum of Lyapunov exponents for quantum systems [59]. In order to treat the out-of-time-order correlator (42) within the Gaussian state approximation, we interpret our thermal Gaussian density matrix as a mixture of pure Gaussian states jX; Pi with randomly distributed classical expectation values X a i and P a i , as discussed in Sec. IV. We further interleave the two commutators in (42) with the identity decompositionÎ ¼ R dX 0 dP 0 jX 0 ; P 0 ihX 0 ; P 0 j in terms of the Gaussian states jX 0 ; P 0 i shifted in coordinate and momentum space by X 0 and P 0 , which leads to Trðρ½X a i ðtÞ;P b j ð0Þ 2 Þ¼ Z dX 0 dP 0 ⟪X;Pj½X a i ðtÞ;P b j ð0ÞjX 0 ;P 0 i ×hX 0 ;P 0 j½X a i ðtÞ;P b j ð0ÞjX;P⟫ c : ð43Þ We then represent each of the commutators in terms of the infinitesimal displacement operators at t ¼ 0, hX; Pj½X a i ðtÞ;P b j ð0ÞjX 0 ; P 0 i and similarly for the second commutator in (43). Calculating the expression for the matrix elements of X between the two Gaussian states and inserting it into (44) and (43), one can show that within the Gaussian state approximation the integral over X 0 and P 0 is saturated by the saddle point at X 0 ¼ Xð0Þ, P 0 ¼ Pð0Þ, at which the integrand is just the square of the expectation value of the single commutator hX; Pj½X a i ðtÞ;P b j ð0ÞjX; Pi. Hence we can read off the quantum corrections to Lyapunov exponents from the X-and P-averaged norms of the Lyapunov distance vector δX a i (45) between the X a i coordinates for the two solutions of Eqs. (9) and (13)  Technically this distance is much easier to calculate than the squared commutator (1). In order to calculate the Lyapunov distance as a function of time, we use the first equation in (45) and consider the distance jδX a i j 2 between two solutions of Eqs. (9) and (13) for which the initial values of the classical expectation values X a i differ by a small random vector ϵ a i with jϵ a i j ¼ 10 −5 . To ensure that the Lyapunov distances grow isotropically in configuration space, we consider at least two different vectors ϵ a i . In all of our simulations, Lyapunov distances were growing at the same rate independently of the choice of ϵ a i . In Fig. 3 we show the time dependence of the Lyapunov distances at different temperatures for N ¼ 5, comparing classical dynamics with the real-time dynamics of both the bosonic matrix model and the full BFSS model.
While the classical dynamics is always chaotic and exhibits a well-defined exponential growth of the Lyapunov distance, quantum corrections from the bosonic sector make the dynamics completely nonchaotic at sufficiently low temperatures, with Lyapunov distances which do not exhibit any growth, but rather oscillate with a characteristic frequency w X given by (30) which can be obtained by linearizing Eqs. (9) in the vicinity of thermal time-independent solution (23). This is an expected nonchaotic behavior in the low-temperature regime of the ungauged bosonic matrix model, where gauging becomes unimportant. Correspondingly, all observables approach their values in the conventional gauged bosonic matrix model, for which the physics in the low-temperature confinement phase does not depend on temperature by virtue of the large-N volume reduction [40,63]. Since Lyapunov exponents vanish at zero temperature, at any temperature in the confinement regime of the bosonic matrix model Lyapunov exponents should also be zero up to 1=N corrections. While there is no strict notion of confinement in the ungauged matrix model, one can still expect an exponential suppression of Lyapunov exponents at low temperatures, by analogy with the low-temperature scaling of energy in (28).
In contrast, for the BFSS model Lyapunov exponents remain finite down to the lowest temperature that we consider, and Lyapunov distances exhibit a clear growth. At temperatures up to T ∼ 1 they follow very closely the classical Lyapunov distances in the early evolution period with t ≲ 2τ 0 L . At later times the growth becomes milder but is still noticeable. This behavior is in qualitative agreement with the absence of a confinement (or confinementlike) regime in the full BFSS model [39,52], which is thus expected to be chaotic at all temperatures.
For higher temperatures Lyapunov distances for the bosonic matrix model and the full BFSS model behave in very similar ways. At early times t ≲ 0.5τ 0 L they follow rather precisely the time dependence of the classical Lyapunov distance. Afterwards they exhibit a rather slow exponential growth at a rate which is noticeably smaller than the classical Lyapunov exponent.
We also note that at sufficiently early times, when the exponential growth of Lyapunov distances has not yet fully developed, jδX a i j 2 exhibits signatures of quasinormal ringing, which we will discuss in more detail in Sec. V C.
To summarize, there are no indications that in either the bosonic matrix model or the BFSS model do quantum corrections lead to faster growth of the Lyapunov distance than in the classical matrix mechanics. At least at sufficiently short evolution times this observation should be qualitatively accurate. To quantify all these observations, we extract the leading quantum Lyapunov exponents λ L by fitting the numerical data for lnðjδX a i ðtÞjÞ with a linear function of the form c 0 þ λ L t. In these fits we disregard the transition regime at early times t=τ 0 L ≲ 1, where Lyapunov distances between the centers of two infinitesimally close wave packets are close to the classical Lyapunov distances but do not yet exhibit a clear exponential growth. Obviously, the inclusion of the early evolution period with t < τ 0 L into the fitting would make our estimates of Lyapunov exponents closer to, but definitely not larger than, the results for the classical matrix mechanics. In particular, for the low-temperature regime of the bosonic matrix model the absence of exponential growth is obviously independent of the fitting range. We thus restrict the fitting range to 1 ≤ t=τ 0 L ≤ 8, when the exponential growth has already set in and the fits have good quality. While the Gaussian approximation should receive corrections here, especially toward larger t, qualitative featurese.g., that the Lyapunov exponent decreases due to the quantum effect-should be correct. The temperature dependence of the leading Lyapunov exponents extracted from these fits is illustrated in Fig. 4 for N ¼ 5 and N ¼ 7.
A comparison of data points for N ¼ 5 and N ¼ 7 indicates a proper 't Hooft scaling of Lyapunov exponents.
The temperature dependence λ 0 L ∼ T 1=4 of the classical Lyapunov exponent immediately suggests that its value becomes incompatible with the MSS bound λ L < 2πT at the temperature T ⋆ ¼ ð 0.292−0.42=N 2 2π Þ 4=3 ¼ 0.015, as one can also see from Fig. 4. As discussed already in the seminal paper [17], this is not a contradiction, since at such low temperatures the classical approximation inevitably breaks We plot the results for a single instance of random initial conditions. The gratings show the period of lowest-frequency oscillations 2π=w X of classical coordinates X a i , P a i around the thermal state, with w X given by (30).
down, and quantum effects become important. Our realtime simulations explicitly illustrate this transition between the classical and the quantum regimes. From Fig. 4 we see that at least for the bosonic matrix model quantum corrections indeed decrease the Lyapunov exponents in such a way that they remain well below the MSS bound at all temperatures.
In the full BFSS model the effect of fermions is to remove the low-temperature confinementlike regime, so that the system remains in the deconfinement phase all the way down to zero temperature independently of gauging [21,39,52]. Correspondingly, the system should remain chaotic at all temperatures, and the Lyapunov exponents should also remain finite. In agreement with these expectations, for the full BFSS Hamiltonian the leading Lyapunov exponent is always finite in our real-time simulations. In particular, in the low-temperature regime it is significantly larger than the corresponding value for the bosonic matrix model and tends to approach the corresponding classical value. As we have already discussed, at low temperatures our description of the full BFSS model is probably not accurate enough due to explicitly broken supersymmetry, and we cannot make any strong statements about the validity of the MSS bound. We can only state that our results are compatible with the possibility that the full BFSS model saturates the MSS bound at very low temperatures, similar to the SYK model.

B. Entanglement entropy generation
Quantum entanglement between different d.o.f. in an interacting system provides a quantitative picture of the "scrambling" and spreading of quantum information. Entanglement can be quantified in terms of the entanglement entropy where jΨi is some pure state characterizing the entire system and A and B are the two complementary sets of d.o.f. which define the decomposition of the Hilbert space H of the system into a direct product Quantum-chaotic systems are expected to "scramble" the information contained in the two subsystems by rapidly entangling the states in H A and H B , whereupon the entanglement entropy quickly reaches some maximal saturation value. For finite-dimensional Hilbert spaces this maximal value is the "Haar-scrambled" entanglement entropy [8,34]. Strictly speaking, in gauge theories (of which the BFSS model is a descendant) the splitting of the physical Hilbert space into a direct product H A ⊗ H B is not completely trivial due to gauge constraints [64,65]. This problem, however, is not relevant for us since we work in the ungauged theory [50] which does not impose the gauge constraints on its Hilbert space by definition.
Numerical calculation of the entanglement entropy is typically a rather nontrivial task, especially for real-time evolution of interacting systems. Since the Gaussian state approximation which we use in this paper evolves pure states into pure states (see Appendix C), it also provides a convenient framework for studying quantum entanglement. Entanglement entropy for Gaussian states can be directly calculated in terms of equal-time correlators (8) of canonical variables [37,[53][54][55][56]. Here the basic observation is that tracing out d.o.f. in subsystem B from the Gaussian density matrix jΨihΨj, as in (46), again yields a Gaussian density matrix which is characterized by the same correlators (8), but restricted to canonical variables in subsystem A which describe N dof ≤ N tot d.o.f. The correlator block matrix of the form (14) which corresponds to the reduced density matrixρ A in (46) is thus obtained from the full correlator matrix by removing the rows and columns which correspond to d.o.f. in subsystem B which are being traced out in (46) and hence has the size 2N dof × 2N dof . Being restricted to subsystem A, the correlators (8) in general describe a mixed state with a nonzero von Neumann entropy (16), which is nothing but the entanglement entropy (46). Since in our setup the bosonic and fermionic d.o.f. communicate only via the classical expectation value X a i , they cannot be entangled quantum mechanically, and we only consider the entanglement between the bosonic d.o.f., thus completely tracing out the fermionic Hilbert space.
We now use the prescription sketched above to calculate entanglement entropy for time-dependent correlators (8) obtained by solving Eqs. (9) and (13). In order to make sure that the time dependence of entanglement entropy exhibits universal features independently of the subsystem size, we consider four different choices of the subsystem A in (46) and X a 1 i can be interpreted as the spatial coordinates of two D0-branes, and X a 2 i and X a 3 i parametrize the excitations of strings which join these branes [3,5]. (iv) Three D0-branes with N dof ¼ 81: this case is analogous to the above case of two D0-branes, but the indices a, b take three values each. This gives nine values in total, out of which three values belong to the Cartan subalgebra, and six values correspond to off-diagonal terms of X i and P i matrices. These ways of separating the bosonic d.o.f. of the BFSS model (3) into subsystems A and B are schematically illustrated in Fig. 5.
The splitting of matrix entries of X i coordinates into diagonal entries which are interpreted as D0-brane coordinates and off-diagonal entries which are interpreted as stringy excitations would be particularly obvious for a uðNÞ Lie algebra, for which the simplest basis of the Cartan subalgebra consists of diagonal matrices with only one unit element. The center of uðNÞ algebra, proportional to the unit matrix, corresponds to the center of mass of the system which we set to zero in order to exclude a trivial flat direction which corresponds to the overall motion of the center of mass of the system. This leads to an suðNÞ algebra, for which this splitting is not so straightforward at finite N, as orthonormal Cartan matrices have several nonzero elements and thus correspond to displacements of many D0-branes simultaneously. Here we simply rely on the fact that at sufficiently large N the difference between suðNÞ and uðNÞ becomes negligible, and all but one of the elements of Cartan matrices are suppressed as 1=N. We will also see that the time dependence of entanglement entropy shows universal features which are independent of a particular splitting of the d.o.f. between A and B.
We calculate the entanglement entropy separately for each random initial condition, that is, separately for each pure Gaussian state jX; Pi in the decomposition (34) of the thermal density matrix. For different typical random initial conditions the time dependence of the entanglement entropy appears to be very similar due to self-averaging at sufficiently large N. Thus if we choose to calculate entanglement entropy for the mixed density matrix (34) averaged over various Gaussian states, we simply obtain a practically time-independent constant contribution to the entanglement entropy which is proportional to N dof . We illustrate the time dependence of the entanglement entropy for N ¼ 5 and three different temperatures in Fig. 6, normalizing it by the corresponding number of d.o.f. N dof . Plots on the left and on the right correspond to the bosonic matrix model and the full BFSS model, respectively. We plot the data for some particular random initial conditions without any averaging. In order to control numerical errors due to time discretization, in Fig. 6 we also plot the overall von Neumann entropy (19) of our pure Gaussian state (also divided by N dof ), which should be zero for continuous time evolution of pure Gaussian states, as discussed in Appendix C. The deviation of the von Neumann entropy from zero thus reflects numerical artifacts due to time discretization. It is considerably smaller than the entanglement entropy for all our simulation parameters.
We observe that the entanglement entropy indeed exhibits an expected universal "scrambling" behavior: a roughly linear growth at early times and saturation at late times. Only for the low-temperature regime of the bosonic matrix model does the growth appear to be so slow that we do not see the onset of saturation up to t ¼ 8τ 0 L . It also turns out that if the number of d.o.f. N dof in the entangled subsystem is much smaller than the total number of d.o.f. N tot , entanglement entropy is approximately proportional to N dof at all times, as could be expected for a thermal entropy.
For the bosonic matrix model the saturation value of the entanglement entropy per d.o.f. S A =N dof at N dof ≪ N tot turns out to be quite close to the von Neumann entropy SðTÞ=N tot per d.o.f. for the thermal Gaussian state given by (16). For illustration, we show the values of SðTÞ=N tot in the left plots of Fig. 6 with cyan horizontal lines.
The observation that entanglement entropy per d.o.f. for a single pure state is close to the thermal von Neumann entropy per d.o.f. illustrates a nontrivial relation between real-time thermalization of pure states of a sufficiently large FIG. 5. Splitting of the matrix entries of X μ ij into the subsystems A and B for the calculation of the entanglement entropy. system and the description of thermal states in terms of mixed density matrices, as discussed in [8]. This relation is a quantum analogue of the equivalence of microcanonical and canonical ensembles for chaotic systems, and provides yet another argument in favor of real-time thermalization in our simulations. In this respect entanglement entropy is a convenient measure of the complexity of a pure state, similar to e.g., Husimi-Wehrl entropy [66].
Approximately linear scaling of entanglement entropy with the number of d.o.f. N dof is also in agreement with the findings of [37], where it was demonstrated that for periodically driven harmonic oscillators the entanglement generation rate is proportional to the sum of N dof largest Lyapunov exponents. For a sufficiently large system with dense Lyapunov spectrum and for a sufficiently small N dof this sum should also scale linearly with N dof . When the number of d.o.f. becomes comparable with the maximal value N tot , the ratio S A =N dof becomes smaller. This behavior is again expected from the identification N dof ↔ N tot − N dof which follows from the equality S A ¼ S B . When N dof is small, the entanglement entropy exhibits noticeable fluctuations which are absent at larger N dof due to self-averaging. For the bosonic matrix model the typical time during which the entanglement entropy reaches saturation is between one and three classical Lyapunov times at T ≳ 1.
In the full BFSS model (plots on the right side of Fig. 6) fermions completely change the dynamics at low temperatures and make the time at which the entanglement entropy saturates significantly shorter than for the bosonic matrix model at the same temperature. For the full BFSS model the thermal entropy per d.o.f. appears to be around 50%-100% larger than the saturation value of the entanglement entropy per d.o.f., which is probably related to the fact that our approximation is in general worse for the BFSS model than for the bosonic matrix model. At higher temperatures the effect of fermions gradually becomes smaller, and at T ¼ 5.0 the time evolution of entanglement entropy is already very similar in both the bosonic matrix model and the full BFSS model.
In order to quantify the timescale for the saturation of the entanglement entropy more precisely, we define the entanglement saturation time τ E by fitting the time-dependent entanglement entropy with a function A tanhðt=τ E Þ. The overall normalization constant A in this fit takes care of the subsystem-dependent late-time saturation value of the entanglement entropy and allows for consistent comparison of τ E between simulation results obtained for different temperatures, N dof and N. τ E also sets the characteristic scale for the entropy production rate, where SðTÞ is the thermal entropy (16). By analogy with Lyapunov exponents, we introduce the inverse entanglement saturation time λ E ≡ τ −1 E . In Fig. 7 we illustrate the temperature dependence of λ E for N ¼ 5 and N ¼ 7 and compare it with the temperature dependence of the classical Lyapunov exponent. Already from Fig. 6 one can see that the entanglement saturation time is practically independent of the number of d.o.f. N dof in subsystem A. For this reason, Fig. 7 only shows the entanglement saturation time for a single D0-brane in order not to clutter the plot. On the other hand, interpreting the entanglement saturation time as a scrambling time in the sense of [17], one could expect a mild logarithmic growth τ E ∼ lnðN dof Þ of τ E with N dof . This observation could mean either that our system is too small for the scaling to be observed or that entanglement saturation time cannot be identified with the scrambling time. The latter interpretation would make sense at least near the classical limit, because the growth of the coarse-grained entropy is characterized by the Kolmogorov-Sinai entropy, which is roughly proportional to the system size.
We also find that in the high-temperature regime for both the bosonic matrix model and the full BFSS model the entanglement saturation time τ E ≡ λ −1 E is very close to the classical Lyapunov time τ 0 L ≡ λ 0 L . This behavior is in sharp contrast to our findings for the quantum Lyapunov exponents, which were extracted from the exponential growth of Lyapunov distances at relatively late times t ≳ τ 0 L . On the other hand, entanglement saturation time probes the earlytime dynamics at t ≲ τ 0 L , where Lyapunov distances for the quantum system follow rather closely their classical counterparts. These results suggest that early-time thermalization at high temperatures takes place before the quantum effect becomes important, so that the classical treatment is justified.
As we approach the low-temperature regime, the entanglement saturation time for the bosonic matrix model quickly decreases and becomes significantly smaller than the classical Lyapunov exponent, in agreement with the nonchaotic nature of the low-temperature regime. On the other hand, for the full BFSS model the temperature dependence of λ E is far less trivial. In particular, around T ∼ 1 it exhibits a rather pronounced growth and becomes almost an order of magnitude larger than the classical Lyapunov exponent, thus indicating an extremely fast generation of entanglement. This behavior is most likely an artifact of our approximation, as the BFSS model is not expected to undergo any finite-temperature phase transition [39,52]. In particular, in the full supersymmetric BFSS model fermions prevent the onset of confinement [39,52] which is observed in the purely bosonic model exactly around T ∼ 1 [40]. Fast entanglement generation which we observe in our simulation of the BFSS model might be interpreted as a kind of "silver blaze" phenomenon. Namely, for the full quantum dynamics of the BFSS model fermions conspire to completely remove the signatures of confinement-deconfinement transition from physical observables. However, since in our approximation supersymmetry is not conserved, the influence of fermions on the entanglement dynamics is most likely overestimated, which leads to the artificial speed-up of entanglement generation around the to-be deconfinement temperature of the bosonic matrix model. In this respect this speed-up would be consistent with the general expectation that the near-critical regime of quantum systems is more chaotic and exhibits faster scrambling [8,67]. At even lower temperatures, λ E for the BFSS model again approaches the classical Lyapunov exponent, but never decreases below it. Finally, in Fig. 8 we compare the temperature dependence of the entanglement saturation time τ E ≡ λ −1 E with that of the characteristic Lyapunov time τ L ≡ λ −1 L , as well as with the MSS bound. The most notable feature is that the entanglement saturation time is always shorter than the characteristic Lyapunov time. These two characteristic timescales become very close to each other and also to the classical Lyapunov time only in the low-temperature regime of the full BFSS model.

C. Quasinormal ringing and quasinormal frequencies
While Lyapunov exponents and entanglement generation define the characteristic times for the onset of chaos and spreading of quantum information, the diffusion-driven approach to thermal equilibrium is characterized by another timescale τ D set by the decay rate of the quasinormal ringing and thus is related to the imaginary parts of quasinormal frequencies.
In our numerical study of Lyapunov distances we drive the system out of the dynamical equilibrium state by introducing a small coordinate shift ϵ a i ; see Eq. (45). When the ringing decays, the system again reaches dynamical equilibrium, but with a slightly different state within the microcanonical ensemble than the initial one. At later times, these states become exponentially far separated, and the exponential growth of the Lyapunov distance sets in. Since the timescale of quasinormal ringing appears to be shorter than the Lyapunov time, we can observe a rather clear decay of quasinormal ringing in the time dependence of Lyapunov distances, as one can see already from Fig. 3. On Fig. 9 we present a better illustration of the quasinormal ringing by showing the time dependence of Lyapunov distances at early times t < τ 0 L only. The ringing behavior is most clear for the bosonic matrix model. For the full BFSS model and especially for the classical dynamics the onset of exponential growth is sometimes so fast that the ringing cannot be clearly identified.
Nevertheless, it is interesting to study the temperature behavior of quasinormal frequencies at least at the qualitative level. To this end we first identify the duration of the ringing as the period of time during which each successive maximum of the Lyapunov distance jδX a i j 2 ðtÞ is smaller than the previous one. The mean distance between successive extrema t ⋆ i of jδX a i j 2 ðtÞ within this time interval is then used to estimate the real part of the quasinormal frequency as ReðwÞ ¼ π=Δt. We use a simple arithmetic mean over all the extrema which lie within the time interval where the ringing can be identified. To obtain the imaginary part of the quasinormal frequency, we consider the differences Δy i ¼ jδX a i j 2 ðt ⋆ iþ1 Þ − jδX a i j 2 ðt ⋆ i Þ and estimate where the mean is again an arithmetic mean over all successive extrema for which Δy i shows the characteristic decaying behavior. This prescription for estimating quasinormal frequencies yields exact results for the quasinormal ringing of the form whereÔ is some physical observable and hÔi 0 is its thermal expectation value. On Fig. 9 we also show the fits of the form (49) as dotted lines, for those datasets for which such fits can be obtained using the above prescription. In practice this means that we can identify a sufficiently large number of consecutively decaying extrema. For a reliable identification of quasinormal frequencies with the smallest imaginary parts (which thus correspond to quasinormal modes with the longest decay time) one needs to analyze the quasinormal ringing at sufficiently late times. For our simulation setup such late times cannot be reached for most We also sketch the MSS bound λ L < 2πT as well as the classical Lyapunov exponent given by (41). datasets due to the onset of exponential Lyapunov growth, and thus our estimates (48) of ImðwÞ are most likely biased toward somewhat larger values.
We show our estimates for the real and imaginary parts of the quasinormal frequency w X associated with the decay of perturbations of X a i in the bosonic matrix model and in the BFSS model with N ¼ 5 in Fig. 10. For the bosonic matrix model we show the mean values of quasinormal frequencies and the corresponding statistical errors, which can be estimated due to relatively good quality of fits. For the full BFSS model the fits are less reliable, and the standard statistical error does not reflect the fitting uncertainty well. For this case, we thus present only scatter plots of our estimates of ReðwÞ and ImðwÞ for different random initial conditions.
The temperature dependence of the real part of w X for the bosonic matrix model is very well approximated by our estimate (30) obtained by linearizing Eqs. (9) in the vicinity of the thermal state (solid line in the left plot in Fig. 10). To be more precise, due to the squaring of δX a i in our definition of Lyapunov distance we observe oscillations at frequency which is exactly 2 times larger than in (30); thus in Fig. 10 we plot 2w X .
For the full BFSS model the temperature dependence of w X cannot be estimated analytically, as linearization of Eqs. (9) and (13) in the vicinity of the maximally symmetric solution (18) with X a i ¼ 0 and P a i ¼ 0 would yield exactly the same result as for the bosonic matrix model. To see the difference we should switch to the interpretation of thermal states in terms of dynamical equilibrium states, where the dynamics cannot be treated analytically. Nevertheless, from the left plot in Fig. 10 we see that the temperature dependence of Reðw X Þ for the BFSS model is quite similar to the one for the bosonic matrix model.
The imaginary part of w X at sufficiently high temperatures turns out to be noticeably larger than the classical Lyapunov exponent. For the bosonic matrix model Imðw X Þ quickly approaches zero at low temperatures, following basically the same trend as for the Lyapunov exponents and the entanglement generation rate. This behavior is in agreement with the general expectation that confining gauge theories are less dissipative, as exemplified, for instance, by the temperature dependence of electric conductivity in QCD [68].
In contrast, for the full BFSS model imaginary parts of quasinormal frequencies are much larger and seem to remain finite even at the lowest temperatures that we consider, which is again in agreement with the absence of a confinement regime in the BFSS model. While at very high temperatures we cannot reliably extract w X from the data, the plot for T ¼ 50.0 in Fig. 9 makes it clear that at this temperature the quasinormal ringing is very similar for the classical matrix mechanics, the bosonic matrix model, and the BFSS model, in agreement with the classicality of the high-temperature limit. The plots in Fig. 9 also suggest that Imðw X Þ for the classical theory is still larger than for the BFSS model, except probably for the lowest temperatures. For the bosonic matrix model Imðw X Þ takes the smallest values.
Finally, let us note that at high temperatures we can use Eqs. (30), (31), (25), and (23) to estimate the real part of the quasinormal frequency of the quasinormal ringing of TrX 2 i as w XX ¼ 4.89T 1=4 , which agrees well with the dependence w XX ¼ 5.152T 1=4 obtained in [60,69] for the bosonic matrix model using more elaborate dedicated simulations. At the highest temperatures where Imðw X Þ can still be estimated in our simulations, we obtain Imðw X Þ=Reðw X Þ ∼ 0.1, which also roughly agrees with the result of [60,69]. According to the dual gravity calculation, near zero temperature the quasinormal frequency is proportional to T and the ratio Imðw X Þ=Reðw X Þ is about 1.7 [70]. While our calculation is too crude to capture the result near T ¼ 0, the observed growth of this ratio follows the right trend.

VI. CONCLUSIONS AND OUTLOOK
Until recently the real-time dynamics of the BFSS model could only be addressed either in the low-energy regime, which is tractable in terms of the dual holographic description, or in the high-energy regime in which the system becomes effectively classical. Our simulations of the corresponding ungauged models [39,50] within the Gaussian state approximation is a step toward bridging the gap between these two regimes.
The Gaussian state approximation appears to be much more accurate for the bosonic matrix model than for the BFSS model, as suggested by the quantitatively good agreement of our equation of state (23) with the results of first-principle Monte Carlo simulations [39]. Our results for the bosonic model are thus likely to be reliable at the quantitative level, while for the BFSS model we can make at most qualitative statements.
We have explicitly studied and confirmed some of the important features of the real-time dynamics of the bosonic matrix model and the BFSS model which fit expectations either based on general grounds [8] or motivated by the dual holographic description [4,5]: (i) Quantum real-time dynamics is characterized by smaller Lyapunov exponents in comparison with the classical system at the same energy. This ensures the validity of the MSS bound [17] for all energies. (ii) The gauged bosonic matrix model becomes confining and nonchaotic at low temperatures [40,63], with Lyapunov exponents being 1=N suppressed. At low temperature the gauged and ungauged model should become exponentially close [39,50], and hence the Lyapunov exponent of the ungauged model should be exponentially small. We indeed observe such a suppression of Lyapunov exponents in our simulations. Decay time of quasinormal modes and the saturation time for entanglement entropy also become very long. (iii) When fermionic d.o.f. are added, the BFSS model exhibits chaotic behavior and fast decay of quasinormal modes at all temperatures, in agreement with the absence of a confining regime all the way down to zero temperature [39,52]. The Lyapunov exponents, however, still appear to be smaller than in the classical theory at all temperatures. (iv) Entanglement entropy shows the expected scrambling behavior [8]: the initial growth followed by  illustrates an equivalence between microcanonical and canonical ensembles for the quantum-chaotic system [8]. (v) The characteristic time τ E at which the entanglement entropy reaches saturation is shorter than the quantum Lyapunov time τ L ≡ λ −1 L (see the discussion below for a possible loophole) and is also much closer to the classical Lyapunov time at high temperatures. In the intermediate-and low-temperature regimes of the BFSS model, τ E appears to be close to or even shorter than the classical Lyapunov time, in agreement with the chaoticity of the BFSS model at all temperatures. Nevertheless, τ E is still longer than the decay time τ D of quasinormal modes. Figures 8 and 10 provide a compact illustration of these findings.
Our results for the time dependence of the entanglement entropy are under the best theoretical control, as they are extracted from the early-time behavior for which the Gaussian state approximation should be quantitatively accurate. At the same time, the relatively smooth time dependence allows for a rather unambiguous definition of the entanglement saturation time τ E . Quasinormal frequencies are also extracted from early-time behavior; however, the extraction of the decay rate of quasinormal modes is not so reliable and unambiguous. On the other hand, Lyapunov distances are extracted from the behavior of Lyapunov distances at relatively late times, where the Gaussian state approximation can be accurate at most qualitatively. Thus while the exponential fits to the time dependence of Lyapunov distances are quite good, our estimates of Lyapunov exponents can still be biased, especially at high temperatures.
Our real-time analysis might be further improved if one considers correlators with more than two canonical variables, e.g., within the n-particle irreducible effective action techniques [71]. While for higher-dimensional gauge theories the application of such methods faces the difficulty of parametrizing and storing higher-order correlators which depend on many momenta, low dimensionality and numerous (super)symmetries of the BFSS model (3) should significantly reduce the number of independent correlation functions which enter the analysis.
In fact, we have already tried to extend Eqs. (9) and (13) by incorporating the connected correlator ⟪ψ a αψ b βX c i ⟫ and requiring the solution to be SUðNÞ and rotationally symmetric. Unfortunately, it turned out that such an extension is not intrinsically consistent in that the correlators ⟪X a iX b j ⟫ cease to be positive definite after a relatively short evolution time. This indicates the necessity to include correlators of even higher orders into the analysis.
It would also be interesting to understand the applicability of the Gaussian state approximation to real-time dynamics of higher-dimensional gauge theories, as well as its interpretation in the context of the conventional scale separation between soft-momenta and hard-momenta gauge fields which allows one to treat the dynamics of soft gauge fields classically. In principle, Gaussian state approximation should extend the range of validity of realtime simulations of the Yang-Mills theory beyond that of the classical dynamics at a numerical cost which scales quadratically with spatial volume, which is thus comparable with the numerical cost of real-time simulations with fermions [14,72]. In particular, in contrast to purely classical dynamics the Gaussian state approximation should be able to incorporate transitions between different topological sectors, and thus might provide an alternative to real-time evolution equations which include Langevin noise to account for the contribution of hard gauge fields [12,13], in particular for problems such as the estimation of the sphaleron rate [73].
Here we give explicit expressions for these conserved quantities: An explicit demonstration of the conservation of these quantities from Eqs. (9) and (13) is a lengthy but straightforward calculation, which we do not present here in order to save space. The conservation of the supersymmetry generatorsQ α in the Gaussian state approximation is a more subtle question. On the one hand, the expectation values hQ α i are zero for our Gaussian state with hψi ¼ 0 and vanishing mixed bosonic-fermionic correlators, and are hence trivially conserved. On the other hand, if the Gaussian state approximation preserved supersymmetry, we should have obtained a zero ground state energy. Since this is not the case, supersymmetry should somehow be violated.
To understand the fate of supersymmetry in more detail, let us first present an outline of the proof of the conservation of supersymmetry generators (6) from the full Heisenberg equations (7), which after some algebra yield Taking into account the vanishing of the gauge constraint (5) on the physical Hilbert space, we can also rewrite (A5) as Now one can use the anticommutativity of the operatorŝ ψ a α and the cyclic symmetry of the structure constants C abc ¼ C bca ¼ C cab to transform this expression into In other words, the indices α, β, and γ are cyclically permuted. But precisely this combination of σ matrices is equal to zero by virtue of the Fierz identity (E3). Hence the time derivative ∂ tQα is equal to zero on the physical Hilbert space. Now we turn to the Gaussian state approximation and take the limit in which the classical expectation values X a i , P a i are so large that their quantum dispersions can be neglected. In this case the evolution of the system is described by the classical equations of motion for X a i and P a i , augmented by the fermionic force, and fermions evolve in the classical background of X a i variables. This corresponds to the CSFT, which is a standard method for addressing the real-time dynamics of fermions interacting with strong gauge fields; see e.g., [14]. Since the fermionic part of the BFSS Hamiltonian (3) is quadratic in the fermionic fields, the fermionic wave function is always Gaussian and is in one-to-one correspondence with the fermionic correlators ⟪ψ a αψ b β ⟫. In this case we can simplify calculations by working directly with the operatorsψ a α which satisfy the Heisenberg equations of motion. In this limit the supersymmetry generators (6) take the form Using equations of motion (9), after some algebra we can represent the time derivative ofQ α as β ⟫ψ c γ ðσ i αβ σ i γδ − δ αβ δ γδ Þ: ðA9Þ The only difference with the expression (A5) for the full quantum evolution is that two out of threeψ operators are now under the vacuum expectation value brackets. This vacuum expectation value is in fact the expectation value of the fermionic force term i 2 C bac σ i αβ ⟪ψ b αψ c β ⟫ on the righthand side of Eq. (9b) which governs the time evolution of P a i . It is now easy to see that due to these brackets, one can no longer cyclically permute differentψ operators and use the Fierz identities (E3) as in (A7) to show the conservation of the SUSY charge. One can construct a similar proof by considering arbitrary Gaussian states with hψ a α i ≠ 0 and hψ a αX b i i ≠ 0, hψ a αP b i i ≠ 0 which mix fermionic and bosonic variables. This makes the derivation significantly more involved, but leads to similar conclusions. We thus conclude that supersymmetry is not preserved in the Gaussian state approximation.
Having fixed the initial conditions for the fermionic correlators ⟪ψ a αψ b β ⟫, we have to solve Eq. (13d) which governs the time evolution of this correlator. This can be done by promoting the basis vectors u a α ðϵÞ and v a α ðϵÞ to time-dependent functions which satisfy the single-particle Schrödinger equations At t ¼ 0, u a α ðϵ k Þ and v a α ðϵ k Þ are defined by (B3). Since the single-particle Hamiltonian (B2) is time dependent, at t > 0, u a α ðϵ k Þ and v a α ðϵ k Þ are in general no longer related to the eigenstates of h αβ ab . It is easy to check that the twopoint function (B11) with time-dependent functions u a α ðϵ k Þ and v a α ðϵ k Þ which satisfy Eqs. (B12) solves equation (13d). In the literature on real-time simulations within the CSFT approximation the functions u a α ðϵ k Þ and v a α ðϵ k Þ are commonly referred to as mode functions [14,72].
This way of solving Eq. (13d) requires at least 2 times less memory and CPU time than a straightforward solution. For illustration, we present the equation of motion for the bosonic momenta P a i in terms of mode functions:

APPENDIX C: SYMPLECTIC STRUCTURE IN THE GAUSSIAN STATE APPROXIMATION
As discussed in Sec. V B, a necessary and sufficient condition for a general Gaussian density matrix to describe a pure state is that the corresponding correlator matrix (14) should have symplectic eigenvalues all equal to f k ¼ 1=2.
In this Appendix we demonstrate that Eqs. (9) and (13) conserve symplectic eigenvalues of the correlator matrix (14) and hence evolve pure states into pure states. For mixed states this property obviously implies the conservation of von Neumann and Rényi entropies.
To begin with, we introduce the condensed index notation A ¼ fa; ig, B ¼ fb; jg and rewrite Eqs. (13) as where V ABCD is the shorthand notation for the coefficients of the quartic term in the Hamiltonian (3), which thus takes the formĤ ¼P 2 A =2 þ V ABCDXAXBXCXD =4. It is now convenient to introduce the symmetric and real matrix V AB ¼ 3V ABCD hX CXD i and to treat the correlators ⟪X AXB ⟫, ⟪X APB ⟫, and ⟪P APB ⟫ as matrices, omitting their indices. This allows one to write equations (C1) in a particularly simple form: We now combine all the correlators in the block matrix Δ given by (14) and introduce the block matrix In terms of the symplectic form Ω defined in (15) and the matrices Δ and ϒ, Eqs. (C2) can be written as The commutator structure on the right-hand side implies that the eigenvalues of the matrix ΔΩ are conserved during the evolution described by Eqs. (C1), which are a compact representation of Eqs. (13). Since the eigenvalues of ΔΩ are the symplectic eigenvalues of Δ which determine whether the state is pure, we have thus proven that pure states are evolved into pure states. While the form of Eqs. (C1) is only valid for models with quartic interactions, our proof can easily be generalized for other Hamiltonians. In particular, the effect of Majorana fermions can easily be incorporated as a time-dependent linear potential in addition to the quartic potential. In order to solve Eqs. (9) numerically, we employ the leap-frog discretization [14,72] with time step δt, which has the advantage of being numerically stable and phase space volume preserving. We enumerate the discrete steps of the numerical evolution by the discrete variable τ ¼ 0; 1; 2; …. In each evolution step we first update the bosonic momenta, where we have used the matrix notation as in (4) for shortness. Next, the fermionic mode functions are updated: Evolution equations for v α take the same form. This evolution scheme corresponds to using the symmetric discretization of the time derivative on the real time axis, which admits also the propagation of lattice doubler modes.
Since these doubler modes correspond to Majorana-Weyl fermions of opposite helicity, together with the physical modes they will form (9 þ 1)-dimensional Dirac fermions. This enhancement of the fermionic Hilbert space could potentially spoil some of the nice properties of the BFSS model; thus it is very important to prevent the doubler modes from being excited. The second equation in (D2) ensures that the doubler modes remain practically unexcited for sufficiently long evolution time [72]. After that we update the bosonic correlators ⟪X a iP b j ⟫ and ⟪P a iP b j ⟫ according to Eqs. (13b): where we have used the shorthand notations of Eq. (C1). Finally, the variables X a i and ⟪X a iX b j ⟫ are updated: The commutator representation (C4) of Eqs. (13) allows one to devise a better discretization which would involve symmetrized time derivatives, similar to discrete equations (D2) for fermionic mode functions. We have not yet implemented this option into our simulations, and we simply achieve a comparable precision by using smaller discrete time step δt.
All the data presented in the main text of the paper were obtained with δt ¼ 2 × 10 −5 =σ xx . Such a rescaling with respect to σ xx is necessary, since at different σ xx the dynamics is characterized by very different timescales which also require different discretization steps to achieve the same accuracy. In order to check the effect of discretization artifacts on our simulations, we have also performed several simulations with 2 times smaller time step δt ¼ 10 −5 =σ and checked that these simulations yield practically the same results for all parameter sets that we have used.

APPENDIX E: PAULI MATRICES IN d = 9 SPATIAL DIMENSIONS
We use the following explicit form of the σ i matrices: where s 1 , s 2 , and s 3 are the conventional 2 × 2 Pauli matrices. We note that despite s 2 being complex and antisymmetric, it always enters the σ matrices twice. Thus, the σ matrices are manifestly real and symmetric.
To demonstrate the conservation of the angular momentum (A1), one also needs the commutation relations between σ matrices and their commutators σ ij ≡ σ i σ j − σ j σ i (which can also be interpreted as the generators of rotations in the space of Majorana-Weyl spinors): For the proof of the conservation of supersymmetry generators (6) outlined in Appendix A one also needs the Fierz identity ½σ i αβ ½σ i γδ þ ½σ i αγ ½σ i βδ þ ½σ i αδ ½σ i γβ ¼ δ αβ δ γδ þ δ αγ δ βδ þ δ αδ δ γβ ; ðE3Þ as well as the following identity for the σ ij matrices: