Holographic Bjorken flow of a hot and dense fluid in the vicinity of a critical point

In this paper we use the gauge/gravity duality to perform the first systematic study of the onset of hydrodynamic behavior in a hot and dense far-from-equilibrium strongly coupled relativistic fluid with a critical point. By employing a top-down holographic construction that stems from string theory, we numerically obtain the full nonlinear evolution of the far-from-equilibrium system undergoing a Bjorken expansion and address the following question: how does hydrodynamic behavior emerge in the vicinity of a critical point in the phase diagram? For the top-down holographic system analyzed in the present work, we find that the approach to hydrodynamics is strongly affected by the presence of the critical point: the closer the ratio between the chemical potential and the temperature is to its critical value, the longer it takes for the system to be well described by the equations of viscous hydrodynamics.


I. INTRODUCTION
Quantum chromodynamics (QCD) [1,2] is the fundamental theory of nature that describes the interactions between quarks and gluons, accounting for the vast majority of the mass of ordinary baryonic matter in the Universe [3]. QCD is also responsible for gluing protons and neutrons together inside atoms allowing for the stability of nuclear matter and a long-standing challenge in the field has been the determination of its equilibrium properties as a function of temperature T and baryon chemical potential µ B [4].
The QCD phase diagram has been the focus of many experimental, observational, and theoretical efforts through the years. In the experimental front, the main tool used to investigate hot and dense QCD matter is heavy ion collisions [5][6][7][8][9], which produce a strongly coupled fluid made of deconfined but highly correlated quarks and gluons, the so-called quark-gluon plasma (QGP) [10][11][12]. In the observational front, the core of compact neutron stars may comprise new phases of quark matter expected to exist when µ B /T 1, such as color superconductivity [13]. While the equation of state of such high baryon dense matter is currently unknown, powerful constraints may be derived at low and asymptotically high densities [14,15] which together with observational data on neutron stars [16,17], including the very recent observation of gravitational waves emitted by a neutron star merger [18,19], provide a much better understanding of QCD under those extreme conditions.
In the theoretical front, first principles lattice QCD simulations can successfully describe many equilibrium * renato.critelli@usp.br † rrougemont@iip.ufrn.br ‡ noronha@if.usp.br properties of the QGP up to µ B /T ∼ 2 [20,21] including the crossover nature of the QCD phase transition at zero chemical potential [22,23]. However, even lattice QCD simulations suffer from severe limitations, as they are currently unable to deal with larger values of baryon chemical potential (due to the fermion sign problem [4]) and also real time out-of-equilibrium phenomena [24], which are of fundamental importance for the spacetime evolution of the fireball produced in heavy ion collisions and also for the cold and dense matter core of neutron stars. In particular, a crucial question for the beam energy scan (BES) program conducted at RHIC [25], the future fixed target experiments at RHIC [26,27], the compressed baryonic matter (CBM) experiment at FAIR [28,29], and also the upcoming experiments at NICA [30], is the location of the QCD critical end point (CEP), which would mark the end of a first order phase transition line expected to exist in the QCD phase diagram at nonzero baryon density [31][32][33][34].
If the QCD CEP is within the reach of low energy heavy ion collisions 1 , the rapidly expanding medium formed in these collisions may generally pass close to this point when it is still far from thermal equilibrium being perhaps not even close to a hydrodynamical regime. Therefore, the investigation of how the presence of a critical point affects the far-from-equilibrium dynamics of the expanding medium, and its subsequent approach to hydrodynamics, is of fundamental importance. In the absence of first principle QCD calculations of the full dynamical evolution of the system, one must resort to effective models. Holographic shockwave-like simulations of a strongly coupled plasma at finite chemical potential (though without critical phenomena) were investigated in [37] while the spinodal instability in a first-order phase transition at strong coupling was studied in [38]. A way to extend hydrodynamics through the inclusion of parametric slowing down and critical fluctuations induced by a critical point has been recently proposed in [39]. Also, recently out-of-equilibrium effects have been shown to produce significant impact on the critical behavior of the cumulants of fluctuations of conserved charges, which has been investigated using the Fokker-Planck equation [40] and also the Kibble-Zurek formalism [41].
In this paper we employ the gauge/gravity correspondence [42][43][44][45] to investigate the effects of a critical point on the hydrodynamization properties of a top-down holographic construction that stems from string theory known as the 1-R charged black hole (1RCBH) model [46][47][48][49][50][51]. This is the holographic dual of a N = 4 Supersymmetric Yang-Mills (SYM) plasma at finite temperature and chemical potential featuring a critical point in its phase diagram. Although not directly related to QCD, the 1RCBH model is a very attractive setup for studying far-from-equilibrium phenomena at finite temperature and density in the presence of a critical point. This well-defined holographic construction allows one to investigate many different physical aspects of a strongly coupled relativistic fluid with chemical potential and critical behavior under well controlled theoretical conditions. As a matter of fact, in previous works the thermodynamics and hydrodynamics [52], the spectra of quasinormal modes [53], and the homogeneous isotropization and thermalization dynamics [54] of the 1RCBH system have been analyzed in detail.
In the present work we perform, for the first time in the literature, a first principles holographic evaluation of the full far-from-equilibrium evolution of a hot and dense relativistic fluid endowed with a critical point -the 1RCBH plasma -, undergoing a boost invariant expansion known as Bjorken flow [55]. The calculations done in the present work are all performed within a single theoretical framework, the gauge/gravity duality, and since the 1RCBH model is a top-down string theory construction, no auxiliary hypotheses or approximations are needed in order to investigate its physical properties. Such a simple and well controlled theoretical setting constitutes the ideal toy model for understanding general aspects of the influence of a critical point on the hydrodynamization dynamics of out-of-equilibrium strongly coupled media. As an example of a general insight coming from the holographic correspondence concerning hydrodynamics, the most robust result obtained from the duality is the nearly perfect fluidity of strongly coupled plasmas, which is encoded in a very small shear viscosity to entropy density ratio, η/s = 1/4π [56][57][58]. This general feature of holographic models is remarkable close to the results of the latest phenomenological simulations describing heavy ion data [59].
The properties of a SYM plasma undergoing Bjorken flow at zero chemical potential has been previously in-vestigated using numerical general relativity techniques [60][61][62], which determined how long it takes for this system to hydrodynamize or, in other words, what is the relevant timescale that characterizes the emergence of hydrodynamic behavior. In the present work, we generalize these studies to the case of the 1RCBH model at finite chemical potential with a critical point. By numerically evolving the full nonlinear partial differential equations of motion of the system in the bulk using a large set of initial conditions that span different values of µ/T , we evaluate the pressure anisotropy, the charge density, the scalar condensate, and the hydrodynamization timescale. The latter describes the time after which a given far-from-equilibrium solution can be reasonably described in terms of the relativistic Navier-Stokes (NS) equations. As the main result of the present work, we show that the onset of hydrodynamics is considerably delayed as the chemical potential is increased towards its critical value.
Notation: We use a mostly plus metric signature and natural units, i.e, = c = k B = 1.

II. THE HOLOGRAPHIC MODEL
The bulk Einstein-Maxwell-Dilaton (EMD) lagrangian that describes the 1RCBH model is where κ 2 5 ≡ 8πG 5 is the five dimensional gravitational constant, φ is the dilaton field with the scalar potential 2 is the Maxwell field associated with the chemical potential in the field theory, and f (φ) = e −2 √ 2 3 φ is the coupling function between φ and A µ .
The Ansatz for the EMD fields in the holographic Bjorken flow, using in-falling Eddington-Finkelstein coordinates, is given by [60] where τ becomes the usual propertime of Bjorken flow at the boundary, r is the radial coordinate of the asymptotically AdS spacetime whose boundary is at r → ∞, ξ denotes the spacetime rapidity, and (x, y) are the transverse spatial directions to the longitudinal flow direction.
The resulting EMD equations of motion for Bjorken flow constitute a set of coupled partial differential equations (PDE's) similar to the homogeneous isotropization case studied in Ref. [54], From Eq. (3f), one obtains, where Φ 2 is a time-dependent coefficient following from the near-boundary expansion of Φ.
To perform the near-boundary expansion of the bulk EMD fields, which is needed to fix the boundary conditions corresponding to the Bjorken symmetry and obtain the one-point functions of the field theory undergoing Bjorken expansion at the boundary, we set as the boundary condition for the metric field the usual expression for the line element in Bjorken flow using Milne coordinates Hence, by imposing that at the boundary the bulk metric field in Eq. (2) is conformally equivalent to Eq. (6), while the dilaton vanishes and the Maxwell field reduces to the chemical potential of the field theory, one works out the following near-boundary expansions where Φ 0 (τ → ∞) = µ is the gauge theory chemical potential and λ(τ ) is an arbitrary function. This function is associated with the residual diffeomorphism invariance of the line element (2) under radial shifts of the form r → r + λ(τ ) [63]. The relevant observables used to probe the hydrodynamization properties of the system, i.e. the energymomentum tensor T µν , the U(1) four-current J µ , and the scalar condensate O φ are obtained from the onepoint functions via holographic renormalization. Their formulas for the 1RCBH model are [54] where ρ 0 is an input corresponding to the initial charge density (see the appendix for details). The nearboundary coefficients a 2 and φ 2 are dynamical quantities and they can be determined once the equations of motion are solved. The energy density, the parallel and longitudinal pressures, and the charge density are, respectively, Due to Bjorken symmetry, all of these quantities depend solely on the propertime variable τ . Because we have a conformal setup, it suffices to look at ε and ∂ τ ε to probe the hydrodynamization properties of the system since

III. VISCOUS RELATIVISTIC HYDRODYNAMICS WITH A CHEMICAL POTENTIAL
In this section we briefly review the hydrodynamics of a Bjorken expanding fluid with nonzero conserved charge [64,65], which describes the late time behavior of our numerical gravity simulations. We restrict ourselves to first order hydrodynamics corresponding to the relativistic Navier-Stokes equations [66]. Higher order derivative corrections to hydrodynamics [67,68] could be implemented once the corresponding transport coefficients (such as the shear viscosity relaxation time) become available for the system under consideration. The improvement in the hydrodynamic description of this system due to the inclusion of higher order gradient corrections, and the subsequent question about the convergence of the series at strong coupling [69,70], is left to a future study.
Due to conformal invariance and the underlying symmetries of Bjorken flow, to first order in the gradient expansion the only contribution to the viscous evolution comes from the shear stress tensor of the system, π µ ν = −ησ µ ν , where the shear tensor is diagonal with σ µ ν ∼ 1/τ . Since the flow velocity is u µ = (−1, 0, 0, 0) and π µ ν is given in terms of the hydrodynamic variables, the NS equations for Bjorken flow with a U(1) global charge reduce to a single equation for the energy density while charge conservation imposes that the charge density associated with the R-charge of the black hole is Thus, the charge density is known once one specifies the initial charge density ρ 0 .
Using that η/s = 1/(4π) and the well-known thermodynamic relation for conformal theories, 4ε = T s + µρ, Eq. (15) reduces to To proceed, we need the equilibrium equation of state of the 1RCBH model [52][53][54], where, Substituting the above results into Eq. (17), and keeping fixed the dimensionless ratio µ/T ≡ x, it follows that Usingε ≡ κ 2 5 ε and multiplying by τ /ε one obtains where we defined the dimensionless time flow Eq. (23) is needed to evaluate the pressure anisotropy in the hydrodynamical regime, according to Eq. (14) (which holds also far-from-equilibrium). The applicability of hydrodynamics in the late time dynamics of planar shockwave collisions in SYM was first investigated in [71] (see [72] for a comprehensive discussion of this problem). Recently, hydrodynamization was also studied in detail in non-conformal shockwave collisions in Ref. [73]. In the case of Bjorken flow, we follow the same convention of previous works [61,62] and estimate the hydrodynamization time w hydro as the timescale after which a given far-from-equilibrium numerical solution begins to satisfy the following inequality, (25) More specifically, this defines the time after which NS hydrodynamics becomes a good description of the evolution.

IV. RESULTS
In this section we present our results for the far-fromequilibrium Bjorken expanding flow in the 1RCBH model (the numerical details can be found in the appendix). In Fig. 1 we show the time evolution of the one-point functions (8) -(12) characterizing the Bjorken expansion for a sample of numerical solutions with many different initial conditions (see the appendix) and four different equilibrium values of µ/T . Subfigure (a) is the most important one, as it shows the evolution of the pressure anisotropy along with the coalescence to hydrodynamics at late times, with the dashed lines indicating the analytical NS result associated with (23). It is clear from subfigure (a) that, as the value of µ/T increases and approaches the critical point, the numerical solution takes longer to coalesce to hydrodynamics. Subfigure (b) shows the evolution of the charge densityρ = κ 2 5 ρ, whose asymptotic behavior is used to extract the equilibrium value of µ/T in each simulation. Finally, subfigure (c) shows the evolution of the scalar condensate (12), which similarly to the equilibrium solution [54], increases with increasing values of the charge density.
Regarding the hydrodynamization time defined in Eq. (25), we plot in Fig. 2 the results for its variation (as a function of µ/T ) with respect to the SYM result at vanishing chemical potential Fig. 2 emphasizes, in a very clear way, the main result of Fig. 1 and of this work: the onset of hydrodynamics is significantly delayed as the chemical potential is increased towards its value at the critical point of the phase diagram.

V. CONCLUSION AND OUTLOOK
In the present work we presented a first principles holographic calculation of the full nonlinear evolution of a hot and dense (i.e., µ = 0) far-from-equilibrium strongly coupled relativistic fluid with a critical point. We investigated how the top-down holographic construction corresponding to the 1RCBH model, which describes a superconformal non-Abelian plasma with a chemical potential, evolves in space and time undergoing Bjorken flow. We found that increasing µ/T towards its critical value considerably delays the emergence of hydrodynamic behavior, as defined by the relativistic Navier-Stokes equations. This feature of the 1RCBH model, if also applicable to QCD, could imply in important differences for correlation functions calculated in and out of equilibrium, with direct impact on the experimental searches for the QCD critical point, since the main signatures of the critical point are usually considered to be the cumulants of fluctuations of conserved charges [34,40,41].
Regarding the far-from-equilibrium properties of the 1RCBH model, it would also be interesting to investigate the presence of hydrodynamic attractor behavior in Bjorken flow [74][75][76][77][78][79][80] and understand how the critical point affects the properties of such an attractor. Previous studies in the case of SYM at µ = 0 and Gauss-Bonnet holography were already done in [78] and [80], respectively.
We finish this paper by stressing that the 1RCBH model has important differences with respect to the QGP. For instance, the 1RCBH model is a conformal field theory, while the QGP is highly non-conformal in the crossover region. Moreover, the dynamic universality classes [81] are different: while QCD is expected to be in the type H dynamic universality class [82], the 1RCBH system is of type B [52]. The fact that the dynamic universality class of the 1RCBH model is of type B is interesting because in this case η/s is finite at the critical point, which means that viscous hydrodynamics is, in principle, well defined everywhere in the phase diagram. Nonetheless, it is also important to stress that phenomenologically realistic non-conformal holographic settings at finite temperature and chemical potential with a critical point may be constructed using the class of bottom-up holographic models first proposed in Ref. [83], which was recently investigated in [35] and [84]. Therefore, the results presented here may be seen as the first steps towards a more realistic phenomenological holographic description of the far-from-equilibrium hot and baryon dense QGP produced in heavy ion collisions.

APPENDIX
This appendix is devoted to clarify the details of the numerical procedure we employed to evaluate the Bjorken flow dynamics of the 1RCBH model. We stress that the steps followed here are similar to the ones described in Ref. [54] which considered the case of homogeneous isotropization.
To solve numerically the PDE's (3a) -(3f), we first redefine the radial domain in such a way that the new domain is finite which means that in this new radial coordinate the boundary lies at u = 0.
The next step is to define subtracted fields X s by removing the trivial information encoded in the leading order terms of the near-boundary expansions of the bulk fields X, 3 We then rewrite Eqs. (3a) -(3f) in terms of the variables We do not write them explicitly here because the expressions are lengthy and not particularly enlightening.
Since we solve the radial problem using the pseudospectral method [85], the radial u−grid is given by the Chebyshev-Gauss-Lobatto grid u k = u 2 1 + cos kπ N − 1 , k = 0, . . . , N − 1, (36) where N is the number of grid points, also known as collocant points. Here, u defines the infrared (IR) limit of the radial grid. The IR limit of the radial u−grid must be chosen in such a way that it covers the entire portion of the bulk geometry causally connected to the boundary, which means that the locus of the event horizon is a good place to set the IR limit. Commonly, this is done by tracking the apparent horizon u h , and then using λ(τ ) to fix the position of u h in the course of the simulation. In our numerical algorithm, we approximated the above condition by using the following relation where δ is a small negative number, typically of the order O(10 −3 ). Since d + Σ is a monotonically decreasing function for increasing values of u, Eq. (38) tells us that u is a little bit behind the apparent horizon u h ; consequently, if the radial grid spans u ∈ [0, u ], we ensure that it covers the whole portion of the bulk geometry causally connected to the boundary. For the initial conditions considered in this work, which are given at the end of this appendix, one can start with a radial grid u ∈ [0, 1] with λ(τ 0 ) = 0, where τ 0 is the initial time used in the simulation, and in all the cases we considered, the condition (38) could be found after the following steps: i) one starts with d + Σ(τ, u = 1), which is generally equal to some negative number (and, therefore, we know that at this point we are behind the horizon, since at the apparent horizon Eq. (37) holds, and beyond it d + Σ > 0); ii) as the simulation proceeds and we march towards the boundary, the value of d + Σ begins to increase and when the condition (38) is met, the corresponding value of u has been found; iii) in order to keep this value of u fixed through the numerical simulation for a given initial condition, we impose that ∂u /∂τ = 0. Bearing this condition in mind, and manipulating the field equations (3a) -(3f) at u = u , the following relation is found where all the functions in this expression are assumed to be evaluated at u = u . In the limit where δ → 0 we recover the usual expression for the stationary apparent horizon [63]. Expressing A in terms of the subtracted field A s in the LHS of Eq. (39), one obtains the following expression for ∂ τ λ Thus, the condition ∂u /∂τ = 0 provides an expression for ∂ τ λ, which is used to update the value of λ(τ ) through the simulation keeping fixed the position of the approximated apparent horizon u .
Evidently, the initial state that we define must contemplate this set as well.
To perform the time evolution of the system one has to use definition (4) of the expansion (7c); and ∂ τ λ from Eq. (40). By doing so, one can evolve in time the fields necessary to start the cascade solution of the nested set of PDE's a 2 (τ + ∆τ ) = a 2 (τ ) + ∆τ (∂ τ a 2 ), (50) λ(τ + ∆τ ) = λ(τ ) + ∆τ (∂ τ λ), where ∆τ is the time step. There are several ways to do the time evolution. In this work, we choose a fourth-order Adams-Bashforth method to evolve in time.
Regarding the initial conditions used in the present work, they were chosen as follows λ(τ 0 ) = 0, 4 If we know B(τ, u) and a 2 (τ ), then we know ∂τ a 2 . This is because the term b 4 (τ ) is proportional to a 2 (τ ) and ∂τ a 2 : b 4 (τ ) = a 2 (τ )+ 3 4 τ ∂τ a 2 (τ )− 1 where Thus, different values for the set {α, β i , ρ 0 } will produce different evolutions and, consequently, different results for the hydrodynamization time (25). In particular, to generate the results presented in this paper, we have randomly selected values for {α, β i , ρ 0 } in the range, 1 ≤ α ≤ 1.05, −0.5 ≤ β i ≤ 0.5 and 0 ≤ ρ 0 3.7. (59) The fact that ρ 0 does not have a clear upper bound, i.e. ρ 0 3.7, is because some initial conditions do not have a well-behaved evolution for e.g. ρ 0 = 3.7. The reason is because some initial profiles tend to reach the vicinity of the critical point with a smaller value of ρ 0 . Indeed, one cannot predict the final value of µ/T by just looking at the initial conditions since this is an equilibrium quantity which is only determined by the late time evolution of the system, constituting, therefore, a posteriori analysis of the numerical data.
Furthermore, just as in the homogeneous isotropization case [54], we also had to use a filter to avoid spurious growth of hard modes [85]. This is done by accessing the spectral coefficients {a i } and damping the high frequency modes using an exponential map of the type a i → e −(n1/N ) n 2 a i . Typical values for the pair (n 1 , n 2 ) used in our work are in the ballpark of (16,18).
To check the numerical consistency of our results, besides standard tests such as the variation of number of collocant points N and the size of the time step ∆τ , we have also analyzed the late time convergence of the numerical solutions to the analytical NS solution (23); we were also able to reproduce the results of Wilke van der Schee's code 5 valid for µ/T = 0. Another non-trivial consistency check is the observation that the value of δ as defined in Eq. (38) does not change in the course of the simulation. Our code that solves the Bjorken expanding far-from-equilibrium 1RCBH model was initially developed in Mathematica, and later translated to C++ using the Eigen linear algebra package [86] and OpenMP to parallelize the code.