Thermodynamics, transport and relaxation in non-conformal theories

We study the equilibrium and near-equilibrium properties of a holographic five-dimensional model consisting of Einstein gravity coupled to a scalar field with a non-trivial potential. The dual four-dimensional gauge theory is not conformal and, at zero temperature, exhibits a renormalisation group flow between two different fixed points. We quantify the deviations from conformality both in terms of thermodynamic observables and in terms of the bulk viscosity of the theory. The ratio of bulk over shear viscosity violates Buchel's bound. We study relaxation of small-amplitude, homogeneous perturbations by computing the quasi-normal modes of the system at zero spatial momentum. In this approximation we identify two different relaxation channels. At high temperatures, the different pressures first become approximately equal to one another, and subsequently this average pressure evolves towards the equilibrium value dictated by the equation of state. At low temperatures, the average pressure first evolves towards the equilibrium pressure, and only later the different pressures become approximately equal to one another.


Introduction
The understanding of the out-of-equilibrium dynamics of matter is an important challenge ubiquitous at all energy scales. A particularly interesting case is the understanding of these dynamics in strongly coupled systems. Examples include strongly correlated electrons, cold atoms and the small drops of Quark-Gluon Plasma (QGP) formed in relativistic colliders such as RHIC or the LHC. The latter case motivates the study of the relaxation process in strongly coupled non-abelian field theories. The gauge/string duality provides a fascinating tool to address this problem in a wide range of theories.
The duality has already provided insights into the dynamics of strongly coupled, deconfined, nonabelian matter of relevance for the heavy ions programme (see e.g. [1] and references therein). The study of the off-equilibrium dynamics of Conformal Field Theories (CFT), most notably N = 4 super Yang-Mills (SYM) theory, has shown that hydrodynamics is a much better approximation to the evolution of this type of matter than ever thought before. Indeed, examples based on flow motions imposed by symmetries [2,3] or by explicit simulations of the collision dynamics [4][5][6][7] have shown that hydrodynamics provides a good approximation to the complete evolution of the system at time and distance scales as small as a fraction of the (local) inverse temperature of the system. This occurs even in situations in which gradient corrections to the hydrodynamic stress tensor are large, extending the applicability beyond a simple gradient expansion (see also [8]). This observation has led to the coining of the term "hydrodynamisation" to refer to the process by which a system comes to be well described by hydrodynamics, in order to differentiate this process from (local) thermalisation. This observation, first made in holographic computations, has now been noted in Boltzmann equation-based analysis of out-of-equilibrium dynamics, when the strength of the coupling in the collision kernel is extrapolated to large values [9]. The success of hydrodynamics to capture the evolution of out-of-equilibrium matter may be at the origin of the strong collective behaviour observed in very small systems, such as Au-Au collisions at RHIC [10][11][12], Pb-Pb [13][14][15][16], p-Pb [17][18][19] and p-p [20] collisions at the LHC. The holographic analysis of collisions of small systems [21,22] supports this viewpoint.
In a CFT the vanishing of the trace of the stress tensor implies that the equation of state, namelȳ is the average pressure, is fixed by symmetry. As a consequence, the equation of state is always obeyed both in and out of equilibrium. We emphasize that the equation of state fixes only the average pressure in terms of the energy density, but not the individual pressures. For this reason the relaxation towards equilibrium in a CFT typically involves "isotropization", namely the process by which the different pressures become approximately equal to one another. The applicability of the gauge/string duality is not restricted to CFTs. By now infinite families of non-conformal examples are known. One of the main new features in these theories as compared to their conformal cousins is that new channels exist for the relaxation of the out-of-equilibrium matter. In particular, in non-conformal theories the equation of state is not fixed by symmetry. As a consequence, out of equilibrium the energy density and the average pressure may fluctuate independently. Therefore, the relaxation towards equilibrium in these theories involves the evolution of the energy density and the average pressure towards asymptotic values related to one another by the equation of state (EoS). When this happens we will say that the system has "EoSized" and we will refer to this process as "EoSization".
Another important motivation for studying non-conformal theories is the connection with hot Quantum Chromodynamics (QCD) and heavy ion collisions. As is well known, QCD is a non-conformal theory even in the limit of vanishing quark masses. State-of-the-art determinations of the QCD equation of state via lattice QCD [23,24] show that, in equilibrium, the trace of the stress tensor normalised by the enthalpy attains values of order one close to the QCD transition. At high temperature this ratio quickly approaches zero, indicating that QCD behaves as an almost-conformal theory in this regime. However, the experimental exploration of the QCD phase diagram via high-energy heavy ion collisions can only reach temperatures a few times larger than the critical temperature. Even though most central, top-energy LHC collisions lead to initial temperatures well into the quasi-conformal regime, the subsequent evolution and cooling of the QGP after production spans all temperature regimes, including those in which non-conformal effects are maximal. In fact, recent attempts for high-precision extraction of the shear viscosity of the QGP have highlighted the need to include the bulk viscosity of the plasma, which is a purely non-conformal effect [25]. Furthermore, off-central collisions both at the LHC and RHIC, as well as lower-energy collisions as those explored at the RHIC energy scan, produce a QGP with a smaller initial temperature. Similarly the apparent success of hydrodynamics in smaller systems such as p-Pb [26] and p-p [27,28] collisions indicate the need to study the properties of deconfined but cooler QCD plasma, where non-conformal effects become significant (see [29] and references therein for a recent review on the hydrodynamic modelling of heavy ion collisions).
In order to study non-conformal theories in a holographic setup we will consider a five-dimensional bottom-up model that nevertheless shares many qualitative features with top-down string models. Specifically, our model is dual to a four-dimensional gauge theory that, at zero temperature, flows from an ultraviolet (UV) fixed point to an infrared (IR) fixed point. This renormalisation group (RG) flow is dual on the gravity side to a domain-wall geometry that interpolates between two AdS spaces. The reason why we require that the flow approaches a fixed point in the UV is that this is the situation in which the holographic duality is best understood. The reason for the IR fixed point is that this guarantees that the zero-temperature solution is smooth in the deep IR. The flow is triggered by a source Λ for a relevant, dimension-three operator in the UV. We will see that this simple model exhibits a rich phenomenology. In particular, we will study the relaxation of small-amplitude, homogeneous perturbations by computing the spectrum of quasi-normal modes (QNM) with zero spatial momentum. We will see that the dominant channel for relaxation in this approximation depends on the value of the ratio T /Λ, with T the temperature of the system. At small T /Λ the system first EoSizes and subsequently isotropises. In contrast, at large T /Λ the order in which these two processes take place is reversed. Although our calculation is done at zero spatial momentum we will argue that, actually, the ordering above is still valid for long-wave-length fluctuations with k T .
Previous analyses addressing the near-equilibrium properties of strongly coupled non-abelian plasmas include [30][31][32][33][34][35][36][37][38]. In particular, the last reference in this list appeared while this paper was being typeset and has some overlap with our observations concerning the different relaxation channels. This paper is organised as follows. In Sec. 2 we introduce the holographic model and discuss its vacuum properties. In Sec. 3 we study black brane solutions and extract from them the equation of state and the viscosities of the model. In Sec. 4 we study the relaxation of small excitations of the system by computing the QNM spectrum of the black branes at different temperatures and zero spatial momentum. Finally, in Sec. 5 we discuss our main findings and place them in the context of the hydrodynamisation of non-abelian plasmas.

A non-conformal holographic model
The holographic model that we will consider consists of five-dimensional Einstein gravity coupled to a scalar field with a non-trivial potential: where κ5 is the five-dimensional Newton constant. For specific forms of V (φ), this action may be viewed as a consistent truncation of five-dimensional N = 8 supergravity. In this paper we will consider a bottom-up model by choosing a potential that is particularly simple and yet shares some of the qualitative properties of these top-down potentials. In particular, we will choose V (φ) to be negative and to possess a maximum at φ = 0 and a minimum at φ = φM > 0. Each of these extrema yields an AdS solution of the equations of motion with constant φ and radius L 2 = −3/V . In the gauge theory each of these solutions is dual to a fixed point of the RG with a number of degrees of freedom N 2 proportional to L 3 /κ 2 5 . 1 We will be interested in domain-wall solutions interpolating between these two AdS solutions. In the gauge theory, these are dual to RG flows from the UV fixed point at φ = 0 to the IR fixed point at φ = φM . The problem of finding those solutions is significantly simplified if the potential can be written globally in terms of a superpotential, W , as In this case, vacuum solutions to the Einstein equations can be easily found. Parametrizing the metric as the solution of the back-reacted gravitational problem is reduced to the first-order equations [39] dA We will choose a simple superpotential characterised by a single parameter, φM , which together with eq. (2.2) yields the potential (2.6) Note that both the superpotential and the potential have a maximum at φ = 0 and a minimum at φ = φM . This choice leads to three important properties of the associated vacuum solution. First, the resulting geometry is asymptotically AdS5 in the UV with radius L, since V (0) = −3/L 2 . Second, the second derivative of the potential at φ = 0 implies that, in this asymptotic region, the scalar field has mass m 2 = −3/L 2 . Following the standard quantisation analysis, this means that, in the UV, this field is dual to an operator in the gauge theory, O, with dimension ∆UV = 3. Third, the solution near φ = φM is again AdS5 with a different radius In this region the effective mass of the scalar field differs from its UV value and it is given by As a consequence, the operator O at the IR fixed point has dimension (2.9) 1 In the case of N = 4 SYM the precise relation would be L 3 /κ 2 5 = N 2 /4π 2 .
To summarize, the vacuum solution describes a RG flow from an UV to an IR fixed point with a smaller number of degrees of freedom, as indicated by the fact that LIR < L. We see that changing φM has two main effects. First, as φM increases the difference in degrees of freedom between the UV and the IR fixed points increases. Second, the dimension of the scalar operator at the IR fixed point decreases with increasing φM , reaching the marginal dimension ∆IR = 4 at φM → ∞. However, in this limiting case the IR fixed point disappears and the background solution becomes singular, as is evident from the fact that the effective AdS radius goes to zero as φM → ∞.
Our simple choice of the superpotential allows us to determine analytically the vacuum solution for arbitrary φM . Solving eq. (2.4), we obtain where Λ is an arbitrary constant that controls the magnitude of the non-normalizable mode of the scalar field. As we will see, in the dual gauge theory side, Λ is identified with the source of the dimension-3 operator O. The presence of this source breaks conformal invariance explicitly.
Noticing that the small field behaviour of the superpotential eq. (2.5) is identical to that of the GPPZ flow [40], we can readily determine the vacuum expectation values (VEV) of the stress tensor and the scalar operator. We begin by expanding the metric and the scalar field in powers of u = Le −r/L in the u → 0 limit. Following [39], we write the 5-dimensional metric in the form and we write the power expansion coefficients of the metric and the scalar field as 14) The expectation values of the field theory operators are then given by To arrive at these expressions we have chosen the superpotential as a counterterm to regularise the on-shell action, which is possible because in our model the superpotential corresponds to a deformation of the gauge theory as opposed to a VEV [41]. We emphasize that these expressions are valid even if the metric gµν does not posses the full Poincaré symmetry but only rotational and translational invariance along the gauge theory directions, as will be the case for the black brane geometries that we will study in the next section. As expected, Eqs. (2.15) and (2.16) imply the Ward identity for the trace of the stress tensor Eqs. (2.10) and (2.11) determine the VEVs in the vacuum of the theory. Let us define the energy density and the pressure p as the diagonal components of the expectation value of the stress tensor, T µν = Diagonal { , p, p, p}. The near boundary behaviour of φ, eq. (2.11), leads to 18) which implies that in the vacuum Note that the explicit breaking of scale invariance means that the trace of the stress tensor is non-zero as an operator. However, the VEV of this operator vanishes in the vacuum state, as implied by trace Ward identity (2.17) together with the fact that O = 0 in the vacuum for our choice of renormalisation scheme. It should be emphasized that even though the trace Ward identity (2.17) is scheme-independent, the individual vacuum expectation values of the trace of the stress tensor and of the scalar operator do depend on the renormalisation scheme. In the model we study here the only scheme ambiguity corresponds to a term of the form Λ 4 ηµν in the expectation value of the stress tensor, accompanied by a term of the form Λ 3 in the expectation value of O, with the relative coefficient such that (2.17) is preserved.
To estimate at which scale non-conformal effects become important, let us perform a change of variables in the holographic direction, which explicitly exploits the relation between the dynamics in the bulk with the physics at different scales in the field theory. Denoting the coordinate by z, we write the metric as with L eff a non-trivial function of z such that L eff (0) = L and L eff (∞) = LIR. In this set of coordinates, at least in the two asymptotic conformal regions, the coordinate z is related to the energy scale, Q, in the gauge theory through z ∼ 1/Q. The relation between z and u is given by 21) and the function L eff is given by In Fig. 1 we show the ratio L eff /L as a function of z for several different values of the parameter φM controlling the physics of the model. We see that the system behaves approximately conformally up to scales of order z ∼ Λ. At this scale, the metric starts to deviate significantly from that of AdS5, and L eff decreases as a function of z. Sufficiently deep in the IR, L eff approaches LIR and the system behaves again as approximately conformal. However, the scale at which this transition occurs depends significantly on the model parameter φM ; as φM increases, the function L eff approaches its asymptotic value more slowly. This different rates at which the IR fixed point is approached have consequences for the finite-temperature behaviour of the dual gauge theory, as we will see in the next section.

Thermodynamics and transport
We will now explore the thermal physics of the gauge theory dual to the gravitational model described in the previous section. 2 To do so, we will search for black brane solutions of the action (2.1). We will follow the method of the master function, introduced in ref. [45], to which we refer the reader for details. 3 Since for the background solution (2.11) the scalar field is a monotonic function of u, we may use the scalar field as a coordinate and express the metric as with h(φ) vanishing at φ = φH , the value of the scalar field at the horizon, i.e. h(φH ) = 0. The region outside the horizon corresponds to 0 < φ < φH . For later convenience, we have expressed the metric in Eddington-Finkelstein form. With this ansatz, Einstein's equations take the form A solution to these equations may be found in terms of a master function G(φ) defined as Manipulating the set of equations (3.2), a non-linear equation for G was found in [45]: Close to the boundary, φ → 0, the solution of this equation behaves as with ∆ the scaling dimension of the dual operator. With our choice of potential (2.6) we have ∆ = 3. Using eq. (3.2), the different metric coefficients are given by . In these expressions, the constants of integration are fixed by requiring that, close to the boundary, the metric and scalar field may be expressed as in eq. (2.13) and eq. (2.14) . At the horizon, the condition h(φH ) = 0 together with the last two equations in (3.2) fix the value of G(φH ). Starting from this fixed value, a power series solution close to the horizon may be found as From these metric coefficients, we can extract the Hawking temperature T and the entropy density s of the black brane: The relation of the different metric coefficients with the master function leads to the following form for the temperature and entropy of the thermal state: These expressions are well suited for the determination of these two quantities from the numerical evaluation of the master equation (3.4).
In Fig. 2 we plot the dimensionless quantity as a function of the inverse temperature for two different values of φM . Since the theory is conformal both at the UV and at the IR, the high and low temperature behaviour of the entropy density must coincide with that of a relativistic conformal theory and scale as T 3 . In particular, for a relativistic CFT, s/T 3 is proportional to the number of degrees of freedom in the theory, which for an SU (N ) gauge theory with matter in the adjoint representation scales as N 2 . For example, for N = 4 SYM but the precise coefficient depends on the specific theory. In terms of the parameters of the dual gravity description this quantity becomes s In our bottom-up setup, the above argument allows us to define the number of degrees of freedom at the fixed points holographically in terms of the effective AdS radius. In particular, the quantity sR should approach 1 at high temperature and (LIR/L) 3 at low temperature, which is confirmed by the plots in Fig. 2.
Using standard thermodynamic relations and the fact that in our renormalisation scheme the vacuum pressure is zero we can determine the pressure and the energy density of the thermal system through Since the theory is not conformal, the trace of the stress tensor in the thermal ensemble does not vanish. Using the Ward identity (2.17), the energy density, the pressure and the scalar condensate at non-zero temperature are related through The thermal expectation value O T may be determined from the normalisable mode of the scalar field in the thermal background via eq. (2.14). Since at T = 0 the scalar VEV vanishes (see eq. (2.19)) this relation implies that = 3p, as expected from the fact that the IR theory is conformal. At T > 0, however, O T = 0, as shown in Fig. 3, and the expectation value of the trace of the stress tensor does not vanish. Note that, unlike at low temperatures, at which O T depends on φM , at high temperatures O T becomes independent of φM . This is easy to understand from the gravitational computation. At high temperatures the value of the scalar field at the horizon is small and, therefore, the physics is sensitive only to the smallfield behaviour of the scalar potential, which is independent of φM . In this limit, the plots in Fig. 3 show that the VEV scales as O T ∼ ΛT 2 . Despite the fact that the trace of the stress tensor at high temperature does not vanish, the theory does behave as a conformal theory. From the gauge theory viewpoint this may be understood from the relative magnitude of the trace of the stress tensor compared to the energy density or the pressure: while at large T the latter quantities scale as T 4 , the trace only grows as T 2 . In Fig. 4 we show the temperature dependence of the ratio of the stress tensor to the enthalpy, which in the thermal-QCD literature is sometimes referred to as the interaction measure. As anticipated, both at low and high temperatures this ratio vanishes, indicating that the theory becomes effectively conformal in these limits. At intermediate temperatures, the value of I is non-zero and depends on φM . As inferred from the behaviour of the entropy, the larger φM the larger the deviations from conformality in the thermodynamic properties of the theory. Because of this behaviour we may use I as a measure of the non-conformality of the theory.
Another way to quantify the non-conformal behaviour of the thermodynamics of the dual theory is the value of the speed of sound. Using thermodynamic identities, the square of the speed of sound may be determined from the inverse of the logarithmic derivative of the entropy, In Fig. 5 we show the temperature behaviour of the deviation of cs from its conformal value, cs = 1/ √ 3, for different values of φM . The qualitative behaviour of this quantity is very similar to that of I. Both at high and low temperatures, the speed of sound approaches its conformal value. At intermediate temperatures we have c 2 s < 1/3 and the deviation from the conformal value grows with φM . The non-conformal behaviour already observed in the equation of state of the system is also reflected in the transport properties of the dual gauge theory plasma. Since this is isotropic, at leading order in gradients transport phenomena are controlled by only two coefficients, the shear viscosity η and the bulk viscosity ζ. Because of the universality of the shear viscosity to entropy ratio [46] in all theories with a two-derivative gravity dual, we have that this ratio in our model takes the same value as in the conformal N = 4 theory, i.e. η/s = 1/4π. In contrast, the bulk viscosity, which would vanish identically in a CFT, is non-zero in our model. Following 4 [47] we determine the bulk viscosity by studying the dependence of the entropy on the value of the scalar field at the horizon, 5 (3.20) The temperature dependence of this ratio is shown in Fig. 6 for different values of φM . The behaviour of this ratio is very similar to that of the interaction measure and the speed of sound: both at low and high temperatures the ratio of the two viscosities vanishes, while at intermediate temperatures attains φM -dependent values that grow with φM . As in the case of − 3p and the interaction measure, the fact that the ratio of viscosities vanishes at high temperatures does not imply that the bulk viscosity itself vanishes. In fact, we have checked numerically that at high temperatures the bulk viscosity scales as ζ ∼ Λ 2 T . Nevertheless, the fact that the ratio of viscosities approaches zero shows that transport is effectively conformal, since the contribution to the hydrodynamic stress tensor of the bulk tensor is suppressed with respect to the shear one. 6 6 Here we are implicitly assuming that the magnitude of the shear tensor is not parametrically suppressed with respect to the bulk one. Should the flow of the system be prepared such that the shear tensor identically vanishes, then transport would be dominated by the bulk tensor. It is interesting to note that the ratio of viscosities at low temperatures violates Buchel's bound as illustrated in Fig. 7. Violations of this bound have been previously encountered in other models such as [45,48,50].

Quasi-normal modes and relaxation
We now turn to the description of the off-equilibrium dynamics of our holographic model. We study the reaction of the system to small perturbations which drive it away from local equilibrium. On the gravity side this problem translates into the study of the relaxation of the black brane solutions constructed above when the different background fields are perturbed. As is well known, this relaxation process is controlled by an infinite set of discrete, damped modes known as QNMs. In this section we will determine the QNM frequencies of the system as a function of the temperature.
Since in our holographic model the scalar field backreacts on the geometry, metric fluctuations couple to fluctuations of the scalar field and they must all be considered simultaneously. Denoting by G (T) the black brane metric in Eddington-Finkelstein coordinates (3.1), we will study fluctuations of the form The dynamics of hMN and ϕ is governed by the linearised Einstein and scalar field equations on the background spacetime G (T) M N . We will use the value of the unperturbed scalar field φ as a coordinate in the holographic direction.
As is well known (see e.g. [51]) not all fluctuations are physical, since reparametrisation invariance leads to a gauge symmetry in the linearised equations of motion. In the presence of a scalar field, the linearised equations of motion are invariant under the transformation In this paper we will study the relaxation of homogeneous disturbances of the plasma. In other words, we will allow for time but not for space dependence of the perturbations. We will consider both isotropic and anisotropic perturbations and we will denote by z the direction of anisotropy. On the gravity side, these perturbations will depend on time and on the holographic radial coordinate. Under these conditions, there are only two independent sets of gauge invariant excitations of the plasma, which may be parametrized by the following combination of fields 8 where haa = (hxx + hyy)/2. The first fluctuation, Zaniso, controls anisotropic perturbations that leave unaffected the expectation value of the scalar operator, the average pressure and the trace of the stress tensor. The non-conformal mode Z bulk controls fluctuations that change the three pressures in an isotropic way and at the same time modify the expectation value of the scalar operator and the trace of the stress tensor. At non-zero spatial momentum these excitations would be coupled to one another and they would include the hydrodynamic modes. Our restriction to the space-independent sector implies that the energy density of the plasma is unchanged by the fluctuations (4.3)-(4.4), since in a homogeneous plasma conservation of the stress tensor reduces to ∂t = 0.
Manipulating the linearised Einstein and Klein-Gordon equations and after a Fourier transform in time, the dynamics of the Zaniso and Z bulk modes are given by the equations where L bulk , L bulk , Raniso, Raniso are linear operators in the holographic direction given by with A, B and h the numerically computed functions which determine the background, given by (3.6)-(3.8). The equation for the anisotropic fluctuation, Zaniso, is that of a massless probe scalar field, while the equation for the bulk fluctuations Z bulk includes an explicit dependence on the potential. The discrete set of normalizable, in-falling solutions of this system of equations are the QNMs. The fact that the equations are linear in the frequency is a consequence of the Eddington-Finkelstein form of the thermal metric eq. (3.1).
Following [33] we use this to determine the QNMs and their associated frequencies by spectral methods, which allow us to reduce the problem of finding the complex-valued spectrum of excitations to an eigenvalue problem. This method is particularly suited for background metrics which are only known numerically. We have also double-checked the results for some representative frequencies with a shooting method. The QNM frequencies depend on the temperature of the plasma. As the temperature changes, each of these complex frequencies follows some trajectory in the complex plane. In Fig. 8 we show these trajectories for the four lowest QNMs of the anisotropic perturbations Zaniso for different values of φM . Each of the points on a given trajectory corresponds to a different value of the temperature. Note that in all panels these trajectories begin and end at the same value, indicated by the crosses, "+". The reason for this is that the Zaniso fluctuations correspond in the gauge theory to fluctuations exclusively of the stress tensor (i.e. with no contribution of the scalar operator). Since the stress tensor is conserved, its dimension is exactly 4 8 Anisotropic fluctuations induced by Z0 = hxy are also possible and independent of the two modes listed in (4.3)-(4.4). However, we will not consider these fluctuations here because at zero spatial momentum the dynamics of Z0 is identical to that of Zaniso.  both at the UV and at the IR fixed points regardless of the value of φM . In a CFT, this information of an operator alone would determine the spectrum of the dual QNMs. Since our theory approaches a CFT in the UV and in the IR, the QNMs associated to the pure-stress-tensor fluctuations Zaniso approach the same limiting conformal values at high and low temperatures. In contrast, at intermediate temperatures all the QNM frequencies possess a smaller imaginary part than their conformal counterparts. The magnitude of this deviation depends on the non-conformality of the theory. For φM = 1, when the non-conformal parameter I, eq. (3.18), is small at all temperatures, the complex plane trajectories of all modes remain close to the conformal value. As φM increases the excursion of all modes in the complex plane deviates more from the conformal values. Note, however, that these paths seem to saturate at very high value of the parameter φM . In particular, even though the change in the number of IR degrees of freedom differs by more than 3 orders of magnitude, the excursion in the complex plane of the simulations with φM = 10 and φM = 100 are very similar. This is in accordance with the small change in non-conformality observed in Fig. 4.
In Fig. 9 we show the complex-plane trajectories of the four lowest QNMs of the bulk mode Z bulk as a function of temperature for different values of φM . In all panels, the blue "+" crosses show the QNMs of a probe scalar field in an AdS black brane background dual to a CFT scalar operator of dimension 3 [53]. Similarly, the red " " squares show the QNMs of a probe scalar field in an AdS black brane background dual to an operator of dimension ∆IR given by eq. (2.9). Since this dimension depends on φM , the position of the red " " squares changes from panel to panel. Based on our discussion of the Zaniso QNMs above, one may expect that in the case of Z bulk the trajectories begin at the blue crosses at high temperature and end at the red crosses at low temperature. However, as we can see from Fig. 9, the trajectories in this case possess a more interesting structure.
In the upper left panel of Fig. 9 we show the trajectories for the φM = 1 potential. For this value,   the effective IR mass of the scalar is such that the first two QNMs of the ultraviolet probe scalar are closer to the real axis than the first QNM of the infrared probe scalar. This ordering determines the trajectories of the QNMs as a function of temperature. As shown in the plot, starting from the IR, the lowest QNM flows towards the closest UV mode in the complex plane, which in this case is the second UV mode. All IR modes follow similar trajectories in such a way that (at least as far as our numerics can resolve) the n-th IR mode flows to the (n + 1)-th ultraviolet mode. 9 As a consequence, there are no available IR modes to which the lowest UV mode can flow into. Therefore, this mode decouples at low temperature, flowing deep into the complex imaginary plane. For the other values of φM displayed in Fig. 9, the positions of the IR and the UV modes alternate in the complex plane, but this does not mean that the flow induced by the temperature is a direct map between these two sets of modes. Even though for the remaining three panels the lowest QNM flows between the lowest modes of the IR and UV theories, in all panels there is always a mode that decouples from the spectrum, although that mode is different for each of the displayed values of φM . The origin of this decoupling is that, after a certain mode, the n-th IR mode flows to the (n + 1)-th UV mode, interrupting the trajectory of the n-th UV mode. When this happens, we observe a phenomenon similar to level anticrossing in quantum mechanics. We have checked that for φM = 1000 (not shown) the complex-plane trajectories are almost identical to the φM = 100 trajectories displayed in the bottom-right panel of Fig. 9. This suggests that the observed structure saturates at large φM and is captured by the φM = 100 plot.
In Fig. 10 we show the temperature dependence of the imaginary (left) and real (right) parts of the first four quasi-normal frequencies for different values of φM . Each plot shows the QNM of the anisotropic (blue) and bulk (red) channels. As already discussed, both of these two sets of modes flow from their values in the UV fixed point to their values in the IR fixed point. As shown in the plots, the effective conformal behaviour of the QNMs at high temperature stops when the temperature becomes of order the source Λ. At higher temperatures, the temperature dependence of both the real and imaginary part of the modes is non-trivial, and it reflects the intricate trajectories in the complex plane displayed in Fig. 9 and Fig. 8. These plots also show explicitly how the disappearance of one bulk QNM occurs at low temperature. The fact that this disappearance seems to be linear in all plots in Fig. 10 clarifies the temperature dependence of this mode. The observed constant slope implies that this quasi-normal frequency becomes temperatureindependent at low temperature (we have explicitly checked this) and therefore it decouples from the IR theory.
As a final remark, we note that the numerical results displayed in Fig. 10 allow us to compare the magnitude of the different modes at the same temperature. As we will discuss in more detail in the next section, the imaginary part of the quasi-normal frequencies is related to the relaxation back to equilibrium of small plasma perturbations. It is interesting to note that the ordering of the imaginary parts of the anisotropic and bulk modes changes with temperature: while at high temperatures the imaginary part of the lowest bulk mode is smaller than that of the anisotropic mode, at low temperatures this order is reversed. This crossing of the imaginary parts of the lowest modes is present for all values of φM . Nevertheless, at φM = 1 this effect is much more prominent, since for this φM the disappearing QNM is the lowest bulk mode at high temperature. In the next section we will discuss the consequences of this behaviour.

Discussion
The behaviour of the QNM with smallest imaginary part, dubbed the lowest QNM, is particularly relevant for understanding the off-equilibrium dynamics of theories with a gravity dual. At non-zero spatial momentum, the lowest QNM of metric perturbations is dual to hydrodynamic excitations of the dual theory. However, in the zero-spatial momentum limit we have considered, the residues of the hydrodynamic poles vanish. 10 In this limit the relaxation back to equilibrium is controlled by the QNM frequencies, with the longest-lived excitation corresponding to the lowest QNM. We will refer to the inverses of the imaginary parts of the frequencies of the lowest QNMs in the different channels as relaxation times. In the non-conformal theory that we have studied, these important time scales have a very interesting behaviour.
As shown in Fig. 10, the relaxation time associated to the anisotropic and bulk channels have a nontrivial temperature dependence, as a consequence of non-conformality. To best understand the origin of this temperature dependence, following [33,34] in Fig. 11 we show the imaginary part of the lowest QNM as a function of δ = 1/3 − c 2 s for the anisotropic channel (left) and the bulk channel (right) for different values of φM . For the anisotropic channel, most of this dependence may be understood as a consequence of the change of the speed of sound, similarly to the holographic constructions analysed in [33,34]. Although the inverse relaxation time is not just a common function of cs for all models, up to small corrections a simple linear dependence of the imaginary part of the lowest anisotropic mode on δ = 1/3 − c 2 s provides a good estimate for the relaxation time in this channel. This simple approximate scaling does not work in the bulk channel, as shown in the right panel of Fig. 11. Unlike the anisotropic channel, the relaxation time is influenced significantly by the change in the scaling dimension of the scalar operator in the high and low temperature phases, which enters only indirectly into thermodynamic properties such as cs. Therefore, the relaxation of strongly coupled gauge theories is, in general, not just controlled by thermodynamic properties, but additional microscopic dynamics of the theory may also be important to understand this complicated process.
The different behaviour of these time scales reflects the fact that the way in which the system relaxes depends on the way it is excited. To focus the discussion, we will restrict ourselves to generic excitations of the stress tensor of the system. Since in a CFT the trace of the stress tensor vanishes by symmetry, in a CFT this trace cannot be affected by fluctuations of the bulk mode. Since in addition the bulk mode is isotropic, it follows that the stress tensor itself in a CFT cannot be affected by the bulk mode. Because of this decoupling, the relaxation of small excitations of the stress tensor is controlled solely by the lowest mode of the anisotropic channel, given by the δ = 0 intercept of Fig. 11 (left). In a non-conformal theory, however, this decoupling does not occur. Because of non-conformality, the fluctuations of the stress tensor and of the operator O mix. As an example, note that small isotropic variations of the pressure of the system, which at finite momentum are part of the sound channel, excite the bulk mode, as it can be easily inferred from the Ward identity eq. (3.17). More generally, the variation of the stress tensor associated to the two fluctuating channels (4.3)-(4.4) is given by where ∆pz and ∆p ⊥ are the diagonal components of the stress tensor along the direction of the anisotropic perturbation and perpendicular to it, and Z bulk and Z (4) aniso are the normalisable modes of the perturbations. These expressions show explicitly how both channels affect the dynamics of the pressure, while only the bulk channel affects the expectation value of the scalar operator. As a consequence, the relaxation of the stress tensor of the system will be dominated by the mode with the smallest imaginary part of the two sets of towers displayed in Fig. 10. As it can be seen in this plot, for all values of φM , relaxation is dominated by different modes at high and low temperatures. The competition between these two channels implies that the relaxation dynamics in our family of holographic models follows different paths at high and low temperatures.
The contributions of the anisotropic and the bulk modes to the stress tensor codify two different physical processes. As explained above, the anisotropic mode controls anisotropic perturbations of the pressure that leave unaffected the energy density, the expectation value of the scalar operator, the average pressure and the trace of the stress tensor. The bulk mode controls fluctuations that change the three pressures in an isotropic way and at the same time modify the expectation value of the scalar operator and the relation between the energy density and pressure given by the equation of state. The relaxation of a generic small stress tensor disturbance therefore requires two distinct process: the "isotropisation" of the system, which amounts to equating the diagonal spatial components of the stress tensor (pressures); and the "EoSization" of the system, with which we only refer to the process by which the trace of the stress tensor attains its equilibrium value. We have carefully defined these two terms to avoid any possible confusion with "thermalization", namely the process by which a system reaches perfect thermal equilibrium.
Consider first the case in which the bulk mode dominates the relaxation process, meaning that its associated lowest QNM decays faster than that associated to the anisotropic mode. In this case the system first relaxes the trace of the stress tensor, such that the pressures of the system no longer fluctuate independently, and only later equates the value of all the pressures to one another. In other words, the system first EoSizes and subsequently isotropizes. This is the behaviour of the holographic models at small values of φM , such as φM = 1, 3, at low temperatures. Since in CFTs the trace of the stress tensor is fixed, this relaxation path is very similar to that in CFTs.
In contrast, consider now the opposite case in which relaxation is dominated by the anisotropic mode, meaning that its associated lowest QNM decays faster than that associated to the bulk mode. In this case the pressure of the system is first isotropized to a value that is not related to the energy density through the equation of state, and only later the subsequent dynamics of this isotropic stress tensor relaxes this value of the pressure to that dictated by the equation of state. At high temperatures, this is the path to equilibration followed by our models, which differs qualitatively from the conformal case. 11 Finally, when the two relaxation times are comparable, as it is the case in the low temperature regime for large values of φM = 10, 100, both of these processes occur simultaneously.
Our calculations are done at zero spatial momentum. At non-zero k the analysis is more complicated because the anisotropic mode splits into the shear, the tensor and the sound modes, and the latter mixes with the bulk mode. Nevertheless, in the coupled bulk-sound system it is still possible to distinguish between those excitations that change the trace of the stress tensor and those that do not. These coupled dynamics 11 Note that the right-hand side of (2.17) may suggest that Λ O must be large in order to cause a significant violation of the equation of state, thus in possible conflict with the linear approximation. It would be interesting to explore this in a non-linear calculation.
will of course modify the EoSization and the isotropization times that we have computed. However, by continuity this modification must be small for small k. Since the QNM frequencies are parametrically of order T , we therefore expect that their ordering will remain the same provided k T . Although the analysis of QNMs can only provide definite answers for the fate of small perturbations offequilibrium, the rich structure exhibited in this relaxation process has implications for the dynamics of initial configurations that are far off-equilibrium. As we mentioned above, the numerical analyses of collisions in N = 4 SYM yield hydrodynamisation times that are comparable to the relaxation times obtained via a QNM analysis. While the microscopic explanation of this observation is not understood, this experience has led the authors of [32] to suggest that the hydrodynamisation of non-conformal theories is basically controlled by the temperature of hydrodynamisation, with small (non-parametric) differences with respect to the conformal case. Following this reasoning, we may estimate how much longer the hydrodynamisation can be in the family of theories that we have studied. Given the mixing of the bulk and anisotropic modes, this longest relaxation is given by the absolute minimum of the (negative) imaginary part of the QNM sets which, as shown in Fig. 10, is always controlled by the bulk mode. Comparing with the relaxation of conformal theories τ conf = 0.73/2πT , this maximal relaxation is τmax/τ conf = 2.1, 2.5, 3.0, 3.15 for φM = 1, 3, 10, 100. These maxima occur at T /Λ = 0.33, 0.19, 0.16, 0.16 for each model respectively. It would be interesting to test explicitly whether the connection with the linearised analysis persists in full numerical simulations of shock collisions in our non-conformal backgrounds. In particular, this would allow the study of the impact of the different relaxation channels on the on-set of hydrodynamic behaviour.