Holographic thermalization in N=4 Super Yang-Mills theory at finite coupling

We investigate the behavior of energy momentum tensor correlators in holographic $\mathcal{N}=4$ super Yang-Mills plasma, taking finite coupling corrections into account. In the thermal limit we determine the flow of quasinormal modes as a function of the 't Hooft coupling. Then we use a specific model of holographic thermalization to study the deviation of the spectral densities from their thermal limit in an out-of-equilibrium situation. The main focus lies on the thermalization pattern with which the plasma constituents approach their thermal distribution as the coupling constant decreases from the infinite coupling limit. All obtained results point towards the weakening of the usual top-down thermalization pattern.


Introduction
Understanding the complicated field dynamics in a heavy ion collision presents a difficult challenge to QCD theorists. Experiments at RHIC and the LHC point towards the conclusion that the quark gluon plasma (QGP) created in heavy ion collisions behaves as a strongly coupled, nearly perfect, liquid [1,2] rather than a weakly interacting gas of quarks and gluons. The strongly coupled nature of the created matter has made the gauge/gravity duality one of the standard tools in describing QGP physics [3,4], supplementing traditional approaches such as perturbation theory or lattice gauge theory.
In its original form the gauge gravity duality relates supergravity on five-dimensional asymptotically anti deSitter space time (AdS) to strongly coupled N = 4 super Yang Mills (SYM) theory living on the boundary of the AdS space. Although SYM is very different from QCD in its vacuum state, it shares many features with QCD in the deconfined phase, such as a finite screening length, Debye screening and broken supersymmetry.
One particularly useful development is the application of the duality to out-of-equilibrium systems by mapping the thermalization process to black hole formation in asymptotically AdS space. This has led to the insight that fluid dynamics becomes a good approximation rather quickly, but this does not mean that the system is isotropic or thermal [5,6,7,8,9].
A particularly important challenge in an out-of-equilibrium system is to identify the thermalization pattern with which the plasma constituents of different energies approach their thermal distribution. On the weakly coupled side classical calculations have shown that the thermalization process is of the bottom-up type, i.e. low energetic modes reach thermal equilibrium first, with inelastic scattering processes being the driving mechanism behind it [10]. In the early stages many soft gluons are emitted which form a thermal bath very quickly and then draw energy from the hard modes. Recently this picture got supported by numerical simulations [11]. In [12] an alternative proposal was made: the thermalization process is driven by instabilities which isotropize the momentum distributions more rapidly than scattering processes 1 . On the contrary, holographic calculations in the infinite coupling limit always point towards top-down thermalization, where the high energetic modes reach equilibrium first, indicating a probable transition between the two behaviours at intermediate coupling [14,15,16,17].
Evaluating non-local observables such as two point functions in a time dependent thermalizing system is an extremely challenging but important task, since they allow to see how different energy/length scales approach thermal equilibrium. One strategy was worked out in [15] where it was shown how fluctuations created near the horizon and dissipation come to a balance to satisfy a generalized fluctuation dissipation theorem. For a different approach to generalize the fluctuation dissipation theorem see [18]. Non equilibrium generalizations of spectral functions and occupation numbers were introduced in [19,20]. Complementary quantities of interest are entanglement entropy and Wilson loops [21,22]. A particularly useful model, due to its simplicity is the collapsing shell model, where the thermalization process is mapped to the collapse of a spherical shell of matter and the subsequent formation of a black hole [23,24,25,26,27,28,29,16,30,31,32,33,34,35,36]. In the limit where the shell's motion is slow compared to the other scales of interest this model was used to study the approach of the spectral density to equilibrium for the components of the energy momentum tensor in [37] and for dileptons and photons in [38,39] 2 . All these studies show the usual top-down thermalization pattern. In addition in [45] the virtuality dependence of the photons was taken into account, showing that onshell photons are last to thermalize, consistent with the conclusions from other models of holographic thermalization [46,47].
Due to its simplicity the collapsing shell model, in the quasi-static limit, even allows one to include the first order string corrections to the supergravity action and leave the infinite coupling limit. In [48] the leading order string corrections to the photon spectral density [49,50] were generalized to an out-of-equilibrium situation, showing indications that the usual top-down thermalization pattern shifts towards bottom-up when finite coupling corrections are included. This observation was strengthened by a quasinormal mode (QNM) analysis at finite coupling in [45]. As the coupling constant is decreased the tower of poles bends towards the real axis, also showing a weakening of the top-down pattern, but this time independent of the thermalization model being used. It still needs to be investigated if the observed change is intrinsic to photons or/and a consequence of the collapsing shell model or of more general validity.
In [51] the AdS-Vaidya solution was used to investigate the thermalization time scale for non-local observables in SU(N) N = 4 SYM theory at finite coupling using geometric probes in the bulk. Interestingly, there the UV modes thermalize faster and the IR modes slower if the coupling constant is decreased from the infinite coupling limit. The authors speculate that the difference between their analysis and the one for photons [48,45] originates from the fact that in order to study current correlators, it is necessary to include the Ramond-Ramond five form field strength in the O(α 3 ) corrections, which produce very large corrections to observables associated with electric charge transport. We will say more about this in the conclusions.
The goal of this paper is to shed light on the above issues by studying energy momentum tensor correlators of a N = 4 SYM plasma and their approach to thermal equilibrium at finite coupling in the collapsing shell model. In the infinite coupling limit the correlators were first studied in [52,53,54] and to next-to-leading order in a strong coupling expansion in [55,56]. The leading order corrections in inverse powers of numbers of colours, N c , was computed in [57,58]. Finite coupling effects on jet quenching were worked out in [59,60]. The out-of-equilibrium dynamics using the collapsing shell model at infinite coupling was considered in [37]. In the paper at hand we fill the missing gap by analyzing the flow of the quasinormal mode spectrum as a function of the coupling constant as well as the approach of the spectral density to its thermal value at finite coupling in the collapsing shell model.
The paper is organized as follows. In section 2 we will review the collapsing shell model and introduce the finite coupling corrections. After that we outline the main parts of the calculation in section 3 and present the results for the quasinormal modes and the spectral densities in section 4. In section 5 we draw our conclusions.

The collapsing shell model
Our aim is to use the collapsing shell model introduced in [23,24] to gain insights into the thermalization process of a strongly coupled N = 4 SYM plasma via the gravitational collapse of a spherically symmetric shell of matter in anti deSitter space. On the field theory this corresponds to the preparation of an exciting state through the injection of energy and the subsequent evolution towards thermal equilibrium.
Following Birkhoff's theorem, outside the shell the background is given by a black hole solution, whereas inside the shell the metric is given by its zero temperature counterpart. The five-dimensional AdS metric is given by where and u ≡ r 2 h /r 2 is a dimensionless coordinate where the boundary is located at u = 0 and the horizon at u = 1. From now on the index '-' denotes the inside and '+' the outside space time of the shell and we set the curvature radius of AdS to L = 1.
The shell can be described by the action for a membrane [37] where g ij is the induced metric on the brane and p is the only parameter that characterizes the membrane. Due to the discontinuity of the time coordinate in the above metric, fields living in the above background have to be matched across the shell using the Israel matching conditions given by where ij is the extrinsic curvature and κ 2 5 = 8πG 5 is Newtons constant in the Einstein frame 3 . From the above equation the trajectory of the shell is determined.
Physical initial conditions for the shell, that could be a good approximation for heavy ion collisions, are determined through the relation of the holographic coordinate with the temperature r h = T π and the saturation scale r s = Q s π together with a vanishing initial velocity [32]. For LHC these values are T ∼ 400 MeV and Q s ∼ 1.23 GeV.
We, however, are not going to treat the dynamical process but work in the quasi-static approximation, where the motion of the shell is slow compared to the other scales of interest and only take snapshots of the shell at certain positions. See appendix A for the exact relation when the quasi-static approximation holds. This condition, however, breaks down at the latest stages of the collapse as can be seen by comparing the Penrose diagram for the black hole space time with the collapsing shell diagram [32].
When the quasi-static approximation is applicable, the matching conditions simplify considerably and explicit calculations in frequency space are possible. In this case the discontinuity of the time coordinate implies that the frequencies measured inside and outside the shell are related through [37] The subscript s denotes the position of the shell at u = u s . The matching conditions at the shell for metric perturbations of the form have been worked out for all the relevant components in [37]. For example for the xy component they are

Finite coupling corrections
In order to leave the strict λ = ∞ limit the leading order string corrections to type IIB supergravity have to be included and this is accounted for by the action [61,62] 3 In our numerical calculations we always set κ 2 5 p = 1.
where γ ≡ 1 8 ζ(3)λ − 3 2 . F 5 is the five form field strength and φ is the dilaton. The C 4 term is proportional to the fourth power of the Weyl tensor The γ-corrected AdS black hole metric derived from the above action can be written as [62,63,64] where the coefficients c i = c i (u) only depend on the dimensionless holographic coordinate u = r 2 h /r 2 and dΩ 2 5 is the metric of a five-dimensional unit sphere. The solution can be written explicitly as with f (u) given in (2) and The γ-corrected relation between r h and the field theory temperature reads r h = πT /(1 + 15γ). Note that the vacuum solution of AdS does not receive γ-corrections [65]. In the next section we will calculate the spectral densities for the different symmetry channels of the energy momentum tensor obtained from the above metric.

Spectral density
In order to see how the plasma constituents approach thermal equilibrium we study the spectral densities of various energy momentum tensor components by considering linearized perturbations of the five dimensional metric where the linear perturbations h µν correspond to the energy momentum tensor of the field theory. Following [54] the metric perturbations can be combined into three gauge invariant fields Z s representing the three symmetry channels, namely spin 0 (sound channel), spin 1 (shear channel) and spin 2 (scalar channel).
Going to momentum space we look at fluctuations of the form In the following it will be convenient to introduce the dimensionless quantitieŝ Deriving the equations of motion for the different symmetry channels is a lengthy and tedious exercise which has been performed in the literature before [56,55,66]. Therefore we shall only describe the main points of the derivation here and guide the interested reader to the relevant references for further details. We will extend the solutions obtained in the hydrodynamic limit [56,55,66] to arbitrary momenta and energies.

Scalar channel
The EoM for the scalar channel is obtained by considering the metric fluctuations h xy . It is convenient to introduce a field By expanding this field to linear order in γ the EoM for the scalar channel takes the compact form [66] where the right hand side only depends on the zeroth order solution. Note that we have an additional term of O(γ) (the last term) in the above equation compared to [66]. This comes from the fact that we defined the dimensionless quantitiesω andq with respect to the γ-corrected temperature T and not T 0 = r h /π as in [66].

Shear channel
Following [54,56], the shear channel is defined by the metric fluctuations Using the gauge condition h xu = 0 and introducing H tx = uh tx /(πT ) 2 and H xz = uh xz /(πT ) 2 one can define the shear channel gauge invariant combination for which one obtains a decoupled second order differential equation for Z 2 to O(γ) upon introducing The equation of motion for the shear channel to O(γ) is given by The source term γJ 2 is of O(γ) and depends only on the zeroth order solution Z 2,0 and derivatives thereof. The explicit form of this lengthy expression can be found in [56].

Sound channel
In order to investigate the sound channel we look at metric perturbations of the form where h = α=x,y h αα is a singlet. After the equations of motion have been derived for the above perturbations the gauge conditions are imposed. In the sound channel there is a subtlety in constructing the gauge invariants due to the non constant warp factor in front of the unit five sphere. See [67,56] for a detailed analysis of this issue. It turns out that after introducingĥ tt =ĉ 1 H tt ,ĥ tz =ĉ 2 H tz ,ĥ =ĉ 2 H andĥ zz =ĉ 2 H zz the gauge invariant for the sound channel is given by Again, we are not showing the lengthy expression for the source term J 3 which only depends on the zeroth order solution and can be found in [56,66]. As before, there is an additional term appearing in the equation of motion due to the different convention of the dimensionless quantitiesω andq.

Solving the EoMs
We are now going to solve the equations of motion and determine the corresponding spectral densities. The equations (19), (23) and (28) have singular points at u = ±1, 0. In the near horizon limit, u → 1, the indicial exponents are given by ∓iω/2, where the minus sign corresponds to the infalling mode and the plus sign to the outgoing mode. In thermal equilibrium, i.e. in the black hole background, one chooses the infalling boundary condition because classically nothing can escape from a black hole. However, in the presence of a shell the solution is a linear combination of the ingoing and outgoing modes where s = 1, 2, 3 and the coefficients c ± are determined by the matching conditions specified below. In order to solve the EoMs (19), (23) and (28) numerically we make the following ansatz for the ingoing and outgoing modes where Z (0,1) s,in,out is regular at the horizon u = 1 and normalized to Z (0,1) s,in,out (u = 1) = 1. We then integrate numerically from the horizon to the boundary.
This solution has to be matched to the zero temperature solution inside the shell. At zero temperature there are no γ-corrections at leading order and therefore we can make use of the results obtained in [37], where it was shown that the inside solution for all channels is given in terms of the Hankel function of the first kind Z s,− (u) = uH The factor f γ s enters by matching the inside frequency to the outside frequency via (5) and is the only source of γ-corrections for the inside solution, which thus takes the form and can be obtained by a simple expansion of the Hankel function given in (31). The matching conditions for all channels have the compact form [37] Z s, and lead to the ratio where all the γ-corrections have to be taken into account. This ratio, in particular the parametric dependence of C 0 and C 1 on the frequency, will play an important role for the thermalization pattern. Having solved for the ingoing and outgoing modes the retarded correlators for the gauge invariants can be calculated via standard AdS/CFT techniques [53,52], which produce where we have dropped the contact terms and χ s,th is the thermal spectral density. All quantities have to be expanded consistently to linear order in γ. The relations between the retarded correlators of transverse stress, momentum density and energy density with the gauge invariant correlators are given by [54] G xy,xy = 1 2 G 1 , and the spectral density is defined as the imaginary part of the retarded correlator χ µν,ρσ (ω,q, u s γ) = −2 ImG µν,ρσ (ω,q, u s , γ).
To see how thermal equilibrium is approached it is instructive to look at the relative deviation of the spectral density from its thermal equilibrium value R s (ω,q, u s , γ) = χ s (ω,q, u s , γ) − χ s,th (ω,q, γ) χ s,th (ω,q, γ) .
This ratio is does not get altered by the relation for the symmetry channels with the retarded correlators (37), therefore we also use the shorthand notation χ s = −2 ImG s .

Results
After describing the main parts of our computation we are now ready to investigate numerically the corresponding results. We will start by analyzing the quasinormal mode spectrum of the different channels obtained from the thermal correlators and investigate the flow of the poles as a function of the coupling. After that we will use the collapsing shell model to analyze the spectral densities and their approach to thermal equilibrium at finite coupling.

Quasinormal mode spectrum
Quasinormal modes characterize the response of the system to infinitesimal external perturbations and are the strong coupling equivalent to quasiparticles and branch cuts at weak coupling [68,69]. They are solutions to linearized fluctuations of some bulk field obeying incoming boundary conditions at the horizon and Dirichlet boundary conditions at the boundary. They appear as poles of the corresponding retarded Green's function and have the generic form where q is the three momentum of the mode, M n and Γ n correspond to the mass and the decay rate of the excitation, respectively. For gravitational perturbations the QNM spectrum was first obtained in the infinite coupling limit in [54] and the diffusion poles in the hydrodynamic limit at finite coupling were worked out in [56,66]. We are extending the existing analysis and study the flow of the tower of QNM obtained in [54] as a function of the 't Hooft coupling λ.
In order to solve for the QNM spectrum we make a Frobenius ansatz for the ingoing Re Ω Im Ω Figure 3: The flow of the QNM in the sound channel for q = 0 (left) and q = 2πT (right). modes and solve recursively for the coefficients a s,n and b s,n by plugging the expansion into the equations of motion (19), (28), (23), while the parameter N is chosen large enough such that the behaviour of both functions is stable. We then make the ansatz for the frequencies, ω = ω 0 + γω 1 and solve numerically for the zeros of Z s,in (0, ω) = Z The results for the flow of the QNMs is displayed in fig. 1 for the scalar, in fig. 2 for the shear and in fig. 3 for the sound channel. In all channels a clear trend is visible. As the coupling constant is lowered from the λ = ∞ limit the imaginary part of ω n increases rapidly, lowering the decay rate of the excitation. There is also a strong dependence on the index n: Higher energetic excitations show a stronger dependence on the coupling with a larger shift towards the real axis. These results point towards a weakening of the usual top-down thermalization pattern, in accordance with the results found in [45], where an equivalent calculation for the R-current correlator was performed. It should be noted, though, that the strong coupling expansion can only be trusted fully if the relative deviation of the poles from the λ = ∞ limit is small, which clearly is not the case for all displayed poles.

Thermalization of the spectral density
Next the behaviour of the spectral density and its deviation from the thermal limit in the collapsing shell model is investigated. To this end we parametrize the momentum of the plasma constituents by q = c ω. For c = 0 the constituents of the plasma are at rest, while for c = 1 they move with the speed of light.
In fig. 4 we show the scalar spectral density and its relative deviation in and offequilibrium in the infinite coupling limit for different values of c. We witness oscillations of the off-equilibrium spectral densities around their equilibrium values and as the shell approaches the horizon the amplitude of the oscillations decreases 4 . From this figure one can also see that high energetic modes are closer to equilibrium than the low energetic ones,  showing the usual top-down thermalization pattern with a dependence on the parameter c. The smaller the value of c the closer the quantity R is to its equilibrium value. Now we are ready to investigate the finite coupling corrections to the relative deviation of the spectral density. In fig. 5 the quantity R is displayed for the scalar channel for two different values of the coupling constant for a fixed position of the shell. For plasma constituents at rest, c = 0, R approaches a constant for large frequencies. But as c is increased, the fluctuation amplitude starts to grow at some critical value of the frequency ω crit , such that the higher energetic modes are further away from their equilibrium value than the low energetic ones, again indicating a weakening of the usual top-down thermalization pattern. This fits nicely into the picture obtained for the QNM where also the higher energetic modes show a stronger dependence on the finite coupling corrections. As the coupling constant is decreased the change of the behaviour shifts to lower ω crit . The results are again in accoreffects of the shell location are minor, we set the shell location in all our plots to the rather arbitrary value us = 1/2. dance with the calculation for the spectral density of the R-current correlator [45]. Since the dependence of the coupling constant is qualitatively the same in the shear and sound channel, they are only shown for one value of the coupling constant in fig. 6.
The parametric change of the relative deviation of the spectral density at finite coupling originates from the behaviour of C 0 and C 1 as defined in (35). As can be seen from fig. 7, C 0 always approaches zero for large frequencies being responsible for the top-down thermalization pattern at infinite coupling. On the other hand, the amplitude of γC 1 is constant for vanishing c and starts growing as c is increased. It is the interplay between C 0 and C 1 that is responsible for the observed pattern in R.
One might be worried about the above results since the quasi-static approximation and the finite coupling expansion employed have finite regions of applicability. However, since we are only taking snapshots, and not treating the dynamical problem, one can few the system as being very close to its initial evolution where the motion of the shell is guaranteed to be slow. After all, we are allowed to inject energy at an arbitrary scale and set the initial velocity of the shell to zero. In addition, from fig. 7 one might conclude that the strong coupling expansion breaks down once γC 1 becomes larger than C 0 . This however is not the case since the finite coupling corrections to the spectral density, which is the physical quantity, are at most of the order of 10%.

Conclusions
In the paper at hand we have studied the thermalization properties of an N =4 SYM plasma at finite 't Hooft coupling. First we analyzed how the plasma reacts to linearized perturbations as a function of the coupling constant through a QNM analysis. Then we studied the approach of the plasma to thermal equilibrium using the collapsing shell model of [24], working in the quasi-static approximation.
The flow of the QNM is depicted in figs. 1, 2, 3 and show a clear trend. As the coupling constant is decreased from the λ = ∞ limit the QNM bend upwards in the complex frequency plane. The increase of the imaginary part shows, according to equ. (40), that finite coupling corrections increase the life time of the excitations. In addition, higher energetic modes show a stronger dependence on the coupling corrections. We interpret this as a weakening of the usual top-down thermalization pattern. This is in accordance with [45] where a similar analysis was performed for virtual photons. This analysis is particularly useful because it is independent of the thermalization model being used and should therefore be of more general validity.
In the collapsing shell model, the deviation of the spectral density from its thermal limit was investigated. The results displayed in figs. 5 and 6 show that outside the limit of infinite coupling, the UV modes are further away from their thermal distribution than the IR modes, indicating a weakening of the top-down thermalization pattern. Both the spectral density and the flow of the quasinormal modes show qualitatively the same pattern as observed for real and virtual photons in [45], suggesting that the change in the thermalization pattern is of more general validity.
The above observations seem, however, to be in direct contradiction with a recent study of thermalization at finite coupling using the Vaidya metric [51], where the UV modes were seen to thermalize even slightly faster than in the infinite coupling limit. One possible explanation for the discrepancy between their work and the calculation for photons presented in [45] is an additional contribution of the Ramond-Ramond five form field strength that has to be added to the action if photons are considered [64]. The analysis presented in this paper shows that the spectral density of the energy momentum tensor, to which this term does not contribute 5 , exhibits the same behaviour as the spectral density for photons. Therefore the additional term for photons is not the source of the discrepancy. Another important difference between the calculation of [51] and the present one, as well as the one of [45], is that the latter two use the quasi-static approximation, i.e. work in the limit of a slowly moving shell, while the first employs the opposite Vaidya limit. In addition, the correlation functions studied in [51] are all so-called geometric probes, meaning that they are only sensitive to the γ-corrected metric and one need not consider fluctuation equa-tions. Which of these differences explain(s) the observed results remains, however, an open question, and a very important topic for further investigation.
For other future directions, it will be important to go beyond the case studied here, where we only take snapshots of the system, and study the time evolution of the correlators within the quasi-static approximation along the lines of [32]. Of course, in the long run the goal is to consider finite coupling corrections to spectral densities in a thermalizing system without using the quasi-static approximation at all.