Continuum extrapolated high order baryon fluctuations

Fluctuations play a key role in the study of QCD phases. Lattice QCD is a valuable tool to calculate them, but going to high orders is challenging. Up to the fourth order, continuum results are available since 2015. We present the first continuum results for sixth order baryon fluctuations for temperatures between $T=130 - 200$ MeV, and eighth order at $T=145$ MeV in a fixed volume. We show that for $T \leq 145$ MeV, relevant for criticality search, finite volume effects are under control. Our results are in sharp contrast with well known results in the literature obtained at finite lattice spacing.

1. Introduction: fluctuations and QCD phases: The main goal of the heavy ion program of many accelerator facilities (e.g., at LHC, RHIC or the upcoming CBM/FAIR) is to create new phases of matter and explore their properties under extreme conditions.Several experimental programs (such as the Beam Energy Scan program at RHIC [1,2]) are designed to search for a hypothetical critical end point in the temperaturebaryon density phase diagram.Some of the most important observables in this quest are fluctuations of conserved charges.In the grand canonical ensemble, they are derivatives of the pressure with respect to the chemical potentials coupled to the charges.In this work, we will calculate such fluctuation observables at zero baryochemical potential.
The physics applications of fluctuation observables are numerous.First, the equation of state of the hot-anddense quark gluon plasma is one of the main inputs of lattice QCD to the phenomenology of heavy ion physics.Fluctuations at zero baryochemical potential µ B are the basis for extrapolations of the QCD equation of state to non-zero µ B , both by means of a truncated Taylor expansion [3], as well as via different resummations of the Taylor series [4][5][6].
Second, fluctuation observables are sensitive to criticality.While first-principle lattice simulations have shown that the chiral transition is a crossover at zero baryon density [7], at larger baryon densities several model calculations predict a critical endpoint in the 3D Ising universality class [8][9][10][11], where the crossover line becomes a line of first order transitions.One of the proposed experimental signatures of such a critical endpoint is a non-monotonic behavior of the fourth-to-second order baryon number fluctuations as a function of µ B [12,13].The extrapolation of this ratio to µ B > 0 is possible, if a sufficient number of Taylor coefficients (fluctuations) are available at µ B = 0. Thus, fluctuations at µ B = 0 are also important for the quest to find the critical endpoint.An important baseline in this search is given by the hadron resonance gas (HRG) model, a non-critical model that describes thermodynamics below the chiral transition at µ B = 0 remarkably well.A reasonable minimum criterium for criticality searches, then, is the presence of solid deviations between HRG predictions and equilibrium QCD.
A different type of criticality -in the O(4) universality class -is also expected to be present in QCD, related to chiral symmetry restoration, one of the most important concepts in heavy ion physics.In the limit of zero light quark masses, the SU (2) × SU (2) chiral symmetry becomes exact, and the crossover transition occurring at physical masses is expected to be replaced by a genuine second order transition [14][15][16].Through the presence of a single scaling variable (a combination of the quark masses, the temperature and the chemical potential), O(4) criticality has imprints on the temperature dependence of higher baryon number fluctuations at µ B = 0 up to physical values of the quark masses (for a model calculation, see e.g.Ref. [17]).The presence of such an O(4) scaling regime is likely the reason why lattice QCD calculation see no sharpening or strengthening of the crossover transition for small chemical potentials [4,18,19].The complicated interplay of O(4) and Ising criticality motivates more involved theoretical approaches to the phase diagram, based on Lee-Yang zeroes [20,21], which have gained popularity in recent years [22][23][24][25][26].These approaches, too, require reliable determinations of the corresponding fluctuations.
Third, fluctuation observables are very sensitive to the degrees of freedom of a thermodynamic system.This fact has been used before to argue, e.g., the existence of further (yet undiscovered) resonances [27][28][29] in the hadron spectrum, which was later confirmed by experiment [30].
Finally, fluctuation observables are central in the study of chemical freezeout in heavy ion collisions.This is an especially interesting avenue of research, since it has the potential of allowing a direct comparison between firstprinciple QCD predictions and experimental data.While the comparison itself has many caveats, due to experimental effects [31][32][33], it is obviously worth pursuing.
2. Current lattice estimates: Conserved charge fluctuations have been a focus of lattice simulations for well over a decade now.Like any other observable, a reliable calculation of these requires a continuum limit extrapolation, via simulations using smaller-and-smaller lattice spacings.Up to second order, they are known in the continuum since 2012 [34].Fourth order fluctuations in the baryon number and strangeness were first continuum extrapolated in 2015 [35].In that case, a large temperature range was considered, showing good agreement with the hadron resonance gas model at low temperatures, as well as good agreement with perturbative calculations [36,37] at high temperatures.Since then, calculations of the fourth order coefficients were pushed to very high precision [38].Thus, up to fourth order, the derivatives in full QCD are known accurately, with the exception of electric charge fluctuations, which suffer from large cut-off effects that make continuum extrapolations difficult [35].
At the sixth and eighth orders, the statistics requirements for the direct determination of the coefficients dramatically blow up, thus, a continuum extrapolation of these coefficients has never been attempted.Because of the high statistics required, all available results on higher fluctuations employ the computationally cheapest discretization: staggered fermions.These suffer from a lattice artefact called taste breaking [39,40], whose effect is to strongly distort the meson spectrum at a finite spacing.One might think that, since baryon number fluctuations directly couple to baryons only, and not mesons, this issue is of no great importance here.Unfortunately, this is only the case at the level of the baryon number B = 1 sector of the Hilbert space.Two-baryon interactions are mediated by mesons, and thus, taste breaking distorting meson physics can have relevant effects on the description of baryon dense matter (or, equivalently, on high order fluctuations at µ B = 0).In fact, the higher the order of the fluctuations, the larger the effect of multibaryon interactions on them.Hence, these fluctuations could be particularly sensitive to taste breaking effects.On a more basic level, higher order µ B -derivatives probe physics at larger µ B , which at a finite spacing will be closer to the cut-off scale of 1/a (a being the lattice spac-ing), making cut-off effects potentially important.
The statistics required for the calculation of higher order fluctuations can be drastically reduced by introducing a purely imginary chemical potential, calculating lower order fluctuations, and fitting their functional dependence on the imaginary chemical potential.The price for this reduction in statistics requirements is that assumptions have to be made on the functional form of the lower order fluctuations, leading to hard-to-control systematic errors.
So far, three collaborations have presented results up to the eighth order with improved lattice actions.In chronological order: First, the Pisa group presented results on lattices with 6 timeslices of 2stout improved fermions [7,41] in Ref. [42].Second, the Wuppertal-Budapest collaboration presented results with 12 timeslices of 4stout improved fermions [43] in Ref. [44].These two calculations took advantage of simulations at imaginary chemical potential.Finally, the HotQCD collaboration presented results with 8 timeslices of HISQ fermions [45], using a direct determination at µ B = 0 (i.e., without imaginary µ B simulations) in Refs.[6,38].For the latter calculation, two orders of magnitude more statistics were collected, compared to the previous two.It is a testament to the efficacy of the imaginary chemical potential method, then, that the error bars on the imaginary chemical potential calculations of the sixth and eighth order coefficients are substantially smaller.
At the current level of precision, discrepancies emerge between the calculations.In particular, the results based on the imaginary chemical potential method are in good agreement with the hadron resonance gas model for low temperatures.On the other hand, the direct calculation shows significant deviations for both observables even at the lowest temperature considered.To shed light on QCD criticality, this discrepancy has to be resolved.
In addition to the different extraction method for the higher order coefficients (direct at µ B = 0 vs indirect from µ 2 B ≤ 0), a potentially more significant difference between the two types of calculation lies in how the chemical potential is defined.Due to the extreme cost of the direct method, for fluctuations of order six and higher, the chemical potential in that case was introduced via the so-called linear prescription.In the other two cases, the exponential definition was used.Derivatives with respect to the chemical potential can be shown to be UV finite by virtue of a U (1) symmetry if the chemical potential is introduced like a constant imaginary gauge field [46].This is the exponential definition.If, instead, a naive linear definition is employed, power-law UV divergences appear in the free energy already in the case of free quarks.By taking enough derivatives with respect to the chemical potential, such power-law divergences disappear.However, this is not the case for logarithmic divergences, the absence of which for the naive linear definition is not proven for the interacting theory.Thus, although the linear definition is computationally cheaper, care should be taken in considering these results.
The linear definition also breaks the exact Roberge-Weiss periodicity [47] of the partition function.Even if one assumes that there are no problems with logarithmic divergences, the loss of Roberge-Weiss periodicity with the linear definition can potentially lead to large cut-off effects, since it effectively means that at a finite spacing, in contrast to the continuum, the free energy gets contributions from Hilbert subpsaces not only at integer, but also at non-integer values of the baryon number, which at the very least, is a non-physical feature at finite spacing.

Lattice calculation of fluctuations up to eighth order
In this letter, we present the first continuum results for baryon fluctuations up to the sixth order for temperatures between T = 130 − 200 MeV, and up to the eighth order for a temperature of T = 145 MeV.Continuum extrapolation is made possible by the introduction of a new discretization, which we call the 4HEX action, that strongly suppresses taste breaking effects compared to all available actions in the literature.Although more costly, we pursue a direct determination at µ B = 0, in order to avoid possible systematic effects due to a choice of fit ansatz, necessary for the imaginary chemical potential method.Moreover, in order to avoid possible issues with the introduction of the chemical potential, we employ the exponential definition to all orders.Due to the extreme statistics cost of the direct method, this endeavour is only feasible in a volume that is smaller than what is typically used in the field, with an aspect ratio LT = 2. Thanks to the availability in the literature of the aforementioned results at finite lattice spacing, but with larger volume, we are able to show that below T = 145 MeV, finite volume effects in our results are under control.Note that this is the relevant temperature range for the search for the elusive critical endpoint of QCD.
The novel lattice action we use for this thermodynamics study, 4HEX, is based on rooted staggered fermions with 4 steps of HEX smearing [48] with physical quark masses, and the DBW2 gauge action [49].This lattice action benefits from dramatically reduced taste breaking effects, compared to all other actions used in the litera-ture.We simulate 16 3 × 8, 20 2 × 10 and 24 3 × 12 lattices to obtain a well-controlled continuum extrapolation.Details on the 4HEX action, the scale setting procedure, and the systematic error estimation can be found in the supplemental material.We calculate fluctuations of the baryon number at zero strangeness chemical potential: . ( We also include results on the strangeness neutral line n s ≡ 0 in the Supplemental Material, which lead to similar conclusions as in the µ S = 0 case.
We use the exponential definition of the chemical potential at all orders in µ B on all our lattices.For the N τ = 8, 10 lattices, we use the reduced matrix formalism to calculate the fluctuations, in the same way as we did in Refs.[50,51].For the N τ = 12 lattice, we use the standard random source method [52].
We show our continuum extrapolated results for χ B 4 /χ B 2 (left) and χ B 6 /χ B 2 (center), together with the corresponding finite lattice spacing results in Fig. 1.The continuum results are obtained together with a spline fit of the temperature dependence.The exact procedure is described in the Supplemental Material.The bands include statistical and systematic uncertainties, consisting of different scale settings and different spline fits of the data.The covered temperature range is 130 MeV ≤ T ≤ 200 MeV.Also shown are the results on the N τ = 8, 10 lattices for χ B 8 /χ B 2 (right).For this observable, we also include the continuum extrapolation at a single temperature of T = 145 MeV.Hadron resonance gas predictions are shown, and they equal 1 in all cases independently from the temperature and the hadron spectrum used.
From Fig. 1, it is apparent how small the cut-off effects of the 4HEX action are, as is the fact that, for T < 145 MeV, the fluctuations in continuum QCD are in very good agreement with the HRG results.ature, with an aspect ratio LT = 4.In the top panel of Fig. 2, we show the continuum χ B 6 (yellow band) alongside previous results by the Wuppertal-Budapest collaboration obtained with the imaginary chemical potential method using 12 timeslices of 4stout improved staggered fermions [44] (black points), and results by the HotQCD collaboration obtained with the direct method and 8 time slices of HISQ fermions [6] (green band).The latter are not direct data, but rather come from a spline interpolation of the direct data.The comparison is fair, since we also use a spline interpolation on our data, to allow for systematic error estimation on the continuum limit.The role of the spline interpolation in thermodynamics studies is to i) reduce the random noise on the direct data, by utilizing continuity in T and to ii) allow one to take temperature derivatives, which are needed for the calculation of certain thermodynamic quantities, such as the speed of sound.The HotQCD collaboration has also used the splined version of their results on high order fluctuations in determining phenomenological quantities (see, e.g., Ref. [6]).Looking at the three results in comparison, a very simple interpretation emerges.First, at temperatures below 145 MeV, our results at N τ = 12 agree with our new continuum extrapolated results, even though the volume is different in the two simulations.On the other hand, the N τ = 8 HotQCD result does not agree with our old N τ = 12 results, even though the physical volumes are the same.This means that, at low temperatures, finite volume effects are small and cut-off effects are large.It appears that these N τ = 8 lattices are too coarse for phenomenological applications.The N τ = 8 results show very significant deviations from the hadron resonance gas.For example, it is noteworthy that the sign of dχ B 6 /dT is opposite to the HRG model prediction at all these low temperatures.Such deviations turned out to be cut-off effects, since the continuum extrapolated results show good quantitative agreement with the HRG.
A similar scenario appears in the bottom panel of Fig. 2, where we compare our new results on 20 3 × 10 lattices with results from the same LT = 4 simulations -Wuppertal-Budapest N τ = 12 [44] (black points) and HotQCD N τ = 8 [6] (green band) -as shown in the case of χ B 6 .We also include our new continuum extrapolation at T = 145 MeV.Besides the markedly smaller errors, our results are in good agreement with each other, showing small volume dependence especially at lower temperature.Moreover, quantitative agreement with HRG model predictions is evident up to T = 145 MeV.Our continuum result at T = 145 MeV confirms these findings.As in the previous case, it appears that cut-off effects for T ≤ 145 MeV are too large to allow for a safe use of results on coarse lattices for phenomenological applications.Finally, we note that the for T > 145 MeV, our old N τ = 12 large volume and new N τ = 10 small volume results do not agree.In particular, the local minimum of χ B 8 is shifted to lower temperatures.This is likely a finite volume effect, due to the crossover transition moving to slightly lower T in a smaller volume.
5. Discussion: In this letter we have reported the first continuum extrapolated results for high order baryon number fluctuations available in the literature.By means of a novel discretization of the QCD action, we were able to carry out a continuum extrapolation using lattice with 8,10 and 12 time slices.We used an aspect ratio of LT = 2.We calculated sixth order fluctuations in the continuum in a temperature range between T = 130 − 200 MeV.We also calculated the eighth order fluctuations in the continuum at a single temperature T = 145 MeV.A comparison of our results with existing results at finite lattice spacing and larger physical volumes showed that, in the temperature regime relevant for the critical point search, volume effects are well under control already for the smaller volume used in our study.In contrast, cut-off effects in previous results in the literature where not always under control, especially at the lower temperatures relevant for constraining the position of the critical endpoint.Thus, phenomenological conclusions based on erroneous coefficients ought to be reexamined in the near future.These include estimates of the radius of convergence and poles of Padé approximants, used to constrain the location of the critical endpoint and the convergence of the Taylor expansion for the equation of state, used as input for hydrodynamic simulations as well as estimates of fluctuation observables at µ B > 0 used to study chemical freeze-out.We plan to revisit all of these points in future publications.
Acknowledgements: The project was supported by the BMBF Grant No. 05P21PXFCA.This work is also supported by the MKW NRW under the funding code NW21-024-A.Further funding was received from the DFG under the Project No. 496127839.This work was also supported by the Hungarian National Research, Development and Innovation Office, NKFIH Grant No. KKP126769.This work was also supported by the NKFIH excellence grant TKP2021 NKTA 64.D.P. is supported by the ÚNKP-23-3 New National Excellence Program of the Ministry for Culture and Innovation from the source of the National Research, Development and Innovation Fund.The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputers Juwels-Booster at Juelich Supercomputer Centre and HAWK at Höchstleistungsrechenzentrum Stuttgart.

SUPPLEMENTAL MATERIAL 4HEX action
In this work we employ simulations with N f = 2+1+1 staggered dynamical quark flavours with physical masses.In the gauge sector we use the DBW2 action [49], and in the fermion sector we use links with 4 levels of HEX smearing [48].
Before quantitative results can be extracted, the action must be parameterized, meaning that the dependence of the bare quark masses on the gauge coupling must be tuned, so that the theory is renormalized in the parameter range of interest.In our simulations we used degenerate up and down quark masses, and fixed the light-strange mass ratio to m s /m ud = 27.63, and the strange-charm mass ratio to m c /m s = 11.85.In order to tune the light quark mass, we used the ratio π /f π = 1.0337,where we used the isospin-averaged pion mass m π = 134.8MeV [53] and the pion decay constant f π = 130.41 MeV.In order to adopt a proper comparison, we corrected both m π and f π for finite volume effects using results from two-flavour chiral perturbation  theory [54].For each value of the coupling, we simulated 3-5 values of the quark masses, then extracted the tuned masses via a linear interpolation in m 2 π /f 2 π .An example of the mass tuning for β = 0.688 is shown in the top panel of Fig. 3.
We parameterized the action in the range β = 0.46 − 0.83, corresponding to lattice spacings in the range a = 0.261 − 0.072 fm.The line of constant physics is shown over the whole range in the center panel of Fig. 3.We set the scale with the pion decay constant f π , as well as with a modified version of the Wilson-flow-based w 0 [55], which we refer to as w 1 Both scale settings are shown in the bottom panel of Fig. 3.The definition for w 0 is the following: where t is the flow time, and W (t) is defined as: where, in turn, ⟨E(t)⟩ is the expectation value of the continuum-like action density G a µν (t)G a µν (t)/4.The quantity w 1 is similarly defined as: and the continuum value of the ratio w q /w 0 is w 1 /w 0 = 1.515 (1).We show the approach to the continuum of this ratio in Fig. 4, where three linear fits differing by the range considered are shown.The final continuum value is obtained combining the different fit results.
The main advantage of the new action we are introducing is the reduced taste breaking one observes in the pion sector.In Fig. 5 we show this as a function of the lattice spacing for our previously used 4stout action (red) and our new 4HEX action (red).It is clear that the deviation of the pion masses from the target value are much smaller, by up to more than one order of magnitude for finer lattices.We will see in the following that this has crucial consequences on the feasibility of a continuum extrapolation with accessibly fine lattices.

Systematic analysis
In order to obtain the continuum extrapolations shown in Figs. 1, 2, we used the following ansatz: where the s i (T ) are a set of basis spline functions.In order to take into account the ambiguity in the temperature definition, we employ settings of the scale with both f π and w 1 , as previously mentioned, and include three different sets of node points to define the basis functions s i (T ).The final results are obtained by combining the 6 = 2×3 analyses to construct a histogram.The width of the histogram defines the systematic error.In the plots we show combined errors, where we assume that statistical and systematic errors add up in quadrature.This systematic error estimation procedure has been used and described in more detail in several of our works, most recently in Ref. [44].We show in Fig. 6 two example continuum extrapolations of χ B 6 , at T = 145, 170 MeV, and the continuum extrapolation of χ B 8 at T = 145 MeV.In all cases, one can appreciate the small cutoff dependence of these quantities.In the case of χ B 6 , we show the continuum extrapolation in pink, together with the result of the spline + continuum fit with the ansatz described in the previous paragraph.
-0.For the different temperatures, the finite Nτ data are shown along with their continuum extrapolation.In the case of χ B 6 , the results of the combined spline and continuum fit (depicted in Fig. 2), are also shown.

Strangeness Neutrality
In the main text we presented results for zero strangeness chemical potential µ S = 0.A phenomenologically more relevant case is that of strangeness neutrality, where the strangeness chemical potential has to be tuned as a function of the baryochemical potential µ S (µ B ) such that the expectation value for the strangeness density is zero n S ∝ ∂p ∂µ S = 0.This is done to every order in the Taylor expansion in µ B .The strangeness chemical potential itself is expanded as follows: where the coefficients s 1 , s 3 and s 5 are determined by solving ∂p ∂µ S = 0 order-by-order in the expansion parameter µ B T .Our results for these coefficients are shown in Fig. 7.We also show the HRG predictions for these coefficients.We observe a good agreement for temperatures below 145MeV.
Once these coefficients are known, we can also calculate the Taylor coefficients of the pressure on the strangeness neutral line: Our results for the Taylor coefficient ratios p 4 /p 2 , p 6 /p 2 and p 8 /p 2 are shown in Fig. 8.We observe mild cut-off effects.Again, the HRG predictions are also shown.We also show our results for p 6 and p 8 in Fig. 9.We show both the finite spacing results 16 3 × 8, 20 3 × 10 and 24 3 × 12, as well as the continuum extrapolated results.Like in the main text, for the eith order coefficients p 8 the 24 3 × 12 result, as well as the continuum extrapolated result is only available at a single temperature of T = 145MeV.Our results are compared to results of the HotQCD collabotion using 32 3 × 8 HISQ lattices from Ref. [6].All conclusions are similar to the case of the coefficients χ B 6 and χ B 8 , with the exception that unlike the case of χ B 6 and χ B 8 , for the coefficients of p 6 and p 8 there are no results from the imaginary chemical potential method using the 4stout 48

FIG. 1 .
FIG. 1.Our lattice results for the ratios χ B 4 /χ B 2 (left), χ B 6 /χ B 2 (center), χ B 8 /χ B 2 (right).For the first two, the continuum extrapolation is shown as a yellow band.HRG model predictions are shown as solid black lines in all cases.

4 .FIG. 2 .
FIG. 2. Our results for sixth (top) and eighth (bottom) order baryon number susceptibilities.The HISQ (splined) results are shown as light blue bands, and the 4stout results are shown are black points.Our new results are shown as a yellow band for the continuum extrapolation of χ B 6 and the continuum extrapolation of χ B 8 (at T = 145 MeV), and as brown points for χ B 8 on the 20 3 × 10 lattice.HRG model predictions are shown as solid black lines.

w 1 FIG. 3 .
FIG.3.Top: Example of quark mass tuning for simulations at β = 0.688.The orange band is a linear fit the data, and the black line is the target value of the ratio m 2 π /f 2 π .Center: Line of constant physics i.e., light quark mass as a function of the inverse bare gauge coupling.The points show the direct determination from the procedure in the top panel, the band is a polynomial fit as a function of β.Bottom: lattice spacing as determined through the pion decay constant (red) and w1 scale (orange).The band is a polynomial fit as a function of the coupling β.

FIG. 5 .
FIG.5.Taste breaking for our previous 4stout action with Symanzik-improved gauge sector (red), HotQCD's HISQ action (black), and our new 4HEX action with DBW2 gauge sector (blue).The improvement in the taste breaking is evident.

TABLE I .
3× 12 lattice action.Number of configurations analyzed on our three lattice geometries.