Non-hydrodynamic quasinormal modes and equilibration of a baryon dense holographic QGP with a critical point

We compute the homogeneous limit of non-hydrodynamic quasinormal modes (QNM's) of a phenomenologically realistic Einstein-Maxwell-Dilaton (EMD) holographic model for the Quark-Gluon Plasma (QGP) that is able to: i) {\it quantitatively} describe state-of-the-art lattice results for the QCD equation of state and higher order baryon susceptibilities with $2+1$ flavors and physical quark masses up to highest values of the baryon chemical potential currently reached in lattice simulations; ii) describe the nearly perfect fluidity of the strongly coupled QGP produced in ultrarelativistic heavy ion collisions; iii) give a very good description of the bulk viscosity extracted via some recent Bayesian analyzes of hydrodynamical descriptions of heavy ion experimental data. This EMD model has been recently used to predict the location of the QCD critical point in the QCD phase diagram, which was found to be within the reach of upcoming low energy heavy ion collisions. The lowest quasinormal modes of the $SO(3)$ rotationally invariant quintuplet, triplet, and singlet channels evaluated in the present work provide upper bounds for characteristic equilibration times describing how fast the dense medium returns to thermal equilibrium after being subjected to small disturbances. We find that the equilibration times in the different channels come closer to each other at high temperatures, although being well separated at the critical point. Moreover, in most cases, these equilibration times decrease with increasing baryon chemical potential while keeping temperature fixed.

Consequently, as discussed e.g. in Ref. [39], the longest-lived non-hydrodynamic QNM's 3 with lowest imaginary part (in absolute value) provide upper bounds for characteristic equilibration times of the late linear stage of the full nonlinear evolution of the medium. In fact, it has been explicitly shown e.g. in Ref. [40] that the lowest QNM's of a top-down conformal holographic plasma with a critical point quantitatively describe the late time behavior of the full nonlinear far-from-equilibrium evolution of the system. On the other hand, earlier nonlinear stages are beyond the scope of these characteristic equilibration times extracted from the QNM analysis and, thus, in general one cannot read e.g. the isotropization and thermalization times of the medium without performing numerical simulations and following the full nonlinear evolution of the system.
Here we initiate the study of how the presence of a critical point affects the non-equilibrium behavior of strongly coupled (non-conformal) baryon dense holographic models. In this work, we investigate the behavior of the lowest, homogeneous non-hydrodynamic QNM's of the bottom-up EMD model of Ref. [23] up to the critical region, from which we extract upper bounds for characteristic equilibration times of the hot and baryon dense QGP specifically in the regime of linearized perturbations. As it will be shown in the course of the next sections, these equilibration times of the plasma in different channels are generally reduced by increasing the baryon chemical potential keeping the temperature fixed.
Definitions: We use a mostly plus metric signature and natural units i.e. = c = k B = 1.

II. EMD HOLOGRAPHY AND THERMODYNAMICS
In this section we give a very brief review of the basic properties of the EMD holographic model of Ref. [23] and its main thermodynamic results. For detailed discussions including technical details on numerics we refer the interested reader to consult the aforementioned work.
The construction of bottom-up dilatonic holographic models seeded with some phenomenological inputs in order to mimic the behavior of the QGP at finite temperature has been originally proposed in the seminal work of Ref. [41] and extended to nonzero baryon densities in [42,43]. The dilaton field is used to dynamically break conformal symmetry in a very specific way dictated by the phenomenological inputs used to engineer the holographic model. Gauge/gravity [44][45][46][47] models constructed in this way may be then used to provide estimates for a large variety of physical observables in the QGP [23,28,29,42,43,[48][49][50][51][52][53][54][55].
The bulk action of the EMD model is given by where κ 2 5 ≡ 8πG 5 is the five dimensional Newton's constant. The action (1) is accompanied by some boundary terms which will play no role for the calculations to be carried out here and, therefore, we omit them. In (1), the dilaton potential V (φ) and the Maxwell-dilaton coupling f (φ) are free functions of the bottom-up EMD construction. There are also two free parameters given by the gravitational constant κ 2 5 and the radius L of the five dimensional asymptotically Anti-de Sitter (AdS) spacetime. Without any loss of generality, we set L = 1 and consider instead a scaling factor Λ corresponding to a fixed energy scale employed to convert observables with mass dimension p evaluated on the gravity side of the holographic duality, which are naturally computed in units of L −p , to physical units of MeV p on the gauge theory side of the holographic correspondence.
The free functions and parameters of the bottom-up EMD model were dynamically fixed in Ref. [23] by matching the holographic equation of state and second order baryon susceptibility evaluated at µ B = 0 to the corresponding QCD results obtained in state-of-the-art lattice simulations [31,56,57], We remark that the results for any other observable, apart from the ones used to fix the EMD parameters in Eq. (2), follow as a legitimate prediction of the EMD holographic model.
Lattice  [30,31]: dimensionless baryon susceptibilities χn ≡ ∂ n (P/T 4 )/∂(µB/T ) n (top), pressure difference (bottom left), and baryon charge density (bottom right) as functions of temperature for different values of the baryon chemical potential. These figures are taken from Ref. [23] (the colored curves in the bottom plots disclose how different Taylor series truncations approximate the full holographic results at different values of µB/T ).
The static and homogeneous charged black hole backgrounds describing isotropic and translationally invariant thermal states at finite density in the gauge theory are obtained by numerically solving the EMD equations of motion, where the Ansatz for the EMD fields may be written as follows [42], The black hole horizon is located at the value of the radial coordinate corresponding to the largest zero of h(r H ) = 0, while the boundary of the asymptotically AdS 5 spacetime lies at r → ∞. The set of coupled second order differential equations of motion for the EMD fields may be numerically solved following the method discussed in Ref. [23] by choosing values for a pair of initial conditions corresponding to the value of the dilaton field at the horizon (φ 0 ) and the value of the derivative of the Maxwell field at the horizon (Φ 1 ). Each value chosen for the pair of initial conditions (φ 0 , Φ 1 ) generates a numerical charged black hole background associated with a definite thermal state at finite baryon density in the gauge theory. By constructing an ensemble of charged black holes, one populates the (T, µ B ) phase diagram and may then proceed to evaluate a wealth of different physical observables.
In Fig. 1 we show the results originally obtained in Ref. [23] for the EMD thermodynamics at finite temperature and baryon density. The result for the second order baryon susceptibility χ 2 strictly at zero chemical potential is not a prediction of the EMD model as it was dynamically fitted to lattice QCD results in order to fix the parameters of the EMD model as mentioned above. On the other hand, all the other results in Fig. 1 are true predictions of the EMD model. One can see that the EMD results are in excellent quantitative agreement with available lattice QCD results for the baryon susceptibilities up to order six, besides also providing a new prediction for the eighth order baryon susceptibility. As far as we know, when it comes to equilibrium quantities, the only other model currently available in the literature which has also achieved this level of quantitative agreement with lattice results is the Cluster Expansion Model (CEM) of Ref. [58] 4 . Also in Fig. 1 we display the EMD predictions for the pressure and baryon charge density at finite temperature and baryon chemical potential. Once again one observes an excellent quantitative agreement with state-of-the-art lattice QCD calculations.
By analyzing the behavior of the baryon susceptibilities in the EMD model beyond the region of the phase diagram currently probed by lattice simulations, in Ref. [23] this model was used to predict the QCD CEP to be located at (T CEP , µ CEP B ) ∼ (89, 724) MeV. As mentioned in Sec. I, a chemical freeze-out analysis was also carried out in Ref. [23] and the following window for the critical region in terms of the center of mass energy of the collisions was found: Though these values of √ s NN are lower than the upcoming Beam Energy Scan II at RHIC, planned fixed-target experiments also at RHIC [35] and future low energy heavy ion collisions at FAIR/GSI [36] should be able to investigate the baryon rich system formed in these low energies. We close this section by remarking that our motivation here to use the bottom-up EMD holographic model of Ref. [23] is that it provides a very practical way to estimate how a strongly coupled QGP at finite temperature behaves at moderately large baryon density in and out of equilibrium. This is still beyond the scope of first principle lattice QCD calculations and, thus, our knowledge about the QGP in this regime must rely on effective models that display some of the QGP's well-known essential features, such as its nearly perfect fluidity in the crossover region. Moreover, in this type of framework, it is mandatory to have at least the thermodynamics of the system (in the regime amenable to the lattice) under investigation properly described, a nontrivial feature displayed by the model of Ref. [23] as discussed above. Also, we note that we are assuming that in the range of the (T, µ B ) plane we are interested in the system is still sufficiently strongly coupled and behaves as a nearly perfect fluid. This is clearly not the case at very large temperatures (where a weakly coupled description must hold [59]) or at very low temperatures in the hadron dominated phase.
As discussed in detail in Ref. [43], for an EMD model at finite temperature and chemical potential in the homogeneous regime of zero wavenumber disturbances, the physically relevant gauge and diffeomorphism invariant linearized perturbations of the system are organized into different representations of the SO(3) rotational symmetry: the quintuplet, triplet, and singlet channels. These perturbations generally possess leading non-normalizable, and subleading normalizable modes at the boundary of the bulk spacetime.
One may obtain different transport coefficients of the medium using the leading modes of these perturbations. More precisely, one solves the linearized equations of motion for the perturbations 5 subjected to the conditions that these disturbances are in-falling at the black hole horizon and normalized to unity at the boundary. In the in-falling Eddington-Finkelstein (EF) coordinates, which we shall use in this section, the causal condition for obtaining in-falling modes at the black hole horizon is simply expressed by requiring that the solutions are regular at the horizon. In this case, one may evaluate the shear viscosity by making use of a Kubo formula that relates this transport coefficient to the retarded 2-point Green's function of the non-normalizable perturbation in the SO(3) quintuplet channel. Similarly, by working with the perturbations of the SO(3) triplet and singlet channels, one may obtain through the use of Kubo formulas the baryon conductivity [51] and the bulk viscosity [28] of the medium, respectively.
On the other hand, by working with the subleading normalizable modes of the perturbations one may obtain the spectra of QNM's of the theory. Specifically, one needs to solve the linearized equations of motion for the perturbations requiring regularity at the horizon while imposing a Dirichlet condition for these perturbations at the boundary. This means that such equations will admit solutions only for a discrete set of complex eigenfrequencies ω, which are the QNM's of the theory. In general, the QNM's ω depend on the wavenumber k of the perturbations and also on the properties of the black hole background, with the latter translating in a dependence on T and µ B . In the present work we are only interested in analyzing the upper bounds for characteristic equilibration times of the medium in the regime of linear perturbations and, for this sake, it suffices to consider only the lowest QNM's in the homogeneous regime with k = 0, which we shall explore in the course of the next sections.
The gauge and diffeomorphism invariant EMD linearized perturbations in the homogeneous regime of disturbances with zero wavenumber, and their corresponding equations of motion, were originally obtained in Ref. [43]. In the SO(3) quintuplet channel the relevant perturbation is given by χ ≡ h ij , where h µν is the perturbation of the metric field and h ij is any of the five traceless spatial components of the graviton. This case will be analyzed in Section III A. In the SO(3) triplet channel the relevant perturbation a ≡ a i represents any of the three spatial components of the perturbation a µ of the Maxwell field -this shall be studied in Section III B. In the SO(3) singlet channel the relevant perturbation is given by where ϕ is the perturbation of the dilaton field φ. We note that the background dilaton field couples the dilaton perturbation to the spatial trace of the graviton through this S-perturbation.
Moreover, the S-perturbation inherits the same ultraviolet asymptotics of the background dilaton field [40,43,65], which is the reason why it has been dubbed in Ref. [40] the "dilaton channel". Furthermore, the SO(3) quintuplet channel may be also called the "external scalar channel" since its equation of motion, to be discussed in what follows, equals that of a massless external scalar field placed on the EMD backgrounds [29,43,70]. Similarly, the SO(3) triplet channel may be also referred as the "vector diffusion" channel, since this channel is related to the charge conductivity and the diffusion coefficient [43,51].
As discussed in Ref. [40], in a more general far-from-equilibrium homogeneous (but highly anisotropic and nonlinear) configuration, each of these channels are associated with physical observables describing the full time evolution of the medium. Namely, the SO(3) quintuplet, triplet, and singlet channels are associated with the pressure anisotropy, the charge density, and the scalar condensate (dual to the bulk dilaton field), respectively. Under fairly general circumstances, the approach of each of these observables toward thermal equilibrium may provide many different characteristic relaxation time scales for the medium [71]. The approach of the pressure anisotropy toward zero gives us the so-called "isotropization time" of the system, while the last equilibration time of the medium is to be naturally identified as the actual "thermalization time" of the system.
As explicitly shown in Ref. [40] for a top-down conformal EMD model with a critical point describing an N = 4 Super Yang-Mills plasma at finite density, the late time behavior of the pressure anisotropy is quantitatively described by linearized oscillations given by the lowest QNM of the SO(3) quintuplet channel [70]. Moreover, in that case, the scalar condensate was found to be always the last observable to equilibrate and, as such, the thermalization time was associated with the approach of the scalar condensate to its equilibrium value, which turned out to be quantitatively described by the lowest QNM of the SO(3) singlet channel. However, one must keep in mind that the isotropization and thermalization times of the system can only be computed by numerically simulating the full nonlinear evolution of the medium, i.e. they cannot be obtained by considering just the characteristic equilibration times associated with the late linear oscillations of the system described in the QNM analysis.
In the present work, even though we are still not investigating far-from-equilibrium properties of the phenomenologically realistic EMD model of Ref. [23], from previous considerations in the literature we expect that the behavior of the lowest QNM's in the SO(3) quintuplet, triplet, and singlet channels will provide correct descriptions of the late time dynamics of the pressure anisotropy, baryon charge density, and the scalar condensate, respectively. In this regard, the characteristic equilibration times of the medium, which we shall obtain here by analyzing the lowest homogeneous non-conformal QNM's of each channel in the EMD model up to the critical region, should be more properly seen as "lower bounds" for the full relaxation times of the medium, since they do not take into account the time spent by the system in earlier nonlinear stages. On the other hand, when considering just small perturbations around thermal equilibrium, the characteristic equilibration times extracted from the QNM analysis indeed provide upper bounds for how fast the disturbed system returns to equilibrium.
Before we proceed to the evaluation of the QNM's in each of the aforementioned channels, let us first specify the in-falling EF coordinates by the relation given below,  where v is the EF time. With this, one may easily translate the equations of motion for the homogeneous EMD perturbations derived in Ref. [43] in domain-wall coordinates to the EF coordinates. This will be done in the following sections. Some technical details on the grids employed to obtain the lowest QNM's in each channel are deferred to Appendix A.

A. SO(3) quintuplet channel
The equation of motion for the χ-perturbation of the quintuplet channel in EF coordinates reads [29], This is a linear second order differential equation which must be solved with appropriate boundary conditions in order to find the corresponding QNM's. As discussed before, the boundary conditions to be imposed are that the solutions must be regular at the black hole horizon and vanish at the boundary.
In order to numerically integrate the equation of motion (5) from the horizon up to the boundary 6 we need to specify the values of the χ-perturbation and its derivative at the horizon. One may simply set χ(r H ) = 1, while the expression for its derivative may be obtained by Taylor-expanding the χ-perturbation and the background around the horizon and plugging these expansions back into Eq. (5), which then becomes an algebraic equation for χ (r H ) whose solution is given by [29] where A 1 is the value of the derivative of the background field A(r) evaluated at the horizon (this is itself a function of the initial conditions (φ 0 , Φ 1 ) [23]). Now that we have the initial conditions (χ(r H ), χ (r H )) required to initialize the numerical integration of Eq. (5), we take the following steps [29]: i. We construct a grid of backgrounds on top of which we shall obtain the QNM's as functions of (T, µ B ). The background grid considered in this work is discussed in Appendix A; ii. We construct on each point within the aforementioned background grid a rough grid of complex frequencies ω which are seeded to the differential equation (5) to be numerically integrated on each point within this ω-grid.
The rough ω-grid we consider in this work has a numerical step-size of 0.03 between consecutive points both in the real and imaginary directions of the complex ω-plane; iii. Once the previous step is accomplished we apply a shooting method upon the ω-grid generated for each point of the background grid. In this shooting method we begin by picking some value for |χ| evaluated at the boundary of the background spacetime and then we lower its value until we can clearly identify isolated clusters of complex eigenfrequencies with small values of the perturbation at the boundary. Since we are only interested in computing the lowest non-hydrodynamic QNM ω 0 , we pay attention to the cluster with lowest imaginary part (in absolute value) and keep decreasing the boundary value of |χ| until only one point remains within that cluster; iv. Next we take these rough estimates of ω 0 for each background point as the centers of finer ω-grids constructed with a numerical step-size of 0.002 between consecutive points both in the real and imaginary directions of the complex ω-plane. We solve again Eq. (5) on top of these finer ω-grids and, for each point on the background grid, we store the eigenfrequency associated with the lowest value of |χ| at the boundary.
In principle, one may repeat this process of ω-grid refinement until the desired numerical accuracy is reached. Moreover, by enlarging the ω-grid size one may also identify excited QNM's besides the lowest one 7 . The same shooting method for obtaining the QNM's, which has been previously employed in Ref. [29], will also be used in the next sections.
The results for the lowest QNM ω 0 of the quintuplet channel up to the critical region are shown in Fig. 2, where we also plot the results for the corresponding equilibration time given by minus the inverse of the imaginary part of ω 0 . As a basic consistency check of the numerics, we note that at high T the real and imaginary parts of ω 0 do return to the correct ultraviolet conformal results [29,63]. One also observes that at fixed temperatures the equilibration time in this channel always decreases with increasing the baryon chemical potential.
From the bottom plot in Fig. 2 one can immediately read some reference values for the equilibration time of the quintuplet channel, e.g. at (T, µ B ) ∼ {(400, 0), (145, 0), (89, 724)} MeV the corresponding equilibration times are, respectively, τ eq ∼ {0.06, 0.26, 0.33} fm/c, with the last point corresponding to the equilibration time at the CEP. By comparing this result with those from the next two sections, we shall see that the equilibration time in the quintuplet channel is the shortest relaxation time of the system.

B. SO(3) triplet channel
The equation of motion for the a-perturbation of the triplet channel in EF coordinates reads  As before, one may simply set a(r H ) = 1 and work out its derivative at the horizon,

T [MeV]
where φ 1 is the value of the derivative of the dilaton field evaluated at the horizon (this is itself a function of the initial conditions (φ 0 , Φ 1 ) [23]). The general steps used to numerically obtain the QNM's are the same discussed in the previous section. In Fig. 3 we show the results for the lowest non-hydrodynamical QNM in the triplet channel, as well as the corresponding equilibration times. As a basic consistency check of the numerics, one notes that at high T the real and imaginary parts of ω 0 go to their correct ultraviolet conformal results, which are known analytically for the vector diffusion channel [37]. In the present case, one notes that for T 215 MeV the equilibration time is reduced by increasing the baryon chemical potential at fixed temperature but at T ∼ 215 MeV all the curves with different fixed µ B meet at the same crossing point, while for T 215 MeV the equilibration time slightly increases with increasing baryon chemical potential at fixed temperatures. This is more clearly illustrated in Fig. 4, where we also As in the previous sections, one may set S(r H ) = 1 and work out its derivative at the horizon, These are the initial conditions used to numerically integrate the equation of motion (9) following the shooting procedure discussed before.
In Fig. 5 we display the results for the lowest non-hydrodynamical QNM of the singlet channel and the corresponding equilibration times. As discussed in Refs. [40,43,65], the S-perturbation inherits the same ultraviolet asymptotics of the background dilaton field and, therefore, the ultraviolet conformal behavior of the QNM's in the dilaton channel depends on the scaling dimension ∆ of the gauge theory operator dual to the dilaton. In the case of the top-down conformal EMD model with finite chemical potential and a critical point investigated in Ref. [40] ∆ = 2, while in the present phenomenological EMD model at finite baryon density ∆ ≈ 2.73 [23]. Consequently, the conformal behavior attained in the ultraviolet by QNM's in the dilaton channels of both EMD models are different.
One notes from Fig. 5 that the equilibration time in the dilaton channel generally decreases with increasing µ B at fixed T . Also from this figure one can read some reference values for the equilibration time in the singlet channel, e.g. at (T, µ B ) ∼ {(400, 0), (145, 0), (89, 724)} MeV the corresponding equilibration times are, respectively, τ eq ∼ {0.11, 0.70, 1.26} fm/c, with the last point corresponding to the CEP. By comparison with the other two channels, the singlet channel gives the longest equilibration time of the medium for the first point above.

IV. CONCLUSIONS AND OUTLOOK
In this work we investigated the homogeneous limit of the lowest non-hydrodynamic QNM's of the EMD model of Ref. [23] at finite temperature and baryon density up to the critical region. We then computed the associated equilibration times in the SO(3) rotationally invariant quintuplet, triplet, and singlet channels. The equilibration times in the different channels come closer to each other at high temperatures, although being well separated at the critical point. In most cases, these equilibration times decrease with increasing baryon chemical potential (while keeping temperature fixed).
We must remark that the homogeneous QNM's evaluated in the present work are not expected to give a precise description of characteristic equilibration times of the QGP under the conditions actually realized in heavy ion collisions since the latter corresponds to a medium that rapidly expands in spatial directions and, in the present analysis, we considered a non-expanding plasma which is slightly out-of-equilibrium.
The first sequel of the present work which we intend to pursue in the near future comprises an investigation of the far-from-equilibrium homogeneous isotropization dynamics of the EMD model. This would allow us to check whether the QNM's obtained here in the SO(3) quintuplet, triplet, and singlet channels indeed describe the late time behavior of the pressure anisotropy, baryon charge density, and the scalar condensate, respectively, of the full nonlinear system of equations. Moreover, in this setting we shall also be able to evaluate the homogeneous isotropization and thermalization times of the EMD system near the critical region.
Future projects also comprise numerical simulations in the EMD model of more realistic expanding scenarions involving for instance a boost-invariant Bjorken flow [72][73][74][75][76][77] and also collisions of holographic shock waves [71,[78][79][80][81][82][83][84]. The latter may provide new insight into the complex interplay between the early time non-equilibrium dynamics and the late time hydrodynamic behavior of the strongly coupled baryon rich QGP near a critical point, which could be formed in upcoming low energy heavy ion collisions at RHIC and FAIR.

T [MeV]
Re[ω 0 ] / 2πT Another interesting question concerns the QNM's of the anisotropic EMD model at finite temperature and magnetic field proposed in Ref. [54], whose thermodynamics is in excellent agreement with lattice QCD results. This setup, in particular, may be of great relevance to high energy peripheral heavy ion collisions where very intense magnetic fields of order eB 0.3 GeV 2 may be produced at the earliest stages of such non-central collisions [85].