Critical behavior of non-hydrodynamic quasinormal modes in a strongly coupled plasma

We study the behavior of quasinormal modes in a top-down holographic dual corresponding to a strongly coupled $\mathcal{N} = 4$ super Yang-Mills plasma charged under a $U(1)$ subgroup of the global $SU(4)$ R-symmetry. In particular, we analyze the spectra of quasinormal modes in the external scalar and vector diffusion channels near the critical point and obtain the behavior of the characteristic equilibration times of the plasma as the system evolves towards the critical point of its phase diagram. Except close to the critical point, we observe that by increasing the chemical potential one generally increases the damping rate of the quasinormal modes, which leads to a reduction of the characteristic equilibration times in the dual strongly coupled plasma. However, as one approaches the critical point the typical equilibration time (as estimated from the lowest non-hydrodynamic quasinormal mode frequency) increases, although remaining finite, while its derivative with respect to the chemical potential diverges with exponent -1/2. We also find a purely imaginary non-hydrodynamical mode in the vector diffusion channel at nonzero chemical potential which dictates the equilibration time in this channel near the critical point.


Introduction
Quasinormal modes (QNM's) are exponentially damped collective excitations [1,2] that define the characteristic behavior of fluctuations of black holes and black branes (for reviews, see [3][4][5][6]). The spectra of QNM's collectively describe the linear part of the decaying fluctuations of a disturbed black hole, a phenomenon known as "quasinormal ringing", which is analogous to the decaying sound emitted by a brass bell when struck by a mallet [7]. For this reason, QNM's are of great interest to astrophysical and cosmological observations since they describe the ringdown of possible black hole remnants of binary stars and black hole mergers, which were pivotal to the direct detection of gravitational waves earlier this year [8,9].
In the context of the holographic gauge/gravity duality [10][11][12][13], the QNM's of asymptotically Anti-de Sitter (AdS) spacetimes carry wealthy information about the near equilibrium behavior of the dual strongly interacting quantum field theory (QFT). In fact, the QNM's associated with the fluctuations of a given bulk field are related to the poles of the retarded Green's function of the dual operator in the QFT [14,15]. These poles describe hydrodynamic and non-hydrodynamic dispersion relations with which one can not only compute hydrodynamic transport coefficients but also derive upper bounds for characteristic equilibration times of the dual QFT plasma [16]. Additionally, non-hydrodynamic modes play an important role in determining the applicability of the hydrodynamic gradient series, as demonstrated by studies in holography [17,18] and also in kinetic theory [19][20][21].
Previous works [22][23][24] have dealt with QNM's in bottom-up Einstein-dilaton constructions [25][26][27][28] exhibiting different kinds of phase transitions at zero chemical potential. Additionally, Ref. [29] investigated the QNM's associated with scalar operators in a topdown N = 2 * non-conformal plasma also at zero chemical potential. On the other hand, in [30] some of us investigated how the QNM's of an external scalar perturbation and, in particular, the equilibration time associated with the imaginary part of the lowest nonhydrodynamical quasinormal frequency, depends on the temperature and baryon chemical potential in a bottom-up, QCD-like Einstein-Maxwell-dilaton model at finite baryon density [31]. In general, through the holographic correspondence, any question regarding the thermalization process in a given strongly coupled gauge theory necessarily involves a study of the QNM's of its gravity dual. These modes describe different timescales in the gauge theory and, close to a critical point, one may expect that the QNM's of the corresponding gravity dual display critical behavior.
Near a critical point, thermodynamical quantities typically display fast variations which enable the definition of critical exponents. Static properties such as single time correlation functions and linear response coefficients to time-independent perturbations display critical behavior which are determined by the underlying equilibrium distribution. However, anomalous behavior is also observed in many dynamical quantities such as the transport coefficients, which depend on the properties of multi-time correlations functions and are not determined by the information contained in the equilibrium distribution. In fact, while static thermodynamical properties of several different physical systems may be grouped into a few different (static) universality classes, dynamical properties associated with slowly varying hydrodynamical fluctuations of a system near criticality do not fit into this static classification scheme, as discussed in detail in [32] nearly 40 years ago. As a matter of fact, the dynamic universality classes reviewed in [32] require the study of hydrodynamic modes, i.e., collective excitations whose frequency vanishes in the case of homogeneous disturbances. While these modes dominate the long time behavior of the system (since they are associated with conserved currents) and can be used to study how transport coefficients (such as the shear viscosity) behave near a critical point, it is conceivable that there is more information about dynamical critical phenomena in multi-time correlation functions that cannot be obtained from their zero frequency limit.
In this paper we initiate the investigation of the critical behavior displayed by nonhydrodynamic modes in strongly coupled gauge theories with gravity duals, i.e., QNM's corresponding to collective excitations in the dual plasma whose frequency does not vanish in the zero wavenumber limit. This novel type of critical phenomena determines the behavior of different characteristic equilibration times of the system at zero wavenumber and, since these QNM's are not directly associated with conserved currents, their behavior at criticality does not follow from the analysis made in Ref. [32]. In the present work, as the first exploration in this new arena, we compute QNM's for an external scalar perturbation and also for the diffusion channel associated with a vector perturbation in the so-called 1-R charge black hole (1RCBH) model [33][34][35][36][37][38]. This is an analytical top-down construction obtained from (4 + 1)-dimensional maximally supersymmetric gauged supergravity, which is holographically dual to a strongly coupled N = 4 super Yang-Mills (SYM) plasma in flat (3 + 1) dimensions with a finite chemical potential under a U (1) subgroup of the global SU (4) symmetry of R-charges. This theory is conformal and its phase diagram is a function of a single dimensionless ratio µ/T , where µ and T are the U (1) R-charge chemical potential and temperature of the black brane background, respectively. The model has a very simple phase diagram with a critical point at µ/T = π/ √ 2 and its static critical exponents were computed in [39] and [40]. Also, the fact that the R-charge conductivity remains finite at the critical point [39] shows that this model belongs to the type B dynamical universality class [32] and the anomalous static critical exponent was found to vanish in [40]. Thus, this model is of mean-field type [40], which was later argued [41] to be a general consequence of the underlying large N c approximation. This simple model provides a useful arena for investigating dynamical phenomena in a strongly coupled plasma at finite temperature and density, even though it does not possess the full set of physical properties (such as chiral symmetry) displayed by the real world quark-gluon plasma (QGP) [42]. In fact, such a model may be useful for discovering new dynamical phenomena associated with critical endpoints in strongly coupled non-Abelian plasmas which could be further investigated in more realistic models of the QGP such as [30,31] , with a view towards applications to the ongoing beam energy scan program at RHIC. Other studies of critical phenomena in holography include Refs. [43][44][45][46][47].
As we are going to show in the next sections, the real and imaginary parts of nonhydrodynamical modes in the external scalar and vector diffusion channels display an infinite slope at the critical point of the phase diagram of the 1RCBH model. This holds true also for higher order QNM's, showing that high frequency modes are also sensitive to the presence of the critical point. In particular, from the imaginary part of the QNM's it is possible to extract the behavior of different characteristic equilibration times in the finite density plasma at criticality (at zero wavenumber) and define a dynamical critical exponent associated with their derivatives with respect to the dimensionless ratio µ/T . We find the same critical exponent 1/2 for all the equilibration times investigated in the different channels. Except close to the critical point, we observe that by increasing the chemical potential one generally increases the damping of the quasinormal black brane oscillations which, consequently, leads to a reduction of the characteristic equilibration times of the dual plasma. However, as one approaches the critical point these equilibration times are enhanced (though they remain finite) and they acquire an infinite slope. We also find a purely imaginary, non-hydrodynamical mode in the vector diffusion channel at nonzero chemical potential and zero wavenumber which dictates the critical behavior of the equilibration time in this channel (this mode was also found in Ref. [48] in the context of a (4 + 1)-dimensional Einstein-Maxwell model).
This paper is organized as follows: in Section 2 we briefly review the 1RCBH model and its thermodynamics and phase diagram. In Section 3 we compute the QNM's of an external scalar perturbation. In Section 4 we compute the QNM's for the vector diffusion channel in the limit of long wavelengths. Finally, in Section 5 we present some final remarks and an outlook. This work is complemented by Appendix A, where we compute the spectral function in the external scalar channel, and by Appendix B, where we present a brief discussion regarding the numerical procedures we used for the computation of QNM's and their critical exponents.
Throughout the paper, we work with natural units = c = k B = 1 and use a mostly plus metric signature.

1-R charge black hole model
For the sake of completeness, in this section we review the thermodynamics of the 1RCBH model [33][34][35][36][37][38]. We closely follow the discussion made in Section 4 of Ref. [49] and supplement it with additional plots to better illustrate the behavior of the thermodynamic properties of the model near the critical point.

Background
The 1RCBH model is described by an Einstein-Maxwell-dilaton (EMD) action, with the dilaton potential and the coupling between the dilaton and Maxwell fields given by, where κ 2 5 is the five dimensional Newton's constant and L is the asymptotic AdS 5 radius. With these profiles for V (φ) and f (φ), one obtains a consistent truncation of maximally supersymmetric gauged supergravity in five dimensions, which is itself a consistent truncation of type IIB superstring theory on AdS 5 × S 5 . This model is a bona fide top-down string theory construction, which is dual to a SYM plasma with a finite chemical potential under a U (1) subgroup of the global SU (4) symmetry of R-charges. For simplicity, in the following we set L = 1.
The 1RCBH solution is defined by , where r is the holographic radial coordinate and the boundary of the asymptotically AdS 5 geometry is located at r → ∞. The radial position of the black brane horizon may be written in terms of the charge Q and mass M of the black brane as follows, The 1RCBH background is, thus, characterized by two nonnegative parameters, (Q, M ) or, alternatively, (Q, r H ). The Hawking temperature of the black brane horizon is given by where denotes a derivative with respect to r. On the other hand, the U (1) R-charge chemical potential reads (2.11)

Phase diagram
The class of solutions corresponding to the 1RCBH model may be parametrized by different values of the dimensionless ratio Q/r H . By dividing Eqs. (2.10) and (2.11), and then solving for Q/r H , one obtains Since Q/r H is nonnegative, (2.12) implies that µ/T ∈ 0, π/ √ 2 . It also follows from (2.12) that for every value of µ/T ∈ 0, π/ √ 2 , there are two different corresponding values of Q/r H , which parametrize two different branches of solutions. As we are going to show in the next subsection, by analyzing the thermodynamics of the 1RCBH solutions one concludes that the point of the phase diagram where these two branches merge, µ/T = π/ √ 2 or, correspondingly, Q/r H = √ 2, is a critical point of a second order phase transition. In order to simplify the equations in this paper we define below a useful variable that smoothly connects the two branches of solutions, where y = 0 parametrizes the critical background geometry with µ/T = π/ √ 2, while y = 1 parametrizes the AdS 5 -Schwarzschild background with zero charge, corresponding to Q = 0 and r H = 0, which implies µ/T = 0. For y = −1 we also have µ/T = 0, but this time r H = 0 and Q = 0, which corresponds to a supersymmetric BPS solution dubbed "superstar" [50] instead of a black hole.
As we are going to see in the next subsection, the thermodynamically stable branch corresponds to the lower sign in Eq. (2.12) with Q/r H ∈ 0, √ 2 or y ∈ 0, 1 , while the thermodynamically unstable branch corresponds to the upper sign in Eq. (2.12) with Q/r H ∈ √ 2, ∞ or y ∈ −1, 0 , with both branches of solutions being smoothly connected at the critical point, Q/r H = √ 2 or y = 0. This is illustrated in Fig. 1.

Thermodynamics: equation of state and susceptibilities
For a SYM plasma, it is known [51] that (2.14) By substituting the expression above in Bekenstein-Hawking's relation [52,53], one can write down the entropy density as follows The R-charge density, ρ = lim r→∞ δS/δΦ , may be written as where the lower/upper signs denote the stable/unstable branches, as in Eq. (2.12) (we note that in this model any dimensionless ratio is always written in terms of µ/T , as expected). From the Gibbs-Duhem relation, dp = sdT + ρdµ, one may compute the pressure, Using Eqs. (2.15), (2.16), and (2.17) one can readily evaluate the internal energy density, ε = T s − p + µρ, obtaining ε = 3p, as expected for a conformal QFT in four dimensions. The heat capacity at fixed chemical potential is given by C µ = (∂s/∂T ) µ , while the nth-order R-charge susceptibility is given by χ n = (∂ n p/∂µ n ) T (note that χ 1 = ρ). At the critical point the heat capacity C µ and the higher order susceptibilities χ n≥2 diverge. In Fig. 2 we show the equation of state and the heat capacity while in Fig. 3 we display the susceptibilities for the stable and unstable branches of the 1RCBH model.
Thermodynamical stability is ensured if the Jacobian J = ∂(s, ρ)/∂(T, µ) is positive. In the 1RCBH model one obtains In the equation above, the quartic term is always positive while the expression in the last parenthesis becomes negative for y ∈ − 1, 0 , justifying the aforementioned classification of stable and unstable branches. 1 Finally, by analyzing the behavior of the R-charge density in (2.16) near the critical point one can see that the static critical exponent δ = 2, as discussed in [39,40,49].

Equation of motion
In the previous section we reviewed the thermodynamical equilibrium properties of the 1RCBH plasma. We now analyze near-equilibrium properties of the system encoded in its quasinormal modes. In this section we calculate the QNM's for an external scalar perturbation ϕ on top of the 1RCBH backgrounds, which is described by the bulk action, The equation of motion following from this action is just the massless Klein-Gordon equation on top of the solution given by Eq. (2.3). We take a plane-wave Ansatz for the Fourier modes of the perturbation, ϕ = e −iωt+i k· xφ (ω, k, r), which for brevity we write simply asφ(ω, k, r) ≡φ(r). The resulting equation of motion then only depends on the frequency ω, the magnitude of the spatial 3-momentum k ≡ | k|, and the background control parameter y.  In what follows we employ the in-falling Eddington-Finkelstein (EF) "time" coordinate defined by, in terms of which the metric (2.3) becomes One of the main advantages of the EF coordinates is that the in-falling wave condition at the horizon, which is associated with the retarded Green's function, becomes automatically satisfied by just requiring regularity of the solutions there. In these coordinates, the equation of motion forφ becomes, We map the radial coordinate r, defined on the interval r H ≤ r < ∞, to a new dimensionless radial coordinate u = r H /r, defined on the interval 0 ≤ u ≤ 1, which is more suitable to be used in the pseudospectral method [54] (to be briefly reviewed in Appendix B). In these new coordinates, the equation of motion for the external scalar perturbation becomes, where the primes now denote derivatives with respect to the new radial coordinate u. From the discussion above, and from the definition of the background control parameter y in Eq. (2.13), one concludes that the dimensionless quasinormal eigenfrequencies, ω/T , will depend only on the dimensionless ratios µ/T and k/T . To completely specify the eigenvalue problem to be solved in order to find the QNM spectra associated to this external scalar perturbation, we still need to impose a Dirichlet boundary condition. From the fact thatφ is a scalar field defined on an asymptotically AdS 5 background, it follows that asymptotically close to the boundary it may be written as ϕ(u) = G(u)+u 4 F (u), with the leading, non-normalizable mode G(u → 0) = J(ω, k) being the source for the QFT operatorÔ dual to the (external) scalar field ϕ, and the subleading, normalizable mode F (u → 0) = Ô (ω, k) being its expectation value. According to the real time holographic dictionary [55], the retarded propagator of the QFT operatorÔ is given the ratio between the normalizable and non-normalizable modes, G R OÔ (ω, k) = − Ô (ω, k) /J(ω, k), therefore, if we impose as a Dirichlet boundary condition the selection of the normalizable mode by setting G(0) = 0 with F (0) = 0, we are left with an eigenvalue problem whose eigenfrequencies correspond to dispersion relations ω/T = ω(k/T ; µ/T )/T describing the poles of G R OÔ , which are the QNM's we are looking for. Then, we set φ(u) = u 4 F (u), with F (0) = 0, from which it follows that, We now have a Generalized Eigenvalue Problem (GEP) for the eigenfunction F (u) and the quasinormal eigenfrequency ω/T , which may be solved as functions of k/T and µ/T . In this work we use the pseudospectral method. Note that Eq. (3.6) also reveals one of the greatest virtues of the EF coordinates, namely, the fact that it reduces the QNM eigenvalue problem from a Quadratic Eigenvalue Problem [56] in the standard set of spacetime coordinates to a GEP, which requires far less computational cost when numerically evaluating the QNM spectra.

QNM spectra and equilibration time
In Fig. 4 we show the evolution of the external scalar QNM spectra for the first 26 poles as we evolve k/T from 0 to 100, both for the AdS 5 -Schwarzschild background at µ/T = 0 and the critical geometry at µ/T = π/ √ 2. We observe the usual non-hydrodynamical QNM structure for the external scalar channel with an infinite series of QNM pairs with Im ω < 0 and Re ω = 0 symmetrically distributed with respect to the imaginary axis [15]. When k/T = 0, by increasing the µ/T ratio one increases the magnitude of the imaginary part of the QNM's, which becomes more appreciable for higher order, faster varying modes. On the other hand, an increase in k/T enhances (suppresses) the magnitude of the real (imaginary) part of the poles. We see that by increasing the chemical potential one generally increases the damping of the quasinormal black hole oscillations, which qualitatively agrees with the result found previously in [30] for a non-conformal, QCD-like bottom-up EMD model describing the physics of the QGP at finite baryon density.
Also, we note that the non-hydrodynamic modes in Fig. 4 remain finite when evaluated at the critical point, even when k = 0. Thus, one can see that the timescales contained in the non-hydrodynamic modes are different than the usual relaxation time quantity τ rel ∼ ξ z , where ξ is the correlation length (which diverges at the critical point) and z is the dynamical critical exponent, which becomes infinitely large at criticality describing the well-known phenomenon of critical slowing down. Nevertheless, in this strongly coupled model the microscopic scales defined by the non-hydrodynamic QNM's still display some critical behavior, as we show below.
In Figs. 5 and 6 we display the imaginary and real parts of the first 4 QNM's as functions of µ/T , for both stable and unstable branches, at k/T = 0 and k/T = 1. We see that at the critical point all the QNM's develop an infinite slope. Moreover, we also note that the effects on the non-hydrodynamic modes due to finite momentum are small for k/T ∼ 1 (especially for the imaginary part), being more pronounced for the lowest QNM's, which seems to be a general holographic property of the dispersion relation of non-hydrodynamics QNM's known as "ultralocality" [57,58].
Following [16], one may define an upper bound for the equilibration time of the plasma, τ eq , using the inverse of minus the imaginary part of the lowest non-hydrodynamical QNM evaluated at zero momentum. This is shown in Fig. 7 (a), from which one can see that far from the critical point the equilibration time of the finite U (1) R-charge density SYM plasma decreases with increasing chemical potential, again in qualitative agreement with what was previously found in [30] in the context of a phenomenologically realistic holographic model for the QGP at finite baryon chemical potential. However, for larger chemical potentials, as one approaches the critical point of the model, this behavior is modified and the equilibration time starts to increase, acquiring an infinite slope at the critical point.
We may associate a new critical exponent with the derivative of the equilibration time, d(T τ eq )/d(µ/T ), since it diverges at the critical point, as shown in Fig. 7 (b). Close to the critical point, where θ is the dynamical critical exponent which we want to calculate. By following the procedure discussed in Section B.3, we performed numerical fits of Eq. (3.7) to the data shown in Fig. 7 (b) for different sets of intervals in µ/T increasingly closer to the critical point at µ/T = π/ √ 2. The final result gives a critical exponent compatible with θ = 1/2. More generally, one may consider different characteristic equilibration times of the medium associated with the different non-hydrodynamic QNM's, where equilibration times associated with higher order modes should be understood as estimates for how fast the system relaxes to equilibrium depending on how rapidly varying are the perturbations to which it is subjected. Operationally, this amounts for computing the inverse of minus the imaginary part of the different non-hydrodynamic QNM's. By doing so, we obtain the same dynamical critical exponent θ = 1/2 associated with all the different characteristic equilibration times of the plasma in the external scalar channel at zero wavenumber.

Equation of motion
In this section we compute the QNM's of the vector diffusion channel in the long wavelength limit. Differently from what was done in the last section where we considered an external scalar perturbation on top of the 1RCBH background, now we need to consider fluctuations of the Maxwell field A µ which is already nonzero in the background and, therefore, we also need to consider disturbances in the background metric g µν and dilaton field φ. At zero spatial momentum the different channels for these fluctuations, at the linearized level, are classified by different representations of the SO(3) rotation group where the first equation is the equation of motion for the background Maxwell field, Φ(r), while the second equation is the equation of motion for the Maxwell perturbation, a(r). We may decouple the perturbations a(r) and s(r) by using Einstein's equations for the metric field, Once again we apply the radial coordinate transformation r → r H /u, which yields, of this more involved case in the present work.
For the vector perturbation the normalizable mode at the boundary corresponds to set a(u) = u 2 F (u), with F (0) = 0, from which one finally obtains,

QNM spectra and equilibration time
With the QNM eigenvalue problem completely specified as discussed above, we can now apply the pseudospectral method (see Appendix B) to numerically solve it. In Fig. 8 we display the QNM spectra for the first 30 symmetric poles in the vector diffusion channel in the limiting cases of µ/T = 0 (AdS 5 -Schwarzschild) and µ/T = π/ √ 2 (critical point). In Figs. 9 and 10 we show the imaginary and real parts of the first 4 complex QNM's as functions of µ/T , for both stable and unstable branches. We see that also in the vector diffusion channel both the real and imaginary parts of the QNM's develop an infinite slope at the critical point.
We remind the reader that at k = µ = 0 the QNM spectra in the vector diffusion channel may be analytically calculated [15], Our numerical calculations at µ/T = 0 agree with this analytical result. The standard hydrodynamical mode ω(k/T = 0)/T = 0 is depicted by the blue diamond in Fig. 8. Since this is a hydrodynamical pole, it does not evolve with the chemical potential if we keep k = 0. This mode determines the R-charge conductivity of the model and the zero frequency limit of this transport coefficient was found in [49] to remain finite at the critical point, as expected for a type B dynamic universality class. However, it was noticed in [49] that the derivative of this quantity near the critical point has infinite slope described by an exponent equal to 1/2, which matches the exponent found in the previous section in the study of non-hydrodynamic modes of different nature corresponding to external scalar perturbations.
On the other hand, the main effect of the chemical potential on the symmetric nonhydrodynamical modes is to increase the magnitude of both the imaginary and real parts of these poles. Therefore, also in the vector diffusion channel one sees that the inclusion of a chemical potential leads to additional damping for the quasinormal black hole oscillations. A novel feature we observe in Fig. 8 is the emergence of a new purely imaginary, non-hydrodynamical pole at finite chemical potential, which comes from ω/T → −i∞ at µ/T = 0 and lies at ω/T ≈ −7.315i at the critical point µ/T = π/ √ 2. For µ/T 2, this new purely imaginary pole becomes the lowest non-hydrodynamical mode, while for lower values of the chemical potential the lowest non-hydrodynamical mode is given by any of the first two symmetric poles with respect to the imaginary axis. Therefore, this new nonhydrodynamical imaginary mode plays a crucial role in the description of the equilibration time of the system in the vector diffusion channel when the chemical potential is large, specially at criticality, when it dominates the physics of the slowest varying perturbations. The appearance of such a purely imaginary mode is an interesting feature of this model that shows that the distinction between transient phenomena at weak and strong coupling, currently understood in terms of their different pattern of non-hydrodynamic modes at zero wavenumber [59] (see also [60]) corresponding to fluctuations around global equilibrium, can become more complicated near a critical point.
We define the upper bound for the equilibration time of the system in the vector diffusion channel as before by taking the inverse of minus the imaginary part of the lowest non-hydrodynamical QNM. The result is shown in Fig. 11 (a). The kink observed in the equilibration time at µ/T ≈ 2 is due to the shift from the regime dominated by the first symmetric poles to the regime dominated by the new purely imaginary mode. This also causes a discontinuity in the derivative of the equilibration time, as seen in Fig. 11 (b). As before, one can calculate the critical exponent associated with this derivative at the critical point and the result is once again compatible with θ = 1/2. This shows that in this model both the hydrodynamic and the non-hydrodynamic modes in this vector diffusion channel have the same critical exponents.

Final Remarks
In the present work we investigated the behavior of the QNM's for the external scalar and vector diffusion channels of the 1RCBH model, which is a top-down gauge/gravity construction dual to a SYM plasma at nonzero U (1) R-charge density. We analyzed the behavior of the QNM's through the phase diagram of the model, which displays a critical point at a second order phase transition. We found that, except close to the critical point, by increasing the chemical potential one generally increases the damping of the quasinormal black hole oscillations, which leads to a reduction of the characteristic equilibration times of the dual plasma. However, as one approaches the critical point these equilibration times are enhanced and they acquire an infinite slope at the criticality. We found that the derivatives of all the characteristic equilibration times of the medium, obtained from the non-hydrodynamic QNM's at zero wavenumber, share the same critical exponent θ = 1/2. Previously, the same value was also found for the critical exponent associated to  the derivative of the DC conductivity extracted from the R-charge diffusive hydrodynamic mode in this model [49]. We also found a purely imaginary, non-hydrodynamical mode in the vector diffusion channel at nonzero chemical potential which dictates the critical behavior of the equilibration time in this channel. The observation that the quasinormal black hole oscillations away from the critical region are additionally damped by a nonzero chemical potential, obtained here for a top-down conformal construction dual to a SYM plasma at finite R-charge density, is consistent with the behavior found previously in [30] for a rather different holographic construction, which involves a bottom-up model with black brane solutions that are engineered to describe the realistic non-conformal physics of the QGP both at zero baryon density [28] and also at finite density [31]. This may indicate that this additional damping in the quasinormal black hole oscillations due to a nonzero chemical potential, and the consequent attenuation of the equilibration time of the dual plasma away from criticality, may be a general holographic property of strongly correlated quantum fluids.
Regarding the purely imaginary, non-hydrodynamical mode found in the diffusion channel at finite density, the common link between the 1RCBH model studied here and the Einstein-Maxwell model investigated in Ref. [48] (where this mode was also found, although with no critical behavior, due to the lack of a phase transition in the model of Ref. [48]) is the presence of bulk electromagnetic fields. Somehow, the Maxwell field changes qualitatively the dynamical response of the system to perturbations. From the point of view of the gravitational dynamics in five dimensions, what we see is that this purely imaginary non-hydrodynamic mode appears for electromagnetically charged asymptotically AdS black holes and that, when the charge of the black hole is large enough, this new mode dominates the dynamics of relaxation towards equilibrium when the black hole is disturbed. From the point of view of the four dimensional dual quantum field theory at the boundary of the asymptotically AdS spacetime background, this new mode dominates the dynamics responsible for the characteristic equilibration time of the plasma at large enough densities in the diffusion channel. Since bulk electric fields map to boundary states at finite density, this may be a general feature of holographic models at finite density. Therefore, one could investigate if the behavior of non-hydrodynamic modes in other holographic models that display critical phenomena possess similar properties to those found in the present study, i.e., the corresponding equilibration times have infinite slope characterized by a single critical exponent θ. In particular, we intend to investigate in the near future these features in two bottom-up constructions of phenomenological relevance for the physics of the QGP: the EMD model at finite baryon chemical potential [30,31] and the anisotropic EMD model at finite magnetic field from Refs. [61,62]. Such studies may be relevant to understand how the presence of a critical endpoint in the QCD phase diagram may lead to new observables associated with, for instance, baryon transport in the baryon rich quark-gluon plasma produced in heavy ion collisions within the beam energy scan program at RHIC.

A Spectral function for the external scalar channel
In this Appendix we study the spectral function associated with the external scalar fluctuation in the bulk. The motivation for this pursuit comes from the fact that QNM's are associated with poles of the retarded Green's function, whose imaginary part defines the spectral function. Therefore, one should expect that the critical behavior found for the QNM's leave somehow an imprint in the spectral function. Here we investigate this issue by considering the spectral function of the external scalar field fluctuation. Also, we note that due to the universal character of the δg y x fluctuation of the metric [63], the same calculation also gives the shear viscosity spectral function.
For the sake of completeness, let us begin by briefly reviewing the holographic computation of the spectral function based on the real time holographic prescription [55] recast using the holographic membrane paradigm [64]. From linear response theory, the expectation value of the QFT operatorÔ dual to the scalar perturbation is associated with the retarded correlator according to (see [65]), where, as before, J(ω, k) denotes the leading mode of the scalar fluctuation at the boundary, which sources the QFT operatorÔ, while the expectation value is associated with the subleading mode, being given by, where Π is the radial canonical momentum conjugate to the scalar perturbationφ. Once more, we impose that the scalar perturbation satisfies an in-falling wave condition at the horizon, which gives the retarded propagator. Our goal here is to compute the spectral function defined by, To do so, we have to first solve the equation of motion that follows from the action (3.1) in the usual coordinates of (2.3), with an in-falling horizon condition at r = r H and lim r→∞φ (r, ω, k) = J(ω, k).
We introduce a bulk response function [64], 3 which allows us to reduce the linear second order differential equation (A.4) to a first order nonlinear Riccati equation, By requiring regularity at the horizon, one obtains the following horizon condition needed to solve the first order flow equation above, where we choose the positive sign, which corresponds to the in-falling wave at the horizon. From the membrane paradigm [64], assuming that the scalar disturbance corresponds to the δg y x fluctuation of the metric, one recognizes this result as the shear viscosity Then, one may write down the following dimensionless ratio, One can now numerically integrate Eq. (A.6) and then use Eq. (A.9) to obtain the normalized spectral function. The results are shown in Fig. 12 for µ/T = 0 (AdS 5 -Schwarzschild) and µ/T = π/ √ 2 (critical point), both evaluated at k = 0. We see that, naively, an increase in µ/T seems to have only a small effect on the spectral function, even as one approaches the critical point. However, this is due to the fact that in the ultraviolet limit, ω/T → ∞, the dimensionless ratio F s /ωη scales as (ω/T ) 3 , which overwhelms any poles or fluctuations in the plot for the spectral function.
Let us now consider a subtraction scheme which removes the scaling (ω/T ) 3 from the spectral function at some ultraviolet cutoff [66][67][68], where a is a constant defined by the asymptotic behavior of the spectral function, We remark that in order to reliably perform the subtraction above, the flow equation (A.6) must be solved with high numerical accuracy. This becomes more difficult for larger values  of the ultraviolet cutoff in ω/T , which also requires one to increase the ultraviolet cutoff used to numerically parametrize the boundary in the radial coordinate. For the numerical evaluation of the first oscillations of the subtracted spectral function, one can safely take as the ultraviolet cutoff in the dimensionless frequency a value around ω/T ∼ 40. In Fig. 13 we display the behavior of ∆F s /ωη as a function of µ/T and ω/T at k/T = 0, while in Fig. 14 we analyze the evolution of the height of its first peak as a function of µ/T , which acquires an infinite slope at the critical point. By computing its derivative and fitting the numerical result close to the critical region using the same functional dependence as in Eq. (3.7), one again obtains a critical exponent compatible with 1/2. Thus, we find that the critical behavior found in the QNM's can also be found, albeit in an indirect manner, in the spectral function.

B Numerical procedures
The main numerical algorithm we employ in this work to find the QNM's is the pseudospectral method [54]. In this Appendix, we briefly review the main steps required to tackle the problem. The main advantage of the pseudospectral method, when compared to other methods used in the literature to calculate QNM's, is the ease of numerical implementation and the accuracy -in general, it requires a modest number of collocation points and basis functions in order to compute several QNM's with high accuracy. The main disadvantage is the requirement of high numerical precision in the intermediate calculations, which still brings a drawback in terms of running time, since typical machine precision calculations result in spurious results for any but the lowest QNM's.

B.1 Brief overview of the pseudospectral method
The general spectral method, of which the pseudospectral (or collocant) method is a particular case, 4 aims to solve the following general non-homogeneous differential equation for a complex function f (u),L f (u) = h(u), whereL is a general differential operator and h(u) is the non-homogeneous term. One may consider, for instance, the basic interval u ∈ [0, 1], which was used to define the radial 4 The main difference between spectral and pseudospectral methods regards the determination of the coefficients ai to be specified in the sequel. As explained in [54], the nomenclature used in the literature is a bit messy: both spectral and pseudospectral methods are known as spectral methods in a broad sense, due to the fact that they use a complete set of orthogonal functions. In the restricted sense, spectral methods (also called non-collocant methods) determine the generalized Fourier coefficients by exploiting the orthogonality of the basis, projecting down the unknown function f (u). On the other hand, pseudospectral methods (also known as collocant methods) use a selected set of points on the function domain in order to build an interpolating polynomial. For the purposes of the present work, we will consider only pseudospectral (or collocant methods) in this restricted sense.
holographic coordinate u in the main text (with the boundary at u = 0 and the horizon at u = 1). In finite difference methods one discretizes the basic interval using a finite grid and introduce finite difference approximations for the derivatives. In both spectral and pseudospectral methods, one instead introduces a subset {φ i (u)} N i=0 of a complete set {φ i (u)} ∞ i=0 of orthogonal basis functions defined on the basic interval, approximating f (u) by its truncated base expansion f N (u), Now we consider the pseudospectral method in order to determine the coefficients a i that best provide an approximation f N (u) for f (u). First, one introduces the residual function R (u; {a i }) defined by, The strategy is to suitably choose a set of points {u j } N j=0 (the so-called collocation points) on the domain of the basis functions φ i , and then fix a i such that the residual function is zero at u j , that is, R(u j ; {a i }) = 0. This generates a system of N linear differential equations with N variables, {a i }, There are several possible choices for the basis functions φ i and collocation points u j , depending on the symmetries of the problem and boundary conditions. A generally proposed basis defined in the interval x ∈ [−1, 1] is given by the Chebyshev polynomials, T n (x). 5 Since u ∈ [0, 1], we can relate these variables by x = 2u−1, then φ i (u) = T i (2u−1). A useful accompanying set of collocant points is given by the so-called Gauss-Lobatto grid, which for u ∈ [0, 1] can be written as,  Table 1: d(T τ eq )/d(µ/T ) intervals used for the fit procedure.

B.2 Implementation details
We used a Chebyshev basis with N basis functions and collocant points, and then employed the Arnoldi method (via the built-in Eigensystem[ ] procedure in Wolfram's Mathematica [69]) to solve the resulting GEP for ω/T . 6 For higher order QNM's, one needs to use high numerical precision from the outset. Denoting the number of floating-point digits used in the calculations by M , the final error estimate of the numerical QNM's is mainly controlled by the number of basis functions N and the numerical precision parameter M . We have checked that using M = 60 and N = 80 yielded the 20 first QNM's with good accuracy -this was the setting used in most of the paper. For higher order QNM's, such as the ones shown in Figs. 4 and 8, we used M = N = 100. We verified the stability of the QNM spectra involving the desired modes by doubling the number of basis points N or the numerical precision parameter M , following the error control procedure discussed in [54].

B.3 Calculation of the critical exponents
Let us now discuss the numerical procedure we followed in order to determine the critical exponent θ of the divergent quantities near the critical point µ/T = π/ √ 2, by performing fits of the asymptotic form (3.7) as a function of µ/T . We used a first order central difference formula with step size h ≡ ∆(µ/T ) in order to compute the required numerical derivatives. By varying h from 10 −5 to 10 −11 , and taking into account that we have an accuracy of several digits in the computed observables, we estimated the error introduced by the numerical differentiation to be of the order of 10 −7 .
Taking as a specific example the calculation of the critical exponent for the equilibration time in the external scalar channel, in order to check that the asymptotic form (3.7) is valid close to the critical point µ/T = π/ √ 2, we split d(T τ eq )/d(µ/T ) into 14 subintervals and then performed a least squares fit of Eq. (3.7) to the resulting data within each subinterval, with each of them populated by 100 points evenly spaced in µ/T . In Table 1 we specify the subintervals used for the determination of the critical exponent of d(T τ eq )/d(µ/T ) in the external scalar channel. In Fig. 15 we display the convergence of the fitted critical exponent θ to the value 1/2. The estimate of the least squares standard error in the value of θ is of the order of 10 −7 for the last interval.