QCD Equation of State at very high temperature: computational strategy, simulations and data analysis

TL;DR

This work delivers a comprehensive non-perturbative determination of the QCD equation of state up to electroweak-scale temperatures using lines of constant physics defined by a Schrödinger-functional coupling and shifted boundary conditions that directly yield the entropy density. The lattice implementation combines Wilson gauge action with O(a) improved Wilson fermions, meticulous renormalization, and a continuum-extrapolated entropy density with about 1% accuracy across nine temperatures from 3 GeV to 165 GeV. By integrating in both the bare gauge coupling and the bare quark mass, and employing perturbative improvements and variance-reduction techniques, the study achieves controlled systematic errors, including finite-volume and topology effects. Comparisons to standard and HTL perturbation theory reveal sizeable higher-order and ultrasoft contributions even at the highest temperatures, while the methodology is sufficiently general to extend to four or five flavors and to other thermal observables, marking a robust platform for non-perturbative QCD thermodynamics at extreme temperatures.

Abstract

We present a detailed account of the theoretical progress and the computational strategy that led to the non-perturbative determination of the QCD Equation of State (EoS) at temperatures up to the electroweak scale reported in [Phys. Rev. Lett. 134, 201904 (2025)]. The two key ingredients that make such a calculation feasible with controlled uncertainties are: (i) the definition of lines of constant physics through the running of a non-perturbatively defined finite-volume coupling across a wide range of energy scales, and (ii) the use of shifted boundary conditions which allow a direct determination of the entropy density thus without the need for a zero-temperature subtraction. Considering the case of QCD with $N_f =3$ massless flavours in the temperature interval between 3 GeV and 165 GeV, we describe the numerical strategy based on integrating in the bare coupling and quark mass, the perturbative improvement of lattice observables, the optimization of numerical simulations, and the continuum extrapolation. Extensive consistency checks, including finite-volume and topological-freezing effects, confirm the robustness of the method. The final results have a relative accuracy of about $1\%$ or better, and the errors are dominated by the statistical fluctuations of the Monte Carlo ensembles. We also compare our non-perturbative results with predictions from standard and hard thermal loop perturbation theory showing that at the level of $\%$-precision contributions beyond those known, including non-perturbative ones due to ultrasoft modes, are relevant up to the highest temperatures explored. The methodological framework is general and readily applicable to QCD with four and five massive quark flavours and to other thermal observables, paving the way for systematic non-perturbative studies of thermal QCD at very high temperatures.

Figures (14)

• Figure 1: Plot of the integrand function in Eq. \ref{['eq:Df_gauge']} as a function of the bare coupling $g_0^2$. Points have been shifted horizontally by $0.03\times(L_0/a-4)$ for better readability.
• Figure 2: Left: plot of the integrand function in Eq. \ref{['eq:Df_quark']} computed at the bare parameters of temperature $T_1$ and at the resolutions $L_0/a=4,6,8,10$, as a function of $m_q/T$. Points have been interpolated with a cubic spline to guide the eye. In most cases, errors are smaller than the markers. Right: the same integrand function is shown at the resolution $L_0/a=6$ and at three temperatures.
• Figure 3: Black dots are the values of the one-loop improved normalized entropy density as a function of $(a/L_0)^2$ at the temperatures $T_0, T_1, ..., T_8$. The red band is our best extrapolation to the continuum limit (fit id3 in Table \ref{['tab:fits_comparison']}). Red crosses are the continuum extrapolated values for $s/T^3$. The horizontal axis is common to all the subplots.
• Figure 4: Left: comparison among the results of different continuum limits. Points have been shifted horizontally by $0.023\times(n/3-1)$ where $n=0,1,...,6$ for fits id0, id1, ..., id6. Right: effect of logarithmic corrections on the continuum extrapolated values of the best fit. For better readability points have been shifted horizontally by $0.01\times\gamma$, with $\gamma$ defined in Eq. \ref{['eq:log_corr']}.
• Figure 5: Comparison of the three parametrizations of the temperature dependence of $s/T^3$, in the temperature interval covered by the non-perturbative data.