Real-time dynamics of matrix quantum mechanics beyond the classical approximation

We describe a numerical method which allows us to go beyond the classical approximation for the real-time dynamics of many-body systems by approximating the many-body Wigner function by the most general Gaussian function with time-dependent mean and dispersion. On a simple example of a classically chaotic system with two degrees of freedom we demonstrate that this Gaussian state approximation is accurate for significantly smaller field strengths and longer times than the classical one. Applying this approximation to matrix quantum mechanics, we demonstrate that the quantum Lyapunov exponents are in general smaller than their classical counterparts, and even seem to vanish below some temperature. This behavior resembles the finite-temperature phase transition which was found for this system in Monte-Carlo simulations, and ensures that the system does not violate the Maldacena-Shenker-Stanford bound $\lambda_L<2 \pi T$, while the classical dynamics inevitably breaks the bound.


Introduction
Thermalization of strongly interacting quantum systems is one of the important problems in different areas of modern physics, ranging from apparent thermalization of the quark-gluon plasma [1][2][3] in heavy-ion collisions to inflation of our Universe and dynamics of ultra-cold quantum gases [4]. Unfortunately, the number of tools which can be used to study non-equilibrium real-time dynamics of generic non-perturbative, non-supersymmetric and non-integrable quantum field theories (in particular, non-Abelian gauge theories) is rather limited. On the one hand, in the regime of large field strengths (large quantum occupation numbers) and small coupling constants one can reliably use the classical equations of motion [5,6], re-summing secular divergences by averaging over quantum fluctuations in initial conditions [5]. On the other hand, in the dilute plasma regime with small quantum occupation numbers one can use the kinetic theory description [7]. The intermediate regime between Speaker. e-mail: pavel.buividovich@physik.uni-regensburg.de. This work was supported by the S. Kowalevskaja award from the A. von Humboldt foundation. the two descriptions with occupation numbers of order of one is very important for matching earlytime strongly non-equilibrium evolution with late-time hydrodynamic behavior [1,2].
In these Proceedings, we aim to go deeper into this intermediate regime of neither large nor small occupation numbers. We outline an extension of the classical dynamics approximation which incorporates sub-leading effects with respect to large occupation numbers by approximately taking into account the quantum dispersion of the wave function of the system. In a few words, the basic idea is to approximate the time-dependent Wigner function by the most general Gaussian function with time-dependent mean and dispersion. This approximation is also used to study short-scale real-time processes in quantum chemistry [8]. However, to our knowledge, it has not been previously used in the context of high-energy physics.
As the simplest prototypical system which features thermalization and (at least classical) chaos and still has non-Abelian structure similar to that of Yang-Mills fields [9], in these Proceedings we consider the matrix mechanics with the Hamiltonian where X i and P i are the canonically conjugate coordinates and momenta which take values in the Lie algebra of S U (N) group. This Hamiltonian arises as a result of compactification of pure continuous Yang-Mills theory from (d + 1) dimensions down to (0 + 1) dimensions. The dynamics described by the Hamiltonian (1) is known to be classically chaotic, so that the distance between initially very close points in phase space grows exponentially with time [9]. The fact that such chaotic classical systems effectively forget about initial conditions after some "thermalization" time can be interpreted as the formation of a black hole state [3,[10][11][12] in the framework of holographic duality between compactified super-Yang-Mills theory and the gravitationally interacting system of D0 branes [13].
However, most of the previous studies of the real-time dynamics of Yang-Mills-type Hamiltonians similar to (1) were based on the classical mechanics approximation, which is only justifiable at sufficiently high temperatures. In this work we use the Gaussian state approximation of [8] to understand how quantum effects might affect the classically chaotic dynamics of the Hamiltonian (1).

Gaussian state approximation: a simple example with classical chaos
Before studying the quantum dynamics for the Hamiltonian (1), in this Section we illustrate the Gaussian state approximation on the example of a simple Hamiltonian This is one of the simplest reductions of the Yang-Mills-type Hamiltonian (1) which still features classical chaos [9]. We start our derivation of quantum corrections to the classical dynamics with the Hamiltonian (2) from the Heisenberg equations for the canonically conjugate operatorsx,ŷ andp x ,p y : We can obtain the equations of motion for the expectation values p x , p y , x , ŷ by averaging equations (3) with respect to some density matrixρ. These equations will contain expectation values of the form xŷ 2 , which should be evolved according to yet another equation, including in turn expectation values of yet larger number of coordinate/momentum operators.
Without any simplifying assumptions we thus get an infinite hierarchy of equations, for which any practical solution is hardly possible. To truncate this infinite set of equations, we approximate the density matrixρ by the most general time-dependent Gaussian function, so that expectation values of products of multiple coordinate and momentum operators can be expressed in terms of the coordinates of the wave packet centers and wave packet dispersion using the Wick's theorem. It is actually more convenient to characterize the most general Gaussian state by the Wigner function where ξ = x, y, p x , p y is the four-component vector of the classical phase space variables, the parametersξ a ≡ ξ a describe the center of the Gaussian wave packet and the matrix characterizes the dispersion of the wave packet in both coordinate and momentum space.
We can now average the Heisenberg equations (3) over our Gaussian state with the Wigner function (4). Using Wick's theorem, we obtain In order to obtain the evolution equations for the dispersions ξ a ξ b , we again use the quantum Heisenberg equations of motion to express the time derivatives of the operator productsξ aξb as Similar equations for other combinations ofx,ŷ andp x ,p y can be obtained by straightforward permutationsx ↔ŷ,p x ↔p y . Averaging the equations (7) over our Gaussian state with the Wigner function (4), we obtain the evolution equations for expectation values ξ aξb of two canonical operators, and, after subtracting the disconnected contributions ξ a ξ b , for the dispersions ξ aξb . The only detail which one has to remember is that the "classical" expectation values Dξρ (ξ) ξ a ξ b ξ c . . . with the Wigner function ρ (ξ) correspond to quantum expectation values Tr ρ ξ aξbξc . . . s of the symmetrized operator products, which can be recursively defined as Applying this procedure to equations (7), we obtain : Equations with other combinations of x, y and p x , p y can be obtained by interchanging x ↔ y. The equations (6) and (9) form the complete system of equations for the time evolution of the coordinatesξ a = x, y, p x , p y of the wave packet center and the dispersion matrix ξ a ξ b . In contrast to the full Schrödinger equation, the number of variables in the equations (6) and (9) grows only quadratically in the number of degrees of freedom.
Let us now compare the numerical solution of the equations (6) and (9) with the solution of the full Schrödinger equation. We consider the pure Gaussian state with initially nonzero expectation values of x and y and the minimal quantum dispersion saturating the uncertainty relation: The variable f controls the dominance of the "classical" expectation values x and y over the quantum dispersions.
To solve the time-dependent Schrödinger equation, we replace the continuum coordinates x and y by a finite lattice with spacing a, with sites labelled by indices for both x and y coordinates. We then perform the leapfrog-type evolution by multiplying the wave function ψ (x, y) → ψ i j by the unitary evolution operator with sufficiently small time step δt. In order to control discretization and finite-volume artifacts, in addition to the solution with lattice spacing a, time step δt and lattice size N s we also consider solutions with parameters {2a, δt, N s }, {a, 2δt, N s } and {a, δt, N s /2}. We then estimate the discretization and finite volume error of the expectation values Ô (t) as the difference between their minimal and maximal values over simulations with these four sets of parameters.
In Fig. 1 we compare the time dependence of the expectation value x (t) of the x coordinate, obtained from the classical equations of motion, from the Gaussian state approximation, and from the numerical solution of the Schrödinger equation. From the topmost plot on the left to the lowest plot on the right we increase the initial expectation values of x and y coordinates (factor f in (10)) as compared to the quantum dispersions which are fixed for all f . We thus move from the quantum regime at small f to the classical regime at large f . We note that as f becomes larger, the classical solution is rescaled as Finally, an important methodological question is how to control the validity of the Gaussian state approximation without having the reference solution of the Schrödinger equation, which is practically unfeasible for sufficiently large number of degrees of freedom. One of the possible intuitive criteria is the importance of terms which contain squares of quantum dispersions in the potential energy V (x, y) = x 2 y 2 /2 = x 2 y 2 + x 2 y 2 + y 2 x 2 + 4xy xy + 2 xy 2 + x 2 y 2 . We have found that the Gaussian state approximation starts being inaccurate when the last two summands in the above expression for V (x, y) become comparable with the kinetic energy p 2 x /2 + p 2 y /2. This replaces the condition of large occupation numbers for the classical dynamics.

Gaussian state approximation for matrix quantum mechanics
In this Section, we apply the Gaussian state approximation to the matrix quantum mechanics with the Hamiltonian (1). To this end we decompose the matrix-valued canonical variables X i and P i in the basis of generators T a of S U (N) Lie algebra, with Tr (T a T b ) = δ ab and [T a , T b ] = iC abc T c : X i = X a i T a , P i = P a i T a . Following the construction of the Gaussian state approximation outlined in Section 2, we  (5). We obtain the following equations of motion: To study how quantum effects influence the classically chaotic dynamics of the Hamiltonian (1), we extend the definition of Lyapunov exponents to the quantum case by considering two solutions of the equations (12) for which the initial Gaussian states differ by a very small shift of the wave packet center X a i → X a i + a i . For the chaotic system, the difference δX a i between the expectation values X a i for these two solutions should first grow exponentially as |δX a i | ∼ | a i | exp (λ L t), where λ L is the leading Lyapunov exponent, and then saturate at a typical scale set by the system size. It is also easy to see that this definition is equivalent to the "out-of-time-order" correlator Our simulations were performed with d = 9 compactified spatial dimensions (motivated by the BFSS matrix model [14]) and S U (5) Lie algebra. The initial values of the coordinates X a i (t = 0) of the wave packet center were randomly drawn from the Gaussian distribution with dispersion σ 2 , the initial values of momenta P a i (t = 0) were set to zero. The quantum dispersions were set to the values X a i X b j = δ ab δ i j σ xx , P a i P b j = δ ab δ i j σ pp , X a i P b j = 0, σ xx = (N (d − 1)) −1/3 /2, σ pp = 1/ (4σ xx ) which minimize the expectation value of the Hamiltonian (1) in the space of pure Gaussian states. We choose a i to be a random vector with | a i | = 0.00001. We illustrate the time dependence of the distance |δX a i | between two close solutions of the equations (12) in the left plot on Fig. 2. For comparison we also show the time dependence of |δX a i | for the classical equations of motion. Since even the classical Lyapunov exponents depend on the scale of X a i variables as λ c L ∼ σ, we use σ t ∼ λ c L t for the time coordinate. The classical dynamics exhibits a clear exponential growth |δX (t) | ∼ exp (λ L t) with time, followed by the expectable saturation of Lyapunov distance at late times for a bound system. For quantum dynamics the growth rate of the Lyapunov distance coincides with the classical one only at some initial period of time, which becomes larger as the dynamics becomes more classical at larger σ. Presumably, this initial period of time can be related to the thermalization of our "artificial" initial state itself, which is then followed by the decay of small perturbations on top of the more universal thermalized state. At later times the leading Lyapunov exponent λ L for the quantum dynamics turn out to be several times smaller.
Let us now consider the temperature dependence of the quantum Lyapunov exponents obtained numerically from the correlator (13). In order to introduce temperature, we rely on the chaotic, "selfaveraging" classical dynamics of the Hamiltonian (1) which implies the equivalence between canonical and micro-canonical ensembles. At sufficiently high energies and temperatures we can use the classical equation of state [15,16] (14) and (15) are used to express temperature in terms of energy. to translate energy E into temperature T . At low temperatures, this equation of state changes due to quantum effects. In particular, in the limit of zero temperature the energy approaches some finite value E 0 of the ground-state energy [17]. To be consistent with the Gaussian state approximation, the ground state energy should be defined by minimizing the Hamiltonian (1) over all pure Gaussian states, which yields E 0 = 3 8 d (d − 1) 1/3 N 1/3 N 2 − 1 . In order to interpolate between the lowtemperature asymptotics E → E 0 (at T → 0) and the classical high-temperature equation of state (14), we use a phenomenological formula Since the first-principle equation of state of [17] would be anyway inconsistent with the Gaussian state approximation, we prefer to use the artifical equation of state (15). The most important conclusions of this work should not dramatically depend on moderate changes in the equation of state. Temperature dependence of quantum Lyapunov exponents is illustrated on the right plot on Fig. 2. For comparison, we also show the temperature dependence of the classical Lyapunov exponents, which is well described by the formula λ L = 0.292 − 0.42/N 2 T 1/4 [16], with the temperature T defined from (14). We see that quantum corrections tend to decrease the Lyapunov exponents at intermediate temperatures. Moreover, at some "critical" temperature T c ≈ 0.6 the leading Lyapunov exponent seem to vanish. Interestingly, first-principle Monte-Carlo simulations indicate that in the large-N limit the matrix quantum mechanics (1) has a finite-temperature confinement-deconfinement phase transition at roughly the same temperature T c ≈ 0.9 [17]. This has to be contrasted with the supersymmetric BFSS model, which is not confining all the way down to the strong-coupling regime at T = 0. By analogy with N = 4 super-Yang-Mills in D = 3 + 1, we expect the MSS bound λ L < 2πT [18] to be saturated in this regime. On the other hand, in the confinement phase of the nonsupersymmetric Hamiltonian (1), it is natural to expect that Lyapunov exponents are 1/N-suppressed, because the large-N limit is temperature independent [17] and hence is equivalent to zero temperature. Our numerical results shown on the right plot on Fig. 2 thus suggest that the growth of quantum corrections at low temperatures and probably the confinement-deconfinement transition prevent the violation of MSS bound, which for classical Lyapunov exponents happens at the temperature T = 0.292−0.42/N 2 2π 4/3 = 0.015 which is much lower than T c .

Discussion and outlook
In these Proceedings we have outlined the Gaussian state approximation for the real-time dynamics of many-body systems, used previously in the quantum chemistry context [8], and demonstrated on a simple example that it is capable of reproducing the essential features of quantum dynamics which are absent for the classical equations of motion. We applied this approximation to study the real-time dynamics of the matrix quantum mechanics, and found that quantum corrections tend to decrease the leading Lyapunov exponent λ L extracted from the simplest out-of-time-order correlator in| P a i (0) ,X a i (t) |in . We have shown that quantum corrections seem to be large enough to ensure the validity of the MSS bound λ L < 2πT at sufficiently low temperatures. Of course, we can make this statement only up to the unquantified systematic error of the Gaussian state approximation, which becomes quantitatively worse at low temperatures (see Fig. 1 for illustration). In particular, it would be interesting to understand in more details whether the vanishing of quantum Lyapunov exponents below some finite temperature is a signature of the finite-temperature phase transition found in [17], or rather an artifact of the Gaussian state approximation.
The generalization of the Gaussian state approximation to Yang-Mills theory should not meet significant conceptual difficulties. The resulting algorithm would require the number of Flops proportional to the square of the number of (spatial) lattice sites, which should allow for simulations on not very large lattices with modest computational resources.