Quasinormal modes and thermalization in Improved Holographic QCD

Following a series of similar calculations in simpler non-conformal holographic setups, we determine the quasinormal mode spectrum for an operator dual to a gauge-invariant scalar field within the Improved Holographic QCD framework. At temperatures somewhat above the critical temperature of the deconfinement transition, we find a small number of clearly separated modes followed by a branch-cut-like structure parallel to the real axis, the presence of which is linked to the form of the IHQCD potential employed. The temperature dependence of the lowest nonzero mode is furthermore used to study the thermalization time of the corresponding correlator, which is found to be of the order of the inverse critical temperature near the phase transition and decrease slightly faster than $1/T$ at higher temperatures.


I. INTRODUCTION
The quantitative description of equilibration dynamics in gauge theories is a notoriously complicated problem, the study of which is motivated by the desire to understand the time evolution of systems such as ultrarelativistic heavy ion collisions or the expanding early universe [1]. In the former context, where the dynamics is governed by the strong nuclear force, it becomes essential to describe thermalization away from the tractable weak coupling regime [2]. For this reason, the gauge/gravity duality has become a standard tool in the field [3], and a considerable amount of work has indeed been devoted to studing equilibration in strongly coupled N = 4 Super Yang-Mills (SYM) theory by mapping the process to the gravitational dynamics of black hole (BH) formation in AdS spacetime [4]. Important milestones in this line of research include e.g. the successful description of shock wave collisions, the observation of a rapid onset of hydrodynamic behavior, and the subsequent understanding that this does not require a full thermalization or even isotropization of the system (see e.g. [5][6][7][8] and references therein).
While a considerable amount of physical insight has been gained from studies of the SYM theory in its infinitely strongly coupled regime, it is clear that the application of the results to real QCD requires an understanding of the quantitative effects of conformal invariance breaking as well as the finiteness of the gauge coupling [9][10][11]. These issues have indeed been studied in many different contexts, in some cases involving even shock wave collisions [12][13][14], but often by restricting to the study of the quasinormal mode (QNM) spectra of BH solutions in different spacetimes [15][16][17]. This is understandable, since the QNM spectrum is a key quantity determining a holographic system's equilibration from small perturbations, as it is dual to the pole locations of the corresponding retarded Green's function on the field theory side. In fact, as argued in [18], the thermalization time of a Wightman function in the field theory dual of a five-dimensional system undergoing gravitational collapse should be inversely proportional to the imaginary part of the lowest QNM determined for the same correlator in thermal equilibrium.
In the paper at hand, we continue work towards understanding equilibration in a QCD-like plasma, but staying on the level of QNMs. In particular, we study the retarded Green's function of an operator dual to a gauge-invariant scalar field in Improved Holographic QCD (IHQCD) -a five-dimensional bottom-up model constructed to mimic large-N c non-supersymmetric Yang-Mills (YM) theory [19][20][21][22][23][24][25]. This work can be considered a direct continuation of similar studies of QNMs in related but often simplified non-conformal models [26][27][28][29] (see also [30]), differing from them through the use of the more realistic IHQCD potential introduced in [25]. Utilizing the fact that the logarithmic running of the QCD gauge coupling is built into IHQCD, we are able to study the QNM spectrum -and thus the thermalization time scale -away from the infinitely strongly coupled regime. As discussed in detail below, our results point towards the emergence of a branch-cut-like structure on the complex frequency plane in the limits of increased temperature and QNM mode number, which both probe the more weakly coupled UV limit of the theory. Our paper is organized as follows. In Section II, we briefly review the model we work in as well as our strategy for determining its bulk thermodynamic properties as well as QNM spectra. The results of this investigation are then presented in Section III, while Section IV is devoted to a detailed discussion of our findings. Some computational details useful in the determination of the QNM spectra are finally relegated to Appendix A.

II. HOLOGRAPHIC MODEL
In this section, we discuss first the construction and basic properties of the holographic model we work with, IHQCD, and after this present some details concerning the determination of different thermodynamic quantities as well as QNM spectra within this framework.

A. Basic equations
Improved holographic QCD is defined by the five-dimensional gravity action together with the metric ansatz [19,20], where z is a radial coordinate chosen in such a way that the UV boundary resides at z = 0. Here, the scalar field φ(z) and the functions b(z) and f (z) are in turn determined from Einstein equations, which in this case take the forms with the dot denoting a rerivative with respect to z. In addition, both the z-coordinate and the function b(z) are closely related to the renormalization scale on the field theory side, so that we may write with the source term for λ scaling like N c g 2 . Finally, the dynamics of the model is required to reproduce the logarithmic running of the coupling at z → 0, where one can choose to work either to one-or two-loop order. This determines the unit of energy Λ, which for the remainder of this section we choose as Λ = 1.
The single most important quantity determining the dynamics of the model is its potential V (λ). Its behavior at small z is dictated by the running of the coupling through Eq. (5), while in the far infrared, i.e. for large z, the behavior of V (φ) is constrained by requiring confinement [20]. The potential we employ reads [25] V (λ) = 12 1 + 88λ 27 which in addition to satisfying the infrared and ultraviolet constraints has been fitted to lattice data for large-N c pure YM theory [31].
As we shall demonstrate in the next section, the different parameters in the potential have been chosen so that the order of the phase transition and the equation of state are reproduced to sufficient accuracy. Varying them allows one to describe phase diagrams with higher order or even crossover transitions [32]. More generally, tuning the functional form of the dilaton potential V (λ) provides holographic realizations of theories with different types of renormalization group flow both in the infrared [26,33] and ultraviolet [34] (see also the discussion in [35]).
As a final remark, we note that the model can be extended by adding another scalar field sourcing theqq operator and allowing for a full treatment of fermion backreaction, resulting in a model commonly referred to as V-QCD [22]. The phase diagram of this extended model has been studied at both at finite temperature [23] and density [24], and it has been extensively used in particular in the description of the dense nuclear and quark matter found inside compact stars [36][37][38].

B. Numerical solution for thermodynamics
To find a numerical solution to Eqs. (4), it is useful to first define the new variable so that the equations are transformed to the first order formṡ that we can approach in a relatively straightforward manner (see e.g. [26] for details).
The horizon values of the fields b(z), f (z) and W (z), with the horizon located at z = z h , are then obtained requiring finiteness and that the leading dependence of λ on the energy scale in the UV be equivalent with the known behavior of the YM gauge coupling at two-loop order.
After generating a family of solutions parametrised by the values λ h ≡ λ(z h ), the thermodynamic behavior of the model can be determined from the relations = T s − p.
The overall scale of thermodynamic quantities is affected by the choice of 4G 5 in the above equations. In principle, this choice should be made by matching with the noninteracting Bose-Einstein limit at asymptotically high temperatures, but to follow standard conventions in the field, we instead choose the scale by optimizing the model's fit to lattice data at temperatures only slightly above the critical temperature T c of the first-order deconfinement transition. This leads to overshooting the expected high-T behavior of thermodynamic quantities by ca. 30%, but as we will see in the following section, produces lower-temperature thermodynamics in excellent agreement with lattice predictions.

C. Quasinormal modes
Being equipped with the solved gravity background as well as with the potential fitted to lattice thermodynamics, we can now proceed to study the QNM spectra predicted by IHQCD. Here, we specialize to a field theory operator dual to a gauge-invariant scalar fluctuation φ(ω, z), noting that the for q = 0 the equation of motion for this operator takes the formφ To aid the forthcoming analysis, we transform this equation into a Schrödinger-type form, with ω 2 playing the role of an eigenvalue. This is achieved by introducing the new variable and defining ψ(u) = √ b 3 φ(u), whereby the equation of motion becomes Solving for the QNMs from this equation is in principle straightforward: the equation is linear, and its solution (for a fixed ω) therefore fully determined by two complex numbers. One of them corresponds to the normalization, while the other can be determined by requiring the solution to be ingoing at the horizon, i.e. ensuring that ψ be proportional to e +iωu at large u. To obtain the QNMs, we then simply need to determine those values of ω, for which the solutions satisfying this boundary condition are normalizable when u → 0.
To implement this boundary condition in the numerical calculation it is convenient to find the solution in the infrared analytically and match it with the numerically determined ultraviolet solution. The details required in implementing the boundary condition are given in Appendix A. A practical complication in the computation arises from the fact that standard spectral methods rely on the ability to expand solutions around u ≈ z = 0 in power series, while in the gravity model we are working with these expansions also contain logarithmic terms. This issue can, however, be resolved through the application of more robust numerical methods to remove the non-normalizable component from the solution. Here, it is crucial to recall that in a numerical calculation the solutions always contain a small part of the non-normalizable solution. However, if the solution is close to the correct one, the divergence due to the non-normalizable solution appears only when u is very small. Then, the method introduced in [26] for finding the QNMs proceeds as follows: for a trial value of ω, find the numerical solution towards the boundary and determine the minimum of |ψ(u)|, denoting its location by u min . The desired QNM is then approximately the value which minimizes u min (ω).

A. Thermodynamic properties
We begin our discussion from thermodynamic quantities, which we however cover only rather briefly, as they have been extensively considered in the literature (see e.g. [23]). Indeed, our focus is mainly on those numerical results that serve as important consistency checks for our model and establish a connection between its finite-temperature phase structure and its spectrum as a function of T .
Specifically, we want to compute the free energy, i.e. the pressure p(T ), as defined in Eqs. (14) in terms of the functions b(λ h ) and T (λ h ). Here, b(λ h ) is a monotonous function, while T (λ h ) first decreases with increasing λ h but then starts to increase; see the illustration in Fig. 1 (left). The behavior in the UV, or small λ h , corresponds to the weak-coupling limit, i.e. to large temperatures, while the domain towards the IR, where temperature increases with λ h , is unstable. This can be seen e.g. by computing the speed of sound squared: it turns out that this quantity is proportional to −T (λ h ), so we must clearly require T (λ h ) < 0.
Evaluating the pressure from high to low temperatures, we next determine the critical  temperature as the point where the pressure becomes negative. From this condition, we find that which allows us to express all dimensionful quantities in units of T c . Our results for the thermodynamics of the model are displayed and compared to the lattice data of [31] in Fig. 1 (right). From here, we observe that our results for the three key quantities -the pressure p, energy density and trace anomaly − 3p -are in very good agreement with the lattice calculation. In particular, the slopes of all three functions are accurately reproduced at temperatures somewhat above the critical one, which is a strong indication that our choice for the dilaton potential in Eq. 7 was indeed a reasonable one. As mentioned in Sec. II B, the price to pay for treating the overall normalization of the pressure as a free parameter is that our results eventually overshoot the expected hightemperature limits of these quantities by some 30%. At the same time, it interestingly turns out that with our normalization convention the UV behavior of energy momentum tensor correlators becomes very accurately reproduced [39].

B. Spectrum of quasinormal modes
Moving next on to the QNM spectra, we first note that at zero temperature, corresponding to λ h → ∞ and f = 1, the potential entering the Schrödinger equation is well approximated by (see [26]) from which one obtains a bound state spectrum with masses m n = 2 √ 2 + n, n = 0, 1, 2, . . .
As λ h is lowered, a numerical calculation performed on the unstable branch of the theoryshows how the potential is perturbed away from the harmonic-oscillator-like form as depicted in the left panel of Fig. 2 (note, however, that the instability analyzed in [27,40] is not present for the potentials we consider). The corresponding spectrum of QNMs, shown in the right panel of Fig. 2, displays a clear and expected pattern: the spectrum moves away from the real axis, with the states becoming broader. It should be highlighted, however, that these results, derived on the unstable branch, can not be straightforwardly related to the physics of the YM theory. Nevertheless, these results connect smoothly with the quasinormal modes determined on the high temperature branch where the system enters the stable high-temperature phase. There the states are observed to become even broader, eventually melting away.
As our physical interest lies in the description of thermalization dynamics at temperatures somewhat (but not excessively) above the critical one, we next specialize to the stable branch of the theory and analyze in more detail the QNM spectrum for temperatures T c < T < 3T c . Of particular interest here is the lowest QNM and specifically its imaginary part, which has been argued to be inversely proportional to the equilibration time for the correlator in question [18]. For temperatures ranging from T c to 3T c , we first display the lowest QNMs on the complex-ω plane in the left panel of Fig. 3, observing that the lowest mode roughly follows the trend Im ω ≈ −iRe ω and ω ≈ 2πT (1 − i). The thermalization time obtained from the lowest mode via the relation [18] is shown as a function of temperature in the right panel of Fig. 3. It is observed to decay slightly faster than 1/T , with the scale being of order τ th 0.5/T c , which is a phenomenologically very reasonable result.
Finally, an interesting observation can be made about the fate of higher QNMs on the stable BH branch as a function of increasing temperature. In our numerical calculation, we find that the string of clearly separated, individual QNMs terminates at a structure that one is tempted to interpret as an emerging branch cut on the complex ω-plane; cf. Fig. 4. This structure always corresponds to a constant value of Im ω, with the number of accessible separate QNMs decreasing with increasing T . This behavior is consistent with the existence of quasistable states as obtained from Eq. (19), with the potential evolving as shown in Fig. (2) (for the unstable branch). We have checked that the existence of this structure is independent of whether we use a leading or next-to-leading order expansion for the analytic solution of the Schrödinger equation in the infrared.

IV. CONCLUSIONS
In the holographic study of gauge theory equilibration, one can identify two somewhat distinct long-term goals. On one hand, there have been frequent attempts to make the physical setting under study more closely resemble realistic heavy ion collisions by considering the collision of shock waves of finite thickness, transverse size, and anisotropy  [5,8,41,42]. At the same time, many groups have been working on modifying the gravity background in the shock-wave setting to have the dual gauge theory feature broken conformality or supersymmetry, or even non-infinite coupling strengths [10,[12][13][14]43]. What is, however, common to both of these lines of work is that one usually works in firmly top-down settings, typically starting from the conformal and supersymmetric N = 4 SYM theory and departing from it one step at a time, which makes progress towards the description of a truly QCD-like theory rather slow.
Independently from the study of shock-wave collisions and other developments in applied top-down holography, systematic efforts towards constructing a bottom-up model to closely resemble non-supersymmetric Yang-Mills theory and QCD have culminated in the development of the Improved Holographic QCD and V-QCD frameworks that have been demostrated to faithfully reproduce the equilibrium thermodynamic properties of the full theories in the large-N c limit [19,20,23,24]. Considering the technical complexity of building a computational framework for studying shock-wave collisions in this type of a setting, it is clearly worthwhile to investigate equilibration within the IHQCD model in a somewhat more modest way: through the study of quasinormal modes that describe the response of the dual field theory system to small departures from equilibrium.
In the paper at hand, we have determined the QNM spectrum for one particular gauge theory operator in the IHQCD framework, dual to a scalar field on the gravity side, for a range of temperatures from T c up to ca. 3T c . In particular, we have studied the temperature dependence of the lowest nonzero QNM to obtain an estimate for the thermalisation time of the corresponding correlator, finding phenomenologically very reasonable results. In addition, we have observed that the number of clearly separate QNMs decreases with increasing temperature -a result we have been able to link to the broadening and melting of the corresponding states. Instead of a clean spectrum of higher individual QNMs separated from each other by the expected 2πT , we have witnessed the emergence of a linear structure parallel to the real axis of the complex frequency plane, which may point towards to the presence of a brach cut. In this context, it is worth pointing out that while we have only shown results for the stable BH branch in Fig. 4, the same qualitative conclusions apply to the unstable branch as well, with the main difference being the somewhat higher number of individual QNMs there.
It is clearly very interesting to contrast our findings with the results of similar exercises carried out in the past. Our present work can be regarded as a rather direct continuation of a series of works by other groups, in which somewhat simpler versions of IHQCD (or other closely related models) have been considered [26][27][28][29], although it should in addition be recalled that some of these works considered correlators different from the one studied by us. Of these papers, Alanen et al. [26] and Janik et al. [27] both employed the IHQCD framework but with simpler choices of the potential, with e.g. [27] not enforcing the condition that the gauge coupling should run logarithmically in the UV. Neither of these references reported the existence of a branch-cut like structure, which is likely related to the IR behaviors of the employed potentials. In this respect, a very interesting point of comparison is the work of Betzios et al. [28], which considered a different model of Einstein-dilaton gravity, dual to a Chamblin-Reall plasma in the IR and having the usual AdS form in the UV. Similarly to us, they witnessed the emergence of a branch cut in the critical limit, albeit exactly on the real axis. It would clearly be very interesting to analyze the reason for this qualitative similarity in a more explicit manner, but we leave this for future work.
Finally, we note that while in the calculation presented above we have worked in the context of pure gauge theory, our results can be extended to QCD in the Veneziano limit by employing the V-QCD model [24]. Performing a similar study of QNM spectra in this framework and thereby analyzing the effect of dynamical quarks on equilibration would clearly be a very interesting, albeit technically more demanding exercise to carry out.

Acknowledgments
We thank Niko Jokela and Matti Järvinen for enlightening discussions and useful comments on the manuscript. The work has been supported by the European Research Council, grant no. 725369, and by the Academy of Finland grants no. 1322507 and 1310310.

Appendix A: Asymptotic solutions
In this Appendix, we briefly discuss the asymptotic expansions of the bulk fields, which turned out to be very useful in the numerical determination of the QNM spectra. In this context, we note that when dealing with the Schrödinger-type equation (19), it is useful to first expand the potential at u → ∞ and solve the equation analytically in this limit. However, this is done more easily after introducing a new variable A such that b(z) = exp(A(z)) and solving the equation near the limit A → A h . For this purpose, we also need to introduce the function so that the fluctuation equation (19) obtains the form To proceed from here, we note the relation wheref (A) ≡ f (A) A and we have set A h = 0 for convenience. Observing thatf is regular and non-zero at A = A h = 0, we can isolate the divergence by integrating by parts where the prime denotes differentiation with respect to A and we have used the fact that q andf grow more slowly than the inverse of log(A)e −A at large A. Next, we define and write where Note that if h(A) is real-analytic, the logarithm appears only in the leading term of the expansion.
In order to write the Schrödinger equation in a more useful form, we have to invert the above relation, i.e. find the function A(u). To this end, we writeû = u − u 0 to simplify the notation. To leading order, we can immediately solve which is very small for largeû, as h 0 < 0.
To get the next-to-leading order term, we attempt to substitute the leading term to Eq. (A7) and expand to order A. The resulting error term is h 1 eû h 0 + O(A 2 ), which leads to the ansatz A(û) = eû h 0 +u 1 eû h 0 (A10) for some number u 1 . Using the fact that eû /h 0 is small, we can expand where the last term is of order A 2 , leading to the ansatz Inserting this again to Eq. (A7), we get which implies It is clear that from here the series would continue further in powers of eû /h 0 . The numerics provide us the Schrödinger potential in terms of A. Writing and inserting the expansion above gives The corresponding equation can be solved analytically, giving where U is Tricomi's confluent hypergeometric function, and L is is the generalized Laguerre polynomial. We have defined the following shorthand notations: Furthermore, V 1 and V 2 are the coefficients appearing in (A16) and C 1 and C 2 are constants of integration. Finally, we note for reference the relations