The QCD crossover at finite chemical potential from lattice simulations

We provide the most accurate results for the QCD transition line so far. We optimize the definition of the crossover temperature $T_c$, allowing for its very precise determination, and extrapolate from imaginary chemical potential up to real $\mu_B \approx 300$ MeV. The definition of $T_c$ adopted in this work is based on the observation that the chiral susceptibility as a function of the condensate is an almost universal curve at zero and imaganiary $\mu_B$. We obtain the parameters $\kappa_2=0.0153(18)$ and $\kappa_4=0.00032(67)$ as a continuum extrapolation based on $N_t=10,12$ and $16$ lattices with physical quark masses. We also extrapolate the peak value of the chiral susceptibility and the width of the chiral transition along the crossover line. In fact, both of these are consistent with a constant function of $\mu_B$. We see no sign of criticality in the explored range.

Introduction-One of the most important open problems in the study of Quantum Chromodynamics (QCD) at finite temperature and density is the determination of the phase diagram of the theory in the temperature (T )-baryo-chemical potential (µ B ) plane. It is now established by first principle lattice QCD calculations that the transition at µ B = 0 is a smooth crossover [1,2] for physical quark masses. Due to the lack of a real phase transition, the crossover temperature is of course ambiguous, since different definitions can lead to different values for it. Observables related to chiral symmetry (i.e. the chiral condensate and its susceptibility) yield a transition temperature around 155 − 160 MeV [3][4][5][6].
Extending our knowledge to the µ B > 0 part of the phase diagram turns out to be very challenging, due to the notorious sign problem. Since this makes direct simulation at finite µ B impossible, the state-of-the-art for finite density QCD on fine lattices is to use one of two extrapolation methods. The first method is the direct calculation of Taylor coefficients [7][8][9][10][11][12][13][14][15][16][17] using simulations at µ B = 0, while the second is to use simulations at imaginary chemical potentials (µ 2 B < 0) where the sign problem is absent, and later perform an extrapolation of different quantities to a real chemical potential (µ 2 B > 0) [18][19][20][21][22][23][24][25][26][27][28][29][30][31]. It is often conjectured that in the (T, µ B ) plane the crossover line, departing from (T c , µ B = 0), eventually turns into a first-order transition line. The point (T CEP , µ CEP ) separating the crossover and the first-order transitions is known as the critical endpoint (CEP), where the transition is expected to be of second order. Though there have been attempts in extracting information about the location of the supposed CEP from lattice simulations [15,26,[32][33][34][35][36][37], these attempts face great difficulties, as extrapolation-type methods have the property that they give reliable results mostly in the immediate vicinity of µ B = 0.
In this letter, we address the problem of calculating the Taylor coefficients of the crossover temperature around µ B = 0, parametrized as: (1) along the phenomenologically relevant strangeness neutrality line. In this work we improve the uncertainty on κ 4 available in the literature [16] by a factor of 6, giving a state-of-the-art determination of the cross-over line in the (T, µ B ) plane. In particular, as we will show, at chemical potentials µ B > 200 MeV the error on the T c extrapolation is dominated by the sub-leading coefficients e.g. κ 4 . The coefficients κ 2 and κ 4 can be calculated with either one of the standard extrapolation methods. A direct evaluation of the µ B derivatives from µ B = 0 ensembles was used in Refs. [38,39]. The current state-of-the art using the µ B = 0 simulation method is Ref. [16], which includes the first continuum extrapolated results for κ 4 . Here we will employ an analytical continuation from imaginary µ B instead, and use lattices as fine as N t = 16. This is motivated by the fact that the signal/noise ratio of higher µ B derivatives is suppressed with powers of the lattice volume, therefore the calculation of higher order derivatives requires very high statistics. Determinations of κ 2 using the imaginary µ B method with continuum extrapolation include Refs. [24,25]. Finally, in Ref. [30] the two methods were compared with a careful check of the systematics, and a very good agreement was found for the coefficient κ 2 .
We also study the strength of the crossover by extrapolating the width of the transition and the value of the chiral susceptibility at the transition to real µ B in the continuum limit. While one always has to be careful not to over-interpret results from extrapolations, we currently do not see any sign of criticality up to µ B ≈ 300 MeV, as the crossover transition does not get narrower or stronger in this region.
On chiral observables in the transition region.-For the lattice simulations we use 4-stout improved staggered fermions with an aspect ratio of LT = 4 and temporal lattice sizes of N t = 10, 12, 16. The details of the simulation setup can be found in the supplemental material.
The main observables in this study are the renormalized dimensionless chiral condensate and susceptibility, respectively defined as: where we assumed isospin symmetry, i.e. m u = m d = m ud . In the above equations, the subscripts T, 0 indicate values at finite-and zero-temperature, respectively. In the following, ψ ψ and χ are always shown after applying the correction to satisfy n s = 0 with zero statistical error (see the supplemental material for details). The peak height of the susceptibility is an indicator for the strength of the transition, while the peak position in temperature serves as a definition for the chiral cross-over temperature. It was pointed out in Refs. [3,4] that different normalizations of the susceptibility, such as using 1/f 4 π or 1/T 4 to define χ in Eq. (2) can shift the peak position by 11 MeV. This difference could be considered as a measure for the broadness of the chiral transition.
Our normalization choice in Eq. (2) was motivated by two observations, shown in Fig. 1 and explained below. These observations (together with the improved statistics and the more accurate tuning of µ S (µ B ) to n S = 0) allow a very precise determination of T c as a function of imaginary chemical potential, which in turn allows a precise determination of the parameters κ 2 and κ 4 . We explored the chiral condensate and susceptibility in a broad range of imaginary baryo-chemical potential. In all panels of Fig. 1, the black curves correspond to µ B = 0. In the left and middle panel we show the chiral condensate and susceptibility as functions of the temperature. By construction, our renormalized condensate is zero at T = 0 and positive at high temperature, because of the explicit vacuum subtraction and the overall negative sign in Eq. (2). In both panels, one can observe the shifting of the transition towards higher temperatures when an imaginary chemical potential is introduced. In the right panel we show the susceptibility as a function of the condensate. Here we converted the statistical error on the condensate into an additional error on the susceptibility, by solving for ψ ψ (T ) = const. and substituting the resulting T into χ(T ) (also taking the correlation of the statistical errors into account). Our first observation on the right panel of Fig. 1 is that the form of the χ( ψ ψ ) curve is simpler than that of χ(T ): a low (e.g. third or fourth) order polynomial can fit the entire transition range with an excellent fit quality. The second observation is that there is virtually no chemical potential dependence in the χ( ψ ψ ) function. This way the susceptibility can be modeled as a low order polynomial of two variables, ψ ψ andμ = µ B /T . Had we used a different normalization for the susceptibility, e.g. χ(T )f 4 π /T 4 as we did in Ref. [5], the peak height would be strongly µ B -dependent and the collapse of the χ( ψ ψ ) curves at different (imaginary) chemical potentials would not happen. The transition line and its analytical continuation.-Keeping the previous observations in mind, one can perform a precise determination of T c , as defined by the peak of χ in Eq. (2) for various values of the imaginary chemical potential. T c (µ 2 B ) can then be fitted for the coefficients κ 2 and κ 4 . This requires the following steps: i) Determine the renormalized condensate ψ ψ and susceptibility χ in a two-dimensional parameter scan in T and Imµ B using lattice simulations. Use these to obtain the susceptibility as a function of the condensate. ii) Search for the peak of χ( ψ ψ ) through a low-order polynomial fit for each N t and Imµ B obtaining ψ ψ c (N t , Imµ B ). iii) Use an interpolation of ψ ψ (T ) to convert the ψ ψ c to T c for each Imµ B /T . iv) Perform a global fit of T c (N t , Imµ B /T c ) to determine the coefficients κ 2 and κ 4 for 1/N 2 t = 0. For this step we use various functions -all containing an independent κ 6 -with coefficients depending linearly on 1/N 2 t . The choice of the fit functions is motivated by the mock data analysis presented in the supplemental material. The total systematic error comes from a pool of 256 analyses: in step i) we have two choices for the scale setting, two choices for the renormalization of ψ ψ and two for the renormalization of χ; in step ii) we use two choices for the fit function used to obtain the maximum of Crossover line from the lattice compared with a prediction from truncated Dyson-Schwinger equations [40] and some estimates of the chemical freezeout parameters in heavy ion collisions [41][42][43][44][45]. Note that the width of the green band is not a representation of the width of the crossover region, it depicts the statistical and systematic errors achievable with the particular definition of the crossover temperature Tc adopted in this work. Note also that the definition of the crossover temperature adopted in Ref. [40] is different from the one used in this work.
χ( ψ ψ ) and two choices for the fit range; in step iii) we use two different interpolations of ψ ψ (T ); in step iv) we use two global fit functions and two choices for the range in Imµ B /T . This leads to a total of 2 8 = 256 ways to analyze our lattice data. These results are combined with a uniform weight. More details on the analyses, the fit qualities and the error estimates can be found in the supplemental material. We finally obtain: We stress that the uncertainties on these two quantities are correlated. We put these results in the context of previous lattice studies in Fig. 2. The extrapolated value of T c (µ B ) is shown in Fig. 3 (green band). Note that the errors on κ 2 and κ 4 are dominated by the statistical errors, as shown in the detailed discussion of the systematic error estimate in the supplemental material. Since Ref. [25] we have more than doubled the statistics and introduced a more precise analysis. The overall  error on κ 2 has reduced slightly. The main result is the extraction of κ 4 . It appears to be a generic feature of deducing Taylor-coefficients from polynomial fits: the increased precision on the input data leads to a sensitivity to a higher order coefficient first, and only later to a reduction of the error of both coefficients. This feature is also clearly seen in the mock data analysis in the supplemental material.
In Fig. 3 we also show the comparison to the leading order Taylor expansion result (using only κ 2 ) and the next to leading order result (using κ 2 and κ 4 ). The latter is very close to our full result (for µ B < 300 MeV), while the leading order result has a much smaller uncertainty. Clearly, κ 2 is precise enough. At intermediate µ B the bottleneck for the precision of T c (µ B ) is the error on κ 4 . We also fitted κ 6 , this turned out to be small enough to be irrelevant for µ B < 300 MeV.
Extrapolation of the transition width and strength.-A natural definition of the width of the susceptibility peak is given by its second derivative at T c as . Unfortunately, evaluating this quantity is numerically difficult, so we introduce a simple width parameter σ as a proxy for ∆T via: with ψ ψ c = 0.285 and ∆ ψ ψ = 0.14. The choice of the range in ψ ψ is such that it is consistent with a linear behavior within our errorbars, meaning that the ratio ∆ ψ ψ /σ can be used as a proxy for d dT ψ ψ | T =Tc as well. The exact range in ψ ψ is chosen such that σ coincides with ∆T at zero and imaginary µ B . A more detailed discussion of the width parameter can be found in the supplemental material.
We can define, for any value of ψ ψ , the temperature function In the left panel of Fig. 4 we show these contours in the (T, µ 2 B ) plane for a selection of ψ ψ values. We show continuum extrapolations and include the systematic errors. For this analysis the chiral susceptibility plays no role, since we use the same interpolations of ψ ψ (T ) as in step iii) of the analysis of the transition line, and two-dimensional fits (continuum and inμ 2 B ) analogous to those in step iv). We conclude that the half width of the transition -shown in the upper panel of Fig. 5 -is consistent with a constant up to µ B ≈ 300 MeV within the uncertainty from the extrapolation (we note that 50% uncertainty is reached at µ B ≈ 280 MeV). We also show ψ ψ as a function of T for several fixed values of real µ B /T in the right panel of Fig. 4, as extrapolated using the contours in the left panel of the same Figure. Finally, as a proxy for the strength of the crossover, we study the value of the chiral susceptibility at the crossover temperature. We get this for each Imµ B and N t as a byproduct of steps i-ii) of the analysis for κ 2 and Top: Half width σ of the transition defined in Eq. (4) using the temperature difference of the contours ψ ψ = 0.31 and ψ ψ = 0.19. In the insert we show a plot of the χ( ψ ψ ) peak, where the shaded region corresponds to ψ ψ c ±∆ ψ ψ /2. Both are extrapolated to real µB. Bottom: Result of a µB-by-µB analysis for the value of the chiral susceptibility at the crossover temperature after continuum extrapolation and including the systematic errors for LTc = 4. The green band shows a linear extrapolation inμ 2 B .
determination of the crossover line, we used the experimentally relevant µ S (µ B ) tuned to keep n S = 0. Based on the observation that the chiral susceptibility as a function of the condensate is a rather simple function, only weakly dependent on the imaginary chemical potential, we were able to obtain the transition temperature as a function of the imaginary chemical potential to very high accuracy. These pure lattice results can be used for further model building, and are summarized in the supplemental material. The high precision data at imaginary µ B in turn allowed us to fit the µ 2 B and µ 4 B Taylor coefficients of the crossover temperature, κ 2 and κ 4 . In particular, while our determination of κ 4 is still consistent with zero, the error is 6 times smaller than the one previ-ously available in the literature, and therefore represents the state-of-the-art in the study of the phase diagram in the (T, µ B ) plane with current lattice techniques. As a byproduct, we also obtain the most precise value for the central temperature of the crossover at µ B = 0 so far, as well as the width of the transition: The present definition was actually included in our earlier list of observables with T c (LT = 3, µ B = 0) = 157 (3)(3) MeV [4]. Recently the HotQCD collaboration has published T c (LT = 4, µ B = 0) = 156.5 ± 1.5 MeV for the peak of the chiral susceptibility in Ref. [16]. Note that all results in this letter were obtained at a fixed aspect ratio of LT = 4. Though this result is consistenct with earlier works, also with LT = 3, finite volume corrections can potentially be relevant with the achieved precision. We also did not take into account isospin breaking and QED effects.
We also studied the strength of the phase transition as a function of µ B by extrapolating our proxy for the transition width and the peak of the chiral susceptibility from imaginary chemical potentials. Even though one has to be careful with extrapolations, we see no sign of the transition getting stronger up to µ B ≈ 300 MeV.

Details of the lattice setup
Throughout this work we use tree-level Symanzik improvement for the gauge action and four levels of stout smearing in the staggered fermion action. We use three lattice spacings in this work, which are given in terms of the (Euclidean) temporal resolution of the isotropic lattices: a = 1/(N t T ), with N t = 10, 12 and 16.
The parameters of our discretization are tuned in such a way that the measured pion and kaon masses are equal to 135 MeV and 495 MeV, respectively, if we use the pion decay constant f π = 130.41 MeV for scale setting [46]. As an alternative scale setting we use the Wilson flow scale w 0 [47]. Continuum extrapolated results do not depend on the choice of the scale setting procedure, while results at finite lattice spacings do. For example, near the transition on our 40 3 ×10 lattice, w 0 f π differs by 2.5% from the continuum value, thus all hadrons appear 2.5% lighter with the w 0 scale setting. On our finer lattices, 48 3 × 12 and 64 3 × 16 this difference reduces to 2% and 1%, respectively, and vanishes in the continuum extrapolated results. Should there be any small deviation between the two f π and w 0 based continuum extrapolations, we consider the difference as part of the systematic error. The actual simulation parameters as well as the bare parameters are given in Ref. [13]. This action has already been used to calculate the equation of state a µ B = 0 [48], fluctuations of conserved charges [13] as well as the cross-correlators [49] and fugacity expansion coefficients [28,29].
For the purpose of renormalization in Eq. (2) in the main text, we calculate the vacuum condensate and susceptibility, on large lattices with Lm π ≈ 4 and N t /N x > ∼ 1.3, for 9 values of the gauge coupling, ranging over β = 3.55 − 4.0126, corresponding to a = 0.19 − 0.063 fm. In this range the bare condensate varies by an order of magnitude. We interpolate between the simulated gauge couplings by fitting the natural logarithm of the T = 0 condensate with a polynomial of order four or five -with χ 2 /n dof = 3.52/4 and 2.6/2 respectively. The bare susceptibility at T = 0, as well as its logarithm, were fitted with a second order polynomial -with χ 2 /n ndof = 8.4/6 and 9.3/6 respectively. The two interpolations of the condensate and the two interpolations of the susceptibility were varied independently in the systematic error analysis.
In this work we generate lattice ensembles for fixed ratios of zero and imaginary µ B /T . With the notation µ B = µ B /T we select eight values where we perform a temperature scan with our three lattice spacings.   Several of these chemical potentials were already used in our earlier work on the transition line [25], where a key point was the use of strangeness neutrality. This means that for each (imaginary) baryo-chemical potential and temperature, the strangeness chemical potential was tuned such that the expectation value of the strangeness density vanishes. In this work, too, a non-zero value of µ B always implies a matching µ S parameter with n S = 0. Thus, our extrapolations in µ B also extrapolate in µ S .
In Table I we list the number of analyzed configurations in the finite-temperature ensembles for the different lattices used in this work.

Mock data analysis
In the main text we calculated the transition temperature with sub-percent precision. This precision is needed when we are attempting to calculate a numerical derivative of T c with respect to the chemical potential. The relevance of the sub-percent errors in the T c determination can be highlighted in a mock analysis that we present below.
We took two model scenarios for the crossover line and generated mock data at imaginary chemical potential. Each mock data set was generated from a fourth order polynomial inμ 2 B with coefficients κ mock 2 , . . . , κ mock 8 . For the first model, we fitted the cross-over line at real µ B from Ref. [40], and obtained the following values for the coefficients: κ mock 2 = 0.0259463, κ mock 4 = −0.0013438, κ mock 6 = 0.000053 and κ mock 8 = −0.00000094. We then generated seven sets of mock data for each value of imaginary chemical potential listed in Eq. (7), assigning to each set a fixed relative error, ranging from 0.1% to 1% (see the left panel of Fig. 6). We fitted each set of mock data T c (μ B )/T c (0) with four fit functions: a 2nd order polynomial, a 1/1 rational function (1 + ax)/(1 + bx), a third order polynomial and its reciprocal. κ 2 and κ 4 were extracted as the leading Taylor coefficients of the fit functions. In the case of the two-parameter fits we dropped the largest imaginary chemical potential, so that all fits had four degrees of freedom.
We repeated the whole procedure with a second model, constructed with a κ mock  Fig. 6).
The purpose of the analysis was to determine whether the fits can reproduce the κ 2 and κ 4 values from the mock data. We found that the two models behave quite differently, even for large errors. In the first one, all the fit functions reproduce the value of κ 4 . In the second model, however, the second order fits yield κ 4 ≈ −0.001, while its actual value is κ mock 4 = 0.001. The third order fits perform much better. For the second model, we find κ 4 = 0.00123 (18) and κ 4 = 0.00099 (16) from the per-mill precision mock data. However, with more precise data or a substantially higher κ 8 value, the third order fits would fail as well. In order for κ 8μ 8 B to be negligible relative to κ 6μ 6 B in a fit whereμ 2 B ∈ [−π 2 , 0], the coefficients κ 8 and κ 6 must be separated by more than an order of magnitude.
Our real lattice data have a relative precision on T c (µ B ) near 2 per-mills. Thus, in our analysis it is essential to take at least third order polynomials to extract the cross-over line up to orderμ 4 B . Note also that, by reducing the relative error on T c (μ B ) from 0.5% to 0.2%, we significantly improved the 2nd order result on κ 2 . However, when switching to 3rd order fit functions the error bars become as large as before. What was gained with the increased precision is an access to κ 4 . This feature of deducing Taylor coefficients from polynomial fits is also pointed out in the main text.

Error estimate
In most results of this work, the statistical errors dominate. They were calculated through the standard jackknife procedure with 48 bins. For each lattice size, temperature and chemical potential we produced several (uncorrelated) Monte Carlo streams. In most cases we calculated the chiral observables after every five doublelength trajectories. The sequence of measurements in the Monte Carlo streams with equal bare parameters are concatenated (after skipping 1000 trajectories for thermalization). The 48 bins are formed by splitting up this sequence evenly.
Systematic errors are introduced every time we make an ambiguous choice in the analysis. Such choices include the scale setting variable, which can be either f π or w 0 , the interpolating function we fit on the zero temperature data and the one we use to fit the renormalized quantities. At several points in the analysis we must make various decisions -e.g. how many data points to include into a fit or what the specific interpolating function should be. Each time we pursue two or more choices, this results in splitting the analysis. The final result is then calculated in many slightly different versions. The width of their distribution yields the systematic error.
In this paper we give combined errors. For this purpose, we consider the cumulative distribution function (CDF) for observable x where we built a weighted sum of a Gaussian CDF corresponding to analysis j with mean m j and variance σ j and j w j = 1. We use uniform weights with the aforementioned cut in the fit quality Q. The upper (lower) edge of the error bar is defined as the value of the observable corresponding to the 84% (16%) level of the cumulative distribution function.
The different sources of systematic error in steps iiv) of the analysis as described in the main text can be summarized as: i) 2 choices of scale setting, 2 choices of the fit used for the renormalization of ψ ψ and 2 choices of the fit used for the renormalization of χ as discussed in the first supplement ii) The function χ( ψ ψ ) is fitted with a quadratic or cubic polynomial on the range given by χ > χ cut . We choose two χ cut values for each polynomial order: we picked the largest and smallest χ cut values that are endorsed by a Kolmogorov-Smirnov test (see Fig. 7). iii) Two different cubic spline interpolations of ψ ψ (T ). The spline node points are selected to be at every other simulation temperature (the separation of our simulation temperatures is 5 MeV or less, see Fig. 1 (left) in the main text). Here we use either every even or every odd simulation point as node points. iv) The global fit ansatz is either T c ( for the first fit and κ 2 = a, κ 4 = b − a 2 for the second. As a further source of systematic error, we either include the largest imaginary chemical potential in this fit, or we drop it. All these choices lead to the 2×2×2×2×2×2×2×2 = 256 analyses mentioned in the main text. The Q values and the obtained κ 2 and κ 4 values for each fit are shown in Fig. 8. The fits are well spread in the Q = [0, 1] interval, and there is no systematic dependence of the coefficients on the order of the used polynomial. There are no too bad or too good fits, either, these were sorted out by the Kolmogorov-Smirnoff tests of Fig. 7. This motivates us to combine the various analyses with a uniform weight. The error is mainly statistical. Among the statistical effects the precision of the peak position in χ( ψ ψ ) is the bottleneck in this computation.
Notice our result on the κ 6 coefficient in the bottom of Fig. 8. We cannot control the systematics of κ 6 , nevertheless, we find it takes a small value, consistent with zero, and is stable for systematic effects. The order of magnitude of κ 6 means that its effect on T c is not significant below µ B ≈ 300 MeV.

Strangeness neutrality
Throughout this work we use a strangeness chemical potential µ S , which is always tuned such that for each simulated µ B , µ S pair we have strangeness density n S = 0. This is the strangeness neutrality condition. We have already calculated the µ S (µ B ) dependence in Ref. [25], where we presented our first continuum extrapolated result on the curvature κ 2 of the transition line.
As an improvement on our analysis from Ref. [25], we correct the renormalized chiral condensate and suscepti- Illustrations of the Kolmogorov-Smirnov test for our χcut fit range selection. For second (left), third (middle) and fourth (right) order fits we selected two fit ranges for the determination of the peak position of the chiral susceptibility. With the three lattice spacing and eight imaginary chemical potentials we had 24 independent data sets for the test. The zig-zag curves show the cumulative distribution function for of the Q value of the maximum fit, the almond shaped band shows the one-sigma expectation of the curves. With higher χcut less points are included, the polynomial model might overfit the data, the corresponding curves go then below the almond. Bad fits, on the other hand, are above the almond, this typically happens if χcut is too low, the fit is too inclusive. In the case of the fourth order polynom always slightly overfit the data, thus we do not fit quartic polynomials. bility for the systematic or statistical deviations from the strangeness neutrality condition. To this end, we calculate the derivatives of the condensate and the susceptibility with respect to the strangeness chemical potential. Let us first write the observables for N f flavors at fixed chemical potential: where: Here H stands for the Hermitian matrix H = M † l M l where M l is the light quark matrix containing the light quark mass m l and the light quark chemical potential µ l = 1 3 µ B . At the same time we introduce the strange quark matrix M s containing the strange quark mass m s and the strange quark chemical potential µ s = 1 3 µ B −µ S . The lattice observable n s is the strange quark density, N f = 2 is the number of light flavors. The factor 1/4 in front of the trace in Eq. (14) is due to the use of staggered quarks.
We then calculate: with We can correct for the small deviations from n s = 0 that we find after performing the simulations, by calculating the change ∆µ s which would restore strangeness neutrality to leading order in a Taylor expansion on a jackknife-by-jackknife basis. We get the Taylor coefficients by averaging Eqs. (14) and (18) We calculate the correction to the condensate and susceptibility to leading order in ∆µ s where ψ ψ and χ are the renormalized observables from Eq. (2) of the main text. Thus, for our corrected ensembles strangeness neutrality is achieved with zero statistical error. The n S = 0 setup may have a systematic error, though, if the leading order expansion in ∆µ s was not satisfactory. However, the ensembles were already tuned to fulfill |n S /N B | < 0.05 even without correction. The correction ∆µ S /µ S < 0.1 introduced here resulted in a relative shift of ∼ 10 −3 or less for both ψ ψ and χ.
The correction that we calculated in this elaborate analysis is found to be smaller than our statistical error in all cases, often by an order of magnitude. Nevertheless, we applied the correction for all lattice spacings that we use in this paper.

Discussion of the width parameter
A natural definition of the width of the transition is given by the curvature at the peak of the susceptibility, i.e.
The second derivative of χ(T ) is numerically difficult to obtain, because it suffers from the systematic errors in the fit of χ(T ). However, we can more easily model χ( ψ ψ ) in a larger range of ψ ψ . Thus, we introduce a proxy δT for the half width ∆T where ψ ψ −1 stands for the inverse function of ψ ψ (T ). If we define ∆ ψ ψ analogously to Eq. (22) as  (10) TABLE II. Continuum extrapolated values of the cross-over temperature Tc, the half width of the susceptibility peak δT , the height of the susceptibility peak χ(Tc), the chiral condensate ψ ψ (Tc) and the width parameter ∆ ψ ψ from our ImμB-by-ImμB analysis. Our final κ2 and κ4 were obtained from a separate, more precise analysis, with correlated Tc values. We give this table to enable other researchers to use our data for further model building.
then using the fact that χ has a maximum at T c , i.e.
we find that δT = ∆T up to higher order corrections in ∆ ψ ψ . The values for this proxy are shown in Table II. The proxy σ used in the main text -given by Eq. (4) -is even simpler. There, we fixed ∆ ψ ψ to its µ B = 0 value for all chemical potentials. This is justified by the last column of Table II: the value ∆ ψ ψ = 0.14 is consistent with all the values in the range. For this reason, at µ B = 0 we have ∆T ≈ δT ≈ σ .

Tabulated continuum extrapolated results
The extrapolation of T c (µ B ) presented in this letter is the result of a global correlated fit. We could have pursued a different strategy: calculate the continuum extrapolated χ( ψ ψ ) curves for each imaginaryμ B , determine T c for each lattice at the givenμ B and then continuum extrapolate T c for fixedμ B . This method has the advantage to yield statistically uncorrelated T c (μ B ) values. Its disadvantage is the lower precision: the simple polynomial continuum extrapolation of the entire χ function has a better fit quality if we restrict the data to a narrow peak region. Thus the best fits have fewer data points. Nevertheless, we tabulate statistically independent transition temperatures in Table II Tc and the width parameter σ at real µB from analytical continuation. The errors are the combined statistical and systematic errors discussed in the main text and the supplemental material.
In Table III we provide the extrapolated cross-over temperature T c and the width parameter σ at real µ B from analytical continuation. They correspond to the green band in Fig. 3 and the blue band in the upper panel of Fig. 5, respectively. Unlike in Table II, the errors here are correlated, since these results are the output of a global analysis.