Holographic QCD phase diagram with critical point from Einstein-Maxwell-dilaton dynamics

Supplementing the holographic Einstein-Maxwell-dilaton model of [O. DeWolfe, S.S. Gubser, C. Rosen, Phys. Rev. D83 (2011) 086005; O. DeWolfe, S.S. Gubser, C. Rosen, Phys. Rev. D84 (2011) 126014] by input of lattice QCD data for 2+1 flavors and physical quark masses for the equation of state and quark number susceptibility at zero baryo-chemical potential we explore the resulting phase diagram over the temperature-chemical potential plane. A first-order phase transition sets in at a temperature of about 112 MeV and a baryo-chemical potential of 612 MeV. We estimate the accuracy of the critical point position in the order of approximately 5-8% by considering parameter variations and different low-temperature asymptotics for the second-order quark number susceptibility. The critical pressure as a function of the temperature has a positive slope, i.e. the entropy per baryon jumps up when crossing the phase border line from larger values of temperature/baryo-chemical potential, thus classifying the phase transition as a gas liquid one. The updated holographic model exhibits in- and outgoing isentropes in the vicinity of the first-order phase transition.


Introduction
The QCD phase diagram exhibits potentially a large variety of structures [3][4][5][6]. Either originating from extrapolations of weak-coupling results or being suggested by models (most notably Nambu-Jona-Lasinio (cf. [7]), linear sigma/quark-meson [8] models in numerous variants), various phases of strongly interacting matter may occur, such as color superconductors (cf. [9][10][11]), or quarkyonic matter (cf. [12]), or chirally restored phases (cf. [13]), or color-flavor locked structures (cf. [14]). While the gas-liquid (GL) first-order phase transition (FOPT) in nuclear matter seems to be well established since some time [15][16][17][18][19][20], the hadronquark (HQ) deconfinement transition still offers a few challenges. At very small or zero net-baryon density corresponding to a small chemical potential (µ)-to-temperature (T ) ratio, µ/T 1, the HQ transition is established as a crossover in 2+1 flavor lattice QCD with physical quark masses [21,22] at a characteristic scale of T c = O(150 MeV). The popular Columbia plot [23] sketches qualitatively the options of the phase structure in dependence of the u, d, s quark masses m u,d,s . For instance, in the chiral limit, m u,d,s → 0, or the opposite infinitely heavy quark-mass limit, m u,d,s → ∞, the deconfinement transition is a FOPT. Due to the sign problem of the fermionic determinant the ab initio lattice QCD evaluations are not yet conclusive with respect to the confinement and chiral restoration transition(s) at non-zero baryochemical potential, in particular for µ/T > 2. Some methods try to avoid or circumvent the sign problem (cf. [24]), e.g. by evaluations at imaginary µ (which need a prescription of iµ → µ) or a Taylor expansion in powers of µ/T with coefficients calculated at µ = 0 (which needs statements on the convergence [25]), or the reweighting method (which needs statements on the density and parameter ranges to incorporate the sign and overlap problem [26]). Other approaches are based on the complex Langevin method [27,28] (see [29] for recent developments) or a recent proposal for a path optimization method [30], which is based on the Lefschetz-thimble pathintegral method [31]. The pertinent uncertainties make the region of larger µ/T interesting. A particularly interesting option is the possibility of a (critical) end point (CEP) of a curve of FOPTs, e.g. T c (µ), setting in at (T CEP , µ CEP ) and running toward the T = 0 axis when imaging the phase diagram in the T −µ plane. The CEP coordinates are yet fairly unconstrained. Plugging model results and QCD-related extrapolations together one arrives at some less conclusive scatter plot (cf. e.g. [24]). Advanced lattice QCD approaches disfavor a CEP position at T /T c (µ = 0) > 0.9 and µ/T ≤ 2 [25]. Experimentally, there are dedicated programs aiming at pinning down the CEP location. For instance, the beam energy scan at RHIC [32] gave hints on some features in the beam energy dependence of selected observables which have been interpreted as CEP signature (cf. [33]). In [34] another view has been launched with the conclusion of having also seen CEP indications. Furthermore, the SHINE (NA61) collaboration at CERN-SPS is also systematically seeking CEP effects [35]. Experiments planned at FAIR and NICA and J-PARC [36] are analogously driven by CEP searches, analoguously as goals by the CBM collaboration [37,38], and the MPD group [39]. Given that challenges from both theory and experiment one can ask whether further theoretical model classes beyond the above mentioned approaches could be useful in exploring the hypothetical FOPT emerging from a CEP. Holographic models, advancing the seminal AdS/CFT correspondence [40][41][42], are thought to mimic essential QCD properties in the strong-coupling regime [43][44][45][46][47] and thus may serve as suitable candidates for such an enterprise. Holographic bottom-up approaches coupled to a self-interacting dilaton with nontrivial potential were particularly successful to describe nonconformal properties of the quark-gluon plasma and QCD [48][49][50][51]. In [1,2] a model formulation has been put forward which displays a critical point in the T−µ plane. While [1,2] focuses on CEP properties and an outline of some transport coefficients, [52,53] employed that holographic model to investigate thermodynamics and further transport quantities at small µ/T , however, the question of the CEP position, based on an adjustment to recent lattice data, and properties of phase diagrams were not addressed. The model rests on the coupled Einstein-Maxwell-dilaton (EMd) dynamics and can be adjusted to QCD thermodynamics, i.e. the equation of state (EoS) and quark number susceptibility at µ = 0. The resulting phase structure is the topic of our present paper. We feel that an update of [1,2] is timely since by now consistent and more precise lattice QCD data are at our disposal. In fact, we find some some qualitatively important modifications in comparison to [1,2] w.r.t. the pattern of isentropes in the phase diagrams as well as the position of the CEP. With respect to the discussion in [54], a FOPT curve is specified by further peculiarities: it can be related either to a GL type or to a HQ type transition. For a discussion contrasting features of GL and HQ phase transitions we refer the interested reader to [54][55][56][57], where the notions of entropic vs. enthalpic transitions as well as congruent and non-congruent material changes are exemplified and representations in other variables than T − µ are exhibited. Such different FOPTs can matter significantly in core-collapse supernova explosions as discussed in some detail in [58]. Motivated by such a relation to astrophysical aspects of the phase structure of strongly interacting matter -not only touching core-collapse dynamics but also neutron (quark core) stars -we unravel here the phase structure of the holographic EMd model. It turns out that the EMd model with adjustments to QCD input belongs to the GL class. That is across the phase boundary both the baryon density n and the entropy density s jump when considering the stable phases. For the GL transition, the entropy per baryon s/n drops down when going into µ or T direction, while at the HQ transition s/n jumps up, according to expectations in [54]. According to the Clausius-Clapeyron equation one finds the critical pressure p(T, µ c (T )) either with positive slope (GL transition) or with negative slope (HQ transition). 1 Our paper is organized as follows. In section 2 we recall the holographic EMd model. The numerical adjustment to lattice QCD data at µ = 0 is described in section 3 and the numerical results for the phase diagrams are presented in section 4, including an analysis of the impact of different assumptions for the susceptibility at small temperatures. We summarize in section 5.

Recalling the holographic EMd model
The holographic model of gravity of a 5-dimensional Riemann space sourced by the coupled Maxwell-dilaton fields is defined in [1,2] by the action where R is the Einstein-Hilbert part, F µν = ∂ µ A ν − ∂ ν A µ with A µ dx µ = Φdt stands for the Abelian gauge fieldà la Maxwell, and φ is a real scalar (dilaton) with self-interaction described by the so called potential V (φ). The Maxwell field and dilaton are coupled by a dynamical strength function f (φ). The Gibbons-Hawking term S GH for a consistent formulation of the variational problem is not needed explicitly in our context. The "Einstein constant" κ 5 is taken as a model parameter. The ansatz for the infinitesimal line element squared highlights that (i) only the dynamics in bulk direction r is considered and (ii) a horizon is admitted at r = r H by a simple zero of the blackness function h. By a gauge choice, one can achieve B = 0 and r H = 0. We solve the field equations following from (1, 2) with the technique described in [1,2]. In a nutshell: One has to numerically integrate from r H + towards the boundary at r → ∞. Requiring regularity of A, h, φ, Φ at the horizon r = r H , defined by h(r H ; r H ) = 0, series solutions for any these functions can be obtained, which yield the initial conditions for the integration. After fixing all gauge redundancies the two remaining independent quantities parametrizing the solutions are φ 0 ≡ φ(r H , r H ) and Φ 1 ≡ ∂Φ ∂r r H . It follows from the horizon By the standard AdS/CFT dictionary, φ A is the source and φ B the expectation value of the boundary theory operator dual to φ. Then one obtains the thermodynamic quantities temperature T , entropy density s, baryo-chemical potential µ and baryon density n as The dimensional scaling factors λ T,s,µ,n are introduced as in [1, 2] to compensate the arbitrary choice κ 5 = L = 1 at intermediate steps and restore afterwards the physical units (here, L is the AdS scale). At the horizon, the radially conserved Gauss charge becomes In such a way one maps out the T −µ plane by suitably chosen pairs (φ 0 , Φ 1 ) which entirely parametrize initial conditions at r H ; (4) and (6) deliver in each point s and n.
The pressure follows from the integration of dp(T, µ) = s(T, µ)dT +n(T, µ)dµ, The given bottom-up approach is to be supplemented by fixing the dilaton potential V (φ) and the dynamical coupling f (φ), e.g. from lattice QCD results at µ = 0. By properly engineering V (φ) one essentially dials the EoS at µ = 0, in particular whether a FOPT is built in (as for pure glue dynamics or QCD in the chiral limit(s)) or a crossover is incorporated (as for 2+1 flavor QCD with physical quark masses), see [64,65] for recent examples and [51,63] for full-fledged pioneering investigations. Adjusting f (φ) at the quark number susceptibility from lattice QCD at µ = 0 completes the model. Beyond p − s − n thermodynamics also fluctuation measures, such as susceptibilities, variance, kurtosis, etc. follow via derivatives of the pressure. The susceptibilities are defined, in general, by By CP invariance, odd susceptibilities χ 3,5,... at µ = 0 vanish. In [1,2], a comfortable formula for χ 2 at µ = 0 is given (see also [52,53]) which allows the matching of f (φ) to lattice data. The EMd model (1, 2) with these input data is then ready to transport the information from µ = 0 to µ > 0, up to T = 0, thus uncovering the T −µ plane. This is very much the spirit of the quasi-particle model [66][67][68], where a flow equation facilitates such a transport.
3. Adjustment to lattice QCD data at µ = 0 Contrary to [52,53] we rely here on a modified previous fit [65] to the 2+1 flavor QCD thermodynamics with physical quark masses [21,22] by the dilaton potential implying ∆ = 2(1 + √ 1 − 3a 1 ). A fit of χ 2 /T 2 to data in [69] by the ansatz in (7) f delivers the parameters Uncharged black hole solutions are numerically generated with initial conditions Φ 1 = 0 and φ 0 ∈ [0.35, 5.0]. The resulting equation of state for µ = 0 is shown in Fig. 1 and the corresponding susceptibilities are exhibited in Fig. 2 (blue curves). 4 We emphasize the consistency of the lattice data in [21] and [22] and select the data of [21] for a comparison. The multi-parameter ansätze (8)-(12) allow in fact a fairly precise description of the available data. 5 3 In [70] we allowed for four independent scale setting parameters λ T,s,µ,n as in [1, 2] and a different ansatz for f (φ), resulting in different values for quantities referring to the µ dependence. Here, we enforce λ T = λ µ and λ s = λ n , as in [52,53], to accommodate the only two scales encoded in κ 5 and L. 4 The fourth-order susceptibility χ 4 is calculated using smoothed spline derivatives w.r.t. n(µ, T = const). We verified the robustness of this numerical procedure for different smoothing conditions. Results for χ 2 /T 2 at finite µ are obtained similarly without smoothing technique. 5 The parameters (9), (11) and (12) represent our best fit values and are used for the following studies. Variations in the order of 0.8 % of the potential parameters (9) and 3 % of the coupling function parameters (11) still allow a good description within the uncertainties of the lattice results.

Phase diagram
Charged black hole solutions with initial conditions φ 0 ∈ [0.35, 4.5] and Φ 1 /Φ max 1 (φ 0 ) ∈ [0, 0.755] result in the thermodynamic phase diagram exhibited in Fig. 3 in various variants over the T −µ plane. Only the stablephase quantities are shown, i.e. in the case of multi-valued solutions at a given T −µ point those with maximum pressure. The CEP coordinates are T CEP = (111.5 ± 0.5) MeV and µ CEP = (611.5 ± 0.5) MeV. These uncertainties are estimated through the numerical calculation of discrete levels of constant temperature and chemical potential. 6 The FOPT curve shows up as kinky behavior of the pressure and jumpy behavior of the entropy density, baryon density and entropy-to-baryon ratio (cf. corresponding panels in Fig. 3). The contour curves in the pressure panel in Fig. 3 are scaled isobars, p/T 4 = const. The pressure increases in µ-direction (e.g. at constant T ). The FOPT curve is steeper than neighboring isobars, as characteristic for the GL FOPT. This implies that the critical pressure p c (T ) = p(T, µ c (T )) increases with temperature, see right panel in Fig. 4. The information on both entropy density and baryon density (see top right and bottom left panel in Fig. 3) can be combined to the contour plot of constant entropy per baryon, see bottom right panel in Fig. 3. The resulting contour curves are isentropes, i.e. paths of gas or fluid elements during an adiabatic expansion (collapse) stage in heavy-ion collisions (stellar core collapse). The scaled pressure, entropy density and baryon density are pushed towards higher values with increasing chemical potential, whereas the entropy-to-baryon ratio is decreasing. For µ ≥ µ CEP the scaled entropy density, scaled baryon density and entropy-tobaryon ratio jump across the FOPT. Comparing the entropy-to-baryon ratios on both sides of a point on the FOPT curve evidences s/n| − > s/n| + , where the label −(+) means approaching the FOPT curve from left (right). In line with the above mentioned Clausius-Clapeyron relation dp c (T )/dT = s n − − s n + 1 n − − 1 n + , this implies in turn that the curve p(T, µ c (T )) as a function of T has positive slope, see Fig. 4 (right panel), as typical for the GL transition. Isentropes meeting the FOPT curve are "incoming" on the deconfined/dense (+) side and "outgoing" on the confined/dilute (−) side. Figure 4  in the T − log n plane. The two-phase coexistence region is depicted by a grey region, where the isentropes (not displayed) are to be constructed by the lever rule. The panel exposes the shape of the isentropes as paths of adiabatically expanding and cooling pieces of matter: In the updated EMd model, isentropes enter and leave the coexistence region. This is in contrast to [1,2] where only incoming isentropes can be found. According to the nomenclature of [57], the updated EMd model is classified as type IA and represents a GL phase transition. In addition to the second-order quark number susceptibility at µ = 0, Fig. 2  exhibits also χ 2 /T 2 for finite µ (green and red curve in the left panel). With increasing chemical potential, χ 2 /T 2 is pushed towards larger values. A maximum is evolving, which transforms into a divergence at the CEP. The updated holographic EMd model is based on a fit to the recent lattice data for µ = 0. Since no lattice results are available for low-temperatures and the critical point is located in a temperature range where the lattice data of χ 2,4 just start, we estimate this uncertainty by allowing parameter variations of the dilaton potential and gauge kinetic function within the above mentioned uncertainties and assuming different generic low-temperature asymptotics for the second-order quark number susceptibility and dilaton potential (i.e. continuing the data exhibited in Fig. 2 to zero with different slopes). The different types just slightly vary T CEP in the order of 5 MeV and µ CEP in the order of 50 MeV, corresponding to a relative uncertainty of approximately 5 % and 8 % respectively.

Summary
In summary we explore here the phase structure of the holographic Einstein-Maxwell-dilaton (EMd) model [1,2] adjusted now at 2 + 1 flavor lattice QCD data with physical quark masses at µ = 0. The EMd model has a first-order phase transition (FOPT) curve setting in at a (critical) endpoint (CEP) with coordinates T CEP ≈ 112 MeV, µ CEP ≈ 612 MeV. By considering parameter variations and different low-temperature asymptotics for the second order quark number susceptibility and equation of state that take the uncertainties of the lattice data into account, we estimate the relative uncertainty of our result for the critical point position in the range up to 8 %. We emphasize that these values are consistent with recent lattice estimates in [25] and the covered range of the phase diagram in [71]. The FOPT curve continues from the CEP towards the T = 0 axis. Recalling the general remarks in [52,53] we refrain from analyzing the region of small temperatures (i.e. we leave the quantum phase transition for separate consideration), which however is relevant w.r.t. neutron (quark) stars and particular scenarios for core-collapse supernova explosions. The EMd model does not include explicitly such QCD relevant aspects as chiral symmetry or confinement. Instead, it accounts implicitly for these fundamental notions by the adjustment at QCD results. In particular, our holographic bottom-up approach does not explicitly include the physics of the chiral condensate. Effects as chiral symmetry breaking and restoration can be studied by taking into account the backreaction of flavored branes in the bulk. A holographic model for such important considerations was put forward in [72] (with further developments in [73][74][75][76][77][78]) for the Veneziano limit of QCD. Instead of accommodating all wanted physics solely in the dilaton potential, the gluon and quark degrees of freedom are shared by the glue sector and a coupled flavor sector being dual to the real-scalar dilaton and a complex scalar for the quark condensate qq . This procedure of adding flavors is appropriate to address also magnetic field effects and enjoys a qualitative agreement with lattice results. The relationship of genuinely non-perturbative properties such as color confinement and spontaneous chiral symmetry breaking poses another issue in this context. The resulting phase diagram in the present (minimalistic) approach with a dilaton solely resembles in many aspects the gas-liquid phase transition. For instance, the critical pressure increases with temperature, contrary to expectations of the hadron-quark transition. In the updated EMd model, isentropes are incoming from the dense phase, enter the coexistence region, run through and leave the critical curve at lower temperature. The updated EMd model exhibits a graceful exit into the pure low-temperature and lowdensity phase. We emphasize the need to supplement (model) phase diagrams by information on isobars or isentropes, for instance, for having access to physics implications. Recently, other ansätze for the dilaton potential and gauge kinetic function were presented in [79], which result in CEP coordinates T CEP = (89 ± 11) MeV and µ CEP = (723 ± 36) MeV, only marginally consistent with our result T CEP = (112 ± 5) MeV and µ CEP = (612 ± 50) MeV, where the uncertainties refer to the above quoted 8 %. Since both setups allow an equally good description of the available lattice data for the equation of state and susceptibilities, we conclude that the underlying holographic model is rather sensitive to both the input data and internal parametrization, which affect the CEP position and phase structure. It is in particular the lacking precision lattice data at T ≈ 100 MeV which seems to hamper a unique determination of CEP coordinates.