Time evolution of a toy semiholographic glasma

We extend our previous study of a toy model for coupling classical Yang-Mills equations for describing overoccupied gluons at the saturation scale with a strongly coupled infrared sector modeled by AdS/CFT. Including propagating modes in the bulk we find that the Yang-Mills sector loses its initial energy to a growing black hole in the gravity dual such that there is a conserved energy-momentum tensor for the total system while entropy grows monotonically. This involves a numerical AdS simulation with a backreacted boundary source far from equilibrium.


Introduction
In the ultrarelativistic heavy-ion collider experiments at RHIC and the LHC, a strongly interacting deconfined state of matter, the so-called quark-gluon plasma (QGP), is produced. On extremely short time-scales, there appears to arise hydrodynamic behavior with a very low value of the specific viscosity [1,2], making the QGP arguably the most perfect fluid ever observed and presenting a challenge to a fundamental description based on perturbative QCD, which should govern the high-energy partonic degrees of freedom of the colliding nuclei with relatively small running coupling α s = g 2 /4π.
In the color-glass-condensate (CGC) framework [3], strong interactions are due to the large occupation numbers of order 1/α s of the gluons liberated in the early stages of the collision, which appear at a semi-hard momentum scale Q s , the so-called saturation scale, with Q s Λ QCD at RHIC and LHC. This out-of-equilibrium gluonic matter with high occupation numbers, dubbed glasma [4], can be approximated by numerical solutions of classical Yang-Mills equations [5][6][7] until occupation numbers become much smaller than 1/α s , where a kinetic description may take over [8]. However, this glasma copiously emits relatively soft gluons with momenta below the scale Q s where the running coupling is much stronger, producing a medium of soft gluons which is expected to play a dominant role in bottom-up thermalization according to [9,10].
In gauge/gravity dual descriptions of heavy-ion collisions [11][12][13], the perturbative QCD part of this evolution is usually ignored or used only to provide initial conditions for a subsequent holographic treatment, with no backreaction of the strongly coupled soft sector of the system on the more weakly coupled hard sector [14,15]. Such a backreaction appears to be particularly interesting given the fact that the gravity-dual description leads JHEP08(2018)074 to the qualitatively different top-down thermalization scenario where high-energy modes are the first to attain a thermal distribution [16]. 1 In [20], a 'semi-holographic' attempt to combine a description of the hard degrees of freedom in terms of glasma field equations and of soft strongly self interacting degrees of freedom in terms of gauge/gravity duality has been presented, where each sector modifies the dynamics of the complement sector through a minimal local coupling of the respective dimension-four operators. In [21] this approach was extended and slightly corrected to permit the construction of a conserved energy-momentum tensor of the total system, which is living in physical Minkowski space, while the gravity-dual sector has a nontrivial boundary metric. (As discussed in [22,23] it is also possible, and ultimately desirable, to have dynamics of the weakly coupled sector take place in yet another nontrivial auxiliary metric.) Moreover, ref. [21] has set up an extremely simple toy model for studying the combined system numerically. This was done by restricting the couplings to one between the two energy-momentum tensors pertaining to the Yang-Mills sector and the holographic sector, but in a spatially homogeneous and isotropic situation where the gravity side did not allow for propagating degrees of freedom. Nevertheless, nontrivial dynamics was observed that consisted of oscillatory energy transfer between the Yang-Mills system and the holographic side, with conservation of the total energy, but no thermalization. While the holographic side included a seed black hole, the latter could not grow absent any propagating modes in the bulk. The toy model of ref. [21] instead permitted controlled studies by having a closed-form solution on the gravitational side and thus a proof-of-concept of the proposed iterative procedure, which indeed showed quick convergence.
In order to study thermalization, which requires propagating degrees of freedom in the bulk that can build up a black hole, one has to either relax the restriction to spatial homogeneity and isotropy or include a coupling of scalar operators as indeed initially proposed. Computationally this means having to deal with numerical AdS calculations with a nontrivial boundary source that gets updated in subsequent iterations of solving the complete system. In this paper we present the successful implementation of this scheme in a variant of the toy model introduced in ref. [21], but reduced by one spatial dimension such that one deals with AdS 4 instead of AdS 5 geometries where it is easier to achieve high numerical accuracy. We expect, however, that this simplification already shows the main features present in the AdS 5 /SYM 4 system to be explored in future work.

The setup
In accordance with the CGC and glasma description of the early stages of the QGP, we model the hard degrees of freedom, i.e. dynamics at the saturation scale Q s , by a classical Yang-Mills theory [3,4]. Following the semiholographic approach to the QGP proposed in [20,21] we extend the glasma picture by including nonperturbative effects of soft degrees of freedom utilizing the gauge-gravity duality. In doing so we replace non-perturbative QCD JHEP08(2018)074 by N = 4 super Yang-Mills theory which at infinite coupling and a large number of colors allows for a dual description in terms of classical supergravity (SUGRA).
The interaction between the hard and soft sectors is established by promoting the marginal couplings of each sector to functions of gauge-invariant operators of their complements. The (coarse-grained) operators at our disposal in the effective description of the soft sector are the energy momentum tensor T µν , the glueball density operator H and the Pontryagin density A. For the scope of this paper we will restrict ourselves to the coupling of the scalar operator H only. This is obtained from the generating functional by varying with respect to the source, e.g.
where W , g µν , and h denote the generating functional, the metric of the space-time, and the source, respectively. The background metric g µν in our context only serves as a computational device and will be set to the Minkowski metric η µν eventually. Note that a non-trivial choice for the source corresponds to a marginal deformation of the theory.
In order to simplify the simulation presented in the following sections we will model the hard sector by a classical Yang-Mills theory in 2 + 1 dimensions and the soft sector by a gravitational theory in 3+1 dimensions. Nevertheless, we are confident that the qualitative behavior of the system will not change in a 3 + 1 dimensional space-time. In the remainder of this section the number of space-time dimensions d will be kept arbitrary.
Starting from the classical Yang-Mills action we can easily deform the hard sector with a scalar operator by adding a source term with the Yang-Mills coupling constant g YM . The non-Abelian field strength in terms of the gauge field is given by First of all we notice that physically this deformation amounts to locally rescaling the Yang-Mills coupling constant g YM . Second, a calculation of the reponse analogous to (2.1) yields In a next step we bring the two deformed sectors into contact by simply adding S YM [A µ , χ] and W [h] supplemented with an interaction term for the two scalar defor-

JHEP08(2018)074
The saturation scale Q s appears for dimensional reasons and is accompanied by a phenomenological dimensionless free parameter β which allows for tuning the interaction of both sectors. 3 Inspecting (2.4) immediately reveals that the two scalar fields, H and h, are auxiliary fields. Their equations of motion, given by indeed connect the deformations of each sector to gauge independent single trace operators of the respective other sector. After integrating out H and h the action becomes This form of the action is discussed in general including the tensor and pseudoscalar coupling channels 4 in [21]. In a recent publication [23] the tensor channel (which is omitted here) was further improved upon by employing a so-called democratic coupling approach.
Let us now turn to the equations of motion for the gauge field A µ and the calculation of the expectation value H. The former read where we introduced the gauge covariant derivative D µ = ∇ µ − iA a µ T a with ∇ µ denoting the Levi-Civita connection of the background metric g µν . For calculating H we employ the holographic dictionary, which maps the generating functional W to the (d + 1) dimensional SUGRA (on-shell) action and operators to fields in the gravity theory satisfying asymptotically anti-de Sitter (AdS) boundary conditions. The relevant terms of the action for our purposes are given by the Einstein-Hilbert action with a cosmological constant coupled to a massless Klein-Gordon scalar field with the AdS radius L set to unity in the remainder of the paper. The bulk equations of motion arising from (2.8) are In d = 4, g 2 YM carries mass dimension and should be understood as g 2 YM = 4παsQ 4−d s . Since in (fourdimensional) CGC applications the energy scale Qs is set by the gluonic degrees of freedom, with αs treated as a small parameter which eventually is set to αs ∼ 0.3 so that initial energy densities are numerically of order Q 4 s , we shall later simply set all dimensionful quantities on the Yang-Mills side to unity in terms of Qs and vary only β. 4 The special case considered here is obtained from the formulation in [21] by omitting the couplings γ and α.

JHEP08(2018)074
In Fefferman-Graham coordinates, where ρ = 0 denotes the location of the conformal boundary of the (d+1) dimensional space-time, the metric and the scalar field have the following asymptotic expansions [24] from which we can read off the expectation values , where X µν and ψ (d) are local functionals of the boundary sources. The leading coefficient in the metric expansion is fixed by the background metric g (0)µν = g µν . The non-normalizeable mode of the scalar field φ (0) is dual to the source in the generating functional W , which in our setup is related to the Lagrange density of the classical Yang-Mills sector The total energy momentum tensor of our semiholographic model is Each of the three contributions is obtained by the variation of the corresponding term in (2.4) with respect to g µν and employing eqs. (2.5). We will refer to the first contribution as the Yang-Mills energy momentum tensor, the second as the holographic energy momentum tensor and to the third as the interaction energy momentum tensor. Note that the expression here differs from the corresponding one in [21] where the terms in t µν YM proportional to H were assigned to the interaction energy momentum tensor. Although the assignment of the contributions to the different sectors is somewhat arbitrary, the advantage of the form presented here is that t µν int only consists of the agents responsible for the deformation, in our case the two simplest gauge independent single trace scalar operators.
The conservation of the energy momentum tensor ∇ µ T µν = 0 is implied by separate Ward identities in the respective sectors of our model The sum of the terms on the right hand side is precisely −∇ µ t µν int . Furthermore, we also want to mention the trace anomalies of the individual sectors, which read where A denotes the holographic conformal anomaly, which is a local functional of the boundary sources and vanishes for the case considered below. Note that in general even if both T µν and t µν YM are tracefree, eg. for g µν = η µν and d = 4, the full system is not conformal due to the contribution g µν t µν int = −d Hh.

Classical Yang-Mills sector
We will now work in a d = 2 + 1 dimensional space-time with g µν = η µν and restrict to isotropic homogeneous SU(2) color gauge fields in temporal gauge, A a 0 = 0, A 3 0 = 0 with a = 1, 2. To further simplify this toy model as far as possible, we make t µν YM diagonal by assuming color-space locking, A a i = δ a i f (t) and A 3 i = 0 with i = 1, 2. The single remaining degree of freedom f (t) satisfies an equation for an anharmonic oscillator with time dependent damping obtained from (2.7) 5 The energy density and the pressure are The source for the dilaton in terms of the YM fields is given by

Holographic sector
To be consistent with the YM sector we also impose homogeneity and isotropy in the spatial field theory directions of the bulk theory of the holographic sector. We make the following ansätze for the metric and the massless scalar field in in-going Eddington-Finkelstein coordinates The equations of motion (2.9) then take the following form 6 has been extensively studied [25,26]. 6 It is interesting to note that these equations are equivalent to those of a homogeneous but anisotropic black brane without scalar matter.

JHEP08(2018)074
where prime denotes radial derivatives f = ∂ r f and the dot-derivative is defined aṡ Near the boundary (r = ∞) solutions to these equations can be expressed as power series in r Fixing the conformal boundary metric to Minkowski ds 2 b = r 2 η µν dx µ dx ν determines the leading coefficients a 0 = 1 and s 0 = 1, and the residual gauge freedom r → r + ξ(v) is fixed by setting the subleading coefficient a 1 = 0. Solving the equations order by order in r gives where the normalizable modes φ 3 (v) and a 3 (v) remain undetermined in this procedure and need to be extracted from the full bulk solution. Furthermore one obtains the relation In order to identify the expectation values of the energy momentum tensor and the scalar operator it is convenient to asymptotically transform the series solutions (2.30) and (2.32) to Fefferman-Graham coordinates (2.10). The relevant coefficients in the Fefferman-Graham expansion in terms of their Eddington-Finkelstein counterparts are given by The expectation values of the energy momentum tensor and the scalar operator are then given by Evaluating the holographic Ward identity (2.14) reproduces the relation (2.33) we find from solving the near boundary expansion

JHEP08(2018)074
3 Solution procedure In this section we describe how we obtain solutions for the time evolution problem of the coupled system (2.18), (2.22)-(2.26) for given values for the couplings β and g YM and given initial energies in the Yang-Mills sector ( ini YM ) and the holographic sector ( ini hol ). Our approach is to solve the coupled system with an iterative procedure that starts with an initial guess for f (t) which we get from a solution to (2.18) for β = 0 Solutions to (3.1) can be written in terms of the Jacobi elliptic function where the integration constant C/g 2 YM = ini YM can be identified via (2.19) with the initial energy in the YM-sector. Without loss of generality we set t 0 = 0, which corresponds to our initial time. From the initial guess (3.2) and a chosen value for β we compute via (2.20) the time dependent boundary source for the gravity system φ (0) (t) = h(t).
For our model defined by (2.6) it is sufficient to consider only positive values of the coupling β. Switching the sign of β amounts to switching the sign of h(t) which we identify with the boundary source φ (0) of the scalar field in the gravitational bulk. Absent any (odd) scalar self interactions the action given in (2.8) is invariant under the transformation φ → −φ. Therefore, the holographic contribution to the energy momentum tensor will not change, while the expectation value H changes sign. However, in the Yang-Mills equation (2.7) as well as in the Yang-Mills energy momentum tensor only the combination βH appears, which are therefore invariant under β → −β. Finally, the product hH in the exchange part of the full energy momentum tensor is likewise invariant.
We numerically evolve the boundary sourced gravity system with the spectral method explained in [12], using 20 Chebyshev grid points in the holographic direction and a 4 th order Adams-Bashforth time stepping algorithm with step size ∆t = 1/800. In order to get a well defined initial value problem resulting in a stable time evolution it is necessary to choose a computational domain in the bulk direction that contains the apparent horizon r ah , defined byṠ(t, r)| r ah = 0, on the initial slice. 7 Initial data for the gravity system are fixed by a 3 (t = 0) = − ini hol /2 and a radial profile for the scalar field which evaluates to φ(r, t = 0) = −β ini YM for the initial guess (3.2) in combination with the Ward identy (2.37). To measure the accuracy of our numerical scheme we monitor in each time step the violation of the constraint equation (2.26) and the Ward identity (2.37) whose absolute values we demand to be smaller than 10 −6 . From the solution of the gravity problem we extract, via (2.35) and (2.36), the time evolution of T µν (t) and H(t) respectively.
As a next step we numerically solve the Yang-Mills equation (2.18) for the new f (t) with the Mathematica routine NDSolve, using f (0) and f (0) from (3.2) as initial conditions and the result for H(t) from the gravity simulation. This completes the first iteration. With the new f (t) we can evaluate the total energy momentum tensor (2.13) and check if it is conserved in time. This will typically not be the case after the first iteration and we have to iterate again. To initialize the next iteration we have to compute the new source for the gravity simulation. Compared to the initial iteration, where we provide an analytic guess for f (t), the updated f (t) in subsequent iterations is known only numerically which introduces via (2.20) some numerical noise in the new source for the gravity simulation. We reduce the numerical noise in the source with a low-pass filter before we feed it to the gravity code. As filtering tool we use the Mathematica routine LowpassFilter and choose a cutoff frequency of 0.1 and filter kernel of length 1. We continue this iterative procedure until the total energy is conserved to O(10 −5 ) or better. The iterative procedure is summarized in figure 1. In appendix A we discuss the procedure in more detail.

Results
Following the procedure outlined above, we compute the gauge field degree of freedom f (t), displayed in the left panel of figure 2 for three different values of β after four iterations. For the other parameters we chose the initial Yang-Mills energy density to be YM /Q 3 s = 1, the initial energy of the strongly coupled sector to be hol /Q 3 s = 0.004 and the Yang-Mills coupling as g YM / √ Q s = 1. As mentioned above, the criterion for an acceptable solution is to yield a constant total energy of the system. The right panel in figure 2 clearly shows that with each iteration the total energy for the choice β = 0.01 stays on its initial value for a longer period of time.
However, with each iteration the numerical errors in the solution of the respective sub-sectors accumulate, see appendix A. We stop the procedure after four iterations in this case, since it provides the optimal trade off between obtaining sufficiently well behaved total energy on the one hand and consistent sub-sectors on the other hand.
Furthermore, we learn from the plot of f (t) in figure 2 that the gauge field decreases in frequency and in amplitude with time. This behavior is more pronounced the higher the coupling β is. As a consequence we find that the energy of the Yang-Mills sector on average is decreasing in time, while the energy of the holographic sector is increasing almost monotonically as the first two panels of figure 3 show.
The third panel of figure 3 displays the interaction energy, which like the gauge field oscillates around zero and decays over time with decreasing frequency. Starting from this observation, one might speculate that eventually all energy from the Yang-Mills sector gets transferred to the strongly coupled sector, with decreasing rate. Note that only when the Yang-Mills sector is empty the source h(t) vanishes and thus the transfer of energy is no longer possible. As long as the source h(t) is varying in time one excites matter fields in the gravitational bulk which fall into the black hole causing its growth.
The choice for the initial values of the energy in the respective sub-sectors is motivated by the CGC picture of heavy-ion collisions, where the YM sector carries essentially all of the initial energy in the form of highly overoccupied gluons at the saturation scale, but the infrared sector to be described by holography is initially empty and thus represented by pure AdS spacetime. Due to numerical issues it is however necessary to start with a small regulator black hole in the gravitational bulk. The two panels in figure 4 compare the gain in the holographic energy for different initial conditions, while the initial Yang-Mills energy is kept fixed. We see that the results are fairly insensitive to the size of this regulator provided that ι := ini hol / ini YM 1 and thus our choice ι = 0.004 used for the plots above is a reasonable one.
In our setup the classical Yang-Mills sector consists only of a single dynamic degree of freedom given by f (t), hence the associated entropy is zero. However, in the holographic sector the area of the apparent horizon provides a commonly used proxy of entropy, which we use as estimate for the lower bound for the entropy in the combined system. In the left plot of figure 5 we show the radial position of the apparent horizon for the case in which β = 0.02. The right plot of figure 5 shows the entropy associated to the apparent horizon for different values of the coupling β. We find that the growth of entropy increases with β. Furthermore, we numerically checked that the effective apparent horizon entropy is monotonically increasing with time in all our simulations.

Conclusion and outlook
In this paper we have extended the toy model used in the first numerical tests [21] of the semiholographic model for heavy-ion collisions proposed in ref. [20] such that the glasma equations are coupled to a holographically described infrared sector which permits the formation of a black hole and thus entropy production and thermalization. The results we discussed present the first successful implementation of a self-consistent numerical AdS/CFT simulation involving a backreacted dynamical boundary source far from equilibrium. As our main result we find that the UV degrees of freedom modelled by a classical Yang-Mills theory lose their energy to the strongly coupled IR degrees of freedom over time (see figure 3), while the total energy is conserved. Motivated by glasma initial conditions one would like to start with all the energy deposited in the Yang-Mills sector and an empty strongly coupled sector. For numerical reasons we have to initialize the gravitational sector with a regulator black hole, but we showed that our results are insensitive to the size of a small regulator. Our numerical findings indicate that eventually the IR degrees of freedom will deplete the energy of the UV degrees of freedom entirely. However, in heavy-ion JHEP08(2018)074 collisions the glasma picture eventually ceases to be applicable, namely when the occupation numbers of the UV degrees of freedom become of order one or less. At that point a quantum effective kinetic theory description of the UV degrees of freedom should replace the description in terms of classical Yang-Mills fields.
We also want to stress that the model we set up can be applied to contexts other than heavy-ion collisions. In principle one can consistently couple any classical field theory to strongly coupled sources following our scheme.
In order to obtain a tractable toy model, we have simplified the semiholographic model by symmetry assumptions and worked in 2 + 1 dimension, which made it easier to obtain high numerical accuracy. However, we expect that in 3 + 1 dimensions our results will not change qualitatively.
In our future studies we plan to work in 3 + 1 dimensions as a next step. To that end we have to improve the numerical stability of our solution procedure. The bottleneck is mainly solving the Yang-Mills equation without introducing too much noise by numerical differentiations. Presently we rely on prebuilt routines of Mathematica when solving the equations, as well as removing some of the noise by filtering techniques. More importantly, we also want to relax the symmetry assumptions, incorporating anisotropies and spatial inhomogeneities by opening up the spin-2 coupling channel in the semiholographic couplings [21,23].
As a first step towards our long-term goal of simulating semi-holographic shock wave collisions with glasma initial conditions we want to study an axisymmetric chromo-electric or chromo-magnetic flux tube in an expanding background geometry along the lines discussed in [30], where our semiholographic model should provide the means to isotropize the system. Moreover, we plan to carry out quantitative comparisons between pure glasma, pure holographic, and semi-holographic glasma calculations.
From what we have learned in the present study we expect our model to be very efficient in converting energy in the UV degrees of freedom to a thermal bath represented by a dual black hole, when the sources vary locally in space time.
Note that although the strongly coupled sector is inherently in a quantum regime, our couplings do not couple fluctuations to the glasma. Thus, in order to improve the model from a more conceptual point of view one should implement couplings to quantum fluctuations. This is the subject of work in progress. A Numerical accuracy of the iterative procedure In this appendix we use an illustrative example, characterized by β = 0.2, g YM / √ Q s = 0.5, ini YM /Q 3 s = 0.1, and ini hol /Q 3 s = 1/250, to demonstrate the numerical feasibility of our iterative procedure. The left plot of figure 6 shows the violation of total energy conservation ∆ tot (t) = ini tot − tot (t), defined as difference between total initial energy ini tot = ini YM + ini hol , and total energy during the time evolution tot (t) = YM (t) + hol (t) + xc (t). Note that initially the exchange energy is zero. Typically, after four iterations our numerical scheme arrives at a solution in which ∆ tot (t)/Q 3 s ≈ O(10 −5 ) or smaller. A characteristic of our algorithm is that at earlier times less iterations are necessary to converge to the true solution than at later times. This effect is shown in the right panel of figure 6, where we plot the Yang-Mills energy in four subsequent iterations. This behaviour is induced by the way we choose our initial guess, which typically requires less improvement at earlier times, because then it differs not so much in amplitude and phase from the true solution. At later times, when already a significant amount of energy is transferred, guess and true solution can have very different amplitude and phase such that several iterations are necessary to achieve sufficient improvement.
In each of these iterations we monitor the Ward identity (2.14) and the constraint (2.26). In the left plot of figure 7 we see that in most time steps, and in all subsequent iterations, the Ward identity is fulfilled to an accuracy better than 10 −12 . In a comparably small number of time steps the accuracy systematically decreases with the number of iterations, but always remains below 10 −7 in this specific example. A similar picture holds for the constraint in the gravity simulation, shown in the right plot of figure 7. Also here, in most time steps the absolute value of the maximum violation (in bulk direction) of the constraint (2.26) remains smaller than 10 −12 in all subsequent iterations, and only for a small number of time steps the error grows with the number of iterations. The origin of this numerical noise, ultimately leading to a break-down of our algorithm, can be traced back to numerical errors introduced when solving the classical Yang-Mills JHEP08(2018)074 equation (2.18). In particular (2.36) shows that derivative of higher order enter the calculation of H, which then in turn complicate the solution of the classical Yang-Mills equation and make the filtering procedure necessary in the first place. Optimizing this part of the simulation is ongoing work.
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.