Hot and dense quark-gluon plasma thermodynamics from holographic black holes
Joaquin Grefa, Jorge Noronha, Jacquelyn Noronha-Hostler, Israel Portillo, Claudia Ratti, Romulo Rougemont
HHot and dense quark-gluon plasma thermodynamics from holographic black holes
Joaquin Grefa, Jorge Noronha, Jacquelyn Noronha-Hostler, Israel Portillo, Claudia Ratti, and Romulo Rougemont Physics Department, University of Houston, Houston TX 77204, USA Illinois Center for Advanced Studies of the Universe, Department of Physics,University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA Departamento de Física Teórica, Universidade do Estado do Rio de Janeiro,Rua São Francisco Xavier 524, 20550-013, Maracanã, Rio de Janeiro, Rio de Janeiro, Brazil (Dated: February 25, 2021)We present new results on the equation of state and transition line of hot and dense stronglyinteracting QCD matter, obtained from a bottom-up Einstein-Maxwell-Dilaton holographic model.We considerably expand the previous coverage in baryon densities in this model by implementingnew numerical methods to map the holographic black hole solutions onto the QCD phase diagram.We are also able to obtain, for the first time, the first-order phase transition line in a wide region ofthe phase diagram. Comparisons with the most recent lattice results for the QCD thermodynamicsare also presented.
I. INTRODUCTION
Significant efforts are underway to search for the quan-tum chromodynamics (QCD) critical point and subse-quent first-order phase transition line at medium- to low-beam energies [1]. Ongoing experiments such as the phaseII of the beam energy scan at the Relativistic Heavy IonCollider (RHIC), including a fixed-target program run-ning at √ s NN = 3 − . GeV [2, 3] and HADES at the GSI,with √ s NN = 1 − GeV [4], are currently looking for theQCD critical point. Additionally, next generation experi-ments such as FAIR at the GSI ( √ s NN = 2 . − . GeV)[5–8] and NICA in Dubna ( √ s NN = 3 − GeV) [9, 10]are being built to precisely determine the QCD equationof state (EOS) and the properties of the strongly interact-ing quark-gluon plasma (QGP) at large baryon densities.Relevant observables in this quest include fluctuations ofconserved charges [11–18], flow [19], and particle yields[20]. For recent reviews see Refs. [1, 21].In order to simulate the evolution of heavy-ion col-lisions at low collision energies, the EOS is needed atlarge baryon chemical potential µ B . First principle latticeQCD calculations provide the EOS at µ B = 0 [22–24].However, due to the fermion sign problem [25], it is notpossible to directly calculate the EOS at finite densities.Nevertheless, one can reconstruct the EOS using suscep-tibilities calculated on the lattice through a Taylor series[26–38], currently limited to µ B /T ≤ (where T is thetemperature). A new expansion has been proposed inRef. [39], which covers a much larger region of µ B withhigh precision. Unfortunately, such an approach cannotcover the whole phase diagram, nor can it accuratelycapture critical behavior. Therefore, one must turn toalternative approaches to describe the matter created inlow-energy collisions, and in the vicinity of a critical point.A promising effective theory should not only reproducelattice QCD thermodynamics results where they are avail-able, but also the QGP’s nearly perfect fluid behavior[40] implied by current extractions of its transport prop-erties from comparisons between model calculations and experimental data [41, 42]. To the best of our knowl-edge, the only effective model available in the literaturethat can simultaneously describe both equilibrium andnear-equilibrium features is the bottom-up non-conformalEinstein-Maxwell-Dilaton (EMD) holographic model pro-posed by some of us in Ref. [43], which is based on thewell-known gauge/gravity duality [44–47] (some previousresults concerning holographic approaches to the physicsof the strongly coupled QGP can be seen in Refs. [48–68]).The EMD model of Ref. [43] predicts a critical point inthe QCD phase diagram at T ∼ MeV and µ B ∼ MeV. However, even though holographic calculations atfinite chemical potentials are not affected by the fermionsign problem, numerical calculations at very large µ B arestill quite challenging in this approach, which was thereason why in Ref. [43] some of us were still unable tolocate the line of first-order phase transition in the regionbeyond the critical point, as we will discuss in detail in thepresent work. Here we considerably expand our previousresults [43] by overcoming most numerical problems andproviding our equation of state over a broad range in tem-perature (2 MeV ≤ T ≤
550 MeV) and baryon chemicalpotential (0 ≤ µ B ≤ ( T, µ B ) plane, we finally locate thefirst-order phase transition line beyond the critical pointof our model. We also present comparisons between ourresults and the latest lattice data from Ref. [39].The paper is organized as follows. In Sections II andIII we review some of the main aspects of the bottom-upEMD model proposed in Ref. [43], which are necessary inthe implementation of our new numerical developmentspresented in detail in Section IV. Also, in Section IV wepresent our results for the thermodynamic quantities ofthe strongly coupled QGP, largely extending the range ofvalues of µ B covered in the phase diagram of the EMDmodel, which allows us to locate the first-order phasetransition past the critical point originally obtained inRef. [43]. In Section V we present our conclusions andfuture perspectives in face of the results discussed here. Inthe present work we employ natural units (cid:126) = c = k B = 1 a r X i v : . [ nu c l - t h ] F e b and a mostly plus metric signature. II. THE HOLOGRAPHIC EMD MODEL
Through the holographic gauge/gravity correspondencedeveloped in string theory, calculations of physical ob-servables in a strongly coupled quantum non-Abeliangauge theory in (flat) four dimensions can be performedby solving the classical equations of motion of a higherdimensional theory of gravity in asymptotically Anti-deSitter (AdS) spacetimes. In the present work we employa five-dimensional bottom-up EMD model defined by thefollowing action [43, 50] S = (cid:90) M d x L = 12 κ (cid:90) M d x √− g ×× (cid:34) R − ( ∂ µ φ ) − V ( φ ) − f ( φ ) F µν (cid:35) , (1)where κ ≡ πG and G is the five-dimensional Newton’sconstant. The EMD action (1) comprises three bulk fieldsin five dimensions: the metric g µν , a real scalar called thedilaton field φ , and a Maxwell field A µ . Additionally, R is the Ricci scalar and F µν = ∇ µ A ν − ∇ ν A µ . We set theasymptotic AdS radius L to unity and introduce as a freeparameter in its place an energy scale Λ , which is goingto be fixed together with κ in Section III C. The singleenergy scale Λ expressed in MeV will be used to writein physical units the gauge theory observables originallycalculated in terms of inverse powers of L on the gravityside of the holographic gauge/gravity correspondence.We note that Eq. (1) is the simplest five-dimensionalaction that can holographically produce a phenomenologi-cally realistic QCD-like effective theory in four dimensionsat finite temperature and chemical potential. In whatfollows, we review some of the main aspects of the EMDmodel already presented in detail in Ref. [43], since theyare important for the new numerical procedure we developin the present work, which shall be discussed in SectionsIV A and IV B.The five-dimensional metric g µν is dual to the stress-energy tensor of the four dimensional quantum gaugetheory and the extra holographic direction may be in-terpreted as a geometrization of the energy scale of therenormalization group flow of the gauge theory [69]. Thedilaton field is used in the present setup to break the con-formal invariance of the theory, with its potential V ( φ ) (and also the free parameters κ and Λ ) being engineeredin a very specific way such as to emulate the behavior ofthe QGP in equilibrium, as inferred from lattice QCD cal-culations at µ B = 0 . The Maxwell field is employed hereto introduce the effects associated with a finite baryonchemical potential, which is done by tuning the couplingfunction f ( φ ) in order to have the holographic baryon sus-ceptibility matching the corresponding lattice QCD resultalso at µ B = 0 . Therefore, as we are going to review inSection III C, all the free parameters of our EMD model are fixed by lattice QCD inputs at zero net baryon density.Consequently, all the observables calculated at nonzero µ B , besides all of those computed at µ B = 0 which werenot used to fix the free parameters of the EMD action,follow as bona fide predictions of our holographic model.We are interested here in five-dimensional, non-rotating,translationally invariant, spatially isotropic, and chargedblack hole backgrounds in thermodynamic equilibrium.In this case, the EMD fields are described by the followinggeneral Ansatz [50] ds = e A ( r ) [ − h ( r ) dt + d(cid:126)x ] + e B ( r ) dr h ( r ) ,φ = φ ( r ) ,A = A µ dx µ = Φ( r ) dt, (2)where r is the holographic coordinate. The radial positionof the black hole event horizon is given by the largest rootof h ( r H ) =0 and the boundary of the asymptotically AdS geometry lies at r → ∞ . The equations of motion (EoM)can be readily obtained φ (cid:48)(cid:48) ( r ) + (cid:20) h (cid:48) ( r ) h ( r ) + 4 A (cid:48) ( r ) − B (cid:48) ( r ) (cid:21) φ (cid:48) ( r ) + (3) − e B ( r ) h ( r ) (cid:20) ∂V ( φ ) ∂φ − e − A ( r )+ B ( r )] Φ (cid:48) ( r ) ∂f ( φ ) ∂φ (cid:21) = 0 , Φ (cid:48)(cid:48) ( r ) + (cid:20) A (cid:48) ( r ) − B (cid:48) ( r ) + d [ln f ( φ )] dφ φ (cid:48) ( r ) (cid:21) Φ (cid:48) ( r ) = 0 , (4) A (cid:48)(cid:48) ( r ) − A (cid:48) ( r ) B (cid:48) ( r ) + φ (cid:48) ( r ) , (5) h (cid:48)(cid:48) ( r ) + [4 A (cid:48) ( r ) − B (cid:48) ( r )] h (cid:48) ( r ) − e − A ( r ) f ( φ )Φ (cid:48) ( r ) = 0 , (6) h ( r )[24 A (cid:48) ( r ) − φ (cid:48) ( r ) ] + 6 A (cid:48) ( r ) h (cid:48) ( r ) ++2 e B ( r ) V ( φ ) + e − A ( r ) f ( φ )Φ (cid:48) ( r ) = 0 , (7)with Eq. (7) being a constraint. Since the backgroundfunction B ( r ) has no dynamics, one may employ a gaugechoice where B ( r ) = 0 in order to simplify the numericalcalculations, as we are going to do in a moment.The equation of motion for Φ( r ) can be integrated toobtain the conserved Gauss charge Q G associated withthe gauge field A µ , Q G ( r ) = f ( φ ) e A ( r ) − B ( r ) Φ (cid:48) ( r ) . (8)From Eq. (6) for the blackening function h ( r ) anotherconserved charge is obtained: the Noether charge Q N = e A ( r ) − B ( r ) (cid:104) e A ( r ) h (cid:48) ( r ) − f ( φ )Φ( r )Φ (cid:48) ( r ) (cid:105) . (9) III. NUMERICAL SOLUTIONS TO THE EOMAND THERMODYNAMIC QUANTITIES
In order to solve the EoM numerically, we need to de-fine a different set of coordinates which we call “numericalcoordinates”, in addition to the so-called “standard coor-dinates”, which will be denoted with a tilde. Both setsof coordinates are defined in the gauge where B ( r ) = 0 .One may calculate the thermodynamic quantities suchas entropy density and temperature from standard holo-graphic formulas using the standard coordinates, in termsof which ˜ h (˜ r → ∞ ) = 1 , as usual. However, to numeri-cally solve the EoM, it is necessary to rescale the standardcoordinates to specify definite values for some of the Tay-lor coefficients in the near-horizon expansions of the EMDfields, as required in order to initialize the numerical in-tegration of the differential equations (3) — (6). Thisrescaling is accomplished using the numerical coordinates,as we discuss next. A. Standard coordinates and thermodynamics
The near-boundary, ultraviolet expansions of the EMDfields in the standard coordinates read [43, 50] ˜ A (˜ r ) = ˜ r + O (cid:0) e − ν ˜ r (cid:1) , ˜ h (˜ r ) = 1 + O (cid:0) e − r (cid:1) , ˜ φ (˜ r ) = e − ν ˜ r + O (cid:0) e − ν ˜ r (cid:1) , ˜Φ(˜ r ) = ˜Φ far + ˜Φ far e − r + O (cid:0) e − (2+ ν )˜ r (cid:1) , (10)where ν ≡ d − ∆ , d = 4 is the number of spacetimedimensions of the boundary gauge theory, ∆ = ( d + √ d + 4 m ) / is the scaling dimension of the gauge fieldtheory operator dual to the dilaton φ ( r ) and m is the massof the dilaton field obtained from the dilaton potential(which will be specified in Section III C).The temperature of the gauge theory fluid correspondsto the Hawking temperature of the black hole solution T = (cid:113) − g (cid:48) ˜ t ˜ t g ˜ r ˜ r (cid:48) π (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˜ r =˜ r H Λ = e ˜ A (˜ r H ) π | ˜ h (cid:48) (˜ r H ) | Λ , (11)where we already introduced the energy scale Λ so thatEq. (11) gives the temperature of the QGP expressed inMeV. The entropy density of the boundary fluid is relatedto the area of the black hole event horizon, A H , via theBekenstein-Hawking formula [70, 71] s = SV Λ = A H G V Λ = 2 πκ e A (˜ r H ) Λ , (12)where V is the 3-dimensional spatial volume. One canalso obtain the baryon chemical potential of the systemfrom the boundary value of the Maxwell field µ B = lim ˜ r →∞ ˜Φ(˜ r )Λ = ˜Φ far Λ , (13) whereas the baryon density is obtained from the boundaryvalue of the radial momentum conjugate to the Maxwellfield ρ B = lim ˜ r →∞ ∂ L ∂ ( ∂ ˜ r ˜Φ) Λ = Q G (˜ r → ∞ )2 κ Λ = − ˜Φ far κ Λ . (14) B. Thermodynamics in the numerical coordinates
For numerically solving the EMD EoM, we considerTaylor expansions of the bulk fields near the black holeevent horizon, (cid:80) ∞ n =0 X n ( r − r H ) n , where X = A, h, φ, Φ .We rescale the holographic coordinate r so that r H = 0 .The fact that the blackening function has a simple zeroat the horizon leads to h = 0 . Also, A = 0 can befixed by rescaling the spacetime coordinates ( t, (cid:126)x ) by acommon factor, while h = 1 can be arranged by rescalingonly the time coordinate t . In addition, one must impose Φ = 0 for Φ dt to be well-defined, since dt has infinitenorm at the horizon. With the Taylor coefficients h , h , A , and Φ determined as aforementioned, the solutionsto the set of equations (3) — (7) via Taylor expansionscan be parametrized by just two coefficients, namely thevalue of the dilaton field calculated at the horizon, φ ,and the derivative of the Maxwell field evaluated at thehorizon, Φ . Indeed, different choices for the pair ofinitial conditions ( φ , Φ ) produce different black holegeometries, each of them corresponding to some definitethermal state of the gauge theory in equilibrium. Thephase diagram of the model can be then populated in the ( T, µ B ) plane by considering a large ensemble of differentblack hole solutions.During the numerical integration of the equations ofmotion, we avoid the singularity at the horizon ( r H = 0 )by starting at a slightly shifted position, e.g r start ≡ − . The boundary can be numerically parameterizedby the value of the holographic coordinate r at whichthe EMD fields have already reached their ultravioletbehavior corresponding to the AdS geometry, which hasa Ricci scalar of R = − . For the set of initial conditionsconsidered here, the corresponding black hole solutionssatisfy this condition for r < ∼ r max = 2 , which is thentaken as the upper bound for the numerical integrationof the EoM.The asymptotics of the EMD fields also imply the follow-ing bound for generating asymptotically AdS solutionsfrom the chosen values of the pair of initial conditions ( φ , Φ ) [43, 50] Φ < (cid:115) − V ( φ ) f ( φ ) ≡ Φ max ( φ ) . (15)In the numerical coordinates, one can show that theultraviolet behavior of the EMD fields is given accordingto [43, 50] A ( r ) = α ( r ) + O (cid:0) e − να ( r ) (cid:1) ,h ( r ) = h far + O (cid:0) e − α ( r ) (cid:1) ,φ ( r ) = φ A e − να ( r ) + O (cid:0) e − να ( r ) (cid:1) , Φ( r ) = Φ far + Φ far e − α ( r ) + O (cid:0) e − (2+ ν ) α ( r ) (cid:1) , (16)where α ( r ) = A far − r + A far . By calculating the constraintEq. (7) at the boundary, one obtains A far − = 1 / (cid:112) h far .Furthermore, by equating the conserved charge (8) evalu-ated at the boundary and at the horizon, one also findsthat Φ far = − (cid:112) h far f (0) f ( φ )Φ . (17)For the kind of calculations we pursue here, we justneed to obtain the behavior of a few ultraviolet expan-sion coefficients of the EMD fields near the boundary,namely h far , Φ far , Φ far , and φ A . As discussed in Ref. [43],one may set h far = h ( r max ) and Φ far = Φ( r max ) , sincethe blackening function and the Maxwell field quicklyreach their respective conformal values. From Eq. (17)one obtains Φ far , while φ A can be reliably estimated byfitting the numerical solution for φ ( r ) using its ultravi-olet asymptotics, φ A e − να ( r ) , within the adaptive range r ∈ [ φ − (10 − ) , φ − (10 − )] .One can show that the thermodynamic variables (11) —(14) can be directly expressed in the numerical coordinatesas follows [43] T = 14 πφ /νA (cid:112) h far Λ , (18) µ B = Φ far φ /νA (cid:112) h far Λ , (19) s = 2 πκ φ /νA Λ , (20) ρ B = − Φ far κ φ /νA (cid:112) h far Λ . (21) C. Fixing the free parameters of the EMD Model
The free parameters of our EMD model, namely, κ , Λ , V ( φ ) , and f ( φ ) , are dynamically fixed by state-of-the-artlattice QCD inputs at µ B = 0 with flavors andphysical values of the quark masses. More specifically, κ , Λ , and V ( φ ) are fixed in order that the holographicequation of state at µ B = 0 closely matches the corre-sponding lattice QCD results from Ref. [23], while f ( φ ) is fixed by requiring that the holographic second orderbaryon susceptibility, also calculated at µ B = 0 , closely
140 160 180 200 220 240 260 280 30000.050.10.150.20.250.3
FIG. 1. Results from the fitting of the holographic suscepti-bility (solid black curve) to the dimensionless second orderbaryon susceptibility χ B ( T, µ B = 0) from lattice QCD [72]. (a) (b)(c) (e) (d) FIG. 2. Thermodynamics at µ B = 0 . Lattice QCD resultsfrom Ref. [23] (red points) are compared to the holographicmodel curves (blue lines): (a) entropy density, (b) speed ofsound squared, (c) energy density (cid:15) , (d) pressure P , and (e)trace anomaly I = (cid:15) − P . matches the corresponding lattice result from Ref. [72]. Inparticular, at vanishing chemical potential it is possible toderive a holographic formula for the dimensionless secondorder baryon susceptibility, χ B ≡ ∂ ( P/T ) /∂ ( µ B /T ) ,which reads as follows [50, 59] χ B ( µ B = 0) = 116 π sT f (0) (cid:82) ∞ r H dre − A ( r ) f − ( φ ( r )) , (22)which is to be evaluated by setting the initial condition Φ to zero. In numerical calculations, we substitute r H → r start and ∞ → r max .In this way, the free parameters of our holographicEMD model are fixed as below, V ( φ ) = −
12 cosh(0 . φ ) + 0 . φ − . φ + 0 . φ ,κ = 8 πG = 8 π (0 . , Λ = 1058 . MeV ,f ( φ ) = sech( c φ + c φ )1+ c + c c sech( c φ ) , (23)where c = − . , c = 0 . , c = 1 . , and c = 100 , withthe corresponding fitting results displayed in Figs. 1 and2. We note that (an approximation for) the pressure canbe easily calculated by integrating the entropy densitywith respect to the temperature, P ( T, µ B = 0) ≈ (cid:90) TT low dT s ( T, µ B = 0) , (24)where we take here T low = 2 MeV (this is the lowest valueof temperature for the black hole solutions generated withthe set of initial conditions considered in the present work,see Section IV A 2).
IV. THERMODYNAMICS AT FINITECHEMICAL POTENTIAL
With the results of Eqs. (18) - (21), we can calculatemany thermodynamic observables at finite temperatureand baryon density. For instance, the internal and freeenergy densities at finite µ B are, respectively, (cid:15) ( s, ρ B ) = T s − P + µ B ρ B , (25) F ( T, µ B ) = − P ( T, µ B ) = (cid:15) ( s, ρ B ) − T s − µ B ρ B (26)from which we can obtain the differential relations d(cid:15) ( s, ρ B ) = T ds + µ B dρ B , (27) dF ( T, µ B ) = − dP ( T, µ B ) = − sdT − ρ B dµ B , (28)so that at fixed µ B , dP ( T, µ B ) | µ B = sdT, (29)and the square of the speed of sound at fixed µ B reads, ˜ c s = dPd(cid:15) (cid:12)(cid:12)(cid:12)(cid:12) µ B = (cid:32) Ts ∂s ( T, µ B ) ∂T (cid:12)(cid:12)(cid:12)(cid:12) µ B + µ B s ∂ρ B ( T, µ B ) ∂T (cid:12)(cid:12)(cid:12)(cid:12) µ B (cid:33) − . (30)Eq. (30) was used in Ref. [43] to calculate the minimum of ˜ c s ( T, µ B ) , which may be used as a “transition line"characterizing the crossover region. However, although(30) is computationally simple to determine along tra-jectories at constant chemical potential, we note that adefinition of the speed of sound that is more relevantto phenomenological applications is the one determinedat constant entropy per particle, which we are going tocompute in this work in IV A. Finally, for completeness,the trace anomaly at finite baryon density is given by I ( T, µ B ) = (cid:15) ( T, µ B ) − P ( T, µ B ) (31) = T s ( T, µ B ) + µ B ρ B ( T, µ B ) − P ( T, µ B ) . A. New numerical procedure
Now we provide some details on the new numericalapproach we developed in the present work, which iscrucial to significantly extend the results originally re-ported in Ref. [43]. With this new numerical procedurewe shall be able to locate the line of first-order phasetransition beyond the critical point of our model and alsoevaluate several thermodynamic observables across the ( T, µ B ) phase diagram, including the phase transition re-gion, where the numerical computations are particularlycomplicated to be performed.
1. Integration of the EMD equations of motion
The equations of motion of the EMD model are solvedwith the MATLab function "ode113". This function im-plements a variable-step, variable-order (VSVO) Adams-Bashforth-Moulton PECE solver of order 13. The pre-cision and stability of this method allow us to explorea wider range of black hole boundary initial conditions( φ , Φ ) than other methods available in MATLab. Theroutine used to integrate the EMD fields and find theQCD thermodynamic observables from Eqs. (18) — (21)checks crucial behavior for the stability and physical con-sistency of the holographic black hole (BH) solutions. ABH-solution is accepted if it satisfies the following require-ments:• The integration of the equations of motion (3) —(6) is finite.• The constraint equation (7) is satisfied.• The dilaton field φ ( r ) tends to zero with the correctultraviolet asymptotics (16) as we approach theboundary.• The near-boundary behavior of all the other EMDfields also respects the correct ultraviolet asymp-totics (16).• The metric coefficient A ( r ) is monotonically increas-ing. FIG. 3. Mapping of an equally spaced, rectangular grid of initial conditions ( φ , Φ ) into an irregular grid of points in the ( T, µ B )plane generated by the corresponding black hole solutions. • The Ricci scalar of the black hole background, R ,is equal to -20 at the ultraviolet radial cutoff r max (meaning that the geometry is already AdS at thispoint).
2. Mapping QCD thermodynamics from the black hole initialconditions
For the results presented in Ref. [43], × black holeswere generated with initial conditions spanning the rect-angle defined by φ ∈ [0 . , and Φ ∈ [0 , . max ( φ ) .Fig. 3 shows how an equally spaced, rectangular grid ofinitial conditions ( φ , Φ ) is mapped into an irregular gridin the ( T, µ B ) plane generated by the associated blackhole solutions.As seen in Fig. 3, a simple rectangular and uniformgrid of initial conditions ( φ , Φ ) produces a wide regionof the ( T, µ B ) plane which is not covered in the QCDphase diagram (shown in white in the figure). In orderto cover the missing section, we introduce here a newway of choosing the black hole initial conditions, which isillustrated in Fig. 4.We first consider Φ = 0 (which implies solutions with µ B = 0 ) and choose the values for φ such that themapping to the solutions in the temperature axis (at µ B = 0 ) is equally spaced in intervals of . MeV from T = 2 MeV to T = 550 MeV. Next, for each chosenvalue for φ , Φ is varied to map the QCD phase diagramcompletely up to µ B = 1100 MeV, leading to the lines ofconstant φ shown in Fig. 4. These lines bend in the QCDphase diagram, giving rise to a region with three layersof competing black hole solutions corresponding to thesame ( T, µ B ) points. In this region, the model is limitedat low T by the end of the lines of constant µ B , were BH-solutions cannot be found using our numerical procedure.It is worth noticing that µ B ∼ MeV is the highestvalue of µ B that can be obtained before the BH-solutions for the more curved lines in Fig. 4 start diverging andbecome unstable, which occurs approximately for valuesof Φ > ∼ . max ( φ ) . In general, a BH-solution cannotbe computed when Φ passes this threshold.For the ensemble of BH-solutions used in the presentwork, each line of constant φ has 3000 BH-solutionsseparated uniformly along these lines and correspondingto different values of Φ , populating the region of the QCDphase diagram within the rectangle defined by T ∈ [2 , MeV and µ B ∈ [0 , MeV, without the holes found inRef. [43] by using a rectangular grid of initial conditions ( φ , Φ ) , as shown in Fig. 3.The precision of the calculations is significantly affectedby numerical noise associated with the fitting of the ul-traviolet coefficients in Eq. (16). The more sensitivecoefficient is φ A , which appears in the holographic ther-modynamic formulas (18) — (21) raised to the powersof − /ν and − /ν . The noise associated with the loss ofnumerical precision is not the same for all lines of constant φ , as shown in the left panels of Figs. 5 — 7.The behavior of the ultraviolet coefficients and thethermodynamic variables, as functions of the BH initialconditions, changes for different lines of constant φ as Φ increases. The value of Φ for the lines close to the QCDphase transition (i.e. lines starting between T = 150 and T = 180 MeV at µ B = 0 ) increases much faster than forthe other lines and its behavior is not as simple as for therest of the lines. Therefore, the treatment of the lines isdifferent depending on their location with respect to theQCD transition line.The strategy to get a smooth mapping is to filter thelines over a large number of BH-solutions. The mappingin Fig. 4 contains 3000 BH-solutions per line of constant φ . Taking a large number of solutions allows us totreat a noisy line with the appropriate filters withoutcompromising its actual behavior.For lines with < T < MeV, Φ / Φ max ( φ ) consid-erably increases (see the color scheme used in Fig. 4, which FIG. 4. Example of how the black hole initial conditionsshould be chosen to map a rectangular region in the QCDphase diagram. allows to identify how the different initial conditions mapinto the ( T, µ B ) plane), and the filtering process consistsin smoothing out these lines using a Cubic SmoothingSpline (CSS) filter which only gets rid of big bumps, andthen filtering the line with a Savitzky-Golay (SG) filter.SG filters are typically used to smooth out a noisy signalwith large noise frequency. For this reason, it is impor-tant to prepare the signal with the CSS filter. The SGfilter employed during this process uses a polynomial ofdegree 3 to interpolate each point with its neighbors. Thenumber of neighbors approximate a range of ± MeV.The rest of the lines are noisy, but the value of Φ / Φ max ( φ ) remains small. In this case, the most noisyultraviolet coefficient is φ /νA , which is corrected by usinga simple polynomial fitting of the form a + bx + cx .The remaining ultraviolet coefficients are filtered withthe SG filter. Notice that the concavity of φ /νA changesfrom positive at small T to negative at large T . Theregion in between is where the BH-solutions can be found with less noise and those lines are the ones that cross thecritical point. Figs. 5 — 7 show the lines of constant φ as functions of µ B for different fixed temperatures, before(blue curves) and after (red curves) the filtering process.Once the lines of constant φ are corrected, they arefitted with a cubic spline to get lines of constant µ B . Thelines of constant µ B are also treated with the SG filter.An example is given in Fig. 8, which shows the baryondensity as a function of the temperature for differentvalues of µ B , before and after the filter.The lines of constant µ B are then fitted with a cubicspline to calculate the pressure, the critical point andthe first-order phase transition line. The next step is tocalculate lines of constant T which, together with thelines of constant µ B , are used to take derivatives of theQCD thermodynamic variables. -4 FIG. 5. Line of Constant φ for T = 60 MeV before and afterthe filtering process. FIG. 6. Line of constant φ for T = 180 MeV before and afterthe filtering process. FIG. 7. Line of constant φ for T = 250 MeV before and afterthe filtering process.
85 87 89 91 930.00.30.60.91.21.5 85 87 89 91 93
FIG. 8. Dependence of the baryon density ( ρ B ) on the temper-ature ( T ) for different values of the baryon chemical potential( µ B ) before and after the filtering process.
3. Finding the transition line and the QCD critical point
With an equally-spaced rectangular grid in the QCDphase diagram, we can start to analyze the region withmultiple solutions. The upper panel of Fig. 9 shows linesof constant φ as Φ increases, where we can distinguishthree types of lines that define a region of overlappingsolutions for the thermodynamics of the holographic EMDmodel. Three different colors have been used to easilyidentify the multi-solution region in the figure. The blackdotted lines are almost parallel and do not cross eachother. Some of the dashed red lines cross each other andalso the black lines. Finally, the solid blue lines on thetop cross the black and red lines and some cross eachother as well. The location where these lines start tointersect can be identified as a candidate point for thecritical end point (CEP) in the QCD phase diagram. Dueto the presence of the first-order phase transition line,the competing phases may appear as solutions of theequations of motion, although only one minimizes thefree energy and represents the true ground state of thesystem. In the crossover region, one expects only onesolution to the equations of motion. However, near the first-order phase transition line, to the right of the criticalpoint, the black hole solutions for the baryon density ρ B and the entropy density s become multivalued functionsof ( T, µ B ) . The first-order phase transition line and themultivalued solutions end precisely at the CEP.In order to find the exact location of the CEP, one cananalyze the second order baryon susceptibility χ B , whichdiverges at the critical point. The behavior of χ B as afunction of T and µ B is shown in the lower panel of Fig.9. With this procedure, we find that the critical pointis located at T CEP ∼ MeV and µ CEP B ∼ MeV, asoriginally reported in Ref. [43].
FIG. 9. The upper panel shows the mapping from the BHinitial conditions ( φ , Φ ) to the QCD phase diagram ( T, µ B ).The plot shows how three kinds of lines of constant φ aremapped into the ( T, µ B ) plane, where the crossing of thelines suggest the location of the CEP. The lower panel showsthe behavior of the second order baryon susceptibility χ B inthe ( T, µ B ) plane. As the chemical potential increases, χ B develops a peak that becomes a divergence at the critical pointlocated at T CEP ∼ MeV and µ CEP B ∼ MeV.
From the highly nonlinear and unequally spaced map-ping showed in Fig. 3, it is possible to obtain the ther-modynamics of QCD on a regular grid in the ( T, µ B ) plane by means of numerical interpolation as done in [43].
65 75 85-2024 10
65 75 85 65 75 85
FIG. 10. Entropy density s (upper panels) and its integral with respect to the temperature, corresponding to the pressure (lowerpanels), for three different values of µ B > µ CEP B ∼ MeV.
Holographic QCD Phase Diagram
FIG. 11. The phase diagram of our EMD model. The inflectionpoint of χ B and the minimum of c s from Eq. (35) are used tocharacterize the crossover region. In particular, the baryon density ρ B was obtained overa regular grid in the interval T = [65 − MeV and µ B = [0 − MeV via numerical interpolation.In this work, however, we obtain an equally spacedgrid in the ( T, µ B ) plane directly from the black holesolutions as described in Section IV A 2, by taking theblack hole initial conditions as shown in Fig. 4. One ofthe advantages of having the thermodynamics over anequally spaced grid in the QCD phase diagram is theopportunity to look at the entropy and baryon density, s and ρ B , respectively, over trajectories of constant T or µ B in the crossover region and near the first-order phasetransition line. For instance, for an isotherm at T > T
CEP or for slices of constant µ B < µ CEP B , the entropy densityand baryon density are single-valued functions, since they do not cross the first-order phase transition line. Onthe other hand, for trajectories of constant T < T
CEP or µ B > µ CEP B , i.e. trajectories that cross the first-orderphase transition line, s and ρ B become multivalued. Sincewe are solving the holographic black hole equations ofmotion, it is reasonable to obtain all extrema of thefree energy which corresponds to the coexistence regionof not only thermodynamically stable minima, but alsothermodynamically metastable and unstable saddle pointsor maxima. In the top panels of Fig. 10, we can observethe characteristic multivalued S-shape for the entropy atthree different slices of µ B > µ CEP B , which means that at agiven T we have three competing BH-solutions. Preciselyat ( T CEP , µ
CEP B ) , the curves for s and similarly for ρ B cease to be multivalued; this characterizes the end of thefirst-order phase transition line at the CEP.Our approach to characterize the first-order phase tran-sition line was to integrate the entropy with respect tothe temperature over the multivalued region, and locatethe point where the resulting curve, corresponding to thepressure or to minus the free energy according to Eqs. (26)and (29), crosses itself. This method is close/analogousto Maxwell’s equal area construction, although computa-tionally easier to implement.Since the QCD deconfinement transition from µ B = 0 up to the critical point is a smooth crossover, there is nounique definition of a transition temperature in this region.However, one may try to characterize this quantity as theinflection point or the extrema of observables sensitive tothe change of degrees of freedom in the transition betweenhadrons and a system of quarks and gluons. In thiswork, we have chosen to characterize the deconfinementtransition in the crossover region by both the inflectionpoint of the second order baryon susceptibility χ B , and0
120 140 160 180 200 220 24000.511.522.533.544.5120 140 160 180 200 220 240024681012141618120 140 160 180 200 220 24000.20.40.60.811.2 B /T=0 B /T=0.5 B /T=2.0 B /T=1.0 B /T=1.5 B /T=2.5 B /T=3.0 B /T=3.5 FIG. 12. Pressure (top panel), entropy density (center panel)and baryon density (bottom panel) as functions of the temper-ature, for different values of µ B /T . Our curves are comparedto the latest lattice QCD results from Ref. [39]. the minimum of the square of the speed of sound c s atconstant entropy per particle. The baryon susceptibilities are generally defined as: χ n = ∂ n ( P/T ) ∂ ( µ B /T ) n , (32)which are basically the coefficients in the Taylor expansionof the pressure P ( T, µ B ) − P ( T, µ B = 0) T = ∞ (cid:88) n =1 n )! χ n ( T ) (cid:16) µ B T (cid:17) n , (33)and the baryon density ρ B ( T, µ B ) T = ∞ (cid:88) n =1 n − χ n − ( T ) (cid:16) µ B T (cid:17) n − . (34)In particular, χ B measures the equilibrium response ofthe baryon density to a change in the chemical potentialof the medium. The square of the speed of sound at con-stant entropy per particle is defined as c s = ( ∂P/∂(cid:15) ) s/ρ B ,but this definition is not practical when one wants tocalculate it on top of a ( T, µ B ) grid of points. For thisreason, it is advantageous to rewrite this state variable interms of derivatives of the pressure along lines of constanttemperature or chemical potential only [35, 73]: c s = ρ B ∂ T P − sρ B ∂ T ∂ µ B P + s ∂ µ B P ( (cid:15) + P )[ ∂ T P ∂ µ B P − ( ∂ T ∂ µ B P ) ] . (35)The lines in the phase diagram corresponding to theminimum of the square of speed of sound computed fromEq. (35) (red, dashed) and to the inflection point of χ B (black, dash-dotted) are shown in Fig. 11. As shown inthe figure, the two lines are separated in temperature at µ B = 0 , while they meet at the critical point. The first-order phase transition line, computed using the schemeimplemented in this paper, is plotted as a blue full linein Fig. 11.The dependence of the transition temperature (definedby the minimum of c s ) on the chemical potential in thecrossover region can be characterized by the followingtruncated series T c ( µ B ) T c (0) = 1 − κ (cid:18) µ B T c (0) (cid:19) − κ (cid:18) µ B T c (0) (cid:19) , (36)with T c (0) = 143 . MeV, κ = 0 . , and κ = − . .In the case of the most recent lattice QCD transition lineobtained from the inflection point of the chiral condensateand its susceptibility in the crossover region, the valuesof κ and κ are . and . , respectively[74]. It should be noted, however, that these expansioncoefficients for the minimum of c s and for the inflectionof the chiral condensate do not need to agree, since thecorresponding transition curves are actually different inthe crossover region.1 FIG. 13. Entropy density (top left), baryon density (top right), pressure (bottom left), and square of the speed of sound atconstant entropy over baryon number as computed via Eq. (35) (bottom right). The CEP is shown on the pressure surface as ared dot.
B. Equation of state
The comparison between the holographic EMD equa-tion of state and the Taylor-expanded lattice QCD equa-tion of state up to µ B /T ≤ [29, 32] was presented inRef. [43]. In that work, we also predicted the location ofthe QCD CEP to lie at ( T, µ B ) ∼ (89 , MeV but, atthat time, due to numerical difficulties, we were unableto identify the location of first-order phase transition linebeyond the CEP and also to calculate the thermodynamicobservables in the phase transition region. These numeri-cal difficulties were solved in the present work through thenew developments discussed in previous sections, namely:i) the new way of choosing the BH initial conditions illus-trated in Fig. 4, which allowed us to cover a much largerregion of the ( T, µ B ) phase diagram, including the loca-tion of the first-order phase transition line; ii) the filteringscheme, which allowed us to obtain smooth results for thephysical observables in the phase transition region, wheretheir computation is plagued by strong numerical noise. In fact, without the filtering process, the result for theentropy density in the phase transition region is so noisythat it becomes impossible to obtain sensible results forthe pressure by integrating the entropy, which in turnmakes it impossible to correctly identify the thermody-namically stable BH-solutions in the multi-solution phasetransition region and the first-order phase transition line.With the aforementioned technical developments, inthe present work we largely extend the coverage of theEMD model on the ( T, µ B ) plane and present our resultsfor the equation of state in a broader region of the phasediagram and also compare our results with the most up-to-date lattice data for the QCD equation of state, whichis now available up to the unprecedentedly high value of µ B /T = 3 . [39].In Fig. 12 we show how the holographic EMD equationof state compares to the lattice data from Ref. [39]. Wesee that the entropy density predicted by the EMD modelis in quantitative agreement with the lattice results for allthe values of T and µ B currently covered by lattice simu-lations. Regarding the pressure, there is also quantitative2agreement for most of the values of T and µ B , althoughthe EMD result starts to deviate from the lattice outcomefor the pressure for µ B /T ≥ . in the high temperatureregion with T > ∼ MeV. With respect to the baryondensity, our results are in quantitative agreement withthe lattice simulations for all the values of µ B /T and tem-peratures up to T ∼ MeV, although the holographicEMD prediction overestimates the lattice results for thebaryon density at high temperatures
T > ∼ MeV when µ B /T > ∼ . . Interestingly enough, for lower temperatures T < ∼ MeV, where the transition from a hadron gas tothe quark-gluon plasma phase takes place, the holographicEMD predictions for the entropy density, the pressure,and the baryon density are in quantitative agreement withthe lattice results all the way up to µ B /T = 3 . , whichsuggests that our prediction for the behavior of the QCDphase transition at nonzero baryon chemical potentials isrobust.In Fig. 13 we show the surface plots of the entropydensity, baryon density, pressure and square of the speedof sound in the ( T, µ B ) plane. We obtained the temper-ature, baryon chemical potential, entropy and baryondensity directly from the holographic dictionary given byEqs. (18) — (21). The pressure is found by integratingthe entropy with respect to the temperature at constantbaryon chemical potential as in Eq. (29), and it can alsobe computed as the integral of the baryon density withrespect to the baryon chemical potential along isothermsas suggested in Eq. (28), which produced the same resultand served as a cross check. The second order baryonsusceptibility χ B shown in Fig. 9 is found as the deriva-tive of the baryon density along the chemical potentialdirection.The critical point manifests itself in the first orderderivatives of the pressure, namely the entropy and baryondensity, where the pronounced gap shown by these statevariables corresponds to the first-order phase transitionline for µ B > µ CEP B . In addition, c s exhibits a dip thatbecomes a zero at the CEP (which is a second-order phasetransition point), and a discontinuity for µ B > µ CEP B along the first-order phase transition line, as expectedfrom thermodynamic considerations. The location of thecritical point is shown on the pressure surface as a redspot.The dependence of the thermodynamic state variableson the temperature along lines of constant µ B is presentedin Fig. 14. There one can see the features mentionedabove more explicitly such as the jump in entropy, baryondensity, and χ B after the first-order phase transition lineis reached. In [43], it was noted that the peak formationin χ B may already indicate that a critical point is presentat larger densities, and here this peak begins to happenaround µ B = 416 MeV. We see that our results appear tobe consistent with this idea.It is worth mentioning that, in order to obtain the c s atconstant entropy per particle we calculate c s from Eq. (35)and it was necessary to remove the noise associated withthe second order derivatives of the pressure. For a region up to µ B ≤ MeV, a SG filter was employed consideringthat the minimum of this observable is not very deep.However, due to the fact that the noise increases nearthe critical point for any observable in addition to theexpected divergences in some derivatives, ∂ T P and ∂ µ B P for example, and the fact that at the critical point weexpected a minimum, it was not possible to remove thenoise without affecting the shape and features of the speedof sound.Therefore, another approach was implemented to obtainthe region of large chemical potential in front of the criticalpoint. By obtaining the isentropic trajectories, one canalso compute the speed of sound by taking the simplederivative ∂P/∂(cid:15) along the isentropic paths. The lines ofconstant entropy over baryon density s/ρ B are relevantin the QCD phase diagram since they approximate thetrajectories that the systems created in relativistic heavyion collisions follow during their evolution when viscouseffects are neglected. In fact, in the ideal case of vanishingviscosity, the quantity s/ρ B is conserved because theentropy generation is only caused by particle generation(although it has been shown that at large baryon densitiesand near a critical point, large deviations from this canbe expected [75]).The regions where the c s computed from Eq. (35) washeavily affected by the noise were removed and replacedwith the information given by computing the same observ-able along the isentropic trajectories through an interpo-lation. Some isentropic lines are shown in the left panel ofFig. 15, along with the dependence of the c s with respectto the temperature for different values of the chemicalpotential (right panel). V. CONCLUSIONS AND FUTUREDIRECTIONS
In this work we presented our most recent results andpredictions on the thermodynamics and phase diagramof strongly interacting QCD matter, obtained througha bottom-up non-conformal holographic approach. Wesignificantly improved our numerical techniques, thus ex-panding our coverage in temperature and baryon densitywell beyond our previous work [43]. We were able toobtain the first-order phase transition line beyond thecritical point for the first time out to µ B ∼ MeV. Agood agreement is found between our equation of stateand the corresponding lattice QCD results available atintermediate densities. The equation of state obtainedhere, covering the first-order phase transition line andthe critical end point in the ( T, µ B ) plane, can be readilyused in hydrodynamic simulations of the matter createdin heavy-ion collisions. Furthermore, these results canbe used to build a bridge to the high-density, low tem-perature region of the QCD phase diagram needed inthe description of neutron star mergers [76]. We pointout that in our model, we do not see c s ever surpassingthe conformal limit of c s = 1 / , despite predictions from3 FIG. 14. Equation of state as a function of the temperature for several values of the baryon chemical potential. The discontinuityexhibited by the entropy density, baryon number density and susceptibility at µ B = 850 MeV corresponds to the line of firstorder phase transition. The divergence of χ B at the critical point is clearly visible in the bottom right panel. FIG. 15. Left: Several isentropic trajectories in the ( T, µ B ) plane. The red dashed line is the minimum of the square of thespeed of sound that terminates at the critical point. The black curve is the first order transition line. Right: Square of thespeed of sound at constant entropy per baryon number for different values of the chemical potential. T = 0 andbaryon densities above nuclear saturation density [77–86].However, the holographic approach employed here is notexpected to be a good guide for the behavior of stronglyinteracting matter at T = 0 . Rather, our model is ex-pected to be useful at finite temperature and chemicalpotential where it is conceivable that QCD matter stilldisplays near perfect fluidity.We expect to report in the near future new results com-ing from the present holographic model, regarding thecalculation of transport coefficients such as the bulk vis-cosity, the baryon and thermal conductivities, the baryondiffusion coefficient, the jet quenching parameter, theheavy quark drag force and the Langevin diffusion co-efficients, all of them evaluated across the entire phasediagram covered in the present work, including the phasetransition region. It would be interesting to see a modelsuch as ours, which is consistent with lattice QCD results, to be applied in real time out-of-equilibrium studies suchas those performed in [62, 87–94]. Such study is, however,beyond the scope of the present work. ACKNOWLEDGEMENTS
This material is based upon work supported by the Na-tional Science Foundation under grant no. PHY-1654219and by the US-DOE Nuclear Science Grant No. DE-SC0020633, US-DOE Office of Science, Office of NuclearPhysics, within the framework of the Beam Energy ScanTopical (BEST) Collaboration. J.N. is partially supportedby the U.S. Department of Energy, Office of Science, Of-fice for Nuclear Physics under Award No. DE-SC0021301.R.R. acknowledges financial support by Universidade doEstado do Rio de Janeiro (UERJ) and Fundação Car-los Chagas de Amparo à Pesquisa do Estado do Rio deJaneiro (FAPERJ). [1] A. Bzdak, S. Esumi, V. Koch, J. Liao, M. Stephanov,and N. Xu, Phys. Rept. , 1 (2020), arXiv:1906.00936[nucl-th].[2] S. collaboration, “Studying the Phase Diagram of QCDMatter at RHIC,” (2014).[3] D. Cebra, S. G. Brovko, C. E. Flores, B. A. Haag, andJ. L. Klay, (2014), arXiv:1408.1369 [nucl-ex].[4] T. Galatyuk (HADES), Nucl. Phys. A , 41 (2014).[5] V. Friese, Nucl. Phys. A , 377 (2006).[6] N. A. Tahir et al. , Phys. Rev. Lett. , 035001 (2005).[7] M. F. M. Lutz et al. (PANDA), (2009), arXiv:0903.3905[hep-ex].[8] M. Durante et al. , Phys. Scripta , 033001 (2019),arXiv:1903.05693 [nucl-th].[9] V. Kekelidze, A. Kovalenko, R. Lednicky, V. Matveev,I. Meshkov, A. Sorin, and G. Trubnikov, Nucl. Phys. A , 884 (2017).[10] V. Kekelidze, A. Kovalenko, R. Lednicky, V. Matveev,I. Meshkov, A. Sorin, and G. Trubnikov, Nucl. Phys. A , 846 (2016).[11] M. A. Stephanov, Phys. Rev. Lett. , 032301 (2009),arXiv:0809.3450 [hep-ph].[12] M. A. Stephanov, Phys. Rev. Lett. , 052301 (2011),arXiv:1104.1627 [hep-ph].[13] R. Bellwied, J. Noronha-Hostler, P. Parotto, I. Por-tillo Vazquez, C. Ratti, and J. M. Stafford, Phys. Rev.C , 034912 (2019), arXiv:1805.00088 [hep-ph].[14] J. Adamczewski-Musch et al. (HADES), Phys. Rev. C , 024914 (2020), arXiv:2002.08701 [nucl-ex].[15] J. Adam et al. (STAR), (2020), arXiv:2001.02852 [nucl-ex].[16] R. Bellwied, S. Borsanyi, Z. Fodor, J. N. Guenther,J. Noronha-Hostler, P. Parotto, A. Pasztor, C. Ratti,and J. M. Stafford, Phys. Rev. D , 034506 (2020),arXiv:1910.14592 [hep-lat].[17] P. Alba, V. M. Sarti, J. Noronha-Hostler, P. Parotto,I. Portillo-Vazquez, C. Ratti, and J. M. Stafford, Phys.Rev. C , 054905 (2020), arXiv:2002.12395 [hep-ph]. [18] D. Mroczek, J. Noronha-Hostler, A. R. N. Acuna,C. Ratti, P. Parotto, and M. A. Stephanov, (2020),arXiv:2008.04022 [nucl-th].[19] B. Kardan (HADES), Nucl. Phys. A , 431 (2019),arXiv:1809.07821 [nucl-ex].[20] J. Adamczewski-Musch et al. (HADES), Nature Phys. ,1040 (2019).[21] C. Ratti, Rept. Prog. Phys. , 084301 (2018),arXiv:1804.07810 [hep-lat].[22] S. Borsanyi, G. Endrodi, Z. Fodor, A. Jakovac, S. D.Katz, S. Krieg, C. Ratti, and K. K. Szabo, JHEP ,077 (2010), arXiv:1007.2580 [hep-lat].[23] S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz,S. Krieg, and K. K. Szabo, Phys. Lett. B , 99 (2014),arXiv:1309.5258 [hep-lat].[24] A. Bazavov et al. (HotQCD), Phys. Rev. D , 094503(2014), arXiv:1407.6387 [hep-lat].[25] O. Philipsen, Prog. Part. Nucl. Phys. , 55 (2013),arXiv:1207.5999 [hep-lat].[26] C. R. Allton, S. Ejiri, S. J. Hands, O. Kaczmarek,F. Karsch, E. Laermann, C. Schmidt, and L. Scorzato,Phys. Rev. D , 074507 (2002), arXiv:hep-lat/0204010.[27] C. R. Allton, M. Doring, S. Ejiri, S. J. Hands, O. Kacz-marek, F. Karsch, E. Laermann, and K. Redlich, Phys.Rev. D , 054508 (2005), arXiv:hep-lat/0501030.[28] S. Borsanyi, G. Endrodi, Z. Fodor, S. D. Katz, S. Krieg,C. Ratti, and K. K. Szabo, JHEP , 053 (2012),arXiv:1204.6710 [hep-lat].[29] A. Bazavov et al. , Phys. Rev. D , 054504 (2017),arXiv:1701.04325 [hep-lat].[30] M. D’Elia and M.-P. Lombardo, Phys. Rev. D , 014505(2003), arXiv:hep-lat/0209146.[31] M. D’Elia, G. Gagliardi, and F. Sanfilippo, Phys. Rev.D , 094503 (2017), arXiv:1611.08285 [hep-lat].[32] J. N. Guenther, R. Bellwied, S. Borsanyi, Z. Fodor, S. D.Katz, A. Pasztor, C. Ratti, and K. K. Szabó, Nucl. Phys.A , 720 (2017), arXiv:1607.02493 [hep-lat].[33] S. Borsanyi, Z. Fodor, J. N. Guenther, S. K. Katz, K. K.Szabo, A. Pasztor, I. Portillo, and C. Ratti, JHEP ,
205 (2018), arXiv:1805.04445 [hep-lat].[34] A. Bazavov et al. , Phys. Rev. D , 074502 (2020),arXiv:2001.08530 [hep-lat].[35] P. Parotto, M. Bluhm, D. Mroczek, M. Nahrgang,J. Noronha-Hostler, K. Rajagopal, C. Ratti, T. Schäfer,and M. Stephanov, Phys. Rev. C , 034901 (2020),arXiv:1805.05249 [hep-ph].[36] J. Noronha-Hostler, P. Parotto, C. Ratti, andJ. M. Stafford, Phys. Rev. C , 064910 (2019),arXiv:1902.06723 [hep-ph].[37] A. Monnai, B. Schenke, and C. Shen, Phys. Rev. C ,024907 (2019), arXiv:1902.05095 [nucl-th].[38] D. Everett et al. (JETSCAPE), (2020), arXiv:2010.03928[hep-ph].[39] S. Borsanyi, Z. Fodor, J. N. Guenther, R. Kara, S. D.Katz, P. Parotto, A. Pasztor, C. Ratti, and K. K. Szabo,(2021), arXiv:2102.06660 [hep-lat].[40] U. Heinz and R. Snellings, Ann. Rev. Nucl. Part. Sci. ,123 (2013), arXiv:1301.2826 [nucl-th].[41] J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu,and U. Heinz, Phys. Rev. C , 024907 (2016),arXiv:1605.03954 [nucl-th].[42] J. E. Bernhard, J. S. Moreland, and S. A. Bass, NaturePhys. , 1113 (2019).[43] R. Critelli, J. Noronha, J. Noronha-Hostler, I. Portillo,C. Ratti, and R. Rougemont, Phys. Rev. D , 096026(2017), arXiv:1706.00455 [nucl-th].[44] J. M. Maldacena, Int. J. Theor. Phys. , 1113 (1999),arXiv:hep-th/9711200.[45] S. S. Gubser, I. R. Klebanov, and A. M. Polyakov, Phys.Lett. B , 105 (1998), arXiv:hep-th/9802109.[46] E. Witten, Adv. Theor. Math. Phys. , 253 (1998),arXiv:hep-th/9802150.[47] E. Witten, Adv. Theor. Math. Phys. , 505 (1998),arXiv:hep-th/9803131.[48] P. Kovtun, D. T. Son, and A. O. Starinets, Phys. Rev.Lett. , 111601 (2005), arXiv:hep-th/0405231.[49] S. S. Gubser, A. Nellore, S. S. Pufu, and F. D. Rocha,Phys. Rev. Lett. , 131601 (2008), arXiv:0804.1950[hep-th].[50] O. DeWolfe, S. S. Gubser, and C. Rosen, Phys. Rev. D , 086005 (2011), arXiv:1012.1864 [hep-th].[51] O. DeWolfe, S. S. Gubser, and C. Rosen, Phys. Rev. D , 126014 (2011), arXiv:1108.2029 [hep-th].[52] J. Casalderrey-Solana, H. Liu, D. Mateos, K. Rajagopal,and U. A. Wiedemann, Gauge/String Duality, Hot QCDand Heavy Ion Collisions (Cambridge University Press,2014) arXiv:1101.0618 [hep-th].[53] A. Ficnar, J. Noronha, and M. Gyulassy, Nucl. Phys. A , 372 (2011), arXiv:1012.0116 [hep-ph].[54] A. Ficnar, J. Noronha, and M. Gyulassy, J. Phys. G ,124176 (2011), arXiv:1106.6303 [hep-ph].[55] A. Ficnar, J. Noronha, and M. Gyulassy, Nucl. Phys. A , 252 (2013), arXiv:1208.0305 [hep-ph].[56] S. I. Finazzo and J. Noronha, Phys. Rev. D , 106008(2014), arXiv:1311.6675 [hep-th].[57] S. I. Finazzo, R. Rougemont, H. Marrochio, andJ. Noronha, JHEP , 051 (2015), arXiv:1412.2968 [hep-ph].[58] R. Rougemont, J. Noronha, and J. Noronha-Hostler,Phys. Rev. Lett. , 202301 (2015), arXiv:1507.06972[hep-ph].[59] R. Rougemont, A. Ficnar, S. Finazzo, and J. Noronha,JHEP , 102 (2016), arXiv:1507.06556 [hep-th]. [60] S. I. Finazzo and R. Rougemont, Phys. Rev. D , 034017(2016), arXiv:1510.03321 [hep-ph].[61] R. Rougemont, R. Critelli, J. Noronha-Hostler,J. Noronha, and C. Ratti, Phys. Rev. D , 014032(2017), arXiv:1704.05558 [hep-ph].[62] R. Rougemont, R. Critelli, and J. Noronha, Phys. Rev.D , 034028 (2018), arXiv:1804.00189 [hep-ph].[63] R. Rougemont, R. Critelli, and J. Noronha, Phys. Rev.D , 045013 (2016), arXiv:1505.07894 [hep-th].[64] S. I. Finazzo, R. Critelli, R. Rougemont, and J. Noronha,Phys. Rev. D , 054020 (2016), [Erratum: Phys.Rev.D96, 019903 (2017)], arXiv:1605.06061 [hep-ph].[65] R. Critelli, R. Rougemont, S. I. Finazzo, and J. Noronha,Phys. Rev. D , 125019 (2016), arXiv:1606.09484 [hep-ph].[66] R. Rougemont, Phys. Rev. D , 034009 (2020),arXiv:2002.06725 [hep-ph].[67] J. Knaute, R. Yaresko, and B. Kämpfer, Phys. Lett. B , 419 (2018), arXiv:1702.06731 [hep-ph].[68] Z. Li, Y. Chen, D. Li, and M. Huang, Chin. Phys. C ,013103 (2018), arXiv:1706.02238 [hep-ph].[69] J. de Boer, E. P. Verlinde, and H. L. Verlinde, JHEP ,003 (2000), arXiv:hep-th/9912012.[70] J. D. Bekenstein, Phys. Rev. D , 2333 (1973).[71] S. W. Hawking, Commun. Math. Phys. , 199 (1975),[Erratum: Commun.Math.Phys. 46, 206 (1976)].[72] R. Bellwied, S. Borsanyi, Z. Fodor, S. D. Katz, A. Pasztor,C. Ratti, and K. K. Szabo, Phys. Rev. D , 114505(2015), arXiv:1507.04627 [hep-lat].[73] S. Floerchinger and M. Martinez, Phys. Rev. C , 064906(2015), arXiv:1507.05569 [nucl-th].[74] S. Borsanyi, Z. Fodor, J. N. Guenther, R. Kara, S. D.Katz, P. Parotto, A. Pasztor, C. Ratti, and K. K. Szabo,Phys. Rev. Lett. , 052001 (2020), arXiv:2002.02821[hep-lat].[75] T. Dore, J. Noronha-Hostler, and E. McLaughlin, Phys.Rev. D , 074017 (2020), arXiv:2007.15083 [nucl-th].[76] V. Dexheimer, J. Noronha, J. Noronha-Hostler, C. Ratti,and N. Yunes, (2020), arXiv:2010.08834 [nucl-th].[77] P. Bedaque and A. W. Steiner, Phys. Rev. Lett. ,031103 (2015), arXiv:1408.5116 [nucl-th].[78] M. G. Alford, G. F. Burgio, S. Han, G. Taranto,and D. Zappalà, Phys. Rev. D , 083002 (2015),arXiv:1501.07902 [nucl-th].[79] I. F. Ranea-Sandoval, S. Han, M. G. Orsaria, G. A. Con-trera, F. Weber, and M. G. Alford, Phys. Rev. C ,045812 (2016), arXiv:1512.09183 [nucl-th].[80] I. Tews, J. Carlson, S. Gandolfi, and S. Reddy, Astrophys.J. , 149 (2018), arXiv:1801.01923 [nucl-th].[81] I. Tews, J. Margueron, and S. Reddy, Phys. Rev. C ,045804 (2018), arXiv:1804.02783 [nucl-th].[82] L. McLerran and S. Reddy, Phys. Rev. Lett. , 122701(2019), arXiv:1811.12503 [nucl-th].[83] P. Jakobus, A. Motornenko, R. O. Gomes, J. Stein-heimer, and H. Stoecker, Eur. Phys. J. C , 41 (2021),arXiv:2004.07026 [nucl-th].[84] E. Annala, T. Gorda, A. Kurkela, J. Nättilä, and A. Vuori-nen, Nature Phys. (2020), 10.1038/s41567-020-0914-9,arXiv:1903.09121 [astro-ph.HE].[85] T. Zhao and J. M. Lattimer, Phys. Rev. D , 023021(2020), arXiv:2004.08293 [astro-ph.HE].[86] H. Tan, J. Noronha-Hostler, and N. Yunes, Phys. Rev.Lett. , 261104 (2020), arXiv:2006.16296 [astro-ph.HE]. [87] J. Casalderrey-Solana, M. P. Heller, D. Mateos, andW. van der Schee, Phys. Rev. Lett. , 181601 (2013),arXiv:1305.4919 [hep-th].[88] J. Casalderrey-Solana, D. Mateos, W. van der Schee, andM. Triana, JHEP , 108 (2016), arXiv:1607.05273 [hep-th].[89] M. Attems, J. Casalderrey-Solana, D. Mateos, D. Santos-Oliván, C. F. Sopuerta, M. Triana, and M. Zilhão, JHEP , 154 (2017), arXiv:1703.09681 [hep-th].[90] M. Attems, Y. Bea, J. Casalderrey-Solana, D. Mateos,M. Triana, and M. Zilhão, Phys. Rev. Lett. , 261601 (2018), arXiv:1807.05175 [hep-th].[91] R. Critelli, R. Rougemont, and J. Noronha, Phys. Rev.D , 066004 (2019), arXiv:1805.00882 [hep-th].[92] M. Attems, Y. Bea, J. Casalderrey-Solana, D. Mateos,and M. Zilhão, JHEP , 106 (2020), arXiv:1905.12544[hep-th].[93] A. Folkestad, S. Grozdanov, K. Rajagopal, andW. van der Schee, JHEP12