The QCD phase diagram from analytic continuation

We present the crossover line between the quark gluon plasma and the hadron gas phases for small real chemical potentials. First we determine the effect of imaginary values of the chemical potential on the transition temperature using lattice QCD simulations. Then we use various formulas to perform an analytic continuation to real values of the baryo-chemical potential. Our data set maintains strangeness neutrality to match the conditions of heavy ion physics. The systematic errors are under control up to $\mu_B\approx 300$ MeV. For the curvature of the transition line we find that there is an approximate agreement between values from three different observables: the chiral susceptibility, chiral condensate and strange quark susceptibility. The continuum extrapolation is based on $N_t=$ 10, 12 and 16 lattices. By combining the analysis for these three observables we find, for the curvature, the value $\kappa = 0.0149 \pm 0.0021$.


Introduction
For heavy ion physics, the most important feature of the phase diagram of Quantum Chromodynamics (QCD) is the line that separates the hadron gas phase from the quark gluon plasma, and the conjectured critical end-point along this line separating cross-over from first order transition [1].
The qualitative form of the phase diagram was sketched four decades ago [2] as a consequence of Hagedorn's exponential spectrum of hadron masses [3]. The order of the transition at zero density has been determined much later, for Nature's selection of quark masses the two high temperature phases are connected through a cross-over [4]. In the absence of a real transition the crossover temperature can be determined but it is ambiguous [5]. Observables that are related to the spontaneous breaking of chiral symmetry (chiral condensate and its susceptibility) give a temperature around 155 MeV [5][6][7][8].
For the first few coefficients in Eq. (1) it is enough to study QCD at small µ B , for which there are several methods. κ can be and has been determined by calculating the µ B -derivative of the chiral condensate using only µ B = 0 ensembles [18,38]. However, the signal/noise ratio of higher µ B derivatives is suppressed with powers of the volume, making this approach impractical beyond µ 2 B order. Lattice calculations are perfectly feasible, though, with imaginary values of the chemical potential [10,19,39]. Setting µ B = iµ I B one avoids the sign problem and the transition line can be studied [40][41][42][43].
In this study we follow the imaginary-µ B approach and go beyond previous studies by a) performing a continuum approximation with lattices up to N t = 16; b) tuning µ S (µ B , T ) such that the strangeness neutrality condition is maintained; c) using several observables: chiral condensate, chiral susceptibility and strange susceptibility; d) comparing the Taylor and the imaginary-µ method for the strange susceptibility; e) calculating the systematic errors from scale setting, fit ranges, analytic formulas, etc. 1 This Letter first gives a brief account of the necessary lattice simulations at zero and finite temperatures. Then the method for setting strangeness neutrality is explained. Finally, we give a detailed description of the analysis and present the continuum results for the curvature.

Simulation setup
This study is part of the 2nd generation staggered thermodynamics program of the Wuppertal-Budapest collaboration [45]. We use a four times stout [46] smeared (ρ = 0.125) staggered fermion action with 2+1+1 flavors, i.e. dynamical up, down, strange and charm quarks. The gauge action uses the tree-level Symanzik improvement. The two light quarks are degenerate, their masses and the strange quark mass are tuned such that the physical pion and kaon mass over pion decay constant are reproduced for every lattice spacing. For the zero temperature runs we kept the volume large Lm π > 4 in the entire lattice spacing range of interest for this study: a = 0.2 . . . 0.063 fm. The charm mass was set to m c /m s = 11.85 [47]. The simulation parameters are detailed in Ref [45]. The overall scale was determined from f π . We used w 0 as an alternative scale setting for the analysis [48].
The chiral susceptibility as well as the chiral condensate require renormalization. The additive divergence is removed by subtracting the vacuum expectation value, the multiplicative divergence canceled by the same factor in the bare quark mass [5,6]. The renormalized condensate and its susceptibility are dimensionful quantities, we use the fourth power of the pion mass to form a dimensionless observable. We do not restrict the chiral susceptibility to the disconnected part. The third observable that we use to identify the µ-dependent transition temperature is the strange susceptibility: thanks to the exact quark number conservation it does not require renormalization.
At finite temperature, we have collected data at zero and at imaginary baryochemical potentials. The µ B = 0 data are used to perform a Taylor expansion on one of our studied observables, and also to obtain a "baseline" for the shifted transition temperatures at other µ B values. The zero density configurations are listed in Ref. [45].
The range of imaginary baryo-chemical potentials is limited by the Roberge-Weiss transition at µ B = iπT [49]. Below a limiting temperature T RW there is no transition as Im [µ B /T ] crosses π, but there is a first order transition above T RW , where the imaginary density is non-vanishing and flips sign at µ I B /T = π. The nature of the transition at T RW depends on the quark masses [50][51][52]. For intermediate masses the system at T = T RW and µ I B /T = π will be critical, and then in the entire range of smaller imaginary chemical potentials we will see a crossover in temperature. Our data suggests that, for physical masses, the latter scenario is realized, namely we are working with a cross-over for all used µ I B /T . We selected six imaginary chemical potential values: We have all six j values for our N t =8, 10 and 12 lattices and only j = 3 . . . 6 for N t =16. The reason for this is the following: j = 0 . . . 5 data are needed to determine the simulation parameters at finite imaginary µ such that the strangeness neutrality condition is fulfilled (see later). The continuum extrapolation for this analysis can be carried out using N t =8, 10 and 12 lattices. For the determination of the curvature κ of the phase diagram we also need the finest N t = 16 lattices. Since j =1 and 2 do not give a statistically very significant contribution to κ we decided not to have these two points in our most expensive N t = 16 ensembles. Therefore, in order to have the same setup for all lattice spacings, the κ determination is based on j = 3, 4 and 5. The j = 6 point is used to estimate higher order effects. This range to find the κ coefficient (µ I B /T 2) is narrower than in earlier studies (e.g. µ I B /T 2.36 in Ref. [42] and µ I B /T 2.6 in Ref. [43]). A broader range of chemical potentials has the advantage that the numerical B has a larger signal/noise ratio. However, more non-linearities appear in a broader range and the results are more prone to systematic errors as the singularity at µ I B ≈ πT is approached. This is the reason (to avoid unwanted systematic uncertainties) why we have taken a smaller µ I B range and we use other methods to increase the signal/noise ratio.
We performed simulations on 32 3 × 8, 40 3 × 10, 48 3 × 12 and 64 3 × 16 lattices, at sixteen temperatures in the temperature range 135. . . 210 MeV. We have generated between 10000-15000 Hybrid Monte Carlo updates, analyzing every 5th of them (every 10th for N t = 16). The configurations have been evaluated for up to fourth order generalized quark number susceptibilities [53] and for the chiral condensate and susceptibility. For µ B = 0 we have 5 . . . 10 times more statistics, this ensures a solid guidance to the fitting procedure.

Strangeness neutrality
The most popular representation of the QCD phase diagram is in the temperature vs. chemical potential plane. The baryo-chemical potential axis leaves room for various interpretations. Ref. [42] used the full baryo-chemical potential including the strange quarks, i.e. µ u = µ d = µ s = µ B /3. In Ref. [43] both However, neither of the recipes µ s = 0 or µ s = µ B /3 maps consistently to the situation that is realized in experiment. In heavy ion collisions, non-strange particles are colliding. Although ss are generated in the collision, the netstrangeness is zero. Therefore we want to tune the chemical potentials to such values which guarantee strangeness neutrality. The light chemical potentials are kept identical (µ u = µ d ), which ensures isospin symmetry also at finite µ B . This corresponds to an experimental situation where Z = 0.5A. Alternatively, one can achieve a different Z/A ratio corresponding to heavy nuclei by tuning µ u and µ d appropriately. This possibility will be discussed later. Requiring strangeness neutrality and fixing the value of Z/A (or alternatively the electric charge/baryon number ratio) uniquely determines all three quark chemical potentials as functions of µ B . For the isospin symmetric case µ u = µ d = µ B /3 so the only non-trivial task is to find the strange quark or strangeness chemical potential.
The strange quark chemical potential (µ s ) is related to the strangeness (µ S ) and baryo-chemical potential (µ B ) as µ s = µ B /3 − µ S . Then µ s = µ B /3 approximates strangeness neutrality at low temperature, and µ s = 0 at high temperature. In this work we determine the strangeness neutral trajectory µ S (µ B , T ) from lattice simulations.
It is relatively straightforward to perform Taylor expansions from µ B = 0 on the trajectory that respects strangeness neutrality. For the equation of state [54,55] and for fluctuations relevant for calculating freeze-out parameters in heavy ion collisions [53,56] this procedure is already standard.
For actual simulations at finite µ I B the strange chemical potential has to be fine tuned for every temperature, baryo-chemical potential and lattice spacing. We solved this challenge by solving the differential equation discretized in µ B with the trivial initial condition ∂ log Z/∂µ S = 0 at µ I B = 0. This equation simply states that the µ I B derivative of strangeness is zero. Using the 2nd order explicit Runge-Kutta scheme, we determine µ I S (µ I B ) using the prescription: ) accurate only, but as an additional correction, we do a small extrapolation for both terms on the RHS of Eq. (4) after each simulation so that the remaining strangeness neutrality violation is not propagated to the next step. For this extrapolation, we need higher order fluctuations [53]. This combination is O(∆µ I B 3 ) accurate in the complete µ I B range. The resulting µ S (µ B , T ) function is interpolated in T and extrapolated in 1/N 2 t and the resulting smooth function is used to start the simulations at µ I B + ∆µ I B . In Fig. 1 we show the resulting strangeness chemical potential.

Analysis details
We calculate the curvature of the phase diagram from three observables. We calculate statistical and systematic errors for all three.
1) Our first observable is the chiral susceptibility χψ ψ /m 4 π . As discussed previously, it requires additive and multiplicative renormalization. For details on this procedure see Ref. [7]. The chiral susceptibility forms a peak at the transition temperature. With increasing imaginary chemical potential this peak is shifted towards higher temperatures, approximately maintaining its height and width. For other normalizations (e.g. χψ ψ /T 4 ) the shape of the function changes more significantly while varying the chemical potential.
We fit χψ ψ (µ I B , T )/m 4 π in a global fit function where for each µ B a different width, height and peak position is allowed, but the other parameters that describe the peak shape are µ B -independent. We use two different modifications to the Lorentzian peak form: and The µ dependent parameters A(µ), W (µ) and T c (µ) describe the change in the height, width and the position of the curve as µ increases. For the zero temperature data which are required for renormalization we use two different interpolations in the inverse gauge coupling: a 6th order polynomial and a simple rational function. We have two options for the scale setting using f π or w 0 and we apply three possible fit windows to select the transition range. In order not to interfere with the shifted temperature dependence the fit windows constrain the value of the susceptibility, not the temperature. The effect of the µ dependent parameters is a shift in T , and a rescaling in T and χ. Applying the inverse transformation to the finite µ I B data points all of them should collapse on the µ B = 0 curve. This is demonstrated in Fig. 2. The advantage of this procedure is that the µ independent parameters can be fitted using the the high-statistics runs at µ B = 0 and the non-vanishing µ I B runs are needed to determine the relative position and rescaling compared to this more complicated functional form. This allows the precise determination of ∆T c (µ I B ) with an error below 0.25 MeV, while T c (µ B ) itself has an error of several MeV. We extract κ by a linear fit of ∆T c vs. µ 2 B in the range 1.2 µ I B /T 2, and extrapolate κ to the continuum. Since the continuum extrapolation of κ had a large χ 2 when all four lattices were used we included only N t = 10, 12 and 16 in the final result, resulting in a good χ 2 for all analyses. In an alternative analysis we made a combined continuum and µ 2 B fit again using only the finest three lattice spacing, and found acceptable χ 2 values again.
2) The chiral condensate ψ ψ r = m q (d log Z/dm q )/m 4 π is a remnant order parameter of the chiral transition. Its inflection point (though it is hard to locate in a finite precision data set) is very close to the peak position of χψ ψ /m 4 π . At finite µ I B the temperature dependence of ψ ψ r is shifted and very slightly stretched.
We find that the data at µ B = 0 (see Ref. [7]) can be very accurately described by simple fit functions. The µ dependence in this case is well described by just two µ B -dependent parameters describing a shifting and rescaling of the renormalized condensate. We use the following parameterizations: and ψ ψ r (µ, Similarly to the chiral susceptibility, we use two possible zero temperature interpolations (6th and 7th order polynomials of the inverse gauge coupling), two scale settings, four fit windows. κ is obtained either from a combined µ 2 B and continuum fit, or separately.
3) The analysis of the strange susceptibility χ S 2 goes along the lines of the chiral condensate. A simplification here is the absence of renormalization. Since this quantity is the most sensitive to the actual value of strangeness, before the analysis we correct for the inaccuracies of the strangeness neutrality condition using the higher µ S fluctuations. Although its inflection point does not have to agree with that of the chiral condensate, we find that the shifting effect of the chemical potential is very similar.
For all three quantities we make a histogram of the results from all analyses. For the chiral susceptibility we have two T > 0 fit forms, two T = 0 interpolations, two scale settings, three fit windows and either separate or combined κ extraction and continuum limit. This results in 2 · 2 · 2 · 3 · 2 = 48 analyses. For the condensate we have the same choices but with four fit windows resulting in 64 analyses. For the strange susceptibility there is no renormalization, thus no T = 0 interpolation is needed which leads to 32 analyses. The central 68% of the histograms estimates our systematic error. The statistical error is obtained from 1000 bootstrap samples. The two errors are of similar magnitude and they are added in quadrature resulting in our final uncertainties.
We summarize our results for the curvature in Table 1.
The histograms of the three quantities can be joined into a single one leading to our combined result based on our three observables with strangeness neutrality: κ = 0.0149 ± 0.0021 .
We also consider the curvature for the case when not only the strangeness neutrality, but also proper charge/baryon density ratio is reproduced (for lead and gold ions: Z ≈ 0.4A). We achieve this by Taylor-extrapolating the strange susceptibility for every finite µ B ensemble to leading order, and fitting as before. We conclude that the difference between Z = 0.4A and Z = 0.5A phase diagrams is negligible for small µ B .
For small enough imaginary chemical potential the analytical and the Taylor method have to give the same curvature at every lattice resolution. In the Taylor method one expands the observables in µ B , the leading coefficients are calculated from µ B = 0 simulations and then used to extrapolate to finite µ B . Fig. 3 shows how this expansion compares to the direct simulations for our j = 5 chemical potential which is the largest one used to extract κ. A comparison in <S>=0 direct simulaton µ s =µ u =µ d extrapolated from µ B =0 µ s =µ u =µ d extrapolated from <S>=0, µ B >0 Figure 3: Comparison of the strange susceptibility obtained from direct simulation and extrapolation for µ I B /T = 5π/8. The blue circles and squares correspond to the strangeness neutral case obtained via extrapolation from µ = 0 and direct simulations, respectively. The green triangles show the full baryo-chemical potential case obtained via extrapolation from µ = 0. There are no direct simulations in this case but one can extrapolate from the strangeness neutral direct point (green dots). As a reference the µ = 0 data are also shown (red crosses).
the case of full baryo-chemical potential is also shown. The extrapolated data are then fitted for κ as if they were simulated at finite µ I B . At N t = 10 we find κ = 0.0131(9) from the direct simulations and κ = 0.0115(10) from the Taylor expansion. The agreement indicates that we are still in the linear regime and the extraction of κ using j = 3, 4, 5 is safe.
Finally we estimate the systematics of the extrapolation to real µ B . We include the j = 6 data points into the analysis and allow for non-linear T c (µ 2 B ) fits. We consider fitting T c (µ 2 B )/T c with the functions 1 + ax, 1 + ax + bx 2 , (1 + ax)/(1 + bx) and (1 + ax + bx) −1 with x = µ 2 B /T 2 . All these functions are analytic in x and they represent various analytic continuations of the T c (µ I B ) imaginary µ phase diagram. The difference between these ansatzes provides a systematic uncertainty for the real µ phase diagram.
Our main results are depicted in Fig. 4. Since the curvature from both the strange susceptibility and from the chiral condensate/susceptibility are consistent with each other we show only one curve. The curvature from the chiral condensate is our most precise result, therefore we present the transition line coming from this observable. The corresponding transition temperature at µ=0 is at 157 MeV. At intermediate real µ B we observe a significant rise in the uncertainty due to the statistical error on the non-linear µ 2 B -dependence and the ambiguity of the analytic ansatz. This sets the range of validity for this study.
The present result indicates a stronger curvature than the one presented in Ref. [38]. There are, however a couple differences between the definitions/ approaches of the curvature of the present analysis and Ref. [38]. Note that the transition is a smooth cross-over, thus different definitions obviously lead to different results. The phase diagram based on the µ-dependent Tc from the chiral condensate, analytically continued from imaginary chemical potential. The blue band indicates the width of the transition. The shaded black region shows the transition line obtained from the chiral condensate. The widening around 300 MeV is coming from the uncertainty of the curvature and from the contribution of higher order terms, thus the application range of the results is restricted for smaller µ values. For completeness, on the right panel we also show some selected non-lattice results: the Dyson-Schwinger result of Ref. [37] and the freeze-out data of Refs. [57][58][59][60][61][62][63].
a. In Ref. [38] we used a vanishing strangeness chemical potential. In the present analysis we use instead vanishing strange density. The reason for this change is to be as close to the experimental situation as possible. In heavy ion collisions the net strangeness is zero.
b. It is emphasized in the discussion of Figure 5 of [38] that only statistical uncertainties were provided. The present analysis estimates systematic uncertainties coming from various aspects of the analysis as discussed earlier. These are comparable to or in some cases even larger than the statistical uncertainties. A similar assumption on the systematics of Ref. [38] would make the tension between the results much weaker.

Acknowledgments
The authors thank G. Endrodi for his valuable comments and suggestions. This project was funded by the DFG grant SFB/TR55. S. D. Katz is funded by the "Lendület" program of the Hungarian Academy of Sciences ((LP2012-44/2012). The work of R. Bellwied is supported through DOE grant DEFG02-07ER41521. C. Ratti is supported by the National Science Foundation through grant number NSF PHY-1513864. An award of computer time was provided by the INCITE program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility sup-ported under Contract DE-AC02-06CH11357. This research also used resources of the PRACE Research Infrastructure resource JUQUEEN at FZ-Jülich, Germany; JUQUEEN as large scale project of the Gauss Centre for Supercomputing (GCS); the QPACE machine supported by the Deutsche Forschungsgesellschaft through the research program SFB TR55; and the GPU cluster at the Wuppertal University.