Energy loss, equilibration, and thermodynamics of a baryon rich strongly coupled quark-gluon plasma

Lattice data for the QCD equation of state and the baryon susceptibility near the crossover phase transition (at zero baryon density) are used to determine the input parameters of a 5-dimensional Einstein-Maxwell-Dilaton holographic model that provides a consistent holographic framework to study both equilibrium and out-of-equilibrium properties of a hot and {\it baryon rich} strongly coupled quark-gluon plasma (QGP). We compare our holographic equation of state computed at nonzero baryon chemical potential, $\mu_B$, with recent lattice calculations and find quantitative agreement for the pressure and the speed of sound for $\mu_B \leq 400$ MeV. This holographic model is used to obtain holographic predictions for the temperature and $\mu_B$ dependence of the drag force and the Langevin diffusion coefficients associated with heavy quark jet propagation as well as the jet quenching parameter $\hat{q}$ and the shooting string energy loss of light quarks in the baryon dense plasma. We find that the energy loss of heavy and light quarks generally displays a nontrivial, fast-varying behavior as a function of the temperature near the crossover. Moreover, energy loss is also found to generally increase due to nonzero baryon density effects even though this strongly coupled liquid cannot be described in terms of well defined quasiparticle excitations. Furthermore, to get a glimpse of how thermalization occurs in a hot and baryon dense QGP, we study how the lowest quasinormal mode of an external massless scalar disturbance in the bulk is affected by a nonzero baryon charge. We find that the equilibration time associated with the lowest quasinormal mode decreases in a dense medium.


Introduction
In high energy heavy ion collisions [1][2][3][4][5] where nonzero baryon density effects are small, it is currently understood that the quark-gluon plasma (QGP) formed in these reactions behaves as a nearly perfect liquid [6] where viscous effects are surprisingly small. For instance, current estimates based on hydrodynamical modeling (which does not include dynamical effects from a nonzero baryonic charge) find η/s between 0.095 [7] and 0.2 [8,9], where η is the shear viscosity and s is the entropy density of the plasma. These values are very close to the uncertainty principle estimate derived three decades ago in [10] and, more importantly, to the celebrated 1/4π result found in strongly coupled holographic non-Abelian plasmas with large number of colors [11,12]. Such a small value for this ratio is an order of magnitude below the result of calculations in transport models involving ordinary hadrons [13] or in QCD at weak coupling [14]. Alternative physical mechanisms involving the effects of additional heavy resonances near the phase transition [15,16] or Polyakov loop degrees of freedom [17] have been devised to obtain a small η/s near the QCD crossover transition [18], though it remains unclear what is the main non-perturbative element that is responsible for the putative near perfect fluidity of the QGP. Much less is known about the perfect fluid nature of the QGP when not only the temperature T but also the baryon chemical potential µ B is large.
In the last years, the RHIC low energy scan has been used to study heavy ion collisions at lower energies with the intention to probe the QCD phase diagram when µ B T while also (hopefully) finding signals that could confirm the existence of a critical point [19,20]. In general, the phase diagram of QCD may be studied as a function of many other variables besides T and µ B such as the other chemical potentials associated with conserved charges, external electric and magnetic fields, quark masses, number of colors N c and flavors N f , and others (for reviews see, for instance, [21][22][23]). An important question that requires further investigation is whether the baryon rich QGP studied experimentally via the RHIC beam energy scan (and in a few years by the future FAIR facility at GSI) is also a strongly coupled liquid that displays nearly perfect fluidity. This seems to be the case given the large value of elliptic flow observed by the low energy scan at RHIC [24], which may indicate that even the hadron gas formed in the last stages of the collisions can exhibit strong collective behavior at large baryon densities [25]. Therefore, one may say that there is now experimental evidence that the significant collective behavior observed in both high and low energies heavy ion collisions is driven by the formation of a strongly coupled, nearly inviscid QGP.
Even though lattice QCD is currently the main non-perturbative tool to understand strongly interacting QCD physics, when µ B = 0 lattice approaches suffer from the famous sign problem of the fermion determinant, which makes the Monte Carlo-based numerical simulations considerably more challenging. While alternative techniques have been developed to circumvent this issue (see, for instance, [26][27][28][29][30][31][32][33][34]) making it possible at least in principle to investigate small to moderate chemical potentials on the lattice, it is certainly desirable to have at one's disposal complementary theoretical tools to study strongly coupled QCD phenomena at nonzero temperature and density both in equilibrium as well as out of equilibrium.
by DeWolfe, Gubser, and Rosen in [42] as a phenomenological gravity dual for the QGP at finite temperature and baryon chemical potential in the strongly coupled regime. This model contains three fields: the bulk metric g µν , dual to the stress-energy tensor of the boundary quantum field theory, an Abelian vector field A µ whose temporal component's boundary value plays the role of the baryon chemical potential, and a real scalar (dilaton) field φ, with a nontrivial profile in the bulk that is responsible for the dynamical breaking of conformal symmetry in the gauge theory, emulating the effects of a dynamically generated infrared scale such as Λ QCD . In Ref. [42], an estimate was obtained for the QCD critical point in the (T, µ B ) plane and the corresponding critical exponents were calculated. Moreover, in Ref. [43], the same holographic setup was used to calculate the bulk viscosity and the baryon charge conductivity for the corresponding hot and dense holographic plasma 2 . We also point out that recently in Ref. [45] an anisotropic Einstein-Maxwell-Dilaton holographic model was constructed to describe a holographic plasma under the influence of a strong external magnetic field and the corresponding holographic equation of state and crossover temperature dependence on the external magnetic field were shown to be in good quantitative agreement with recent lattice data [46] for values of the magnetic field relevant for noncentral ultrarelativistic heavy ion collisions.
These bottom-up models, though not directly derivable from top-down string theory constructions, rely on the conjecture regarding the validity of the holographic dictionary under more general circumstances in a way that closely resembles the majority of the ongoing efforts to apply holographic approaches to strongly interacting condensed matter systems [47]. The key idea here is that these 5-dimensional, holographic effective models for the strong coupling limit of many-body QCD phenomena are able to describe the nearly perfect fluidity properties of strongly coupled non-Abelian gauge theories through the physics of black holes. Given that nearly perfect fluidity is a fundamental property of holography [11,12] and that these bottom-up models can be consistently tuned to describe the known equilibrium properties of the QGP calculated on the lattice [48], their current application in the calculation of non-equilibrium phenomena may not only be natural but also necessary in order to obtain quantitative predictions near the crossover transition that can be later used in heavy ion phenomenology, such as done in [38].
In the present work, we use the general holographic setup proposed in [42,43] to obtain some equilibrium and out-of-equilibrium properties of a hot and baryon rich strongly coupled QGP. We compare our holographic equation of state computed at nonzero baryon chemical potential, µ B , with recent lattice calculations [49] and find surprising quantitative agreement for the pressure and the speed of sound for µ B ≤ 400 MeV. Moreover, we find that the crossover temperature decreases with µ B in our model in a way that quantitatively matches the behavior found on the lattice [50,51]. This shows that this holographic construction provides a consistent framework to study both finite temperature and nonzero baryon density effects at strong coupling. We obtain holographic predictions for the temperature and µ B dependence of the drag force and the Langevin diffusion coefficients associated with heavy quark jets as well as the jet quenching parameterq and the shooting string energy loss of light quarks in this baryon dense holographic plasma. We find that all these quantities generally display a nontrivial, fast-varying behavior as a function of the temperature near the crossover. We also find that the energy loss of heavy and light quarks in the plasma generally increases in the dense QGP as one increases µ B . Thus, even though this strongly coupled liquid cannot be described in terms of well defined quasiparticle excitations, the general expectation that jet probes should lose more energy as they plow through a baryon dense medium remains valid at strong coupling. Moreover, we initiate a study on how thermalization occurs in a hot and baryon dense QGP by computing the lowest quasinormal mode of an external massless scalar disturbance in the bulk (which gives information about the linearized approach to equilibrium of different hydrodynamic channels in the spatially uniform limit) taking into account a nonzero µ B . We find that the equilibration time associated with this lowest quasinormal mode decreases in a baryon dense medium, which may be relevant for the phenomenology of low energy collisions currently under investigation at RHIC. This paper is organized as follows. In Section 2 we review in detail the bottom-up holographic setup proposed in Refs. [42,43]. In Section 3 we fix the unknown functions and parameters of the model using (2 + 1)-flavor lattice QCD data for the equation of state and baryon susceptibility at µ B = 0 from Refs. [49,52], compute the equilibrium properties of the plasma when µ B = 0, and finally compare it with lattice results at nonzero baryon density [49]. In Section 4 we compute the energy loss of heavy quarks and light quarks in this model. Our investigation regarding the equilibration of the plasma at nonzero µ B is done in Section 5. Some rather technical discussions on the structure of our holographic model are deferred to Appendix A. We finish this paper in Section 6 with a discussion of our main results and we also point out some related studies to be pursued in the future.

The holographic model
We begin this section by reviewing the main points in the derivation of the classical backgrounds proposed in the approach of Refs. [42,43] to holographically model QCD thermodynamics at finite temperature and baryon chemical potential. The slight modifications to this approach which shall be implemented here, regarding the determination of the free functions and parameters of the holographic model, will be discussed in Section 3.1.
The Einstein-Maxwell-Dilaton holographic model of [42,43] is given by the following where S GHY is the Gibbons-Hawking-York boundary term [53,54] necessary to establish a well-defined variational problem with Dirichlet boundary condition for the metric field, and S CT is the counterterm action needed to regularize the ultraviolet divergences of the on-shell action, which is done through the holographic renormalization procedure [55][56][57][58][59]. 4 Since these two boundary terms only contribute to the on-shell action, which will not be needed in this work, we do not write down their explicit form. The holographic model in (2.1) contains two unknown functions: the dilaton potential V (φ) and the Maxwell-dilaton gauge coupling, f (φ). These functions, together with the value of the 5-dimensional gravitational constant G 5 and the temperature scale in MeV, will be fixed in Section 3.1 using the µ B = 0 lattice data.
In order to holographically model a field theory at finite temperature and nonzero chemical potential, we take the following charged (and spatially isotropic) black brane Ansatz for the bulk gravity fields where the radial location of the black hole horizon, r H , is given by the largest simple root of the equation h(r H ) = 0 (we consider coordinates where the boundary of the asymptotically AdS 5 space is located at r → ∞). Also, for simplicity, we set to unity the AdS radius.
With Ansatz (2.2), the equation of motion for the dilaton field following from the action (2.1) is where the prime denotes derivative with respect to the radial direction. The relevant equation of motion for the Maxwell field is Combining the independent components of Einstein equations, we arrive at the following set of equations of motion for the metric One could also consider adding to the action an Abelian Chern-Simons term of the form A ∧ F ∧ F .
However, as discussed in [42], this term does not contribute to the classical equations of motion and it vanishes on-shell for the solutions considered here, which have only electric but no magnetic charge.
We will use the set of equations of motion (2.3)-(2.6) and the constraint (2.7) to numerically solve for the unknown functions A(r), h(r), φ(r), and Φ(r). Due to reparametrization invariance of the radial coordinate, B(r) may be conveniently chosen in order to simplify the calculations. We will specify a convenient metric gauge in Sections 2.1 and 2.2.
There are two conserved charges related to the equations of motions: the Gauss charge Q G , and the Noether charge Q N , which are given by (2.8) The equation of motion for the Maxwell field, Eq. (2.4), may be written as dQ G /dr = 0, while the equation of motion (2.6) for the blackening function, h(r), may be written as dQ N /dr = 0. Since these are conserved charges in the radial direction, one may evaluate (2.8) at any value of r.
We should also remark the limitations of the holographic model discussed here. This model cannot describe phenomena directly related to the breaking of chiral symmetry and its restoration at finite temperature and density, for which one would need to include flavored branes in the bulk [44]. Also, the holographic setup discussed here cannot describe hadron thermodynamics: the pressure of the holographic plasma in the deconfined phase goes like ∼ N 2 c , while the pressure in the confined hadronic phase goes like ∼ N 0 c , requiring quantum string corrections in the bulk in order to be properly taken into account 5 . Furthermore, our holographic dual defined in the classical gravity limit asymptotes to a strongly coupled ultraviolet fixed point where conformal invariance is restored and, therefore, the gauge theory investigated here is not asymptotically free at very high temperatures or densities (as it is the case in QCD). Therefore, this model may be applicable in the description of the strongly coupled QGP with T ∼ 150 − 400 MeV and, as the results in Section 3.2 will show, µ B = 0 − 400 MeV, 6 which is perhaps the regime in which there is still great theoretical uncertainty in QGP modeling. 5 Note also that this model does not display confinement in the strict sense of an area law for the Wilson loop at zero temperature and zero density, as expected for a boundary non-Abelian gauge theory with dynamical fermions. By suppressing the fermion dynamics one obtains a pure glue theory which does confine probe charges according to the Wilson criterion, which may be achieved by adequately choosing the functional form of the dilaton potential (see, for instance, [60][61][62][63]). 6 This defines some kind of "holographic Goldilocks regime": the plasma cannot be too hot or too cold for the approach described here to be useful.

Thermodynamics
For convenience, we now choose the gaugeB(r) = 0 with the AdS 5 boundary located atr → ∞ and the horizon atr =r H . 7 In this gauge, we have (2.9) and the near-boundary, far from the horizon asymptotics for the bulk fields are given by [42]Ã Here ν ≡ d − ∆, where d = 4 is the number of dimensions at the boundary and ∆ = (d + √ d 2 + 4m 2 )/2 is the scaling dimension of the boundary gauge theory operator dual to the dilaton field, with m being the effective dilaton mass obtained from the expansion of the dilaton potential close to the boundary. For the potential we consider here (which will be specified in Section 3.1), ∆ ≈ 3, and the bulk dilaton field is dual to a relevant operator in the boundary field theory that creates a renormalization group flow from an ultraviolet fixed point towards a nonconformal state in the infrared.
The temperature in the dual field theory is given by the Hawking temperature of the black hole solution 8T The entropy density is given by the Bekenstein-Hawking formula [64,65] where κ 2 ≡ 8πG 5 . The chemical potential is given by the boundary value of the Maxwell fieldμ The tilde is used here to denote the standard coordinates associated with the conditionB(r) = 0 where the blackening function goes to unity at the boundary. In these coordinates, we will obtain standard expressions for thermodynamical quantities such as the temperature and entropy density. However, in order to numerically solve the equations of motion for the bulk fields, we will need to rescale these coordinates and we reserve the variables without the tilde to denote the corresponding rescaled "numerical" coordinates, which we shall specify in Section 2.2. 8 We use a "hat" to designate thermodynamical observables expressed in terms of the quantities on the gravity side. The counterparts of these observables without the hat will be used to refer to thermodynamical quantities expressed in powers of MeV, as will be discussed in Section 3.1.
while the charge density is associated with the boundary value of the radial momentum conjugate to the gauge field 9 where we used Eqs. (2.8) and (2.10).

Numerical procedure
In order to numerically solve the equations of motion (2.3)-(2.6) we follow [42] and adopt numerical coordinates where the near-horizon Taylor expansions for the bulk fields X = {A, h, φ, Φ} are given by X(r) = ∞ n=0 X n (r − r H ) n with The location of the horizon fixed at r H = 0 may be obtained by rescaling the radial coordinate. Clearly, h 0 = 0 comes from the fact that, by definition, h(r) has a simple zero at the horizon while h 1 = 1 may be obtained by rescaling t. Also, A 0 = 0 may be obtained by rescaling (t, x) by a common factor. Furthermore, Φ 0 = 0 is required in order to have a well defined gauge field A = Φ(r)dt at the horizon 10 . The reason for using this set of numerical coordinates comes from the fact that in order to numerically solve the equations of motion we need to specify numerical values for each parameter in the near-horizon Taylor expansions for the bulk fields, which is possible after implementing the coordinate rescalings discussed above. With Eq. (2.15) at hand, all the remaining coefficients in the near-horizon Taylor expansions for the bulk fields may be determined from a two-parameter initial condition, (φ 0 , Φ 1 ), by solving the equations of motion order by order in these expansions.
The numerical strategy we adopt here to solve the equations of motion (2.3)-(2.6) proceeds as follows 11 . In order to avoid the singular point of the equations of motion at the horizon r H = 0, we start our numerical integration slightly above it at r start = 10 −8 , and use second order near-horizon expansions, X(r start ) = X 0 + X 1 r start + X 2 r 2 start + O(r 3 start ), to initialize the numerical integration going up to r max = 2, which we choose to be the ultraviolet cutoff of our calculations 12 . Using Eq. (2.15) we then determine the remaining coefficients of the second order near-horizon expansions, A 1 , A 2 , h 2 , φ 1 , φ 2 , and Φ 2 as functions of the initial conditions, (φ 0 , Φ 1 ), by substituting the second order 9 In Eq. (2.14) the metric determinant is included in the definition of the Lagrangian density, i.e. S = M 5 d 5x L. 10 As mentioned in [66], dt has infinite norm at the horizon such that if Φ(rH ) = Φ0 = 0 one would find that A = Φ(r)dt would be ill-defined at the horizon. 11 Of course, to numerically solve the equations of motion we need to specify the dilaton potential, V (φ), and the Maxwell-dilaton gauge coupling, f (φ), which will be done in Section 3.1. 12 At this value of the radial coordinate our numerical solutions already reach the ultraviolet fixed point corresponding to the AdS5 geometry. . The boundary conditions specified near the horizon, which are necessary to initialize the numerical integration of the second order differential equations (2.3)-(2.6), are then given by X(r start ) and X (r start ).
Any asymptotically AdS 5 geometry must satisfy R(r max ) = −20. For each value of φ 0 , there is a bound on the maximum value of Φ 1 above which the solutions are not asymptotically AdS 5 . In order to determine this bound, we note that in the gauge B(r) = 0 the equation of motion (2.5) gives A (r) = −φ (r) 2 /6 ≤ 0 so that A(r) is a concave function of the radial coordinate. Since for asymptotically AdS 5 geometries in the gauge B(r) = 0 the function A(r) must increase for large r and, since it is a concave function of r, this means that A(r) must increase monotonically from the horizon towards the boundary. This, in turn, means that the derivative of A(r) at the horizon is positive, i.e., for asymptotically AdS 5 geometries one must have A 1 > 0. By plugging the near-horizon expansions into the constraint equation (2.7) and evaluating it at the horizon using the conditions (2.15) we obtain Since the dilaton potential we will work with is negative definite and the Maxwell-dilaton coupling is positive definite, and since A 1 > 0 for asymptotically AdS 5 geometries, Eq. (2.16) gives us the bound 13 [42]

Coordinate transformations
Here we will express various relevant quantities in the standard coordinates in terms of the numerical coordinates of the gauge B(r) = 0 discussed above. In order to do so, we first need to inspect the far-from-horizon (r → ∞) asymptotics of the bulk fields in the numerical coordinates, which are known to have the following form [42] After evaluating the constraint (2.7) at the boundary and using (2.18), we obtain where we also used that φ → 0 at the boundary where the potential V → −12. The Gauss charge in Eq. (2.8) evaluated at the horizon is given by Q G (r H ) = f (φ 0 )Φ 1 , using (2.15).
When evaluated at the boundary, one finds Q G (r → ∞) = −2f (0)Φ far 2 / h far 0 using (2.18) and (2.19). Since this is a conserved charge in the radial direction, it follows that As we shall see in a moment, besides Φ far 2 that may be calculated using Eq. (2.20) once we determine h far 0 , we also need Φ far 0 and φ A to compute thermodynamical quantities such as the temperature (2.11), the entropy density (2.12), the chemical potential (2.13), and the charge density (2.14) directly in terms of quantities extracted from the numerical solutions for A(r), h(r), φ(r), and Φ(r). The numerical solutions for h(r) and Φ(r) are found to quickly converge to their asymptotic values at large r and, thus, we can set h far 0 = h(r max ) and Φ far 0 = Φ(r max ). On the other hand, we fix φ A through the following near-boundary fit evaluated in the interval r ∈ [r max − 1, r max ]: φ(r) = φ A e −νA(r) . 14 In order to expressT ,μ,ŝ, andρ in terms of h far 0 , Φ far 0 , Φ far 2 , and φ A , we need to derive some relations between the standard and numerical coordinates. Takingφ(r) = φ(r), ds 2 = ds 2 , andΦ(r)dt = Φ(r)dt we obtain 15 by comparing the near-boundary asymptotics (2.10) and (2.18) for r → ∞r (2.27) 14 We checked that A(r) ≈ α(r) within this interval. 15 As mentioned in [42], the following relations are valid for φA > 0. If φA < 0 one must replace φA → |φA|.
Using the above relations and Eqs. (2.11)-(2.15) we can express the temperature, baryon chemical potential, entropy density, and baryon charge density as functions of the coefficients of the near-boundary asymptotics in the numerical coordinateŝ In the next section we will express these thermodynamical quantities in physical units (which shall be then denoted without the hat) using lattice data for the thermodynamic functions and the quark susceptibility at vanishing baryon chemical potential.

Fixing the model parameters and comparison with lattice data at nonzero baryon density
We use lattice data from [49,52] at vanishing baryon chemical potential 16 to determine the unknown functions in our effective holographic model: the dilaton potential and the Maxwell-dilaton gauge coupling. Once this is done, the model is fully specified and any quantity computed at µ B = 0 can be interpreted as a direct prediction from the holographic setup. In fact, later in this section we will compare our results for the thermodynamic functions with µ B = 0 to the corresponding lattice data.
Some technical aspects regarding the form of our dilaton potential and Maxwell-dilaton gauge coupling and their relations with the thermodynamic stability of our black hole solutions will be discussed in Appendix A.

Extracting
The gauge field in the bulk is associated with a conserved charge at the boundary and in our bottom-up model we ensure that this conserved charge plays the role of a baryon charge by matching the susceptibility computed in the model (at zero chemical potential) with the corresponding baryon susceptibility χ B 2 (T, µ B = 0) computed on the lattice. A holographic formula for the susceptibility at zero chemical potential was obtained in Ref. [42] and we will review this calculation below since the final result in (3.8) is needed to fix the Maxwell-dilaton coupling function f (φ). 16 We remark that there is now agreement between different lattice collaborations when it comes to the thermodynamic properties of the QGP at µB = 0 [68,69].
The baryon susceptibility is given by the second derivative of the Helmholtz free energy density F with respect to the baryon chemical potential and, using standard thermodynamical relations, one can express it as [42] which is evaluated at a fixed temperature. At zero chemical potential (and, consequently, zero charge density), we have Since we want to obtain a holographic formula to compute (3.2) using black hole physics, we approximate Φ(r) as a linear perturbation in the equations of motion such that Φ (r) 2 ∼ 0. In this limit, the equations of motion for A(r), h(r), and φ(r) decouple from the equation for Φ(r) and they may be numerically solved yielding uncharged black hole geometries dual to a field theory at finite temperature and zero chemical potential [48].
Using the near-horizon expansions (2.15) and the near-boundary asymptotics (2.18) one can express Φ far 0 as follows Now one needs to find a suitable formula for Q G . Since the Gauss charge (2.8) is conserved in the radial direction, we can take it outside the following integral From the definition of the Gauss charge (2.8) and fixing the metric gauge B(r) = 0, we see that the left hand side of Eqs. (3.4) and (3.5) are equal and, thus .
Also, using (2.28) and (2.30) we obtain 17 It is convenient to consider the dimensionless ratio χ B 2 /T 2 since it asymptotes to a constant at high temperatures.

Finally, plugging (3.6) and (3.7) into (3.3) we arrive at
. (3.8) Note that in (3.8) there is no dependence on the linear perturbation Φ(r) so that χ B 2 (µ B = 0)/T 2 can be evaluated using the solutions corresponding to the uncharged black holes with zero chemical potential without having to explicitly solve the equation of motion for Φ(r).
In this work we have numerically generated ∼ 10 3 uncharged black hole solutions by setting Φ 1 = 0 and varying φ 0 in the holographic setup discussed in the previous sections. Guided by the functional forms for the dilaton potential V (φ) considered in [42] and by the recent lattice data [49] for the speed of sound squared, c 2 s , and the normalized pressure, p/T 4 , in (2 + 1)-flavor QCD at zero baryon chemical potential, we fixed, respectively, the dilaton potential and the gravitational constant to be [38,45] From this potential, we see that the effective dilaton mass is m 2 ≈ −3 (therefore, satisfying the Breitenlohner-Freedman bound [70,71] for massive scalar fields on asympotitcally AdS 5 geometries) and that the scaling dimension of the relevant gauge field theory operator dual to the dilaton is ∆ ≈ 3, as anticipated in Section 2.1.
Some remarks are in order at this point. The speed of sound squared at zero baryon chemical potential was calculated using is the internal energy density. Also, since we have not renormalized the Einstein-Maxwell-Dilaton action in the present work, we did not calculate the pressure directly from the renormalized free energy density. Rather, here we followed [38] and computed pressure differences with respect to a very small reference pressure p 0 evaluated at some low temperature T low by integrating the entropy density calculated via the Bekenstein-Hawking formula (2.12) where we took T low = 22 MeV [38,45]. Then, the trace anomaly was calculated at zero chemical potential using its definition I(T, µ B = 0) ≡ (T, µ B = 0) − 3p(T, µ B = 0). In Fig. 1 we show a comparison between the model calculations for c 2 s , p/T 4 , and the trace anomaly at µ B = 0 and the corresponding lattice data [49]. One can see that the choice of parameters in (3.9) gives an overall good description of the lattice data (though the trace anomaly is underpredicted by the model at high temperatures). s , normalized pressure p/T 4 , and normalized trace anomaly I/T 4 at zero baryon chemical potential as functions of the temperature. The data points correspond to lattice calculations for QCD with physical quark masses from [49].
Guided by the functional forms for the Maxwell-dilaton gauge coupling f (φ) used in [42] and by the recent lattice data from Ref. [52] (which is in agreement with data from other lattice groups [72]) for the baryon susceptibility of (2 + 1)-flavor QCD at zero baryon chemical potential and physical quark masses, we fixed this coupling in this work to be In Fig. 2 we compare the results for χ B 2 from our holographic model with potential V (φ) and coupling f (φ) given by (3.9) and (3.12) with the lattice data from Ref. [52]. Our choice for the Maxwell-dilaton coupling gives a very good description of the data in the crossover transition though the agreement becomes worse at high temperatures.
In order to perform a meaningful comparison of our holographic results with the lattice data, one needs to express the quantities computed using the black hole physics (which are defined in units of the asymptotic AdS radius) in terms of physical units. This can be done Holographic result for the baryon susceptibility at zero baryon chemical potential as a function of the temperature. The data points correspond to lattice calculations for QCD with physical quark masses from [52]. by matching the temperature at which c 2 s displays a minimum in our calculations with the corresponding temperature extracted from the lattice data 18 , which gives the energy scale We can use this to express any dimensionful quantity on the gravity sideX in terms of its counterpart X with mass dimension [MeV p ] through X = Λ pX .
Our procedure for fixing the free parameters of the holographic model, and our results for them discussed above, differ from those presented in Ref. [42]. Here, we fitted the lattice data at zero baryon chemical potential from [49,52] while in Ref. [42] the data used was from an older calculation performed in [73]. Also, in [42] four independent scaling parameters were introduced to express T , µ, s, and ρ in physical units while here we employed a single scaling parameter given by (3.13).
With the dilaton potential and the Maxwell-dilaton gauge coupling fixed by lattice data at µ B = 0, the holographic model is now fully specified and the numerical charged black hole geometries obtained from it may be used to derive holographic predictions for the physics of the dual field theory at nonzero baryon chemical potential. In the next section, we are going to compare our results for the thermodynamics in the presence of nonzero baryon density with the lattice data [49] for baryon chemical potentials up to µ B = 400 MeV.
We close this section by remarking that the forms of the dilaton potential V (φ) and Maxwell-Dilaton gauge coupling f (φ) fixed in Eqs. (3.9) and (3.12), respectively, were only partially constrained by lattice data at µ B = 0 in the range of temperatures T ∼ 100 − 400 MeV. In the deep infrared, for large values of φ (which generically translate into small values of T ), one could in principle deform V (φ) in infinitely many different ways such as to obtain very different results for low T physics, while still keeping basically the same agreement with lattice data in the range of temperatures spanned by T ∼ 100 − 400 MeV. This is the reason why we shall not consider values of T outside this region in the present work. We also remark that, since V (φ) and f (φ) were fixed by lattice data at zero chemical potential, all of our results for nonzero µ B physics constitute genuine holographic predictions. At µ B = 0, the results we shall present for transport observables are also predictions, although the results for the equation of state and baryon susceptibility shown in Figs. 1 and 2, respectively, are not predictions, but instead the results of fits to lattice data used to fix V (φ) and f (φ), as explained before.

Model predictions for the equation of state at µ B = 0
Let us begin this section by briefly reviewing the relevant thermodynamic relations that will be used to obtain the holographic equation of state at nonzero baryon chemical potential. The internal and free energy densities of the system at finite µ B are given by, respectively from which follow the differential relations Therefore, at fixed µ B , one finds dp(T, fixed µ B ) = sdT , (3.18) such that the speed of sound squared at any fixed value of the chemical potential reads Also, for the trace anomaly at nonzero chemical potential one obtains In order to determine the holographic equation of state at finite baryon density, we numerically generated a large number of ∼ 10 5 charged black holes on a rectangular grid We first compute the entropy density using the Bekenstein-Hawking relation (2.12) and then integrate it with respect to the temperature, while keeping the chemical potential fixed to obtain the pressure differences. One can see in Fig. 4 that the holographic result for the normalized pressure at finite chemical potential agrees quantitatively very well with the lattice data. We find also a good quantitative agreement between our holographic results and the lattice data for the speed of sound squared at nonzero baryon chemical potential while for the trace anomaly the disagreement found at large temperatures and zero density remains.
In Fig. 3, one notes from the behavior of the normalized entropy density that the change of the degrees of freedom of the system from a hadron gas to a QGP remains a smooth analytical crossover for the values of T and µ B considered here, as also seen in current lattice simulations. Therefore, as in [42], a putative holographic critical point in our model could only appear at larger values of µ B , which are beyond the values probed experimentally via the low energy collisions at RHIC.
In Fig. 4, we also note that the minimum of c 2 s and the peak of the normalized trace anomaly are pushed towards lower temperatures with increasing chemical potential, which indicates a reduction in the crossover temperature as one increases the baryon chemical potential. In Fig. 5 we defined the crossover region at nonzero baryon chemical potential in our holographic model using the inflection point of the normalized entropy density and the minimum of the speed of sound squared. One can see in Fig. 5 that the crossover temperature decreases with increasing baryon chemical potential, which is the behavior also found in other (non-holographic) calculations (see, for instance, [50,51,[74][75][76]).
We follow [50,51] and fit our numerical results for the crossover lines displayed in Fig.  5 to the series expansion Combining the results for the curvature κ of the crossover lines T c (µ B ) calculated from the inflection point of s/T 3 and the minimum of c 2 s , we obtain the holographic estimate κ ≈ 0.013, which is in remarkable agreement with very recent results obtained by different lattice collaborations: κ (I) lattice = 0.0135(20) [50] and κ (II) lattice = 0.0149(21) [51]. Furthermore, we can also give a holographic prediction for the fourth order coefficient in (3.21), which could be tested in future lattice simulations: λ ≈ 5.4 × 10 −5 .
We should also point out that for high temperatures (above 300 MeV) recent results for thermodynamical quantities obtained from three-loop Hard Thermal Loop perturbation theory are in good agreement with the lattice data for baryon chemical potentials up to 400 MeV [77]. In the near T c region, however, the weak coupling calculations start to deviate from the lattice data. Presumably, this region involves a somewhat large coupling and the holographic results obtained here provide a good quantitative description of the Within the region delimited by these curves, the degrees of freedom of the system are changing from a hadron gas to a QGP phase. One can see that the crossover temperature decreases with increasing baryon chemical potential. lattice data. This has motivated the detailed study of some non-equilibrium properties of the plasma presented in the next sections.

Energy loss of heavy and light quarks
The energy loss experienced by fast moving probes in the QGP [78,79] is of great interest to the heavy ion community (see [80] for a review and [81] for a recent summary of the results from theoretical approaches defined at weak coupling). In this section we present our predictions for the heavy and light quark energy loss in the strongly coupled QGP near the crossover transition both at zero and nonzero baryon density. While the phenomenological relevance of hard probes in low energy heavy ion collisions is only marginal, it is interesting from a theoretical point of view to understand how nonzero baryon density effects may change the energy loss of heavy and light quarks moving through the baryon dense strongly interacting plasma.
The results to be discussed in this section make use of calculations involving classical strings in asymptotically AdS 5 backgrounds. Since our backgrounds support a nontrivial Maxwell field, one could consider in principle minimally coupling the string end points to the gauge field. In the case of infinitely heavy quarks, to be discussed in the calculations of the drag force and the Langevin diffusion coefficients in Sections 4.1 and 4.2, respectively, the minimal coupling between the string and the gauge field is suppressed in the t' Hooft coupling relatively to the Nambu-Goto action, therefore, it does not contribute at leading order in the t' Hooft coupling. In the case of the jet quenching parameter to be discussed in Section 4.3, the minimal coupling does not play any role because the contributions coming from each string endpoint (located at the boundary) cancel each other out. In the case of the light quark energy loss in the shooting string setting to be discussed in Section 4.4, the minimal coupling may be relevant if one considers a flavor brane setup as the one discussed in [82], for instance. However, we postpone to a future work the generalization of the shooting string scenario with a minimal coupling to a Maxwell field.

Heavy quark energy loss
Here we compute the drag force experienced by heavy quarks in the plasma as a function of T and µ B using the holographic model constructed in the previous sections.
In order to calculate the energy loss that a heavy quark experiences in a strongly coupled medium, we will use the trailing string model of [83,84] (see also Refs. [85][86][87]). In this model, a heavy quark probe moving with a constant velocity v (say, in the spatial x direction) in the strongly coupled thermal plasma is represented holographically by an open string with an endpoint moving with the same velocity at the boundary, while the rest of the string trails behind it in the bulk with the other endpoint being located at a black hole horizon formed on top of the string worldsheet. As the quark moves through a hot and dense medium, such as the one described by our holographic model, it loses energy and momentum through drag force dp x /dt which is holographically computed using the energy flow dE/dx from the endpoint of the string at the boundary towards the worldsheet horizon.
A classical trailing string is described by the Nambu-Goto action, which in the Einstein frame has the following form where we used the 5D non-critical string theory-inspired result to relate the metrics in the string and Einstein frames [60,61] 19 . In (4.1), α = l 2 s is the square of the fundamental string length and is the induced worldsheet metric with X µ (τ, σ) describing the worldsheet embedding with the background metric g µν , which is considered here to be expressed in the standard coordinate system in Eq. (2.9). 19 We adopt the same conventions for the action (2.1) as in [42,43], which differ from the conventions in [60,61] by a simple rescaling of the scalar field,φ = 3 8 φ, whereφ is the dilaton field from [60,61] and φ is the dilaton field used here and in [42,43]. The background metrics in string and Einstein frames are then related by g We choose the static gauge (σ, τ ) = (r,t) for the worldsheet coordinates so that the worldsheet embedding for an endpoint moving in thex-direction is given by X µ (r,t) = (r,t,x(r,t), 0, 0). Following [83,84], we choose the stationary trailing string Ansatzx(r,t) = vt + ξ(r), with ξ(r) describing the string radial profile. The Nambu-Goto action with this Ansatz has the following form The equation of motion for ξ(r) implies that its conjugate radial momentum, Π ξ = ∂L NG /∂ξ , is a constant of motion. Solving this relation for ξ (r) one obtains where we chose the positive root which corresponds to the string trailing behind the endpoint. The numerator in Eq. (4.4) changes sign at somer =r given by the solution of the following equation As discussed in Refs. [88][89][90], the induced worldsheet metric has the form of a 2D black hole geometry with a horizon atr =r . The only way for ξ(r) to remain real for allr is that the denominator in Eq. (4.4) also changes sign at the worldsheet horizonr . This requirement fixes the value of the conjugate momentum The heavy quark drag force is given by the flux of the momentum from the endpoint into the bulk of the string, which turns out to be given precisely by Π ξ [83,84] Plugging in the Ansatz (2.9) for the background metric into Eq. (4.7) and expressing everything in terms of the numerical coordinates (as discussed in Section 2.3), one finally obtains the following dimensionless ratio (which is negative-definite) where we defined the 't Hooft coupling λ t = (α ) −2 . In the numerical coordinates, the worldsheet horizon defined by Eq. (4.5) is obtained by numerically solving the following equation We show in Fig. 6 our numerical results for the heavy quark drag force normalized by its conformal limit (which is only attained at high temperatures) as a function of T for several values of µ B . The conformal result at µ B = 0 is obtained by using the AdS 5 -Schwarzschild background, which gives the analytical result F conformal [83,84]. In the left panel we set v = 0.6 while in the right panel we take v = 0.99. One can see that the value of the normalized drag force is very sensitive to the quark velocity and a nontrivial µ B dependence and a non-monotonic T dependence only appear at high velocities. This is because, for low velocities, r from (4.9) is close to the background horizon where Φ(r) is small and the drag force is less sensitive to the finite density effects 20 . At v = 0.99, one can see that the normalized drag force peaks in the phase transition region T ∼ 150 − 200 MeV. These results illustrate that the heavy quark energy loss is sensitive to the details of the phase transition and this should be taken into account in phenomenological strongly coupled models used in the calculation of the heavy quark nuclear modification factor [91,92]. Moreover, the fact that a nontrivial baryon chemical potential dependence and a non-monotonic temperature dependence are only observed at large velocities suggests that charm quarks should be generally more sensitive to the medium properties in comparison to bottom quarks. Furthermore, one can see in Fig. 6 that the peak produced at high velocities increases significantly when µ B is raised to 400 MeV and so does the overall energy loss. It is also interesting to note that the peak in the energy loss is pushed towards lower temperatures with increasing chemical potentials, which is in agreement with the fact that in our model the crossover temperature decreases with increasing µ B .

Heavy quark Langevin diffusion coefficients
The holographic prescription for treating Langevin processes was studied in Refs. [88,89] in the context of a strongly coupled N = 4 Supersymmetric Yang-Mills (SYM) plasma and the generalization to non-conformal holographic duals was put forward in Ref. [90]. Here we use the general formulas derived in [90] to compute the Langevin diffusion coefficients for heavy quarks in the zero-frequency limit, which corresponds to the long time behavior of the stochastic diffusion processes in the hot and dense medium described by our holographic model.
In the trailing string model the effects of thermal fluctuations on the heavy quark propagation are not taken into account. Following Refs. [88][89][90], we consider the Brownian motion of the heavy quark as it moves through the medium, which may be approximately modeled by a linearized local Langevin equation describing the fluctuations of the heavy quark trajectory with a constant velocity (see, e.g., [90] for a derivation and discussion of the approximations). Following the standard holographic dictionary, from the minimal coupling interaction for an external heavy quark, dtδX µ (r → ∞, t) F µ (t), one identifies the boundary value of the fluctuation of the trailing string worldsheet embedding, δX µ (r → ∞, t), as the source for a boundary gauge invariant operator, F µ (t), which plays the role of a random force acting on the heavy quark. The momentum dependent Langevin diffusion coefficients can be calculated from the symmetrized real time 2-point correlation function of this boundary force operator and the Langevin diffusion constants are then obtained from the zero frequency limit of these coefficients. The standard way [93][94][95][96] to obtain these correlation functions holographically is to solve the linearized bulk equations of motion for the trailing string fluctuations δX µ (r, t) and compute the on-shell action. However, in the zero momentum limit of these correlation functions, an equivalent and simpler approach involves the membrane paradigm [97] since in this case it is not necessary to explicitly solve the equations of motion. This was used in [90] to compute the parallel and perpendicular Langevin diffusion constants for a general 5-dimensional nonconformal holographic background, where the Einstein-frame background metric is written in the conformal gauge (with boundary at u = 0) In (4.10), T is the Hawking temperature of the worldsheet horizon defined by u computed via f (u ) = v 2 , and b 2 (S) (u) = e √ 2/3φ(u) b 2 (u) is the warping factor in the string frame 21 .
In order to relate our background (2.9) in theB(r) = 0 gauge to the background (4.11) in the conformal gauge, we equate (2.9) and (4.11) and require thatφ(u) =φ(r) and (t, x) = (t, x), which gives We can now use this to write (4.10) in our standard coordinates and then express those equations in the numerical coordinates by following the discussion in Section 2.3, (4.15) Our numerical results for the perpendicular (upper panel) and parallel (lower panel) Langevin diffusion coefficients normalized by their respective conformal limits are shown in Fig. 7. The conformal values for these coefficients are κ conformal ⊥ /(T 3 √ λ t ) = π γ(v) and κ conformal /(T 3 √ λ t ) = πγ 5/2 (v) [88,89]. In the left column we set v = 0.6 while on the right we take v = 0.99. Similarly to the curves for the heavy quark drag force in Fig.  6, one can see that these transport coefficients display a more appreciable µ B dependence and a non-monotonic T dependence when the heavy quark velocity is close to unity. This non-monotonic dependence with the temperature at µ B = 0 may play an important role in phenomenological setups involving strongly coupled modeling of thermal fluctuations in heavy quark energy loss. For instance, it would be interesting to see how an improved treatment of the strongly coupled medium properties near the phase transition affects the study performed in Ref. [98], which employed the conformal limit of the Langevin diffusion coefficients. Moreover, we find that these coefficients increase with µ B , which indicates that the thermal fluctuations of the trailing string become more relevant in a baryon dense medium.

Light quark energy loss: the jet quenching parameter
Even though weak coupling approaches to jet quenching phenomena seem to be consistent with heavy ion data [81], one should keep in mind that the weak coupling description of the medium used in these approaches does not incorporate the defining property of the     QGP, i.e., its nearly perfect fluidity. Energy loss in a hot and dense medium involves several different scales [80] and, as discussed in [99,100], strong coupling approaches to the interaction between the jet and the medium at scales of order ∼ T may be called for. In this section we compute the jet quenching parameterq [99,101] (see also [87,100]) associated with light quark probes in our nonconformal backgrounds with nonzero baryon chemical potential. For recent lattice results at zero baryon chemical potential regarding this quantity see [102].
The original formula forq in holography, derived in [103], was inspired by the dipole approximation [104] for the perturbative QCD jet quenching parameterq [105] where L and L − (with L L − ) are the sides of an adjoint light-like rectangular Wilson loop described by a quark-antiquark pair. The authors of [99,101] used this to define a non-perturbative version ofq that would be valid at strong couplinĝ (4.17) According to the standard holographic dictionary [106] (see also Refs. [107][108][109]) the expectation value of a Wilson loop in the fundamental representation in a strongly coupled SU (N c ) gauge theory in the large N c limit is given by the exponential of minus the on-shell Nambu-Goto action for a string worldsheet whose boundary coincides with that Wilson loop. Hence, we have from (4.17)q where the extra factor of 2 comes from the fact that at large N c the adjoint Wilson loop is twice the value of the fundamental loop.
Since we are interested in a light-like Wilson loop, we define light-cone coordinates leavingr,ỹ, andz unchanged. In these coordinates, the background metric (2.9) has the following form The geometric configuration considered here corresponds to a string attached to a quarkantiquark pair at the boundary, separated by a distance L and moving through a much longer (light-like) distance L − , tracing out a rectangle at the boundary. The string sags into the bulk describing an U-shaped string worldsheet whose tip turns out to lie at the radial position of the bulk horizon, as we will see in a moment. It is convenient to use the static gauge (τ, σ) = (x − ,ỹ) such that the worldsheet embedding reads X µ (x − ,ỹ) = (r(x − ,ỹ), 0,x − ,ỹ, 0). Since we are taking the L − L limit, there is translational invariance in thex − -direction and the worldsheet embedding reduces to X µ (x − ,ỹ) = (r(ỹ), 0,x − ,ỹ, 0). (4.21) The Nambu-Goto action for this string configuration is given by Since the integrand in Eq. (4.22) does not depend explicitly onỹ, we can identify the following constant of motion Solving Eq. (4.23) forr (ỹ), we obtain the following on-shell relation Note that for sufficiently small C 0 the term in the brackets in (4.24) is always positive andr (ỹ) vanishes only whenh(r) vanishes, which happens at the bulk horizonr =r H and it corresponds to the tip of the string worldsheet. Since we place the origin of thẽ y-axis at the center of the L-side of the rectangular Wilson loop, the U-shaped profile implies the symmetryr(ỹ) =r(−ỹ), and, consequently, the turning point must correspond tor H =r(ỹ = 0). Then, after integrating Eq. (4.24) we have where, formally,r max → ∞. For small values of C 0 one may expand the integrand of (4. 25) in powers of C 0 to obtain at leading order .

(4.26)
This shows that a small value of C 0 corresponds to a small L (when compared to L − ).
Plugging the solution (4.24) in Eq. (4.22) we obtain the on-shell Nambu-Goto action where we used (4.25) andr(ỹ) =r(−ỹ) symmetry. Following [101], we now compute the disconnected contribution for this correlator corresponding to the action of two straight strings in the bulk. In order to compute these disconnected string configurations, it is convenient to define another static gauge for the worldsheet parameters, namely, (τ, σ) = (x − ,r), such that the corresponding worldsheet embeddings are of the form X µ (x − ,r) = (r, 0,x − , 0, 0). The sum of the actions for these two disconnected string configurations is then given by . .

(4.29)
Finally, plugging (4.29) into (4.18) and going back to the numerical coordinates, as discussed in Section 2.3, we arrive at the formula for the jet quenching parameter in our backgroundq . (4.30) Our numerical results for theq parameter (4.30) normalized by its conformal limit, q CFT / √ λ t T 3 = π 3/2 Γ(3/4)/Γ(5/4) [99,101], are shown in Fig. 8. The jet quenching parameter displays a peak in the crossover region and its overall value increases with increasing baryon chemical potential. As observed in Ref. [110], the behavior of the curves for the jet quenching parameter are qualitatively similar to the ones found for the normalized trace anomaly (Fig. 4) and, consequently, the peak of this observable may be also employed to characterize the crossover temperature. The phenomenological relevance of a peak in the energy loss near the phase transition was pointed out in [111] and, more recently, in [112].
Our results show that the energy loss of light probes in this plasma is not only sensitive to µ B but also to the rapid change in degrees of freedom that occurs at the phase transition. The qualitative temperature dependence ofq for µ B = 0 found here is consistent with the analysis performed in [81], though the overall magnitude of our result is larger than the values extracted via weak coupling calculations in [81]. Also, note that Fig. 8 is consistent with our calculations for the energy loss of rapidly moving heavy quarks presented in the previous sections. More importantly, even though this strongly coupled medium cannot be understood in terms of quasiparticles,q is overall enhanced in a baryon dense system. Let us make a brief comparison between our result illustrated in Fig. 8 at µ B = 0 (solid curve) and the lattice result recently obtained in Ref. [102]. There, by taking the gauge coupling g 2 QCD ≈ 2.6 at T ≈ 398 MeV, it was estimated thatq is around 6 GeV 2 /fm for QCD with two light flavors at this temperature (we note that this value is larger than the phenomenologically extracted result of [81]). From our result computed at µ B = 0 one obtainsq ≈ 9.7 √ λ t T 3 at T ≈ 398 MeV and, taking λ t = g 2 QCD N c ≈ 2.6 × 3 = 7.8, one finds the holographic estimateq(T ≈ 398 MeV, µ B = 0) ≈ 8.6 GeV 2 /fm. This is quite close to the lattice result even though our holographic model was constructed to mimic the thermodynamics of (2 + 1)-flavor QCD instead of QCD with only two light flavors considered in Ref. [102].
Finally, let us note that it is also possible (as discussed, for instance, in [90]) to define a Langevin jet quenching parameterq ⊥ (associated with the transverse momentum broadening of heavy quarks in the Langevin processes) from the perpendicular Langevin diffusion coefficient calculated in Section 4.2:q ⊥ = 2κ ⊥ /v. This parameter is different than the jet quenching parameterq [101] since the latter is strictly valid for light-like trajectories (see, for instance, the discussion in [89]). However, as one can see in Fig. 7, this Langevinq ⊥ (which is trivially related to κ ⊥ ) at high velocities displays a qualitative behavior similar to the one seen in the jet quenching parameterq in Fig. 8.

Light quark energy loss: the shooting string
In this section we will use the holographic model of "shooting strings" proposed in [113] to calculate the light quark energy loss in our nonconformal backgrounds. The shooting string is the simplest phenomenological implementation of the finite endpoint momentum approach [114] to express the energy loss of light quarks in a way that can be used to calculate observables such as the nuclear modification factor R AA at RHIC and LHC 22 . As shown in [114], adding finite momentum at the endpoints of classical strings allows them to travel further in the AdS 5 -Schwarzschild background than the previous falling string configurations [118,119]. At the same time, it was argued that strings with finite momentum at their endpoints may provide a more natural description of energetic light quarks plowing through a hot, strongly coupled plasma. These results were reassuring as they had the potential to reconcile the apparent quantitative inconsistencies of the first attempts at holographic light quark energy loss of [120,121] and the LHC suppression data for light hadrons [122].
When finite momentum is added to the endpoint, it turns out that its motion simplifies [114]: the endpoint must follow null geodesics, while the rest of the string sags behind it. In this approach, the quark energy loss is naturally described as the flow of energy from the endpoint into the bulk of the string, as the latter represents the color field generated 22 Other holographic models of light quark energy loss include the holographic version of the brick problem of [115] and its hybrid phenomenological implementation in [116], as well as the recent work of [117] where different holographic definitions of light quark jets were explored. by the quark. This quantity has a simple analytical expression [113]: where we assumed that the endpoint is moving in thex-direction, and that the string frame metricG (S) µν is diagonal in (t, x,r) coordinates and depends only on the radial coordinatẽ r, as is the case in our background (2.9). The constant R parametrizes the geodesic that the endpoint is following, and is equal to the ratio of energy and momentum that a free particle following this geodesic would have conserved.
In order to cast (4.31) into a phenomenologically usable form, we need to relate R to some observable quantity and expressr in terms ofx, which involves inverting the geodesic trajectoryx geo (r). A simple way to do this was proposed in [113] where it was shown that, for phenomenological reasons, it makes sense to consider the endpoints that are initially close to the horizon. In this case, the endpoint is "shot" towards the boundary, and, for moderately large maximumr that the endpoint can reach, the energy loss can be well approximated by using the critical geodesic with R = 1. This is called the shooting string limit and for an N = 4 SYM plasma (corresponding to the AdS 5 -Schwarzschild background) it leads to the following expression for the energy loss To obtain the analog of (4.32) in our numerical backgrounds, we first need to find the critical geodesics numerically. To do that, we will parametrize the geodesics with timet and seek to obtainr geo (t) andx geo (t). Along any geodesic in background (2.9), the energy and the momentum (of a test particle) are conserved: where η(t) is an auxiliary field (the "metric" on the worldline). Defining R ≡Ẽ geo /p x,geo , we can solve for dx geo /dt dx geo dt =h (r geo (t)) R . (4.35) Plugging this in the null condition ds 2 = 0 allows us to solve for dr geo /dt: We can now in principle solve for null geodesics by dividing (4.35) and (4.36) and integrating the resulting dx geo (r)/dr. However, close to the maximumr that the geodesic reaches (which we callr * and is given by R 2 =h(r * )), the quantity dx geo /dr diverges and the numerical procedure is unstable. This is precisely the reason why we chose to parametrize the geodesics with coordinatet. The first order equation (4.36) still vanishes atr =r * , resulting in "stiffness" of the numerical system: in order to avoid this, we will use the second-order equation forr geo (t), which is just the geodesic equation (with an additional term that corrects for the fact thatt is not an affine parameter).
Finally, in order to take the shooting limit, we choose as larger * as the numerical code allows, solve for the null geodesic trajectory with initial conditionsr geo (0) =r H and x geo (0) = 0, plug the solutionr geo (t) in (4.31), which explicitly looks like and check that changingr * around that value does not affect dẼ/dx within some phenomenologically relevant range ofx. The results for a range of different chemical potentials and temperatures are shown in Fig. 9.   One can see in Fig. 9 that the greatest deviations from the conformal limit occur at early times (smallx), which is a novel feature of the shooting string setup: at small x the endpoint is close to the horizon and hence sees the strongest effects of conformal symmetry breaking. We also see that, at early times, the energy loss in this model is larger than its conformal limit, while at later times, it becomes smaller and asymptotes to the same functionalx-dependence as the conformal one (which is, again, expected, since in the shooting limit, largex means that the endpoint is closer to the boundary). Finally, the higher the temperature, the closer the ratios in Fig. 9 are to unity, correctly recovering the conformal limit. Furthermore, the energy loss increases with increasing µ B , which is consistent with expectations based on simple density arguments (though this strongly coupled dense medium does not admit a quasiparticle description).
Another interesting observation about Fig. 9 is that, in the phenomenologically relevant range ofx ∼ 0 − 1 fm, the shooting string energy loss is larger relatively to its conformal limit at lower temperatures than it is at higher temperatures. This has potentially interesting phenomenological implications for the problematic simultaneous match of the nuclear suppression factors R AA of light hadrons at RHIC and LHC [123]. As shown in [113], if one tries to fix the parameters of a conformal model by matching, say, the R AA data at the LHC and then use the same set of parameters at RHIC, the RHIC R AA ends up being above the data indicating that the predicted energy loss is simply too small. However, in this more realistic nonconformal model that can assess the rapid change in degrees of freedom in the phase transition, we see that the energy loss at RHIC (which is, on average, colder than LHC) is higher than the conformal expectation and, thus, this will result in a better match with the data and perhaps even completely close the "gap" between RHIC and LHC. We hope to come back to this issue in the near future.
It is also interesting to consider energy loss as a function of temperature at some fixed distancex, as shown in Fig. 10. We see that, at early times (left plot), the light quark energy loss (normalized by its conformal limit) approaches unity from above and looks in fact more similar to the heavy quark drag force at low velocities (Fig. 6) than toq (Fig.  8). At late times (right plot), this behavior changes qualitatively and this ratio approaches unity from below, developing features qualitatively similar toq. This nontrivial behavior is quite specific to the shooting string setup and it would be interesting to see how this is translated into observables such as R AA and especially the elliptic flow parameter v 2 , which is very sensitive to the specific path dependence of energy loss.
One should also note how, at early times, the energy loss is almost insensitive to the baryon chemical potential while at later times it becomes more sensitive to it. This is again a reflection of the shooting string setup in which, at early times, the endpoint is close to the horizon whereΦ(r) = 0. A similar behavior has been observed before in the drag force (Fig. 6) where at low velocities one is more sensitive to the metric close to the horizon.   In summary, our holographic model allows for the calculation of quantities that characterize the energy loss of heavy quarks (see Figs. 6 and 7) and light quarks (see Figs. 8, 9 and 10) in a hot and baryon dense strongly interacting plasma that not only behaves as a nearly perfect fluid but also displays thermodynamic properties that are in quantitative agreement with current lattice QCD calculations in the crossover region. Energy loss is found to overall increase with µ B while also displaying a nontrivial, fast-varying behavior as a function of the temperature near the crossover. These results can be readily implemented in phenomenological studies of the nuclear modification factor in the QGP such as, for instance, [124].

Equilibration time of a baryon rich strongly coupled plasma
In this section we will explore some generic aspects concerning the equilibration properties of our strongly coupled plasma at finite temperature and baryon chemical potential by computing the lowest quasinormal mode (QNM) associated with an external massless scalar perturbation in the bulk. The main idea behind this is that non-hydrodynamical degrees of freedom 23 play an important role in the transition from a non-equilibrium state to an equilibrated plasma 24 and these modes in a strongly coupled gauge theory correspond, via the holographic duality, to the quasinormal modes of a black hole background in the bulk [129,130]. The motivation for studying this subject comes from the interest in the thermalization properties of the strongly coupled quark-gluon plasma produced in heavy ion collisions. Here we will focus on the lowest quasinormal modes as these represent the dominant, longest-lived non-hydrodynamical degrees of freedom that can be used to provide an upper bound on the equilibration time in the dual plasma.
Quasinormal modes for nonconformal holographic duals have been recently investigated in [131] and [132,133], motivated by the potential effects that the breaking of conformal invariance might have on the equilibration of strongly coupled plasmas that are qualitatively more similar to the actual quark-gluon plasma formed in ultrarelativistic heavy ion collisions. In [131] the authors looked at a class of bottom-up nonconformal Einstein-Dilaton models with the dilaton potential tuned to reproduce the QCD equation of state at zero baryon chemical potential 25 , while in [132] and [133] the authors studied the rate of equilibration in two nonconformal top-down models, the N = 2 * plasma [134] and the Klebanov-Strassler model [135], respectively. These studies have found that even when the deviation from conformal behavior was large, as it happens for instance around the crossover phase transition, the imaginary parts of the lowest quasinormal mode frequencies 23 A non-hydrodynamical mode is defined by the following condition for its dispersion relation: ω(k = 0) = 0. 24 This can be easily understood in kinetic theory where the timescale associated with the lowest quasinormal mode enters in the effective hydrodynamic description via the relaxation time coefficient associated with the shear stress tensor [125]. This detail is important in studies of the applicability of different hydrodynamic theories when compared to exact kinetic theory solutions, such as in [126][127][128]. 25 One of the dilaton potentials used in [131] was the same we used here in (3.9) and, previously, in [38].
differed by roughly a factor of two from their conformal limit. Although this does not represent a qualitative change, a factor of two difference in the exponential damping of these degrees of freedom close to the phase transition may be of phenomenological importance.
In fact, as originally pointed out in [125,136], the fact that the non-hydrodynamical modes in spatially isotropic and uniform strongly coupled models possess a nonzero real part [137] implies that the way these systems approach their universal hydrodynamic limit is radically different than what is observed in the case of dilute gases described by the relativistic Boltzmann equation. In the latter, spatially isotropic and uniform dilute gases have purely imaginary non-hydrodynamic modes, which lead to a characteristic exponential decay of these excitations. However, the nonzero real part found for these modes at strong coupling leads to an oscillatory component for the dissipative currents (such as the shear stress tensor) that seems to be a novel feature revealed by the holographic correspondence. This property has not yet been incorporated in the hydrodynamical modeling of the QGP [9] though it is clearly seen in real time, holographic calculations [138].
Here we investigate the effect of a nonzero baryon chemical potential on the lowest QNM of an external massless scalar field in the limit of zero spatial momentum, k → 0, as in [131][132][133]. In this limit, as discussed in [131], all the relevant information about the QNM's of the model may be extracted from the study of an external massless scalar perturbation. We will also check how the QNM frequency varies with k and see that it shows only a very mild dependence up to k ≈ 2πT , in accordance with the so called "ultralocality" property first pointed out in [138] for a conformal plasma and later found to hold also in nonconformal plasmas at zero chemical potential [131].
We start by writing our background (2.2) in the in-falling Eddington-Finkelstein coordinates since in those coordinates the condition that the mode is ingoing at the horizon translates into simple regularity there. We define the Eddington-Finkelstein time v, as follows dv = dt + e −A dr h . (5.1) With this transformation, the metric becomes The equation of motion of an external, minimally coupled massless scalar field 26 in this geometry is just the Klein-Gordon equation, and we are considering the (spatially isotropic) harmonic solution in the boundary directions with a nontrivial dependence on the radial coordinate: Ψ = e −iωv+ikx ψ(r). With this Ansatz, the equation of motion has the following form This equation needs to be solved with the condition that ψ is regular at the horizon while it obeys a Dirichlet condition at the boundary, which means that this equation for a given background will admit solutions only for discretely many complex frequencies ω. We searched for these frequencies numerically using a simple shooting method where one solves (5.4) on a grid of complex ω, specifying some initial value for ψ at the horizon and its derivative as given by the near horizon expansion of (5.4): where we set k = 0. We note that here we are working with the numerical coordinates and we used the near-horizon expansions of the metric from Section 2.2 (in particular, A 1 is given by (2.16)). Once the numerical solution is obtained, we look for the QNM frequencies by demanding that the value of ψ at some large r near the boundary vanishes 27 .  Our results for the lowest QNM mode ω 0 with vanishing momentum are shown in Fig.  11, where one can see how the temperature dependence of the real and imaginary parts of ω 0 change as we vary the baryon chemical potential. Our result for the imaginary part of ω 0 at zero chemical potential is very similar to the one found in [131] as it exhibits a factor of two decrease (in absolute value) close to the phase transition while returning to the conformal limit of

T [MeV]
= −1.373 at high temperatures. This indicates that, generically, the typical equilibration times close to the phase transition may in fact be a factor of two larger than what one would expect from the conformal results. In the same plot we can also see the effects of the finite chemical potential on Im[ω 0 ], which show that, by increasing µ B , one generally gets smaller equilibration times. The real part of ω 0 , shown in the right plot, is showing very little (relative) deviation from the conformal limit 27 In practice, one could start with some rough grid in complex ω centered around the conformal value and then find clusters of small values of |ψ(r → ∞)|. These points can then be used as centers for a finer grid and one can keep repeating this until the required accuracy is reached. as one varies both the temperature and the baryon chemical potential. If a strong coupling description based on holography is appropriate for the QGP formed in low energy heavy ion collisions at RHIC, our results show that this plasma should still be a nearly perfect liquid 28 in which the equilibration process towards hydrodynamic behavior is similar to that found at zero baryon density [138], which is not included in the current numerical hydrodynamic codes. In Fig. 12 we show the k-dependence of the imaginary and real parts of the frequency of the lowest QNM for four different backgrounds that span our phase space with temperatures between 140 and 400 MeV and baryon chemical potentials between 0 and 400 MeV. We see that the k-dependence of the imaginary parts is quite mild while the real parts quickly attain the conformal linear dispersion.

Concluding remarks and perspectives
In this work, we used a 5-dimensional Einstein-Maxwell-Dilaton bottom-up holographic model [42,43] that provides a comprehensive framework to investigate both the equilibrium and out-of-equilibrium properties of a baryon rich strongly coupled QGP. The parameters of the model were determined using current lattice data for (2 + 1)-flavor QCD thermodynamics at zero baryon chemical potential. Once all the parameters were fixed by µ B = 0 lattice data, we proceeded to compare the thermodynamics of our holographic model with µ B = 0 lattice results obtaining a good quantitative agreement for the pressure and the speed of sound for baryon chemical potentials up to 400 MeV. We also determined the crossover region in our holographic model for this range of chemical potentials and we found that the crossover temperature is lowered as one increases the value of the baryon chemical potential. The curvature of the crossover band in our model is in remarkable agreement with recent lattice results. Therefore, this holographic model was shown to adequately describe many of the equilibrium properties of the hot and dense QGP near the crossover transition, besides also displaying the nearly perfect fluidity property inherent to holographic approaches [11,12].
This model was used to investigate two important aspects involving non-equilibrium processes that can take place in the QGP: the energy loss of light and heavy quarks and the equilibration time associated with non-hydrodynamic modes. We found that the energy loss of both light and heavy flavors near the phase transition is enhanced once finite baryon density effects are taken into account. This justifies the expectation that highly energetic probes should lose more energy in denser media at a given temperature even in the absence of a quasiparticle description. We also found that the inclusion of nonzero baryon density effects in the nonconformal plasma does not change the qualitative behavior displayed by the QNM's of a massless scalar excitation: the spectrum of QNM's still possesses both real and imaginary parts at zero spatial momentum and their values are reduced by a factor of nearly two at the crossover transition [131], though we found that the equilibration time associated with the lowest lying QNM becomes shorter in a dense medium.
The holographic model used in the present work may be also employed to obtain predictions for many other physical observables that are relevant for the study of the QGP formed in heavy ion collisions. Some projects we hope to pursue in the near future include the generalization of the recent calculations performed in Ref. [38] involving the second order hydrodynamic transport coefficients by taking into account a nonzero µ B , which would entail in the calculation of new transport coefficients associated with baryon diffusion in the dense plasma. In this regard, by using the very same Einstein-Maxwell-Dilaton holographic model constructed in the present work, some of us have recently obtained in Ref. [140] holographic predictions for the T and µ B dependence of some of the transport coefficients associated with the baryon charge sector, namely, the baryon and thermal conductivities, and also the baryon diffusion constant. In Ref. [141], the same model was used to obtain predictions for the electric charge sector of the hot and baryon dense QGP, with calculations of the T and µ B dependence of the electric conductivity and electric diffusion constant, and also the thermal photon and dilepton production rates. It would also be interesting to see how the holographic Polyakov loop [142,143] and the imaginary part of the rectangular, time-like Wilson loop [144,145]

A Some technical aspects of the model
In this Appendix we discuss some technical aspects of our holographic model concerning mainly the forms of V (φ) and f (φ) fixed in Eqs. (3.9) and (3.12), respectively, and their relations with the thermodynamic stability of our black hole solutions.
In the present work, as done originally in [42,43], we restricted our calculations to positive values of the initial condition φ 0 corresponding to the value of the dilaton field at the horizon. In this case, we checked that all the numerical solutions we produced have φ(r) ≥ 0. We show in Fig. 13 the plots for the dilaton potential V (φ) and the Maxwell-Dilaton coupling f (φ) fixed in Eqs. (3.9) and (3.12).
By following [146,147], we identify the fixed points of our model: • Ultraviolet fixed points: these are given by the maxima of the potential at φ = 0 and φ ≈ ±10.72; • Standard infrared fixed points: these are given by the local minima of the potential at φ ≈ ±7.81; • Hyperscaling violating infrared fixed points: these are given by φ → ±∞; • AdS 2 infrared fixed points: these are given by the vanishing of V /V + f /f and are located at φ ≈ 0.06, 0.70, 10.98, 13.45.
We remark that the ultraviolet fixed point reached by our solutions corresponds to AdS 5 geometries with a vanishing dilaton field at the boundary.
Note that we may replace the non-monotonic part of our potential by a monotonic piece smoothly connecting the deep infrared at large φ to the part of the potential we effectively used in our present calculations and the results obtained in the region of the phase diagram explored in our work, 100 MeV ≤ T ≤ 400 MeV and 0 ≤ µ B ≤ 400 MeV would not change since the replaced piece of the potential has been always hidden behind the horizon in our current calculations. Such replacement, however, would modify the results in other regions of the phase diagram, such as the physics at very low T , and it would also modify (part of) the structure of fixed points of the model. However, as discussed before, the infrared region of our potential, corresponding to large values of φ or low values of T , was not constrained by lattice data and, therefore, we have no interest in this region in the present work. In fact, as discussed in [42], by adequately adjusting the large φ piece of V (φ), one could in principle have a variety of different behaviors at low T .
We also note that the non-monotonicity of our dilaton potential poses some technical complications if one intends to explore larger regions of the plane of initial conditions (φ 0 , Φ 1 ) and the corresponding generated region in the (T, µ B ) phase diagram. This is because each local extrema of the dilaton potential corresponds to a singular point of the Einstein-Maxwell-Dilaton equations of motion. The numerical integration of these equations of motion in cases when extra singular points between the boundary and the horizon are present is complicated and we have not considered such cases in the present work 29 . This, in turn, limits the range of initial conditions one considers in practice, because if one wants to avoid the complicated task of numerically integrating the equations of motion with extra singular points between the boundary and the horizon, then the local extrema of V (φ) must be always hidden behind the horizon (as we have done in all of our current calculations). In fact, if we take φ 0 above 7.81, then we stumble in an extra singular point for any value of Φ 1 . The same happens for φ 0 < 7.81 if we take high enough values of Φ 1 (even if it satisfies the physical bound Φ 1 < Φ max 1 (φ 0 )) and the dilaton field φ becomes a non-monotonic function of r eventually probing the region where φ(r) > 7.81. We have excluded all such initial conditions from our calculations by restricting the ranges of φ 0 and Φ 1 explored in the present work. However, we remark that the range of initial 29 The singular points corresponding to the boundary and the horizon have been properly dealt with in our numerical calculations, as discussed in Section 2.2.
conditions we employed is enough to cover the relevant region of the phase diagram which we are interested in, 100 MeV ≤ T ≤ 400 MeV and 0 ≤ µ B ≤ 400 MeV, and in fact, a fairly broader region than that.
Next, we give evidence that the black hole solutions we generated in the region of interest of the phase diagram are unique. In order to do this, we examine the mapping between the initial conditions and the relevant thermodynamic state variables, i.e., the mapping (φ 0 , Φ 1 ) → (T, µ B ), exploring the (extended) range of initial conditions comprised by φ 0 ∈ [0.05, 7.3] and Φ 1 ∈ [0, 0.5 Φ max 1 (φ 0 )]. These solutions are displayed in Fig. 14.
As one can see in the bottom plot in Fig. 14, for the range of initial conditions analyzed here, there are no competing black hole solutions in the region of the phase diagram which we are interested in our work since no two solutions with different (φ 0 , Φ 1 ) give the same (T, µ B ) inside this region (note that each point is only crossed by one blue line and one golden line). Moreover, as one can see from the top plots in Fig. 14, extending the region covered in the plane of initial conditions just add points in the (T, µ B ) plane which are more and more distant from the phenomenologically relevant region studied in our work (for very high values of T and µ B , which were not explored in the present work, one notes that there are indeed competing black hole solutions). This analysis constitutes evidence that in the region of the phase diagram we worked with in our manuscript the black hole solutions are, in fact, unique.
In spite of the above argument, let us now assume the unlikely possibility that, if we consider some initial conditions outside the region of the plane of initial conditions we covered, then a hypothetical competing branch of black hole solutions could fall within the region of the phase diagram which we are interested in the present work. This could only be the case for solutions where φ(r) probes values above 7.81, since in our scanning we have excluded this possibility for φ(r) < 7.81. In this case, since the region where φ(r) probes values above 7.81 was always hidden behind the horizon in our calculations, then we could in principle just modify the form of the potential in that region such as to remove any possible competing branch of black hole solutions hypothetically falling inside the relevant region of the phase diagram investigated in the present work, while maintaining the solutions we used in our calculations completely unmodified. Consequently, all of our current results would remain exactly the same (which is indeed desirable, given the quantitative agreement between our holographic thermodynamics and corresponding lattice data even at nonzero µ B ).
The last question we address in this Appendix regards the stability of the black hole solutions used in our calculations. Since we argued before that it is unlikely that there are competing branches of black hole solutions in the relevant region of the phase diagram covered in the present work 30 , one still has to show that the black hole solutions we derived have larger pressure than the thermal gas solution, which is just the vacuum solution 30 And that, even if a hypothetical extra competing branch were indeed present, we could in principle simply remove it by deforming the large φ part of V (φ), which has been always hidden behind the horizon in the present work without modifying our current results. with the imaginary Euclidean time direction compactified over a circle with circumference β = 1/T . In the present model, it is numerically difficult to obtain very low temperature geometries, and in particular the vacuum solution, because one would need to integrate the equations of motion with the extra singular points corresponding to the local minimum and the local maximum of V (φ). Here we present an argument of plausibility for why the black hole solutions we obtained should be expected to be thermodynamically preferred, i.e. to possess a larger pressure than the thermal gas solution. We divide our argument in two parts: 1. First, we calculated the discrete version of the Jacobian J = ∂(s, ρ)/∂(T, µ B ) defined in [42] for all the black hole solutions generated within the region of interest of the phase diagram and checked that J > 0 for all of them, which shows that there are no thermodynamically unstable black hole solutions inside the region of the phase diagram investigated in our work; 2. Second, since there are no unstable black hole solutions, these solutions could be at least metastable, if they had smaller pressure than the thermal gas solution. But this is not plausible, because we can go as low as T ≈ 0.005 MeV at µ B = 0 (this geometry is generated by the initial conditions (φ 0 , Φ 1 ) = (7.8, 0)) and check that the pressure of the black hole solutions monotonically increases with increasing T . Since the pressure of the thermal gas solution must be just a constant, which is usually matched with the black hole pressure in the extremal limit of vanishing horizon [148], it cannot be larger than the pressures we calculated using black hole solutions in the region of the phase diagram explored in our work, where T 100 MeV.
Therefore, the black hole solutions we derived should not only be thermodynamically stable but they should also be thermodynamically preferred over the thermal gas solution in the entire region of the phase diagram studied in our work. Furthermore, at nonzero µ B this conclusion also holds since, as we have shown before, the pressure of the black hole solutions increases with increasing µ B (what is quantitatively confirmed by current lattice data up to µ B = 400 MeV), while the thermal gas solution does not depend on µ B . Then, the numerical black hole solutions we used to compute the physical observables in the present paper are indeed expected to be the dominant saddle points of the model.
In the future, it would be interesting to use a purely monotonic dilaton potential that can describe the lattice data for the QCD equation of state at µ B = 0 since in this case the technical complications we commented upon in this Appendix would not be present. An example of such a potential can be found in Ref. [149].