Real time quantum gravity dynamics from classical statistical Yang-Mills simulations

We perform microcanonical classical statistical lattice simulations of SU(N) Yang-Mills theory with eight scalars on a circle. Measuring the eigenvalue distribution of the spatial Wilson loop we find two distinct phases depending on the total energy and circle radius, which we tentatively interpret as corresponding to black hole and black string phases in a dual gravity picture. We proceed to study quenches by first preparing the system in one phase, rapidly changing the total energy, and monitoring the real-time system response. We observe that the system relaxes to the equilibrium phase corresponding to the new energy, in the process exhibiting characteristic damped oscillations. We interpret this as the topology change from black hole to black string configurations, with damped oscillations corresponding to quasi-normal mode ringing of the black hole/black string final state. This would suggest that α′ corrections alone can resolve the singularity associated with the topology change. We extract the real and imaginary part of the lowest-lying presumptive quasinormal mode as a function of energy and N.


Introduction
General relativity breaks down when curvature singularities appear. How these singularities are resolved in a consistent extension of general relativity is a very important issue.
The topology change from black hole (BH) to black string (BS) [1] is an interesting physical process for which the singularity resolution is essential: although rich dynamics beyond linear perturbation theory is expected based on results in numerical relativity (see e.g. [2]), the change in topology requires physics beyond general relativity.
In string theory, gauge/gravity duality [3] provides us with a well-defined description of this BH/BS transition in terms of a gauge theory dual [4,5] (see also [6][7][8][9]). To review the main aspects of this description, let us consider the example of two-dimensional maximally supersymmetric SU(N ) Yang-Mills (2d maximal SYM) compactified on spatial circle of circumference r x , which is conjectured to possess a dual string theory description. In fact, this gauge theory has two dual description, one being type IIB string theory on R 1,8 ×S 1 with N D1-branes wrapped on the circle, and the other one being type IIA string theory on a circle of circumference r x = (2π) 2 α rx with N D0-branes (known as T-dual). In terms of the 2d maximal SYM gauge theory description, the information about the individual D0-branes on the spatial circle is encoded in the phases of the eigenvalue distribution of the Wilson line winding on circle (referred to as Wilson loop in the following). Depending on the distribution of D0-branes along the circle, there can be various phases, such as a black hole phase (corresponding to a localized, or "gapped", distribution of Wilson loop phases), and a black string phase (corresponding to an "ungapped" distribution of Wilson loop phases). In the gravity dual of 2d maximal SYM gauge theory, one can further distinguish between a uniform black string phase (corresponding to a uniform distribution JHEP01(2019)201 Figure 1. The conjectured correspondence between the distribution of the phases of Wilson loop eigenvalues (top row) and the topology of the black hole/black string dual configurations (bottom row), with black hole, uniform black string and wavy string configurations corresponding to localized ("gapped"), uniform and nonuniform ("ungapped") phase distributions, respectively. Dotted lines indicate periodic boundary conditions on the circle. Figure adapted from ref. [11].
of Wilson loop phases), and a wavy string phase (corresponding to a nonuniform distribution of Wilson loop phases); see figure 1. Note that these phases have been explicitly been constructed in ref. [10].
Equilibrium properties of the 2d SYM gauge theory can be studied by using lattice Monte Carlo simulations, which are non-trivial to set up [12][13][14][15][16][17][18][19][20][21][22][23], but seem to be able to offer fully non-perturbative insights [24][25][26][27][28]. As a result of numerical studies in 2d maximal SYM [28][29][30] combined with numerical and analytical studies of the thermodynamic properties of D-branes [4,5,11,31,32], the conjectured equilibrium phase diagram sketched in figure 2 has started to emerge. Specifically, at low temperature T (large temporal radius β = T −1 ), two phases of localized and uniform eigenvalue phase distributions at small and large spatial radius r x , respectively, are separated by a first order phase transition line. Above a critical temperature a third phase corresponding to a non-uniform eigenvalue distribution arises, while the phase transitions soften to be of second and third order, respectively.
This tremendous progress non-withstanding, the question of topology change unfortunately cannot be answered using simulations of equilibrium quantities, because it requires information about real-time evolution.
Real-time information in quantum field theories are extremely difficult to obtain, because there is no known way to study the real-time dynamics of quantum systems at a reasonable computational cost. Absent any breakthrough development in the numerical treatment of real-time quantum systems, we take a more pedestrian approach in the present work by studying the real-time evolution in the classical statistical approximation, the non-perturbative lattice technology of which has been well developed in the context of relativistic heavy-ion collisions [33][34][35][36][37] (see also [38,39]). From a gauge theory perspective, the real-time dynamics of fully quantum 2d SYM is expected to be well approximated by its classical dynamics for modes which are highly occupied. For bosons in equilibrium, this is the case for the modes with energy ω obeying βω 1, or the low-energy modes. At small temporal circle radius β or high temperature, many bosonic modes can be well approximated by their classical dynamics. Additionally, at high temperature fermionic modes develop a large thermal mass, so they can be expected to decouple and thus not contribute to observables that are insensitive to the number of degrees of freedom of the theory. Thus at high temperature, one can expect the dynamics of 2d SYM to be at least qualitatively well approximated by the classical dynamics of purely bosonic Yang-Mills theory. This expectation has indeed be confirmed in previous numerical studies in other dimensions [40][41][42][43][44] and a comparison of Monte Carlo results between 2d SYM and 2d SU(N) pure gauge theory [11]. (Note that generalizations of the classical statistical framework to include quantum effects to some extent may be possible, cf. [45][46][47][48].)

JHEP01(2019)201
Not surprisingly, some aspects of the 2d SYM dynamics cannot be captured by the purely bosonic classical statistical approximation. For instance, by neglecting the fermions one is studying the phase diagram of the purely bosonic theory shown in figure 2, cf. refs. [5,11,49,50]. However, at high temperature, the 2d SYM and pure Yang-Mills phase diagrams become indistinguishable, such that the classical statistical simulations performed as part of this work could reasonably be expected to offer qualitative insights into the black string phase, and potentially also the deconfined black hole phase.
From string theory point of view, the high temperature (weak coupling) regime should describe a highly stringy system, where the α corrections are large but one can still tune the g s correction by dialing N . By performing direct numerical simulations of the high temperature gauge theory dynamics, one can study whether the α correction is sufficient

JHEP01(2019)201
for resolving the singularity associated with the topology change, and how the g s correction affects the result. A key asset of our classical statistical simulations is that we should be able to see how a perturbation of a black hole/black string rings down as the system (re-)approaches equilibrium. In the context of classical gravity, this ring-down process is encoded in the quasi-normal mode spectrum [51], which we attempt to measure numerically in purely bosonic Yang-Mills.
This paper is organized as follows: section 2 contains background on the simulation method, a particular scaling symmetry and the equilibrium phase diagram. Section 3 discusses rapid quenches of the system energy and corresponding signals of topology change, including quasinormal mode ringing. We summarize and conclude in section 4.

Equilibrium phase diagram
We study SU(N ) Yang-Mills theory with 8 scalars on a spatial circle, in Minkowski signature such that the dynamics is 1+1 dimensional. The theory will be set up using the tools from lattice gauge theory, by discretizing 9-dimensional classical Yang-Mills on a lattice where eight dimensions are toroidally compactified (see ref. [11] for more discussion about this point).

Simulation method
In order to solve the equation of motion while preserving gauge invariance it is useful to make use of standard lattice formulations of gauge fields. Our setup is based on our earlier work in ref. [11] on quantum Monte-Carlo simulations, which we briefly review here in order to keep this work self-contained. 1 We replace continuum 9-dimensional Euclidean space by an isotropic cubic lattice such that x i = ax i , i = 1, 2, . . . 9 withx taking on integer values and a being the (spatial) lattice spacing. With g denoting the strong coupling constant, the continuum gauge field variables A i (x) are replaced by link variables U i (x) = e aA i (x) = e −igaA a i (x)T a which are elements of the SU(N ) Lie group and live on links between lattice sitesx i and obeying U † . This allows to use the standard single plaquette definition for the Hamiltonian density and U is defined through [11,39]: (Note that denotes the sum over all spatial loops on the lattice starting from sitex i with only one orientation, e.g. ≡ 1≤i<j≤d .) The lattice equations of motion are calculated from the Hamiltonian equation and one finds

JHEP01(2019)201
when using the definition E i (x) = −ga 2 T a E a i (x) and discretizing time in units of t = a∆tt witht an integer. The update rule for the electric field is found by requiringḢ ≡ dH(t) dt = 0 [11,36,39]: where the gauge staple S ij is defined as in [39] as Note that for negative values of j, a gauge link is traversed in the opposite direction, e.g.
The Hamiltonian for the system is given by Similarly, the Gauss law constraint on the lattice is given by To prepare initial conditions which satisfy G(t) = 0, we simply take E i = 0 and start with link variables are Gaussian-random with predefined magnitude. Note that this is different than the initial conditions chosen for standard (quantum, Euclidean) lattice Monte Carlo simulations. Once initial conditions are specified, we determine the total system energy from measuring (2.5) and real time evolution on the lattice is then performed by using set of equations (2.3), (2.4) to time-step forward the fields U i , E i . Our evolution scheme is accurate up to O(a 2 ) corrections in time.
One of the main observables in this work will be the Wilson loop defined as a product over link matrices (2.7) We will study the absolute value of the normalized trace of the Wilson loop, It should be pointed out that for spatial directions with only one site and periodic boundary conditions, the corresponding gauge field becomes a scalar, e.g. A i → X i . In the following we consider the situation where eight spatial directions are compactified on a point, while the ninth direction is allowed to be large, so that A i → {X I , A x } with I = 1, 2, . . . 8.

Relation to matrix models and scaling symmetry
In continuum, the lattice-discretized theory described above corresponds to the classical Lagrangian where µ = {t, x} and F µν , D µ are the two-dimensional field strength and gauge covariant derivative, respectively. In the 't Hooft large-N limit, the 't Hooft coupling λ = g 2 N , the compactification radius r x and the energy per d.o.f E/N 2 are fixed to be of order N 0 . Since in two dimension λ is dimensionful, in the following we employ dimensionless units such as In temporal gauge A t = 0, the classical equations of motion in continuum become which obey the following scaling symmetry With this rescaling, we can set r x = 1. Therefore, the energy E and the number of colors N are the only relevant simulation parameters.

UV problem in classical statistical lattice simulations
Classical statistical simulations suffer from a well-known ultraviolet instability (this is the same as the Rayleigh-Jeans law for black body radiation which is cured by quantum mechanics). Given a finite lattice discretization scale a (the lattice spacing), classical dynamics will eventually start to populate modes close to the lattice UV cutoff scale a −1 at late times. The dynamics of these high momentum modes, however, should involve quantum effects which are absent in the classical statistical simulations. Therefore, once a simulation starts to become sensitive to the lattice UV scale, the resulting dynamics can no longer be trusted, and the simulation has to be stopped. Nevertheless, classical lattice statistical simulations with a fixed UV cutoff a −1 can successfully be used for system properties dominated by modes in the IR. In practice, we prepare initial conditions for the simulations which are well localized in the IR, and then run the simulations for sufficiently short times before modes start to pile-up at the UV lattice cutoff scale. To check if pile-up has occurred, we monitor the Fourier transform of the magnetic component of Hamiltonian, which has successfully been used as an indicator for this purpose in the past [35,53]. be seen that the system starts out dominated by IR physics, but eventually flows to the UV. As long as the system is dominated by IR physics, observables such as the expectation value of the Wilson loop (2.7) are unaffected by the classical UV instability, but become sensitive as soon as pile-up of modes in the UV occurs, as shown in the right panel of figure 3. For this reason, the results reported below have been obtained by only using simulation data for times t < t UV , where t UV denotes the time when UV pile-up has first occurred.

Results for equilibrium phase diagram in classical statistical simulations
As a reminder, classical statistical simulations are performed in the microcanonical ensemble (at fixed energy E) instead of the grand-canonical ensemble (fixed temperature T ) employed in quantum simulations of lattice field theory. Real-time information on observables may be obtained by averaging time-dependent results over classical ensembles. In order to connect to equilibrium properties of the system, one would want to measure observables after the system has become time-independent (thermalized) at time t ≥ t therm , but before the simulation becomes contaminated by UV artefacts t ≤ t UV . We find that both t therm and t UV depend on the normalized energy E/N 2 in simulations.
In practice, for fixed E/N 2 we consider the time-evolution of observables such as the ensemble-averaged Wilson loop expectation value shown in figure 3, and perform an additional time-average over t therm ≤ t ≤ t UV to increase statistics. Results from this procedure for the Wilson loop expectation value as a function of energy for various N are shown in figure 4. We find that results for different N ≥ 16 cluster around an apparently universal curve, indicating that |W | (E/N 2 ) is approximately independent of N .  At large E/N 2 , where our classical-statistical simulations can reasonably be expected to be a good approximation of the full quantum results, we find that |W | → 0, in accordance with bosonic quantum theory simulations [11]. Conversely, bosonic quantum theory simulations indicate that |W | (E/N 2 ) = 1 for low temperature, whereas from figure 4 it can be seen that classical statistical simulations result in |W | (E/N 2 ) → 1 for E → 0, because all link variables are equal to the identity matrix in this limit. Furthermore, full quantum theory simulations exhibit confinement for low temperatures [11], whereas classical statistical simulations do not show confinement. Therefore, as expected, classical statistical simulations can not be used to study properties of, or the transition to, the confined center broken phase of the full bosonic quantum theory (cf. figure 2).
However, one may ask if classical statistical simulations can be used to probe the properties of the quantum theory in the deconfined phase, notably the region close to the center symmetry phase transition shown in figure 2. In order to study the phase structure more precisely, the distribution of the Wilson line phases ρ(θ) are shown in figure 5 as a function of energy.
At high energy, one finds that the phase distribution ρ(θ) in the classical statistical simulations is ungapped, see figure 5. This is consistent with the full bosonic quantum simulations [11], and as such is consistent with the conjectured gravity picture of a black string phase at high temperature that was eluded to in the introduction.
Conversely, at low energy E/N 2 → 0, we see clear evidence for a gapped distribution of Wilson loop phases in the classical statistical simulations that can be well-fit by a semicircle distribution semi-circle GWW model As indicated in figure 2, a phase transition from ungapped to gapped Wilson line phase is expected as the temperature is lowered at fixed circle length in the full bosonic quantum theory. (In SYM, the transition in the microcanonical ensemble has been found to be first order in ref. [10]). A priori, it is not obvious that such a transition should be accessible using classical statistical simulations. However, the fact that the behavior of the Wilson loop eigenvalues changes qualitatively between low and high energy suggests that classical statistical simulations are able to access both the ungapped and gapped deconfined phase separated by either a phase transition or rapid analytic crossover transition, respectively. Thus, while classical statistical simulations cannot be expected to be quantitatively accurate approximations in the gapped phase, our results strongly suggest that they can be The gapped Wilson loop phase distribution observed at low energy in the classical statistical simulation is consistent with the conjectured gravity dual picture of a localized black hole configuration. As the energy in the classical statistical simulations is increased, the fit of ρ(θ) in terms of the semi-circle distribution (2.13) leads to increasing fit parameter c (see figure 6), where c = c crit = π having the natural interpretation 2 of separating blackhole from black string configurations at a critical normalized energy of where the uncertainty is primarily coming from the residual N -dependence of c, cf. figure 6. Note that at for c crit = π, the expectation value for the Wilson loop becomes We find that for E/N 2 ≥ E/N 2 critical , the Wilson loop distribution is rather well fit by another one-parameter form given by which is reminiscent of the analytically known distribution for the non-uniform string [31,55] in the Gross-Witten-Wadia model [56,57]. Best-fit values for k are indicated in figure 5.

JHEP01(2019)201
Our current results are not precise enough to say anything about the nature of the transition from gapped to ungapped Wilson loop phase in classical statistical simulations, nor can we confirm or rule out the presence of a third phase at high temperature expected from gravity, namely uniform Wilson loop phase distributions (cf. figure 1).

Real-time response to quenches
The results from the previous section strongly suggests that classical statistical simulations may be used to probe the properties of the black-hole/black-string transition that are in qualitative, if not quantitative, agreement with full quantum theory simulations in equilibrium.
Unlike current quantum simulations, no obstacle prevents the application of classical statistical simulations to non-equilibrium problems, suggesting that our method can be used to probe the real-time dynamics of the topology change from black hole to black string phases. Assuming that the real-time evolution of the Wilson loop phase distribution ρ(θ, t) possesses a good large-N limit, then such topology changes may take place within a finite time even at large-N .

Black-hole/black-string topology change
In order to study topology changes in a controlled manner in classical statistical simulations, we employ the following protocol: 1. Generate a classical statistical gauge field configuration with initial energy E/N 2 < 750 (E/N 2 > 750) expected to be in the black hole (black string) phase 2. Evolve the gauge field configuration until t ≥ t therm so that early-time transients have disappeared 3. Rapidly quench the system energy to a fixed final value of E/N 2 > 750 (E/N 2 < 750) without changing the Wilson line phases distribution or violating the Gauss law constraint (2.6). Note that the new energy indicates an equilibrium configuration in the respective other phase. This quench can be achieved by multiplying the electric field by a constant q, 4. Measure the real-time response of observables as the system tries to attain the new equilibrium black string (black hole) configuration Repeating the above protocol for many configuration with fixed initial energy and averaging over this classical ensemble of configurations leads to real-time results for observables shown in the following. We have repeated the above procedure to study the real-time response of other quenches, verifying in particular that it is also possible to observe the inverse process of gapped phase to ungapped phase, which we associate with a topology change process from black-hole to black-string configurations.

Possible observation of quasinormal modes
The real-time evolution of Wilson loop phase distribution ensemble-averages show oscillations around the new equilibrium configuration after a rapid quench, cf. figure 7. These oscillations are particularly visible in the ensemble-averaged Wilson loop expectation value shown in figure 8, and are apparently displaying little sensitivity to the choice of N ≥ 16. Similar oscillations in classical statistical simulations are ubiquitous, with frequency and damping rates associated with the mass and width of quasi-particles that are in good agreement with perturbative quantum field theory, cf. refs. [35,58].
By contrast, in classical gravity, oscillations of excited geometries for instance of black holes are characterized in terms of their quasinormal mode ringdown behavior, cf. ref. [51]. If the conjectured relation between Wilson loop phase distribution and geometry holds, this . Real and imaginary part of lowest-lying presumptive quasinormal mode frequency ν as a function of final system energy for N = 32. Multiple entries correspond to quenching protocols with different initial, but the same final energy, e.g. E/N 2 = 500 → 700 and E/N 2 = 900 → 700, and have been slightly displaced on the plot to increase visibility.
naturally leads to the interpretation of linking quasiparticle oscillations in gauge theory with quasinormal mode oscillations of black holes in string theory.
We are able to obtain estimates for both the real and imaginary part of the lowest-lying mode frequency ν by fitting the location and height of the extrema of |W |(t) to a form |W |(t) ∝ Re e iνt . Results from this fitting procedure for various final energies E/N 2 are shown in figure 9. Also shown in figure 9 are results from the fitting procedure when changing the initial energy but leaving the final energy after the quenched unchanged. Our results seem to imply that the extracted results for ν show little sensitivity to initial system energies, instead only depending on the final E/N 2 . The behavior of ν seems to be smooth and continuous with energy, and as such is apparently insensitive to the phase change near (2.14). Furthermore, our results for ν are consistent with those from matrix model simulations which do not suffer from UV problems [59], to which our simulations reduce to in the zero volume limit.
The smooth analytic behavior of ν on E/N 2 in classical statistical simulations may be obtained from the classical scaling symmetry outlined in section 2.2. Under the scaling symmetry, frequencies are expected to scale as ν → αν, and the energy density scales as → α 4 . At fixed temperature, assuming the energy density in classical statistical simulations to be approximately independent of volume, we therefore expect ν ∝ E N 2 1/4 at fixed volume. Note that this scaling argument would apply to any time-dependent quantity, such as for example for the Lyapunov exponent λ L ∝ (E/N 2 ) 1/4 , cf. ref. [44].
Results shown in figure 9 indicate that ν increases with energy, qualitatively consistent with the finding that the quasinormal mode frequency calculated in type II supergravity on black p-brane background, increases with temperature regardless of the value of p [60]. Quantitatively, the quasinormal mode frequency in type II supergravity for p = 1 differs from the classical statistical results shown in figure 9, which is expected because supergravity results are applicable for small temperatures whereas classical statistical approximations are quantitatively accurate at high temperatures.

JHEP01(2019)201 4 Summary and conclusions
In the present work, we have performed classical statistical simulations of microcanonical ensembles in SU(N ) Yang-Mills theory with eight scalars on a circle. Depending on the energy, we have identified two distinct equilibrium phases of the Wilson loop eigenvalue distribution which qualitatively correspond to those expected for black holes and black strings in the conjectured dual gravity picture. We found that gapped Wilson loop equilibrium phase distributions occurred for energies E < E crit , while ungapped distributions occurred for E > E crit , with our estimate for E crit given in (2.14). Our present results were not precise enough to decide if there is an actual phase transition at E = E crit as opposed to an analytic cross-over in classical statistical simulations, nor can we confirm or rule out the presence of a second transition to a phase of uniform Wilson loop distributions at very high energies.
We were able to perform real-time measurements of Wilson loop distributions following a quench in system energy from one phase to another. We found characteristic oscillations in the phase distribution and the Wilson loop expectation value, which showed very little sensitivity on the number of colors for N ≥ 16 or the state the system was in before the quench. Interpreting these oscillations as the string-theory analogue of quasinormal mode ringdown of black holes in classical gravity, we were able to extract estimates for the real and imaginary part of the lowest-lying presumptive quasinormal mode as a function of energy. Our results were found to be in qualitative, but not quantitative, agreement with analytic calculations of quasi-normal modes in the supergravity approximation.
There are several natural extensions to this work. For instance, one might be able to study Lyapunov exponents in classical statistical simulations similar to refs. [43,44,61], including 1/N effects.
Another natural extension would be to consider simulations of Yang-Mills on a 2dimensional torus rather than a circle, where a much richer phase diagram is expected based on gravity calculations [62]. Within the present simulation environment based on ref. [11], such change is operationally trivial and just corresponds to extending the number of lattice sites along a second direction.
To conclude, based on our results obtained in this work we expect that classical statistical simulations of SU(N ) Yang-Mills theory could become a valuable tool to study real-time phenomena in quantum gravity that are otherwise hard to access.

JHEP01(2019)201
Since then we had several conversations regarding the real-time dynamics of quantum gravitational systems. The classical Yang-Mills simulation is one of the options Joe encouraged us to try. Perhaps we are still at a very primitive stage, but we hope that we can make a steady progress toward Joe's dream! Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.