HotQCD collaboration

# Equation of state and speed of sound of (2+1)-flavor QCD in strangeness-
neutral matter at non-vanishing net baryon-number density.

D. Bollweg  Physics Department, Brookhaven National Laboratory, Upton, New
York 11973, USA  RIKEN-BNL Research Center, Brookhaven National Laboratory,
Upton, New York 11973  D. A. Clarke  Fakultät für Physik, Universität
Bielefeld, D-33615 Bielefeld, Germany  Department of Physics and Astronomy,
University of Utah, Salt Lake City, Utah, United States  J. Goswami  RIKEN
Center for Computational Science, Kobe 650-0047, Japan  O. Kaczmarek  Fakultät
für Physik, Universität Bielefeld, D-33615 Bielefeld, Germany  F. Karsch
Fakultät für Physik, Universität Bielefeld, D-33615 Bielefeld, Germany
Swagato Mukherjee  Physics Department, Brookhaven National Laboratory, Upton,
New York 11973, USA  P. Petreczky  Physics Department, Brookhaven National
Laboratory, Upton, New York 11973, USA  C. Schmidt  Fakultät für Physik,
Universität Bielefeld, D-33615 Bielefeld, Germany  Sipaz Sharma  Fakultät für
Physik, Universität Bielefeld, D-33615 Bielefeld, Germany

###### Abstract

We update results on the QCD equation of state in (2+1)-flavor QCD with non-
zero conserved charge chemical potentials obtained from an eighth-order Taylor
series. We present results for basic bulk thermodynamic observables of
strangeness-neutral strong-interaction matter, i.e. pressure, number
densities, energy and entropy density, and resum Taylor series results using
Padé approximants. Furthermore, we calculate the speed of sound as well as the
adiabatic compression factor of strangeness-neutral matter on lines of
constant entropy per net baryon number. We show that the equation of state
($P(n_{B}),\epsilon(n_{B})$) is already well described by the $4^{\rm
th}$-order Taylor series in almost the entire range of temperatures accessible
with the beam energy scan in collider mode at the Relativistic Heavy Ion
Collider.

###### pacs:

11.10.Wx, 11.15.Ha, 12.38.Aw, 12.38.Gc, 12.38.Mh, 24.60.Ky, 25.75.Gz, 25.75.Nq

##  I Introduction

Fundamental properties of hot and dense matter controlled by the strong force
are encoded in the equation of state (EoS). The EoS is sensitive to the change
of degrees of freedom that dominate properties of strong-interaction matter at
low temperatures and densities on the one hand and high temperature and
densities on the other hand. It describes fluctuations of e.g. energy and
particle densities that signal the occurrence of phase transitions. The EoS
finds application in the hydrodynamic modelling of hot and dense matter
created in heavy ion collisions [1] and in constraining the “cosmic
trajectory” of matter in the expanding early universe [2].

The EoS of strong-interaction matter is obtained from quantum chromodynamics
(QCD) at finite temperature and non-zero values of the conserved charge
chemical potentials. It is well studied at vanishing values of the chemical
potentials [3, 4, 5]. For non-vanishing values of the chemical potentials that
couple to net baryon number ($\mu_{B}$), net electric charge ($\mu_{Q}$) and
net strangeness ($\mu_{S}$), lattice QCD calculations have to face the well-
known sign problem. Currently, this renders direct numerical calculations with
non-zero, real-valued chemical potentials impossible. For this reason
numerical calculations at non-zero $(\mu_{B},\mu_{Q},\mu_{S})$ have been
performed by either using Taylor expansions in terms of the chemical
potentials [6, 7] or by performing numerical simulations at imaginary values
of the chemical potential [8, 9]. Using the latter and analyzing results with
various ansätze for an analytic continuation to real values of the chemical
potential, the EoS of (2+1)-flavor QCD has been determined [10, 11]. Following
the Taylor expansion approach, results for the QCD EoS have been obtained from
up to sixth-order expansions [12, 13].

We extend here the Taylor expansion approach for the EoS in (2+1)-flavor QCD
using recent high-statistics results for Taylor expansion coefficients up to
eighth-order [14] in ($\mu_{B},\mu_{Q},\mu_{S})$. We will make use of the
Padé-resummation of Taylor series discussed in [14] and determine direct
Taylor series results as well as Padé-resummed approximants for pressure, net
baryon-number, energy and entropy densities. Making use of our high-statistics
data at temperatures in the vicinity of the pseudo-critical temperature we
will, for the first time, present contributions to higher-order derivatives of
bulk thermodynamic observables, i.e. the speed of sound and adiabatic
compressibility of strangeness-neutral strong-interaction matter. These
quantities have been calculated previously at vanishing values of the chemical
potentials [15, 4, 5].

We focus on a situation most relevant in the context of heavy ion collision
experiments, in which overall net strangeness number density ($n_{S}$)
vanishes and the ratio of net electric charge ($n_{Q}$) to net baryon-number
($n_{B}$) densities is close to the isospin-symmetric limit of strangeness
neutral matter, $n_{Q}=n_{B}/2$. In collisions of heavy nuclei, e.g. gold or
lead, the ratio $n_{Q}/n_{B}$ is approximately 0.4. Earlier studies of the QCD
EoS with a strange quark mass tuned to its physical value, and degenerate
light quark masses tuned to reproduce the physical, light hadron spectrum,
have been performed for different values of $r\equiv n_{Q}/n_{B}$. They have
shown that differences in bulk thermodynamic observables, arising from the
deviation of $r$ from the isospin-symmetric value, $r=0.5$, are small [13]. We
therefore stick here to an analysis of the thermodynamics of strangeness-
neutral, isospin-symmetric systems, which is equivalent to setting the
electric charge chemical potential $\mu_{Q}$ to zero.

This paper is organized as follows. In Section II we summarize basic
thermodynamic relations and outline calculations in QCD with non-vanishing
chemical potentials using the Taylor series approach and the resummation of
Taylor series using Padé approximants. In Section III we present our results
on Taylor expansion coefficients for bulk thermodynamic observables and
discuss their properties at low and high temperature by comparing with HRG
model calculations and $\order{g^{2}}$ high-temperature perturbation theory,
respectively. We furthermore discuss the structure of expansion coefficients
in the vicinity of the pseudo-critical temperature. In Section IV we analyze
thermodynamic observables as function of the baryon chemical potential. We
determine the baryon chemical potential as function of net baryon-number
density to arrive at the QCD equation of state, $P(n_{B})$. We furthermore
present results for basic thermodynamic observables on the pseudo-critical
line as well as on lines of constant ratio of entropy to net baryon number.
The latter are discussed in Section V. We finally give our conclusions in
Section VI. In three appendices we give some details on directional partial
derivatives (A) and discuss $\order{g^{2}}$ perturbative results for the EoS
of strangeness-neutral matter (B) as well as the calculation of the isentropic
speed of sound (C).

##  II Thermodynamics in (2+1)-flavor QCD and global constraints

We present here the basic thermodynamic observables that will be analyzed by
us in later sections using results from up to eighth-order Taylor expansions
of the pressure of QCD with a strange quark mass tuned to its physical value
and two degenerate light quark masses that correspond in the continuum limit
to a physical pion mass value of about $135$ MeV.

We also briefly summarize basic relations used in the Taylor expansion
approach to the thermodynamics of QCD at non-zero values of the conserved
charge chemical potentials and discuss our approach to the utilization of
Padé-resummed Taylor series for thermodynamic observables. Further details on
these three topics can be found in [14].

###  II.1 Thermodynamic observables

In a grand canonical ensemble the pressure, $p$, is a function of temperature
$T$ and a set of chemical potentials $\vec{\mu}=(\mu_{B},\mu_{Q},\mu_{S})$,
which couple to the currents of conserved charges for net baryon number ($B$),
electric charge ($Q$) and strangeness ($S$). The pressure is given in terms of
the logarithm of the grand canonical partition function

$\displaystyle\frac{p}{T^{4}}$ $\displaystyle=$
$\displaystyle\frac{1}{VT^{3}}\ln\mathcal{Z}(T,V,\hat{\mu}_{B},\hat{\mu}_{Q},\hat{\mu}_{S})\;.$
(1)

We will often use dimensionless variables and observables obtained by
rescaling the dimensionful observables with appropriate powers of the
temperature, e.g. $\hat{\mu}_{X}\equiv\mu_{X}/T$ or $\hat{p}\equiv p/T^{4}$.
With this we obtain for the number densities,

$\hat{n}_{X}=\frac{1}{VT^{3}}\frac{\partial\ln\mathcal{Z}(T,V,\hat{\mu}_{B},\hat{\mu}_{Q},\hat{\mu}_{S})}{\partial\hat{\mu}_{X}}\;,\;X=B,\
Q,\ S\;,$ (2)

and the energy density,

$\displaystyle\hat{\epsilon}$ $\displaystyle=$
$\displaystyle\frac{1}{VT^{2}}\frac{\partial\ln\mathcal{Z}(T,V,\hat{\mu}_{B},\hat{\mu}_{Q},\hat{\mu}_{S})}{\partial
T}$ (3) $\displaystyle=$ $\displaystyle
3\hat{p}+T\frac{\partial\hat{p}}{\partial T}\;.$

The entropy density reads

$\displaystyle\hat{s}$ $\displaystyle=$
$\displaystyle\hat{\epsilon}+\hat{p}-\hat{\mu}_{B}\hat{n}_{B}-\hat{\mu}_{Q}\hat{n}_{Q}-\hat{\mu}_{S}\hat{n}_{S}\;.$
(4)

Enforcing the constraint for isospin-symmetric systems, $n_{Q}=n_{B}/2$ in
strangeness-neutral matter ($n_{S}=0$) is equivalent to demanding $\mu_{Q}=0$
[13]. Thermodynamic observables thus become functions of $T$ and
$\hat{\mu}_{B}$ only, i.e. $\hat{p}\equiv p(T,\hat{\mu}_{B})/T^{4}$ and
$\hat{s}=\hat{\epsilon}+\hat{p}-\hat{\mu}_{B}\hat{n}_{B}$.

In addition to the basic bulk thermodynamic observables introduced above, we
also will calculate their temperature derivatives, imposing external
constraints. We calculate the isentropic speed of sound in a strangeness-
neutral medium with fixed $n_{Q}/n_{B}$, [16],

$\displaystyle c_{s}^{2}\equiv c^{2}_{\vec{X}}$ $\displaystyle=$
$\displaystyle\left(\frac{\partial
p}{\partial\epsilon}\right)_{\vec{X}}=\frac{\left(\partial p/\partial
T\right)_{\vec{X}}}{\left(\partial\epsilon/\partial T\right)_{\vec{X}}}\;.$
(5)

with $\vec{X}=(s/n_{B},n_{Q}/n_{B},n_{S})$. This also gives the adiabatic
compressibility

$\displaystyle\kappa_{s}$ $\displaystyle=$
$\displaystyle\frac{1}{n_{B}}\left(\frac{\partial n_{B}}{\partial
p}\right)_{\vec{X}}=\frac{1}{c_{s}^{2}(\epsilon+p-\mu_{S}n_{S})}\;.$ (6)

The calculation of $c_{s}^{2}$ and $\kappa_{s}$ requires taking directional
derivatives, which can be done numerically on lines of fixed constraints
$\vec{X}$ in the parameter space $(T,\vec{\mu})$, or using analytic relations.
Some details on partial derivatives at fixed constraints are given in Appendix
A.

###  II.2 Taylor expansions

We start with the expansion of the pressure, $\hat{p}=p/T^{4}$, in terms of
the three conserved charge chemical potentials
$(\hat{\mu}_{B},\hat{\mu}_{Q},\hat{\mu}_{S})$,

$\displaystyle\hat{p}$ $\displaystyle=$
$\displaystyle\frac{1}{VT^{3}}\ln\mathcal{Z}(T,V,\hat{\mu}_{B},\hat{\mu}_{Q},\hat{\mu}_{S})$
(7) $\displaystyle=$
$\displaystyle\sum_{i,j,k=0}^{\infty}\frac{\chi_{ijk}^{BQS}}{i!j!\,k!}\hat{\mu}_{B}^{i}\hat{\mu}_{Q}^{j}\hat{\mu}_{S}^{k}\;,$

with $\chi_{000}^{BQS}\equiv p(T,0)/T^{4}$ and expansion coefficients
$\chi_{ijk}^{BQS}$,

$\chi_{ijk}^{BQS}\equiv\chi_{ijk}^{BQS}(T)=\left.\frac{\partial\hat{p}}{\partial\hat{\mu}_{B}^{i}\partial\hat{\mu}_{Q}^{j}\partial\hat{\mu}_{S}^{k}}\right|_{\vec{\mu}=0}\;.$
(8)

These expansion coefficients are cumulants of conserved charge fluctuations
evaluated at vanishing chemical potential111We often suppress the argument
($T$) of these generalized susceptibilities and also suppress superscripts and
subscripts of $\chi_{ijk}^{BQS}$ whenever one of the subscripts vanishes, e.g.
$\chi_{i0k}^{BQS}\equiv\chi_{ik}^{BS}$..

From this expansion and the constraints, we obtain Taylor expansions for
pressure and number densities,

$\displaystyle\Delta\hat{p}$ $\displaystyle\equiv$
$\displaystyle\frac{p(T,\mu_{B})}{T^{4}}-\frac{p(T,0)}{T^{4}}=\sum_{k=1}^{\infty}P_{2k}(T)\hat{\mu}_{B}^{2k}\;,$
(9) $\displaystyle\hat{n}_{X}$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{\infty}N_{2k-1}^{X}(T)\hat{\mu}_{B}^{2k-1}\;\;,\;\;X=B,\
Q,\ S\;\;.$ (10)

The construction of the expansion coefficients222See [17] for this notation
and explicit definitions of $\bar{\chi}_{0}^{B,n}$.,
$P_{2k}\equiv\bar{\chi}^{B,2k}_{0}/(2k)!$ and $N_{2k-1}$ has been discussed
previously for the general case of strangeness-neutral matter and a fixed
ratio $n_{Q}/n_{B}=r$ [13, 14]. This gives the chemical potentials
$\hat{\mu}_{Q}$ and $\hat{\mu}_{S}$ as Taylor series in terms of
$\hat{\mu}_{B}$,

$\displaystyle\hat{\mu}_{Q}(T,\mu_{B})$ $\displaystyle=$ $\displaystyle
q_{1}(T)\hat{\mu}_{B}+q_{3}(T)\hat{\mu}_{B}^{3}+q_{5}(T)\hat{\mu}_{B}^{5}+\dots$
(11) $\displaystyle\hat{\mu}_{S}(T,\mu_{B})$ $\displaystyle=$ $\displaystyle
s_{1}(T)\hat{\mu}_{B}+s_{3}(T)\hat{\mu}_{B}^{3}+s_{5}(T)\hat{\mu}_{B}^{5}+\dots\;.$
(12)

Here we only consider the isospin-symmetric case, $n_{Q}/n_{B}=1/2$, which
fixes $\mu_{Q}=0$. In that case one finds, for instance, a simple relation
between the expansion coefficients of the net baryon-number density and the
pressure, $N_{2k-1}^{B}=2kP_{2k}$.

Taking derivatives of the expansion coefficients with respect to temperature
we also obtain Taylor series for the energy and entropy densities,

$\displaystyle\Delta\hat{\epsilon}$ $\displaystyle\equiv$
$\displaystyle\frac{\epsilon(T,\hat{\mu}_{B})}{T^{4}}-\frac{\epsilon(T,0)}{T^{4}}=\sum_{k=1}^{\infty}\epsilon_{2k}(T)\hat{\mu}_{B}^{2k}\;,$
(13) $\displaystyle\Delta\hat{s}$ $\displaystyle\equiv$
$\displaystyle\frac{s(T,\hat{\mu}_{B})}{T^{3}}-\frac{s(T,0)}{T^{3}}=\sum_{k=1}^{\infty}\sigma_{2k}(T)\hat{\mu}_{B}^{2k}\;.$
(14)

We note that to obtain Eqs. 13 and 14 one must take the temperature
derivatives before applying the constraints. In the case $\mu_{Q}=0,\
n_{S}=0$, as well as for $\mu_{Q}=0,\ \mu_{S}=0$ the expansion coefficients
are directly related to the expansion coefficients of the Taylor series of the
pressure,

$\displaystyle N_{2k-1}^{B}(T)$ $\displaystyle=$ $\displaystyle
2kP_{2k}(T)\;\;,$ (15) $\displaystyle\epsilon_{2k}(T)$ $\displaystyle=$
$\displaystyle 3P_{2k}(T)+T\frac{{\rm d}P_{2k}(T)}{{\rm d}T}\;\;,$ (16)
$\displaystyle\sigma_{2k}(T)$ $\displaystyle=$
$\displaystyle(4-2k)P_{2k}(T)+T\frac{{\rm d}P_{2k}(T)}{{\rm d}T}\;\;.$ (17)

We also introduce the notation $\hat{\epsilon}_{0}\equiv\epsilon(T,0)/T^{4}$
and $\hat{\sigma}_{0}\equiv s(T,0)/T^{3}$.

###  II.3 Resummation with Padé approximants

In [14] we analyzed the resummation of Taylor series for the pressure and net
baryon-number density using $[n,m]$ Padé approximants. We showed that in the
temperature range [135 MeV:175 MeV], the [4,4] Padé approximants for the
pressure have poles only in the complex plane. The location of these poles
provide an estimate for the range of validity of the Taylor series. For
chemical potentials in this range, one finds good agreement between the Padé
approximants and the eighth-order Taylor series for pressure and net baryon-
number density333We note that the recently developed multi-Padé approach [18]
provides similar information on the location of poles in the complex plane and
may, in the future, also be used to calculate thermodynamic observables.. We
will use here the Padé approximants for the Taylor series of the pressure also
to determine energy and entropy densities as well as other thermodynamic
observables introduced in Sec. II.1.

As has been discussed in [14], the Taylor series for the pressure of isospin-
symmetric matter with either $\mu_{S}=0$ or $n_{S}=0$ may conveniently be
rewritten as

$\displaystyle\frac{\Delta p(T,\mu_{B})}{T^{4}}$ $\displaystyle=$
$\displaystyle\frac{P_{2}^{2}}{P_{4}}\sum_{k=1}^{\infty}c_{2k,2}\bar{x}^{2k}\;,$
$\displaystyle=$
$\displaystyle\frac{P_{2}^{2}}{P_{4}}\left(\bar{x}^{2}+\bar{x}^{4}+c_{6,2}\bar{x}^{6}+c_{8,2}\bar{x}^{8}+...\right)\;\;,$

with $\bar{x}=\sqrt{P_{4}/P_{2}}\ \hat{\mu}_{B}$ and

$\displaystyle
c_{6,2}=\frac{P_{6}P_{2}}{P_{4}^{2}}\;\;,\;\;c_{8,2}=\frac{P_{8}P_{2}^{2}}{P_{4}^{3}}\;\;.$
(19)

Figure 1: The $n{}^{\rm th}$-order Taylor expansion coefficients, $P_{n}(T)$,
for the Taylor series of the pressure of (2+1)-flavor QCD as function of
$\hat{\mu}_{B}=\mu_{B}/T$ versus temperature. Shown are the expansion
coefficients $n!P_{n}$ for isospin-symmetric, strangeness-neutral matter
($\mu_{Q}=0$, $n_{S}=0$). The solid magenta lines show results from a
corresponding HRG model calculation using the QMHRG2020 list of hadron
resonances [14]. Yellow bands show the location of the pseudo-critical
temperature $T_{\rm pc}(0)=156.5(1.5)$ MeV [19]. The green bands at high
temperature show the $\order{g^{2}}$ perturbative result using a
renormalization scale in the range $[4\pi T,8\pi T]$.

Using expansion coefficients from fourth, sixth and eighth-order Taylor series
for the pressure we construct several $[n,m]$ Padé approximants,

$\displaystyle P[2,4]$ $\displaystyle=$
$\displaystyle\frac{\bar{x}^{2}}{1-\bar{x}^{2}+(1-c_{6,2})\bar{x}^{4}}\;\;,$
(20) $\displaystyle P[4,2]$ $\displaystyle=$
$\displaystyle\frac{\bar{x}^{2}+(1-c_{6,2})\bar{x}^{4}}{1-c_{6,2}\bar{x}^{2}}\;\;,$
(21) $\displaystyle P[4,4]$ $\displaystyle=$
$\displaystyle\frac{(1-c_{6,2})\bar{x}^{2}+\left(1-2c_{6,2}+c_{8,2}\right)\bar{x}^{4}}{(1-c_{6,2})+(c_{8,2}-c_{6,2})\bar{x}^{2}+(c_{6,2}^{2}-c_{8,2})\bar{x}^{4}}.$
(22)

The resulting approximants for the pressure are then given by

$\left(\frac{\Delta
p(T,\mu_{B})}{T^{4}}\right)_{[n,m]}=\frac{P_{2}^{2}}{P_{4}}P_{[n,m]}\;.$ (23)

For the number density we construct the $[3,4]$ Padé approximant.

Similar Padé approximants can be obtained for the energy density and entropy
densities by replacing the expansion coefficients $P_{2k}$ with the
corresponding expansion coefficients $\epsilon_{2k}$ and $\sigma_{2k}$. Rather
then constructing in this way Padé approximants for various thermodynamic
observables, we also directly determine Padé-based approximants from Eq. 23,
using the relevant thermodynamic relations, e.g. Eqs. 3, 4 and 10, for the
energy, entropy and net baryon-number densities. The Padé-based approximants,
called P-Padé in the following, ensure thermodynamic consistency among
different observables and, for instance, ensure that the singularities in the
approximations for, $n_{B}/T^{3}$, $\epsilon/T^{4}$ and $s/T^{3}$ coincide
with those of the Padé approximation for the pressure. The Padé-based (P-Padé
) result for the energy and entropy densities are obtained by taking
appropriate partial derivatives of the $[n,m]$ Padé approximation for the
pressure with respect to temperature and chemical potential, and using Eqs. 3
and 4,

$\displaystyle\left(\frac{\Delta\epsilon}{T^{4}}\right)_{[n,m]}$
$\displaystyle=$ $\displaystyle 3\left(\frac{\Delta
p}{T^{4}}\right)_{[n,m]}+T\frac{\rm d}{{\rm
d}T}\left(\frac{P_{2}^{2}}{P_{4}}P_{[n,m]}\right)\;,$ (24)
$\displaystyle\left(\frac{\Delta s}{T^{3}}\right)_{[n,m]}$ $\displaystyle=$
$\displaystyle\left(\frac{\Delta\epsilon}{T^{4}}\right)_{[n,m]}+\left(\frac{\Delta
p}{T^{4}}\right)_{[n,m]}$ (25)
$\displaystyle-\hat{\mu}_{B}\frac{P_{2}^{2}}{P_{4}}\frac{{\rm
d}P_{[n,m]}}{{\rm d}\hat{\mu}_{B}}\;.$

##  III Taylor expansion coefficients of the pressure in (2+1)-flavor QCD

We update here our previous analysis of the EoS of (2+1)-flavor QCD performed
in lattice QCD calculations with the highly improved staggered quark (HISQ)
action and a tree-level improved gauge action [13].

Figure 2: Derivatives of Taylor expansion coefficients $P_{2}$, $P_{4}$,
$P_{6}$, and $P_{8}$ with respect to temperature. The solid lines at low
temperatures show results from a corresponding HRG model calculation using the
QMHRG2020 list of hadron resonances [14]. Yellow bands show the location of
the pseudo-critical temperature $T_{\rm pc}$ [19] and the green bands at high
temperature show the $\order{g^{2}}$ perturbative result using a
renormalization scale in the range $[4\pi T,8\pi T]$.

###  III.1 Data sets and analysis details

For our calculation of bulk thermodynamic observables we use high-statistics
data sets for $(2+1)$-flavor QCD with light and strange quark masses tuned to
physical values ($m_{\ell}/m_{s}=1/27$). These data sets have been generated
using the HISQ action with tree-level coefficients and a tree-level improved
gauge action [13]. Gauge configurations were generated using SIMULATeQCD [20,
21]. New data sets exist for the temperature interval $[125~{}{\rm
MeV}:175~{}{\rm MeV}]$ and have been presented in [14]. Additional updates
exist for the lowest temperature, where we doubled the number of gauge field
configuration to 2.2M.

Compared to our previous calculation of the EoS [13], the data sets for
$N_{\tau}=8$ and $12$ contain more than an order of magnitude larger
statistics and, moreover, include entirely new data sets on lattices with
temporal extent $N_{\tau}=16$.

At temperatures $T>180$ MeV, we added to our analysis data from calculations
with a slightly larger light quark mass ($m_{\ell}/m_{s}=1/20$), which had
been used by HotQCD previously [13]. In all cases results have been obtained
on lattices with spatial extent $N_{\sigma}=4N_{\tau}$.

Continuum extrapolations of $P_{2}$ and $P_{4}$ have been performed in three
different temperature intervals, using different ansätze. In the temperature
interval $[135~{}{\rm MeV}:175~{}{\rm MeV}]$ we use rational polynomial
functions as described in [14] and [22]. For the high-temperature part $T\geq
175~{}{\rm MeV}$, we use a polynomial ansatz to interpolate the data to the
continuum limit. Moreover, for low temperatures in the range $[125~{}{\rm
MeV}:145~{}{\rm MeV}]$ we consider an HRG-motivated ansatz, which takes into
account the exponential drop of thermodynamic observables at low temperatures
and incorporates cutoff-dependent corrections,

$\displaystyle f(T)=A\left(\frac{B}{T}\right)^{3/2}{\rm
e}^{-B/T}\left(1+d_{1}T+d_{2}T^{2}+\frac{c}{N_{\tau}^{2}}\right)\;.$ (26)

Finally, we match all the three functions and their first derivative at
$T=135~{}{\rm MeV}$ and at $T=175~{}{\rm MeV}$. For the interpolation of
$P_{6}$ and $P_{8}$ we only used spline interpolations of the $N_{\tau}=8$
data set. Much of the data analysis in this study was facilitateurd through
the AnalysisToolbox [23]. Statistical uncertainty represented by bands in all
figures is calculated through bootstrap resampling. Central values are
returned as the median with lower and upper error bounds given by the 32% and
68% quantiles, respectively. Spline interpolations, when needed, are cubic
with evenly spaced knots. They are calculated using the LSQUnivariateSpline
method of SciPy [24]. Central values and error bands shown for the lattice QCD
results are smoothed using splines. Temperature derivatives of lattice QCD
data are calculated by fitting the temperature dependence with a spline, then
calculating the derivative of the spline numerically.

###  III.2 Taylor expansion coefficients for the pressure and the trace
anomaly

The Taylor series in terms of the conserved charge chemical potentials
$\vec{\mu}=(\mu_{B},\mu_{Q},\mu_{S})$ have been reorganized to obtain the
Taylor series for the pressure in strangeness-neutral, isospin-symmetric
matter given in Eq. 9. Results for the expansion coefficients
$(2k)!P_{2k}\equiv\bar{\chi}_{0}^{B,n}$ are shown in Fig. 1. For the fits and
interpolations of different order expansion coefficients shown in this figure
we use the same approach as discussed in [14]. I.e. for the
$\order{\hat{\mu}_{B}^{2}}$ expansion coefficient the band shows a continuum
extrapolation based on $N_{\tau}=6,\ 8,\ 12$ and $16$ datasets; for the
$\order{\hat{\mu}_{B}^{4}}$ expansion coefficient we show a continuum
extrapolation based on $N_{\tau}=6,\ 8$ and $12$ datasets and for the higher-
order expansion coefficients we only use results from our high-statistics
calculations on lattices with temporal extent $N_{\tau}=8$, where more than
1.5 million gauge field configurations have been generated at each temperature
value. In this case the curves show spline interpolations of the $N_{\tau}=8$
data.

In the following we discuss several features of the expansion coefficients in
the high, low and intermediate temperatures ranges, respectively.

####  III.2.1 High-temperature region

It is apparent that the expansion coefficients $P_{2k}$ approach the high-
temperature ideal gas limit rapidly. In fact, for temperatures $T\gtrsim 250$
MeV, deviations from the ideal gas values are only about $12$% for $P_{2}$ and
$P_{4}$. Moreover, at these temperature values the expansion coefficients
$P_{2k}$ start being consistent with zero for $k\geq 3$, as is the case in a
massless, ideal gas. For $k=1$ and 2, deviations from the ideal gas limit are
well described by $\order{g^{2}}$ high-temperature perturbation theory,
although for a detailed quantitative comparison, refined resummed approaches
such as hard thermal loop (HTL) [25] or 3-$d$ effective theory (EQCD) [26] may
be necessary, as discussed in Refs. [27, 28, 29]. In Fig. 1 we show the ideal
gas limit result (red line) and $\order{g^{2}}$ corrections (green band),
using a 2-loop running coupling [30] $g^{2}(T)$ with a renormalization scale
$k_{T}\pi T$ with $4\leq k_{T}\leq 8$ as discussed in Appendix B.

For the higher-order expansion coefficients, $k\geq 3$, the ideal gas limit as
well as $\order{g^{2}}$ corrections vanish, which also is consistent with the
results shown in Fig. 1.

We also note that for vanishing quark masses, the strange quark sector in
strangeness-neutral, isospin-symmetric matter does not contribute to the
perturbative expansion of the pressure up to $\order{g^{2}}$, i.e. although in
this case the strangeness chemical potential is non-zero, $\mu_{S}=\mu_{B}/3$,
the flavor chemical potential for strange quarks vanishes, $\mu_{s}=0$. The
rapid approach to perturbative behavior of the Taylor expansion coefficients
of pressure, and as such also to the $\hat{\mu}_{B}$-dependent contribution to
the pressure, is in contrast to the behavior seen at $\hat{\mu}_{B}=0$, where
it was found that even at $T\simeq 300$ MeV deviations from the ideal gas
limit amount to almost 50% for the pressure and 25% for the energy density
[5].

####  III.2.2 Low-temperature region

At low temperatures lattice QCD results for $P_{2}$ and $P_{4}$ approach
results obtained in hadron resonance gas (HRG) model calculations using non-
interacting, point-like resonances. In Fig. 1 QCD results for the pressure
coefficients are compared to HRG model calculations that use the QMHRG2020
list of hadron resonances [22]. It is apparent that differences between QCD
and HRG model calculations show up earlier and are more pronounced with
increasing order of the expansion coefficients. While $P_{2}/P_{2}^{{\rm
HRG}}$ deviates from unity by about 10% at $T\simeq T_{\rm pc}$, this
deviation reaches already 50% for $P_{4}/P_{4}^{{\rm HRG}}$. Moreover, it is
evident from the comparison of the temperature dependence of $P_{4}$
calculated in QCD and the QMHRG model, respectively, that the slope of
$P_{4}(T)$ differs already significantly for $T>140$ MeV.

Deviations from HRG model calculations thus are even more apparent in the
slope of $P_{2k}(T)$. The $T$-derivatives, $P^{\prime}_{2k}=T{\rm
d}P_{2k}/{\rm d}T$, of the expansion coefficients are shown in Fig. 2. As
expected $P^{\prime}_{2}(T)$ starts deviating from HRG model calculations
already at $T\simeq 140$ MeV and differences are about $30$% in the vicinity
of the pseudo-critical temperature, $T_{\rm pc}$. For $P^{\prime}_{4}(T)$
differences are significant in the entire low-temperature region, $T\geq 125$
MeV.

The functions $P^{\prime}_{2k}(T)$ are the Taylor expansion coefficients of
the $\hat{\mu}_{B}$-dependent part of the trace anomaly,

$\frac{\Delta(\epsilon-3p)}{T^{4}}=\sum_{k}P^{\prime}_{2k}\hat{\mu}_{B}^{2k}\;.$
(27)

We thus can expect that differences between HRG and QCD results will be larger
for energy and entropy densities than for pressure and net baryon-number
density as the former two observables receive contributions from
$P^{\prime}_{2k}(T)$. This will be discussed in Sec. III.3.

####  III.2.3 Intermediate temperature range: Pseudo-critical region

The expansion coefficients shown in Fig. 1, as well as their $T$-derivatives
shown in Fig. 2, make it clear that at temperatures in the vicinity of $T_{\rm
pc}$, deviations from the asymptotic behavior at low and high temperature are
large. As has been noted before, the temperature dependence and the relation
between subsequent expansion coefficients in the vicinity of $T_{\rm pc}$
resemble many features expected from universal scaling in the vicinity of a
second-order phase transition. The maximum in $T{\rm d}P_{2}/{\rm d}T$ is
close to $T_{\rm pc}$. In fact, as $P_{2}(T)$ is an energy-like observable
[31], its derivative with respect to $T$,

$\displaystyle P^{\prime}_{2}=T\frac{{\rm d}P_{2}}{{\rm
d}T}\equiv\left.\frac{T}{2}\frac{\partial^{3}\hat{p}}{\partial
T\partial\hat{\mu}_{B}^{2}}\right|_{\hat{\mu}_{B}=0}\;,$ (28)

behaves like a specific heat, which in the chiral limit will develop a
pronounced peak at the chiral phase transition temperature $T_{\rm c}$. This
can be used as one of the definitions of a pseudo-critical temperature.
However, unlike the magnetization-like susceptibilities, which diverge at
$T_{\rm c}$, the energy-like susceptibility will only lead to a local maximum
at $T_{\rm c}$. Using the maximum of $P^{\prime}_{2}$ as a definition for
$T_{\rm pc}$ thus may be more strongly affected by regular contributions to
the partition function. Its maximum, on the other hand, characterizes the
temperature at which $P^{\prime\prime}_{2}(T)$ vanishes. As the maximum and
minimum of $P^{\prime\prime}_{2}(T)$ on the left and right of this crossing
point approach each other and will diverge in the chiral limit, the crossing
point itself is a well defined estimator for a pseudo-critical temperature.
The temperature at this crossing point will converge to the uniquely defined
$T_{\rm c}$ in the chiral limit. We find for the location of the maximum of
$P^{\prime}_{2}$,

$\displaystyle{\rm d}P^{\prime}_{2}/{\rm d}T=0\;\;\Rightarrow\;\;T_{{\rm
pc},P_{2}}=153.8(7)(5)\ {\rm MeV}\;,\;~{}$ (29)

where the first error is statistical and the second error reflects the
uncertainty in defining the $T$-scale used in all our calculations.

Figure 3: Comparison of the fourth-order Taylor expansion coefficient,
$P_{4}$, of the pressure and the temperature derivative of the second-order
expansion coefficient, $P_{2}$. The red and green bands reflect the errors on
$P^{\prime}_{2}$ and $P_{4}$, respectively. The inset shows results for
$\kappa_{2}^{B}$ obtained from an appropriately scaled ratio of
$P^{\prime}_{2}$ and $P_{4}$ as discussed in the text. From this we find
$0.012\leq\kappa_{2}^{B}\leq 0.022$.

Similarly, the maximum and minimum, visible in $P^{\prime}_{4}$ (Fig. 2 (top,
right)) will diverge in the chiral limit. The temperature at which
$P^{\prime}_{4}$ changes sign thus also is a good observable to define a
pseudo-critical temperature. From the continuum-extrapolated fits shown in
Fig. 2 we thus deduce another estimator for $T_{\rm pc}$,

$\displaystyle P^{\prime}_{4}=0\;\;\Rightarrow\;\;T_{{\rm
pc},P_{4}}=155.3(2)(5)\ {\rm MeV}\;.\;~{}$ (30)

Both estimators for $T_{\rm pc}$, based on temperature-like derivatives of
$p/T$, are in good agreement with results obtained previously by the HotQCD
Collaboration using magnetization-like observables only, i.e. $T_{\rm
pc}=156.5(1.5)$ MeV [19].

The fact that the temperature derivative of $P_{4}$, on the one hand, and the
temperature derivative of the mixed observable $P^{\prime}_{2}$ (involving
$T$\- and $\hat{\mu}_{B}$ derivatives), on the other hand, lead to quite
similar results for the value of a pseudo-critical temperature is naturally
understood in terms of universal behavior in the vicinity of a critical point.
Close to $T_{\rm c}$ higher-order derivatives of the pressure are dominated by
the contribution of the so-called singular part of $p/T$,

$p/T\sim f_{f}(z),$ (31)

where $f_{f}$ is the universal scaling function corresponding to the
universality class of the chiral phase transition444The chiral phase
transition generally is expected to belong to the $3$-$d$, $O(4)$ universality
class, although larger symmetry groups may become relevant, if the anomalous
$U(1)$ symmetry of QCD gets effectively resorted [32].. The argument $z$ of
the scaling function $f_{f}(z)$ and its prefactor generally depend on the
quark masses. However, as we are interested here only in derivatives with
respect to $T$ and $\hat{\mu}_{B}$, keeping the quark masses fixed, we may
think of the argument of the scaling function as just being proportional to
the reduced temperature $t$, i.e. $z\equiv c_{m}t$, with $t$ given by a
leading-order Taylor expansion,

$\displaystyle t=\frac{1}{t_{0}}$ $\displaystyle\left(\frac{T-T_{\rm
c}}{T_{\rm c}}+\kappa_{2}^{B}\hat{\mu}_{B}^{2}\right.$ (32)
$\displaystyle\left.+\order{\hat{\mu}_{B}^{4},\hat{\mu}_{B}^{2}(T-T_{\rm
c}),(T-T_{\rm c})^{2}}\right)\;.$

Here $t_{0}$, as well as the chiral phase transition temperature $T_{\rm c}$
and the curvature coefficient $\kappa_{2}^{B}$, are non-universal parameters.
This suggests that up to a proportionality factor, taking two derivatives with
respect to $\hat{\mu}_{B}$ is equivalent to taking one derivative with respect
to $T$ [33, 34].

The anticipated relation between derivatives with respect to $T$ and
$\hat{\mu}_{B}$, respectively, is apparent in the structure of the expansion
coefficients shown in Figs. 1 and 2; $P^{\prime}_{2}$ and $P_{4}$ receive
singular contributions from $f^{\prime\prime}_{f}(z)$,

$\displaystyle{P_{2}^{s}}^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{2T}{T_{\rm c}}\kappa_{2}^{B}Af_{f}^{\prime\prime}(z)+{\rm
subleading}$ $\displaystyle P_{4}^{s}$ $\displaystyle=$
$\displaystyle-(\kappa_{2}^{B})^{2}Af^{\prime\prime}_{f}(z)+{\rm
subleading}\;,$ (33)

where we summarized common factors in the factor
$A=(h_{0}/2t_{0}^{2}T^{3})~{}h^{-\alpha/\beta\delta}$.

At temperatures in the region between $T_{\rm c}$ and the pseudo-critical
temperature for physical quark mass values, the relative magnitude of
$P^{\prime}_{2}$ and $P_{4}$ is well described by scaling relations. Matching
values of $P^{\prime}_{2}$ and $P_{4}$ in the vicinity of $T_{\rm c}$ thus
ensures that their $T$-derivatives, which will diverge in the chiral limit,
will become identical when approaching this limit.

In Fig. 3 we show the expansion coefficient $P_{4}$ and compare it with an
appropriately rescaled temperature derivative, $P^{\prime}_{2}$, of the
expansion coefficient $P_{2}$. To determine the rescaling factor
$\kappa_{2}^{B}$ we calculate the ratio $2TP_{4}/T_{\rm c}P^{\prime}_{2}$ in
the temperature interval $[130~{}{\rm MeV}:156.5~{}{\rm MeV}]$ and average the
results in this interval. As can be seen in the inset of Fig. 3 these
estimates increase with increasing $T$, indicating the relevance of regular
contributions to the expansion coefficients $P_{2k}$, which lead to deviations
from a $T$-independent result for $\kappa_{2}$ that one would expect to find
in the scaling regime. For the curvature coefficient we thus find
$0.012\leq\kappa_{2}^{B}\leq 0.022$. This result is in good agreement with
other determinations of the curvature coefficient $\kappa_{2}^{B}$ [19, 35,
36] that are based on properties of magnetization-like observables.

###  III.3 Taylor expansion coefficients for energy and entropy densities

Figure 4: Second-order (top, left), fourth-order (top, right), sixth-order
(bottom, left) and eighth-order (bottom, right) expansion coefficients of
pressure ($P_{2k}$), energy density ($\epsilon_{2k}$) and entropy density
($\sigma_{2k}$) for a isospin-symmetric, strangeness-neutral medium
($\mu_{Q}=0,n_{S}=0$). Bands at high temperature show the $\order{g^{2}}$
perturbative result using a renormalization scale in the range $[4\pi T,8\pi
T]$. The yellow band indicates the location of the pseudo-critical temperature
at $T_{pc}$ [19].

As discussed in Sec. II.2 in the case of isospin-symmetric, strangeness-
neutral matter, Taylor expansion coefficients of net baryon-number, energy and
entropy densities are all given in terms of the expansion coefficients
$P_{2k}$ and their temperature derivatives $P^{\prime}_{2k}$. In particular,
the Taylor series for the $\hat{\mu}_{B}$-dependent contribution to the trace
anomaly $(\epsilon-3p)/T^{4}$ is given in terms of the $T$-derivatives,
$P^{\prime}_{2k}$, only. Combining results for $P_{2k}$ and $P^{\prime}_{2k}$,
shown in Figs. 1 and 2, we then obtain the expansion coefficients for Taylor
series of the energy and entropy densities. We show the results for these
expansion coefficients, together with those of the pressure series555Note that
the expansion coefficients for the number density are just proportional to
those of the pressure series and thus are not shown separately in Fig. 4., in
Fig. 4. As discussed for the expansion coefficients $P_{2k}$, also
$\epsilon_{2k}$ and $\sigma_{2k}$ are continuum extrapolated for $k=1,\ 2$ and
we give spline interpolations of the $N_{\tau}=8$ results for $k=3,\ 4$. We
show a comparison with HRG model calculations for the expansion coefficients
of the pressure. For energy and entropy we do so only for the second-order
expansion coefficients. As discussed in connection with Fig. 2, differences
between HRG model and QCD calculations are large for higher-order expansion
coefficients in the entire temperature range .

At high temperatures we again show results from $\order{g^{2}}$ perturbation
theory (PT). In the case of the energy and entropy densities this also
includes contributions arising from $T$-derivatives of the coupling
$g^{2}(T)$, although this might be considered as a part of the $\order{g^{4}}$
contribution in a perturbative expansion. We thus consider our $\order{g^{2}}$
PT-ansatz as a model for describing the high-$T$ behavior of the Taylor
expansion coefficients and subsequently also for the bulk thermodynamic
observables.

Figure 5: Pressure (left) and net baryon-number density (right) versus
temperature for several values of the baryon chemical potential
$\hat{\mu}_{B}\equiv\mu_{B}/T$. Shown are results obtained in different orders
of the Taylor series of pressure and net baryon-number density of isospin-
symmetric, strangeness-neutral matter. The yellow bands highlight the
variation of $\Delta p/T^{4}$ and $n_{B}/T^{3}$ with $\hat{\mu}_{B}$ at
$T_{\rm pc}(\hat{\mu}_{B})$.

##  IV Bulk thermodynamics of (2+1)-flavor QCD

We will use here the Taylor series of the pressure and related observables to
calculate thermodynamic observables in an isospin-symmetric, strangeness-
neutral medium as function of temperature and baryon chemical potential.
Further on, we will eliminate the baryon chemical potential in favor of other
external control parameters like the net baryon-number density or constant
ratio of entropy over net baryon number.

###  IV.1 Density dependent contribution to bulk thermodynamic observables

In [14] we already used the new data obtained in the temperature interval
$T\in[125{\rm MeV}:175{\rm MeV}]$ to construct $8^{\rm th}$-order Taylor
series for the pressure and the related $7^{\rm th}$-order series for the net
baryon-number density. We had shown there that these expansions are reliable
for $\hat{\mu}_{B}\lesssim 2.5$ for the pressure and $\hat{\mu}_{B}\lesssim
2.0$ for the number density, respectively. It also could be shown that these
expansions agree well with corresponding Padé approximants. Here we extend
this discussion up to $T=280$ MeV using existing data from our previous
analysis of Taylor series in QCD [13], where a somewhat larger light quark
mass, $m_{\ell}/m_{s}=1/20$, has been used in the high-temperature region.
Using such a somewhat larger light quark mass at high temperature is of no
concern as $m_{\ell}/T\ll 1$ and has only little influence on thermodynamic
observables in this region. We furthermore present results for the energy and
entropy densities, which require the temperature derivatives of the expansion
coefficients of the pressure.

####  IV.1.1 EoS: Pressure versus net baryon-number density

In Fig. 5 we show the number density (right) and the density dependent
($\mu_{B}\neq 0$) contribution to the pressure (left) for the case
$(\mu_{Q},n_{S})=(0,0)$ and for several values of $\hat{\mu}_{B}$. Similar
results for the case $(\mu_{Q},\mu_{S})=(0,0)$ have been shown in [14]. At
high temperatures, $T\gtrsim 200$ MeV, the Taylor series converges rapidly. As
can be seen from the expansion coefficients presented in Fig. 1, in this
temperature range the $4^{\rm th}$-order expansion coefficient, $P_{4}$, is
more than two orders of magnitude smaller than the $2^{\rm nd}$-order
coefficient, $P_{2}$. Furthermore, within our current statistical accuracy the
$6^{\rm th}$-order coefficient, $P_{6}$, is consistent with zero. Taking into
account the statistical error on $P_{6}$, also its magnitudes is at least two
orders of magnitude smaller than the $4^{\rm th}$-order expansion coefficient.
For $T\gtrsim 200$ MeV we thus show results for the pressure and net baryon-
number density also for larger chemical potentials, e.g. $\hat{\mu}_{B}\simeq
3$. Using even larger values for $\hat{\mu}_{B}$ in this high-temperature
regime seems to be possible.

For $T\leq 200$ MeV contributions from higher-order expansion coefficients
become important. The $6^{\rm th}$ and $8^{\rm th}$-order expansion
coefficients, which are clearly non-zero, show a lot of structure.
Nonetheless, for $T\leq 150$ MeV their contribution to the Taylor series of
pressure and net baryon-number density is small666Note that also in a hadron
resonance gas, where the $\hat{\mu}_{B}$-dependent contribution is
proportional to $\cosh(\hat{\mu}_{B})-1$, the $\order{\hat{\mu}_{B}^{4}}$
Taylor series provides more than 90% of the exact result. Including also the
$\order{\hat{\mu}_{B}^{6}}$ contribution results in agreement with the exact
HRG result to better than 1% for $\hat{\mu}_{B}\leq 2.5$. for
$\hat{\mu}_{B}\leq 2.5$. This is in accordance with estimates for the radius
of convergence $(\mu_{B}/T)_{conv}$ of the Taylor series [14] that suggest
$(\mu_{B}/T)_{conv}>2.5$ also in this temperature range. We also note that in
this temperature range and for $\hat{\mu}_{B}<2.5$ Taylor series results are
consistent within errors with those obtained from the $4^{\rm th}$-order
series. For temperatures in the vicinity of $T_{\rm pc}$, i.e. for $150~{}{\rm
MeV}<T<200$ MeV, contributions from higher-order Taylor expansion
coefficients, however, need to be taken into account already for
$\hat{\mu}_{B}>2$.

This suggests that in the three different temperature intervals $4^{\rm th}$
and $3^{\rm rd}$-order Taylor series provide a good approximation to the
$8^{\rm th}$-order Taylor series for pressure and net baryon-number density,
as long as the baryon chemical potential is smaller than (i) $2.5T$ for $T\leq
150$ MeV, (ii) $2T$ for $150~{}{\rm MeV}<T<200$ MeV and (iii) $3T$ for $T\geq
200$ MeV. In these cases we may invert $p(\mu_{B})$ and $n_{B}(T,\mu_{B})$
exactly and provide a simple analytic expression for $p(T,n_{B})$, which
provides useful insight into the structure of the QCD EoS in a parameter range
that covers almost the entire parameter range accessible with the beam energy
scan (BES-II) at the Relativistic Heavy Ion Collider (RHIC) in collider mode.
In particular, at the two lowest beam energies $\sqrt{s_{{}_{NN}}}=7.7$ GeV
and 11.5 GeV the baryon chemical potentials at the freeze-out temperature
$T_{f}(\mu_{B})$ has been estimated to be about $\mu_{B}/T_{f}=2.8$ and $2.0$,
respectively [37].

Inverting the $\order{\hat{\mu}_{B}^{3}}$ Taylor series for the net baryon-
number density we obtain the chemical potential as function of $T$ and
$n_{B}$,

$\displaystyle\hat{\mu}_{B}(T,n_{B})$ $\displaystyle=$ $\displaystyle
y(T,n_{B})^{1/3}-\frac{N_{1}^{B}(T)}{3N_{3}^{B}(T)}y(T,n_{B})^{-1/3}\;,$ (34)
$\displaystyle y(T,n_{B})$ $\displaystyle=$
$\displaystyle\frac{N_{1}^{B}(T)}{N_{3}^{B}(T)}\left(\frac{\hat{n}_{B}}{2N_{1}^{B}(T)}+\sqrt{\frac{N_{1}^{B}(T)}{27N_{3}^{B}(T)}+\left(\frac{\hat{n}_{B}}{2N_{1}^{B}(T)}\right)^{2}}\right)\;,$

where $N_{1}^{B}(T)=2P_{2}(T)$ and $N_{3}^{B}(T)=4P_{4}(T)$ are the
coefficients of $\order{\hat{\mu}_{B}^{2}}$ and $\order{\hat{\mu}_{B}^{4}}$
terms in the net baryon-number density, which are obtained either from the
Taylor series for net baryon-number density or in $\order{g^{2}}$ high-
temperature perturbation theory.

Inserting these results for the chemical potential in the $4{}^{\rm th}$-order
Taylor series for the pressure we obtain a low-order approximation for the EoS
of strangeness-neutral matter, i.e. $p(T,n_{B})$,

$\displaystyle\frac{\Delta p(T,n_{B})}{T^{4}}$ $\displaystyle=$ $\displaystyle
P_{2}(T)\hat{\mu}_{B}^{2}(T,n_{B})+P_{4}(T)\hat{\mu}_{B}^{4}(T,n_{B})\;,$ (35)
$\displaystyle\;.$

Figure 6: Lines of constant net baryon-number density in the
$\hat{\mu}_{B}$-$T$ plane. Solid bands indicate results obtained by
numerically solving the $\order{\hat{\mu}_{B}^{5}}$ series of $\hat{n}_{B}$
for $\hat{\mu}_{B}$. Dashed lines indicate HRG results while short bands at
high temperature indicate the $\order{g^{2}}$ perturbative result with $4\leq
k_{T}\leq 8$. The yellow band indicates $T_{\rm pc}(\hat{\mu}_{B})$. In the
inset, we show the result from numerically inverting the Taylor series divided
by the analytic inversion of the $\order{\hat{\mu}_{B}^{3}}$ result.

Going beyond the low-order approximation for the EoS we can solve
$\hat{n}_{B}(T,\hat{\mu}_{B})$ numerically for $\hat{\mu}_{B}$ at fixed $T$.
We then obtain $\hat{\mu}_{B}(T,n_{B})$, which can be inserted in the $8^{\rm
th}$-order Taylor series for $p(T,\hat{\mu}_{B})$ and thus allows us to
determine $p(T,n_{B})$.

In Fig. 6 we show results for lines of constant net baryon-number density in
the $\hat{\mu}_{B}$-$T$ plane. Results are shown for several values of
$n_{B}/n_{0}$, with $n_{0}=0.16/{\rm fm}^{3}$ denoting nuclear matter density.
They are compared to QMHRG model calculations at low $T$ and $\order{g^{2}}$
perturbation theory at high $T$, respectively. The inset shows a comparison of
the numerically inverted Taylor series result and the exact inversion of the
$3^{\rm rd}$-order Taylor series for $\hat{n}_{B}(\hat{\mu}_{B})$ given in Eq.
34. Here we show results only in $T$-intervals where $\hat{\mu}_{B}$ stays
below the maximal values given above. As can be seen, the temperature
dependence of $\hat{\mu}_{B}$ characterizing lines of constant $n_{B}$ are
well determined from the $\order{\hat{\mu}_{B}^{3}}$ Taylor series for the net
baryon number-density for $\hat{\mu}_{B}\leq 2.5$.

Inserting $\hat{\mu}_{B}(T,n_{B})$ in the Taylor series for the pressure we
obtain the contribution to the QCD EoS of strangeness-neutral, isospin-
symmetric matter that depends on non-zero net baryon number density. Results
for $\Delta p/(n_{B}T)$ as function of $n_{B}/T^{3}$ for several values of $T$
are shown in Fig. 7 (left). Again we only show results for $\Delta p/(n_{B}T)$
for values of $n_{B}/T^{3}$ which correspond to $\hat{\mu}_{B}\leq 2.5$ for
$T\leq 200$ MeV and $\hat{\mu}_{B}\leq 3$ for $T>200$ MeV. In Fig. 7 (left) we
show corresponding results as function of $n_{B}/n_{0}$.

For $T\gtrsim 240$ MeV it is evident that $\Delta p(T,n_{B})/(Tn_{B})$, rises
almost linearly in $n_{B}/T^{3}$ and is close to the ideal gas result. I.e.
$\Delta p(T,n_{B})$ is proportional to $(n_{B}/T)^{2}$. In this temperature
range the density-dependent part of the QCD EoS is not only close to the ideal
gas result, but also agrees rather well with the $\order{g^{2}}$ perturbative
EoS already for $T\geq 220$ MeV. This is shown by the grey bands in Fig. 7
(right) for which we used the $\order{g^{2}}$ perturbative EoS for the case
$(\mu_{Q},n_{S})=(0,0)$, given in Appendix B, with a renormalization scale
$k_{T}\pi T$, using $4\leq k_{T}\leq 8$.

Figure 7: Pressure divided by net baryon-number density versus $n_{B}/T^{3}$
(left) and $n_{B}/n_{0}$ (right), respectively. Shown are results for
strangeness-neutral, isospin-symmetric matter at several values of $T$. In the
left hand figure we compare results obtained from the full Taylor series for
the pressure with those obtained in $\order{\hat{\mu}_{B}^{4}}$ only (dashed
lines). In the right hand side the grey bands show a comparison with
$\order{g^{2}}$ high-T perturbation theory. The bands shown in both figure are
shown up to values of $n_{B}/T^{3}$ or $n_{B}/n_{0}$ corresponding to
$\hat{\mu}_{B}=2.5$ for $T<200$ MeV and $\hat{\mu}_{B}=3$ otherwise.

From Eqs. 34 and 35 we expect at low densities to find $\Delta
p\sim(n_{B}/T)^{2}$. The notion of low density, however, is temperature
dependent. As can be seen from Eq. 34, the leading order coefficient
$N_{1}^{B}(T)$ sets the scale for $n_{B}$. As $N_{1}^{B}(T)$ drops
exponentially at low temperature one quickly enters a region of “high density”
already at densities smaller than nuclear matter density, $n_{0}\simeq
0.16{\rm fm}^{-3}$, whereas at high temperature the low density regime
($\Delta p\sim(n_{B}/T)^{2}$) persists into a region well above $n_{0}$.

At low densities we thus expect that the density dependent contribution
($\Delta p$) to the total pressure is proportional to $(n_{B}/T)^{2}$ while
the $n_{B}$-dependence weakens. At high densities and, within the
approximation given by Eq. 35, one finds $\Delta p\sim n_{B}^{4/3}$.

Results for $\Delta p(T,n_{B})/n_{B}$, obtained with this
$\order{\hat{\mu}_{B}^{4}}$ approximation (dashed lines), are compared with
the full Taylor series results (bands) in Fig 7 (left). The bands are shown up
to a maximal density that can be reliably reached with our current $8^{\rm
th}$-order Taylor series in different temperature regions, i.e. we demand that
$\hat{\mu}_{B}\leq 2.5$ for $T\leq 200$ MeV and $\hat{\mu}_{B}\leq 3$ at
higher temperatures. As can be seen the dashed lines indeed provide a good
approximation for the EoS at temperatures $135~{}{\rm MeV}<T\leq 150$ MeV and
for $T\geq 200$ MeV. Only for $150~{}{\rm MeV}<T<200~{}{\rm MeV}$ we find that
an EoS based on the $4^{\rm th}$-order Taylor series is not sufficient. Here
differences become visible in a region corresponding to $\hat{\mu}_{B}>2$.

####  IV.1.2 Energy and entropy densities

The calculation of Taylor series for the energy and entropy densities requires
the derivatives of Taylor expansion coefficients of the pressure with respect
to temperature, $P^{\prime}_{2k}\equiv T{\rm d}P_{2k}/{\rm d}T$. These
expansion coefficients are shown in Fig. 2. They directly give the Taylor
series for the $\hat{\mu}_{B}$-dependent contribution to the trace anomaly, as
introduced in Eq. 27. In Fig. 8 we show $\Delta(\epsilon-3p)/T^{4}$ for
several values of the baryon chemical potential. This should be compared to
the trace anomaly at $\hat{\mu}_{B}=0$ ( Fig. 5 (right) of [5]). We note that
at $T_{\rm pc}$ and for $\hat{\mu}_{B}=2.5$ the fermionic contribution to the
trace anomaly is as large as the contribution at $\hat{\mu}_{B}=0$ but rapidly
drops at larger values of the temperature. For $T\geq 200$ MeV the
$\hat{\mu}_{B}$-dependent contribution to the trace anomaly stays below 10% of
the value of the trace anomaly at $\hat{\mu}_{B}=0$.

A consequence of this is that the maximum in the trace anomaly, which at
$\hat{\mu}_{B}=0$ is reached at $T\simeq 200$ MeV [5], gets shifted to smaller
temperatures and comes closer to the pseudo-critical temperature. In fact, as
can be seen in Fig. 8, the location of the maximum of
$\Delta(\epsilon-3p)/T^{4}$ is close to the pseudo-critical temperature
determined from chiral observables,

$T_{\rm pc}(\hat{\mu}_{B})=T_{\rm
pc}^{0}\left(1-\kappa_{2}^{B,f}\hat{\mu}_{B}^{2}+\order{\hat{\mu}_{B}^{4}}\right)\;,$
(36)

with $T_{pc}^{0}=156.5(1.5)$ MeV and $\kappa_{2}^{B,f}=0.012(4)$ denoting the
curvature coefficient of the pseudo-critical line in strangeness-neutral
matter777For a summary of different determinations of the curvature
coefficient $\kappa_{2}^{B}$, using various external constraints, see for
instance [38, 39]. [19]. For large values of the chemical potential, where
$\Delta(\epsilon-3p)/T^{4}$ is dominated by the $\hat{\mu}_{B}$-dependent
contribution, violations of the conformal relation $\epsilon=3p$ thus are
maximal on the pseudo-critical line.

Figure 8: The $\hat{\mu}_{B}$-dependent contribution to the trace anomaly in
$(2+1)$-flavor QCD for several values of $\hat{\mu}_{B}$. The yellow band
shows the line $\Delta((\epsilon-3p)/T^{4})(T_{\rm pc}(\hat{\mu}_{B}))$.

Figure 9: $\hat{\mu}_{B}$-dependent part of the energy (left) and entropy
(right) densities versus temperature. Shown are results from $2^{\rm nd},\
4^{\rm th},\ 6^{\rm th}$ and $8^{\rm th}$-order Taylor expansions. The yellow
bands indicate how $\Delta\epsilon/T^{4}$ and $\Delta s/T^{3}$ change with
$\hat{\mu}_{B}$ at $T_{\rm pc}$.

Combining the results for $\Delta(\epsilon-3p)/T^{4}$ (Fig. 8) and $\Delta
p/T^{4}$ (Fig. 5) we obtain the Taylor series results for the
$\hat{\mu}_{B}$-dependent contributions to $\epsilon/T^{4}$ and $s/T^{3}$,
respectively. These observables are shown in Fig. 9. As discussed in the case
of $\Delta p/T^{4}$ and $n_{B}/T^{3}$, we also find for
$\Delta(\epsilon-3p)/T^{4}$, as well as for $\Delta\epsilon/T^{4}$ and $\Delta
s/T^{3}$, that the $4^{\rm th}$-order Taylor expansions provide good
approximations to the full $6^{\rm th}$ and $8^{\rm th}$-order series results
for $T\geq 200$ MeV and $T\leq 150$ MeV whenever $\hat{\mu}_{B}\leq 2.5$ in
the low-$T$ region and $\hat{\mu}_{B}\leq 3$ at high $T$. Higher-order
corrections become more important in the calculation of the entropy density at
high temperature. As can be seen in Fig. 9 (right) the $6^{\rm th}$ and
$4^{\rm th}$-order correction to $s/T^{3}$ agree with each other in the region
$T>200$ MeV only up to $\hat{\mu}_{B}\simeq 2.5$ but start to differ for
larger $\hat{\mu}_{B}$.

As can be seen in Fig. 9 for large $\hat{\mu}_{B}$ the $4^{\rm th}$-order
Taylor expansion results for $\Delta\epsilon/T^{4}$ and $\Delta s/T^{3}$ show
a wiggly behavior in the intermediate temperature interval, $150~{}{\rm
MeV}<T<200~{}{\rm MeV}$, which results from the pronounced minimum of
$P^{\prime}_{4}$ in this temperature interval. This structure gets smoothed
out by higher-order corrections. We will show in the next subsection that a
Padé-resummed Taylor series shows a much smoother behavior in this temperature
range.

Figure 10: Energy density in units of $T^{4}$ divided by net baryon-number
density versus $n_{B}/T^{3}$ (left) and the corresponding figure for entropy
density (right). Shown are results for strangeness-neutral, isospin-symmetric
matter at several values of $T$ using $6^{\rm th}$-order Taylor expansions for
$\epsilon/T^{4}$ and $s/T^{3}$ and a $5^{\rm th}$-order Taylor expansions for
$n_{B}/T^{3}$. Dashed lines show a comparison with results based on $4^{\rm
th}$ and $3^{\rm rd}$-order expansions, respectively, using Eq. 34 to obtain
$\epsilon(n_{B})$ and $s(n_{B})$.

We use the $6^{\rm th}$-order Taylor expansion results for $\Delta p/T^{4}$,
$\Delta\epsilon/T^{4}$ and $\Delta s/T^{3}$ to obtain the latter two as
functions of $T$ and $n_{B}$. We again compare with corresponding $4^{\rm
th}$-order approximations using,

$\displaystyle\frac{\Delta\epsilon(T,\hat{n}_{B})}{T^{4}}$ $\displaystyle=$
$\displaystyle\epsilon_{2}(T)\hat{\mu}_{B}^{2}(T,n_{B})+\epsilon_{4}(T)\hat{\mu}_{B}^{4}(T,n_{B})\;,$
(37) $\displaystyle\frac{\Delta s(T,\hat{n}_{B})}{T^{3}}$ $\displaystyle=$
$\displaystyle\sigma_{2}(T)\hat{\mu}_{B}^{2}(T,n_{B})+\sigma_{4}(T)\hat{\mu}_{B}^{4}(T,n_{B})\;,$
(38)

with $\hat{\mu}_{B}(T,\hat{n}_{B})$ taken from Eq. 34. Results for
$\Delta\epsilon(T,\hat{n}_{B})/n_{B}$ and $\Delta s(T,\hat{n}_{B})/n_{B}$ as
functions of $\hat{n}_{B}$ are shown in Fig 10.

####  IV.1.3 Comparison of Taylor series and their resummation using Padé
approximants

As seen already in the analysis of the Taylor expansion for the pressure and
net baryon-number density, Padé approximants agree well with the Taylor series
themselves at low $\hat{\mu}_{B}$ up to the region where we estimated the
latter to provide reliable results [14]. We will extend this approach here to
the analysis of Taylor series for the energy and entropy densities.

We use Padé approximants for thermodynamic observables derived from the Taylor
series of the pressure in two ways. On the one hand we construct Padé
approximants based on the Taylor series for a given observable, e.g. the
energy and entropy density series given in Eqs. 13 and 14 can be resummed
using Padé approximants similar to that of the pressure series given in Eq. 23
by just replacing the expansion coefficients $P_{2k}$ by $\epsilon_{2k}$ or
$\sigma_{2k}$, respectively. On the other hand we use the P-Padé , i.e.
appropriate derivatives of the Padé approximants for the pressure, for the
energy and entropy densities as given in Eqs. 24 and 25.

Figure 11: Comparison of $6^{\rm th}$ (top) and $8^{\rm th}$ (bottom) order
Taylor expansion results for $\Delta p/T^{4}$ and $\Delta\epsilon/T^{4}$ with
corresponding [2,4] and [4,4] Padé and P-Padé approximants. The inset in the
bottom panels shows differences in the [2,4] and [4,4] as well as in [4,2] and
[4,4] Padé approximations, which are all based on $8^{\rm th}$-order Taylor
series. The yellow bands highlight the variation of $\Delta p/T^{4}$ and
$\Delta\epsilon/T^{4}$ with $\hat{\mu}_{B}$ at $T_{\rm pc}(\hat{\mu}_{B})$.

In Fig. 11 (left) we compare $6^{\rm th}$ and $8^{\rm th}$-order Taylor series
for $\Delta\hat{p}$ with corresponding [4,2], [2,4] and [4,4] Padé
approximants introduced in Eqs. 20-23. Corresponding results for
$\Delta\hat{\epsilon}$ are shown in Fig. 11 (right). In the figure for the
pressure (top, left) we compare the $6^{\rm th}$-order Taylor expansion
results with the two possible $[n,m]$ Padé approximants that use up to
$6{}^{\rm th}$-order expansion coefficients. As can be seen the [4,2] Padé
agrees with the Taylor series result while the [2,4] Padé differs from these
two in the temperature interval $150~{}{\rm MeV}\lesssim T\lesssim 180~{}{\rm
MeV}$. In the (bottom, left) figure shows a comparison of the $8^{\rm
th}$-order Taylor expansion results with the [4,4] Padé approximant. They are
in good agreement with each other. Moreover, as can be seen from the inset,
for large values of $\hat{\mu}_{B}$ the [2,4] and [4,4] Padé approximants stay
in much better agreement with each other than the [4,2] and [4,4] Padé
approximants or, equivalently, the $6^{\rm th}$ and $8^{\rm th}$-order Taylor
series. Similar conclusions can be drawn for the energy density shown in the
right hand part of the figure.

In Fig. 12 we compare at fixed values of the temperature $8^{\rm th}$-order
Taylor expansion results with the [4,4] Padé approximant as well as the P-Padé
results for pressure, energy and entropy densities that we have discussed
above and in Sec. II.3. We generally find that the P-Padé results for bulk
thermodynamic observables are in better agreement with the Taylor series
results than the [4,4] or [3,4] Padé approximants. This may not be too
surprising, as both approaches are based on a thermodynamically self-
consistent set of expressions.

We consider results for bulk thermodynamic observables based on the $8^{\rm
th}$-order Taylor series for the pressure reliable when the Taylor series
results agree with the [4,4] Padé approximants as well as the P-Padé
approximants shown in Fig. 12. These parameter ranges are similar for all four
bulk thermodynamic observables shown in Fig. 12 and also agree with current
estimates for a radius of convergence of the Taylor series given in [14]. We
thus conclude that the region of reliability of current EoS results is
slightly temperature dependent, increasing from $\hat{\mu}_{B}\simeq 2.2$ at
$T=135~{}{\rm MeV}$ to $\hat{\mu}_{B}\simeq 3.2$ at $T=180~{}{\rm MeV}$. The
latter also holds for higher values of the temperature where our results only
are based on $6^{\rm th}$-order Taylor series expansions.

Figure 12: Top to bottom: Comparison of Taylor series and Padé approximants
for pressure, net baryon-number, energy and entropy densities, at three values
of the temperature, $T=135$ MeV (left), 155 MeV (middle) and 175 MeV (right).
Taylor expansion results are shown for $8^{\rm th}$ and $7^{\rm th}$-order
(number density), respectively. The corresponding Padé approximants are [4,4]
and [3,4]. Also shown are P-Padé approximants (P-Pade) which are appropriate
$T$\- and $\mu_{B}$-derivatives of the [4,4] Padé approximant for the pressure
as discussed in Sec. II.3.

###  IV.2 Bulk thermodynamic observables

Having discussed the $\mu_{B}$-dependent contribution to bulk thermodynamic
observables in the previous subsection we can now combine these results with
those obtained at vanishing values of the chemical potential, for which we use
the continuum-extrapolated results obtained for $(2+1)$-flavor QCD by the
HotQCD collaboration [5].

In Fig. 13 we show the total pressure (left) as well as energy (middle) and
entropy (right) densities for $\hat{\mu}_{B}\in[0:2.5]$ in the entire
temperature range analyzed by us. For $\hat{\mu}_{B}=3$ we show results in the
region $T\geq 180$ MeV only.

Also shown in Fig. 13 are results from HRG model calculations for
$\hat{\mu}_{B}=0$ and $2.5$. It is apparent that the HRG model results show a
stronger $\hat{\mu}_{B}$-dependence than the QCD data. Already in the vicinity
of $T_{\rm pc}$ differences between HRG model calculations and QCD increase
with increasing values of $\hat{\mu}_{B}$. This is to be expected as higher-
order Taylor expansion coefficients start to differ from HRG model results
already at rather low temperatures, $T<140$ MeV, and differences get large in
the transition region (see Fig. 4). These differences are larger for the
energy and entropy densities than for pressure and number density. In
particular, we find that the $\order{\hat{\mu}_{B}^{2}}$ coefficients in QCD
are about 30% smaller than in HRG model calculations (see Fig. 4) and the
ratio of $\order{\hat{\mu}_{B}^{4}}$ and $\order{\hat{\mu}_{B}^{2}}$ expansion
coefficients of the energy density calculated in lattice QCD and HRG models,
respectively, differs by almost a factor 5, i.e. $\epsilon_{4}(T_{\rm
pc})/\epsilon_{2}(T_{\rm pc})=0.013(2)(6)$ in QCD compared to
$(\epsilon_{4}(T_{\rm pc})/\epsilon_{2}(T_{\rm pc}))_{\rm HRG}=0.065(1)$ in
the HRG model, using the QMHRG2020 hadron list. As a consequence, with
increasing $\hat{\mu}_{B}$ the difference between the energy density
calculated in QCD and HRG models increases on the pseudo-critical line,
$T_{\rm pc}(\hat{\mu}_{B})$.

As discussed in [13], the curvature coefficients $\kappa_{2}^{\epsilon}$ and
$\kappa_{2}^{s}$ for lines of constant energy and entropy densities,
respectively, are related through the second-order expansion coefficients
$\sigma_{2}$ and $s_{2}$,

$\displaystyle\frac{\kappa_{2}^{s}}{\kappa_{2}^{\epsilon}}$ $\displaystyle=$
$\displaystyle\frac{\sigma_{2}}{\epsilon_{2}}=1-\frac{P_{2}}{\epsilon_{2}}=0.872(3)(5)\;.$
(39)

Lines of constant energy density thus are flatter than those for constant
entropy density, i.e. the entropy decreases on a line of constant energy
density. Using the coefficients $\epsilon_{2}$ and $\sigma_{2}$, shown in Fig.
4, we get

$\displaystyle\kappa_{2}^{\epsilon}$ $\displaystyle=$ $\displaystyle
0.0104(12)\;\;,\;\;\kappa_{2}^{s}=0.0091(11)\;,$ (40)

which should be compared to the curvature coefficient $\kappa_{2}^{B,f}$ of
the pseudo-critical line given in Eq. 36. These results thus suggest that,
within current statistical errors, energy and entropy densities stay constant
on the pseudo-critical line or drop slightly.

Figure 13: Pressure (left) as well as energy (middle) and entropy (right)
densities versus temperature for several values of the baryon chemical
potential. Figures show results for the case $n_{S}=0$, $n_{Q}/n_{B}=0.5$ in
the temperature interval [130 MeV:280 MeV]. Figures on the top show results
from $6^{\rm th}$ order Taylor series and the figures on the bottom have been
used using Padé and P-Padé approximants, respectively. The Taylor expansions
are based on the continuum and spline interpolated data shown in Fig. 1 for
the temperature interval [130 MeV:175 MeV]. At larger temperatures we used
data from our earlier analysis of bulk thermodynamics in a sixth-order Taylor
expansion [13]. Results for $\hat{\mu}_{B}=0$ have been taken from [5].

Using the results for $\epsilon/T^{4}$ shown in Fig. 13 (middle) we find for
the energy density on the pseudo-critical line,

$\epsilon(T_{\rm pc}(\hat{\mu}_{B}))=\begin{cases}370(40)(30)~{}{\rm
MeV/fm}^{3}\;,\;\hat{\mu}_{B}=0\\\ 330(28)(53)~{}{\rm
MeV/fm}^{3}\;,\;\hat{\mu}_{B}=2.5\end{cases}$ (41)

Here the first error is the statistical error on $\epsilon$ at $T_{\rm
pc}(\hat{\mu}_{B})$, while the second error reflects the systematic
uncertainty arising from the uncertainty on $T_{\rm pc}(\hat{\mu}_{B})$, which
rises from 1.5 MeV at $\hat{\mu}_{B}=0$ to 4 MeV at $\hat{\mu}_{B}=2.5$.

Figure 14:  Lines of constant entropy per baryon number in the
$T$-$\hat{\mu}_{B}$ plane. Solid bands indicate results obtained by
numerically solving $s/n_{B}$ derived from the $\order{\hat{\mu}_{B}^{6}}$
pressure series for $\hat{\mu}_{B}$. Dashed lines indicate HRG results while
the (almost straight vertical bands indicate the $\order{g^{2}}$ perturbative
result with $4\leq k_{T}\leq 8$. The yellow band indicates $T_{\rm
pc}(\hat{\mu}_{B})$. In the inset, we show the result from this numerical
inversion divided by the corresponding inversion coming from the
$\order{\hat{\mu}_{B}^{4}}$ pressure series.

##  V Thermodynamics at constant ratio of entropy to net baryon number

Figure 15:  Total pressure (left), energy density (middle) and entropy density
(right) versus temperature for several values of the ratio $s/n_{B}$.

Equilibrated strong interaction matter created, for instance in heavy ion
experiments, follows lines of constant entropy per net baryon number while
expanding and cooling after the initial collision of nuclei. The composition
of the initially colliding nuclei thus fixes the conserved charge content of
the matter created in this collisions, i.e. $n_{S}=0$ and $n_{Q}/n_{B}$ fixed.

Also for the calculation of thermodynamic observables such as the speed of
sound or the adiabatic compressibility of matter, one needs to determine the
lines of constant $s/n_{B}$ in QCD [12, 11]. We have determined the lines of
constant $s/n_{B}$ in a strangeness-neutral medium with fixed $n_{Q}/n_{B}$.

In Fig. 14 we show such lines in the $\hat{\mu}_{B}$-$T$ plane for $s/n_{B}$
and fixed ($n_{S},n_{Q}/n_{B})=(0,0.5)$ in the range $25\leq s/n_{B}\leq 400$.
This roughly corresponds to the range covered by BES-II at RHIC in the range
of beam energies $7.7~{}{\rm GeV}\leq\sqrt{s_{{}_{NN}}}\leq 200~{}{\rm GeV}$.
Also shown in Fig. 14 is the asymptotic behavior of lines of constant $n_{B}$
at high temperature and the approach to the HRG model at low temperature. As
can be seen, QCD results for lines of constant $s/n_{B}$ differ much more from
high-$T$ perturbation theory then the lines of constant $n_{B}$. This, of
course, arises from the $\mu_{B}=0$ contribution to the total entropy density,
which approaches the ideal gas limit slowly. The inset in Fig. 14 shows that
these lines are already well determined using $4^{\rm th}$ order Taylor series
only.

In Fig. 15 we show the temperature dependence of pressure, energy and entropy
densities on lines of constant $s/n_{B}$. Obviously, the $T$-dependence of
these observables is similar to that obtained for fixed $\hat{\mu}_{B}$. Above
the pseudo-critical temperature, however, these observables rise somewhat
faster than at $\mu_{B}=0$, i.e. at $s/n_{B}=\infty$. This is expected as Fig.
14 shows that the value of the chemical potential on lines of constant
$s/n_{B}$ increases above $T_{\rm pc}$ with increasing temperature.

The speed of sound in strangeness-neutral matter with fixed $n_{Q}/n_{B}$
reflects the temperature dependence of $p$ and $\epsilon$, shown in Fig. 15.
As defined in Eq. 5 we have,

$\displaystyle c_{s}^{2}$ $\displaystyle=$ $\displaystyle\left(\frac{\partial
p}{\partial\epsilon}\right)_{s/n_{B},n_{Q}/n_{B},n_{S}}$ (42) $\displaystyle=$
$\displaystyle\frac{(\partial p/\partial
T)_{s/n_{B},n_{Q}/n_{B},n_{S}}}{(\partial\epsilon/\partial
T)_{s/n_{B},n_{Q}/n_{B},n_{S}}}\;.$

The numerator is related to the entropy density evaluated at fixed
$\vec{X}=(s/n_{B},n_{Q}/n_{B},n_{S})$, while the denominator is a specific
heat defined at fixed $\vec{X}$. We note, however, that the temperature
derivatives taken at fixed $\vec{X}$ differ, of course, from those taken at
fixed $\vec{\hat{\mu}}$. Some relations for constrained partial derivatives
are given in Appendix A, and in Appendix C we give an explicit expression for
$c_{s}^{2}$.

For $n_{S}=0$ we obtain from Eqs. 5 and 6 for the speed of sound ($c_{s}^{2}$)
and the closely related adiabatic compressibility ($\kappa_{s}$)

$\displaystyle c^{2}_{s}$ $\displaystyle=$
$\displaystyle\frac{N_{\vec{X}}}{(\epsilon+p)\ D_{\vec{X}}}\;,\;$
$\displaystyle\kappa_{s}$ $\displaystyle=$
$\displaystyle\frac{1}{c_{s}^{2}(\epsilon+p)}=\frac{D_{\vec{X}}}{N_{\vec{X}}}\;.$
(43)

The functions $N_{\vec{X}}(T,\vec{\mu})$ and $D_{\vec{X}}(T,\vec{\mu})$ are
given in terms of second-order cumulants of the conserved charges $(B,\ Q,\
S)$, the entropy density and its derivatives with respect to $T$ as well as
the three chemical potentials. All these observables themselves are functions
of $T$ and $\vec{\mu}$. One thus may evaluate $c_{s}^{2}$ by either taking
numerically the $T$-derivatives appearing in Eq. 42 or use directly Eqs.
62-65. The results shown in Fig. 16 are based on the former approach.

As can be seen in Fig. 15 (left) and (middle) the shape of $p(T)$ and
$\epsilon(T)$ varies little when changing $s/n_{B}$ and follows a similar
pattern in both quantities. In the temperature range currently accessible to
us we thus do not expect to observe a strong dependence of the speed of sound
and the adiabatic compressibility on $s/n_{B}$. This is indeed apparent from
the results for both observables, calculated for several fixed ratios
$s/n_{B}$, as shown in Fig. 16. At low temperatures we show results for
$c_{s}^{2}$ only down to $T$-values that can be reached at fixed $s/n_{B}$ for
$\hat{\mu}_{B}<2.5$. As can be seen in Fig. 14 this allows to explore the
entire available temperature range ($T>135$ MeV) for $s/n_{B}\geq 30$, but
limits the temperature range for smaller values of $s/n_{B}$.

For temperatures $T\gtrsim T_{\rm pc}$ the speed of sound slightly increases
with decreasing $s/n_{B}$ and steadily increases with $T$. At temperatures
just below $T_{\rm pc}$ the opposite trend is expected to show up in HRG model
calculations: $c_{s}^{2}$ decreases with increasing $s/n_{B}$ as well as $T$.
At low temperatures, the speed of sound is known to have a “bump” arising from
the interplay of baryon and light meson contributions to $c_{s}^{2}$. This
leads to the minimum in $c_{s}^{2}$ observed close to $T_{\rm pc}$. As can be
seen in the inset of Fig. 16 (left) the bump at low $T$ becomes more shallow
for smaller values of $s/n_{B}$. In fact, depending on details of the baryon
component of the HRG model, model calculations suggest that $c_{s}^{2}$ will
stop developing a bump at low temperatures at $s/n_{B}\simeq 15$. This is
apparent from the inset in Fig. 16 (left) where we show some results on the
speed of sound calculated in a HRG model using the QMHRG2020 list of hadrons.

As $c_{s}^{2}$ as well as $\epsilon$ and $p$ increase with decreasing
$s/n_{B}$ it is clear that the adiabatic compressibility shown in Fig. 16
(right) will decrease with decreasing $s/n_{B}$. Moreover we find that it
monotonically decreases with increasing temperature.

Figure 16: Speed of sound (left) and adiabatic compressibility (right) in
strangeness-neutral, isospin-symmetric matter versus temperature. Shown are
results for several values of $s/n_{B}$. The limit $s/n_{B}=\infty$
corresponds to the case of vanishing chemical potentials. Dashed lines at low
temperatures indicate QMHRG2020 model calculations, at high temperatures they
show the non-interacting quark-gluon gas results. In the inset HRG model
calculations at lower temperatures are shown. The yellow band indicates
$T_{\rm pc}$.

##  VI Conclusions

We presented an update on the Taylor expansion of basic thermodynamic
observable (pressure, energy and entropy densities) of strangeness-neutral,
isospin-symmetric strong-interaction matter. We showed that the
$\hat{\mu}_{B}$-dependent parts of these observables approach the high-
temperature perturbative regime much more rapidly than the contribution to
these observables at $\hat{\mu}_{B}=0$. The $\hat{\mu}_{B}$-dependent
contributions are well described by $\order{g^{2}}$ perturbation theory
already for $T\gtrsim 250$ MeV.

We furthermore constructed the QCD equation of state ($P(n_{B})$,
$\epsilon(n_{B})$). We showed that the Taylor series can be resummed using
Padé approximants. The regime where Taylor expansions and Padé approximants
agree with each other gives further confidence in the validity range of
current Taylor expansion results varying from $\hat{\mu}_{B}\simeq 2.5$ at low
temperatures to $\hat{\mu}_{B}\gtrsim 3$ at $T>200$ MeV.

We showed that $\order{\hat{\mu}_{B}^{4}}$ Taylor series results for the
pressure provide a good description of the EoS in the entire parameter range
currently accessible in the beam energy scan performed at RHIC in collider
mode. This allows for an analytic ansatz for the $\hat{\mu}_{B}$-dependent
part of the EoS of strangeness-neutral, isospin-symmetric matter, valid at
temperatures $T\gtrsim 135$ MeV and $\mu_{B}/T\lesssim 2.5$, which requires as
input only the two leading-order Taylor expansion coefficients ($P_{2}(T),\
P_{4}(T))$ of the pressure series.

For the first time, we presented results for the speed of sound in
strangeness-neutral matter at non-zero net baryon-number density. We showed
that the minimum in $c_{s}^{2}$, present at $\hat{\mu}_{B}=0$ and for
$s/n_{B}>100$ at temperatures close to the pseudo-critical line, will
disappear at smaller values of $s/n_{B}$, i.e for large $n_{B}$ in the
vicinity of $T_{\rm pc}$.

_Acknowledgments.—_ This work was supported by: (i) The U.S. Department of
Energy, Office of Science, Office of Nuclear Physics through the Contract No.
DE-SC0012704; (ii) The U.S. Department of Energy, Office of Science, Office of
Nuclear Physics within the framework of Scientific Discovery through Advance
Computing (SciDAC) award Fundamental Nuclear Physics at the Exascale and
Beyond; (iii) The Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) - Project numbers 315477589-TRR 211 and “NFDI 39/1” for the
PUNCH4NFDI consortium; (iv) The grant 05P2018 (ErUM-FSP T01) of the German
Bundesministerium für Bildung und Forschung; (v) The grant 283286 of the
European Union. D.B. was supported by the Intel Corporation.

This research used awards of computer time provided by: (i) The INCITE program
at Oak Ridge Leadership Computing Facility, a DOE Office of Science User
Facility operated under Contract No. DE-AC05-00OR22725; (ii) The ALCC program
at National Energy Research Scientific Computing Center, a U.S. Department of
Energy Office of Science User Facility operated under Contract No. DE-
AC02-05CH11231; (iii) The INCITE program at Argonne Leadership Computing
Facility, a U.S. Department of Energy Office of Science User Facility operated
under Contract No. DE-AC02-06CH11357; (iv) The USQCD resources at the Thomas
Jefferson National Accelerator Facility.

This research also used computing resources made available through: (i) a
PRACE grant at CINECA, Italy; (ii) the Gauss Center at NIC-Jülich, Germany;
(iii) the GPU-cluster at Bielefeld University, Germany.

We thank the Bielefeld HPC.NRW team for their support and Dietrich Boedeker
for very helpful discussions.

##  Appendix A Constrained partial derivatives

We summarize here relations for partial derivatives of thermodynamic
observables with respect to temperature, keeping specific external conditions
$(x,y,z)$ fixed,

For any thermodynamic function $f(T,\mu_{B},\mu_{Q},\mu_{S})$ we have

$\displaystyle\left(\frac{\partial f}{\partial T}\right)_{(x,y,z)}$
$\displaystyle=$ $\displaystyle\left(\frac{\partial f}{\partial
T}\right)_{(\mu_{B},\mu_{Q},\mu_{S})}$ $\displaystyle+\left(\frac{\partial
f}{\partial\mu_{B}}\right)_{(T,\mu_{Q},\mu_{S})}\left(\frac{\partial\mu_{B}}{\partial
T}\right)_{(x,y,z)}$ $\displaystyle+\left(\frac{\partial
f}{\partial\mu_{Q}}\right)_{(T,\mu_{B},\mu_{S})}\left(\frac{\partial\mu_{Q}}{\partial
T}\right)_{(x,y,z)}$ $\displaystyle+\left(\frac{\partial
f}{\partial\mu_{S}}\right)_{(T,\mu_{B},\mu_{Q})}\left(\frac{\partial\mu_{S}}{\partial
T}\right)_{(x,y,z)}\;.$

Similarly one has for two thermodynamic functions
$f(T,\mu_{B},\mu_{Q},\mu_{S})$ and $g(T,\mu_{B},\mu_{Q},\mu_{S})$ the relation

$\left(\frac{\partial f}{\partial g}\right)_{(x,y,z)}=\frac{\left(\partial
f/\partial T\right)_{(x,y,z)}}{\left(\partial g/\partial T\right)_{(x,y,z)}}$
(45)

In Eqs. A and 45 the derivatives of the chemical potentials are taken on lines
of constant $x(T,\mu_{B},\mu_{Q},\mu_{S})$, $y(T,\mu_{B},\mu_{Q},\mu_{S})$ and
$z(T,\mu_{B},\mu_{Q},\mu_{S})$ in the space of external parameters
$(T,\mu_{B},\mu_{Q},\mu_{S})$. In the lattice QCD context we usually work in
the parameter space ($T,\hat{\mu}\equiv\mu/T$). Moreover, we conveniently work
with reduced, i.e. dimensionless, thermodynamic observables, i.e. we want to
replace e.g. $\hat{\epsilon}=\epsilon/T^{4}$, etc.

Changing the partial derivatives $\partial_{\mu_{B}}$ to
$\partial_{\mu_{B}/T}$ and introducing reduced observables is straightforward,
as these derivatives are taken at fixed $T$. We have for an observable that
has dimension of $T^{n}$ the relation,

$\left.\frac{\partial
f}{\partial\mu_{B}}\right|_{T}=\left.T^{n-1}\frac{\partial\hat{f}}{\partial\hat{\mu}_{B}}\right|_{T}\;.$
(46)

Rewriting the temperature derivatives one has to be a bit more careful,

$\displaystyle\left.\frac{\partial f}{\partial T}\right|_{\vec{\mu}}$
$\displaystyle=$ $\displaystyle\left.\frac{\partial f}{\partial
T}\right|_{\vec{\hat{\mu}}}-\sum_{i=B,Q,S}\frac{\mu_{i}}{T^{2}}\left.\frac{\partial
f}{\partial\mu_{i}/T}\right|_{T}$ $\displaystyle=$
$\displaystyle\left.\frac{\partial T^{n}\hat{f}}{\partial
T}\right|_{\vec{\hat{\mu}}}-\sum_{i}\frac{\mu_{i}}{T^{2}}T^{n}\left.\frac{\partial\hat{f}}{\partial\hat{\mu}_{i}}\right|_{T}$
$\displaystyle=$ $\displaystyle
T^{n-1}\left(n\hat{f}+T\frac{\partial\hat{f}}{\partial
T}\Bigg{|}_{\vec{\hat{\mu}}}-\sum_{i}\hat{\mu}_{i}\frac{\partial\hat{f}}{\partial\hat{\mu}_{i}}\Bigg{|}_{T}\right)\;.$

##  Appendix B Ideal gas and high-temperature perturbation theory for isospin-
symmetric, strangeness-neutral matter

At large values of the temperature, the pressure approaches that of a massless
ideal gas of quarks and gluons. In the case of massless quarks, the leading
order, ideal gas result as well as the $\order{g^{2}}$ correction are second-
order polynomials in $\hat{\mu}_{f}^{2}$ [40, 41],

$p/T^{4}=\hat{p}_{id}-g^{2}\hat{p}_{2}\;,$ (48)

with the ideal gas term

$\displaystyle\hat{p}_{id}$ $\displaystyle=$
$\displaystyle\frac{8\pi^{2}}{45}\left(1+\frac{21n_{f}}{32}\right)$ (49)
$\displaystyle\hskip
14.22636pt+\sum_{f=u,d,s}\left[\frac{1}{2}\left(\frac{\mu_{f}}{T}\right)^{2}+\frac{1}{4\pi^{2}}\left(\frac{\mu_{f}}{T}\right)^{4}\right]\;,$

and the $\order{g^{2}}$ correction

$\displaystyle\hat{p}_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{6}\left(1+\frac{5n_{f}}{12}\right)$ (50)
$\displaystyle\hskip
2.84544pt+\frac{1}{2\pi^{2}}\sum_{f=u,d,s}\left[\frac{1}{2}\left(\frac{\mu_{f}}{T}\right)^{2}+\frac{1}{4\pi^{2}}\left(\frac{\mu_{f}}{T}\right)^{4}\right]\;,$

Expressing the flavor chemical potentials in terms of conserved charge
chemical potentials,

$\displaystyle\hat{\mu}_{u}$ $\displaystyle=$
$\displaystyle(\hat{\mu}_{B}+2\hat{\mu}_{Q})/3\;,$
$\displaystyle\hat{\mu}_{d}$ $\displaystyle=$
$\displaystyle(\hat{\mu}_{B}-\hat{\mu}_{Q})/3\;,$ $\displaystyle\hat{\mu}_{s}$
$\displaystyle=$
$\displaystyle(\hat{\mu}_{B}-\hat{\mu}_{Q})/3-\hat{\mu}_{S}\;,$ (51)

the strangeness number density can be written as

$\displaystyle\hat{n}_{S}$ $\displaystyle=$
$\displaystyle-\frac{1-q_{1}-3s_{1}}{3}\left(1+\frac{1}{2\pi^{2}}g^{2}\right)$
(52)
$\displaystyle\times\left(\hat{\mu}_{B}+-\frac{(1-q_{1}-3s_{1})^{2}}{18\pi^{2}}\hat{\mu}_{B}^{3}\right)\;$

Here we introduced the ratio of chemical potentials,
$s_{1}=\hat{\mu}_{S}/\hat{\mu}_{B}$ and $q_{1}=\hat{\mu}_{Q}/\hat{\mu}_{B}$.
From Eq. 52 we find the constraint for a strangeness-neutral ideal gas,

$s_{1}=\frac{1-q_{1}}{3}\;.$ (53)

Inserting this constraint in Eq. 51 one finds that up to $\order{g^{2}}$ in
perturbation theory the strange quark chemical potential vanishes in
strangeness-neutral matter, $\mu_{s}=0$. Moreover, the electric charge
chemical potential vanishes ($\mu_{Q}=0$) in the isospin-symmetric case
$n_{u}=n_{d}$. Up to $\order{g^{2}}$ we thus may write for the
$\hat{\mu}_{B}$-dependent contribution to the pressure of an isospin-symmetric
medium and the two

Figure 17: Baryon chemical potential (left) and the $\mu_{B}$-dependent part
of the pressure (right) of a 3-flavor, ideal gas of massless quarks versus net
baryon-number density ($n_{B}$). Shown are results for the cases
$(\mu_{Q}=0,n_{S}=0)$ ($A_{f}=2)$ and $(\mu_{Q}=0,\mu_{S}=0)$ ($A_{f}=3)$.
Results are rescaled with powers of $n_{B}$ appropriate in the low density
(main panel) and high density (insets) limits, respectively. Dashed lines show
the respective low and high density limits.

cases888 Because of isospin symmetry and strangeness-neutrality conditions
that we impose, one has $\mu_{s}=0$ at high temperature. For this reason,
strange quarks do not contribute to the density dependent part, in particular
the net baryon-number density of the ideal gas and LO perturbative limits.
This is reflected in the parameter $A_{f}$. of interest to us (i) $\mu_{S}=0$
or (ii) $n_{S}=0$,

$\displaystyle\Delta p/T^{4}$ $\displaystyle=$
$\displaystyle\frac{A_{f}}{18}\left(1-\frac{g^{2}}{2\pi^{2}}\right)\left(\hat{\mu}_{B}^{2}+\frac{1}{18\pi^{2}}\hat{\mu}_{B}^{4}\right)\;,$
(54) $\displaystyle n_{B}/T^{3}$ $\displaystyle=$
$\displaystyle\frac{A_{f}}{9}\left(1-\frac{g^{2}}{2\pi^{2}}\right)\left(\hat{\mu}_{B}+\frac{1}{9\pi^{2}}\hat{\mu}_{B}^{3}\right)\;.$
(55)

where $A_{f}=3$ for $\mu_{S}=0$ and $A_{f}=2$ for $n_{S}=0$.

For the coupling $g^{2}$ we use the two-loop running coupling, allowing for a
free parameter in the renormalization scale, $k_{T}\pi T$, as it also is used
in more refined, resummed approaches, e.g. hard thermal loop (HTL)
calculations [26],

$\displaystyle\beta_{0}$ $\displaystyle=$
$\displaystyle\frac{11-2n_{f}/3}{4\pi}$ $\displaystyle\beta_{1}$
$\displaystyle=$ $\displaystyle\frac{102-38n_{f}/3}{16\pi^{2}}$ $\displaystyle
g^{2}_{1loop}$ $\displaystyle=$
$\displaystyle\frac{2\pi}{\beta_{0}\ln\left(k_{T}\pi
T/\Lambda_{\overline{MS}}\right)}\;,$ (56) $\displaystyle g^{2}$
$\displaystyle=$ $\displaystyle
g^{2}_{1loop}\left[1-\frac{\beta_{1}}{\beta_{0}^{2}}\frac{\ln(2\ln(k_{T}\pi
T/\Lambda_{\overline{MS}}))}{2\ln(k_{T}\pi
T/\Lambda_{\overline{MS}})}\right]\;,$ (57)

with $n_{f}=3$. As a scale factor we use the average of recent values
summarized in the last FLAG report $\Lambda_{\overline{MS}}=339(12)$ MeV [42].

From Eqs. 54 and 55 we easily find the low and high density limits of the EoS
of a massless, strangeness-neutral, isospin-symmetric quark gas to
$\order{g^{2}}$ in high-$T$ perturbation theory,

$\frac{\Delta
p}{T^{4}}=\begin{cases}\frac{9}{2A_{f}(1-g^{2}/2\pi^{2})}\left(\frac{n_{B}}{T^{3}}\right)^{2}&\;,\;\frac{n_{B}}{T^{3}}\ll
1\;,\\\
\frac{1}{4}\left(\frac{81\pi^{2}}{A_{f}(1-g^{2}/2\pi^{2})}\right)^{1/3}\left(\frac{n_{B}}{T^{3}}\right)^{4/3}&\;,\;\frac{n_{B}}{T^{3}}\gg
1\;.\end{cases}$ (58)

For arbitrary net baryon-number densities we invert Eq. 55 to obtain the
baryon chemical potential as function of $n_{B}/T^{3}$ in $\order{g^{2}}$
perturbation theory,

$\displaystyle\hat{\mu}_{B}(\hat{n}_{B})$ $\displaystyle=$ $\displaystyle
y^{1/3}-\frac{3\pi^{2}}{y^{1/3}}\;,$ (59)

with

$\displaystyle y$ $\displaystyle=$
$\displaystyle\frac{3^{4}\pi^{2}}{2A_{f}(1-g^{2}/2\pi^{2})}$
$\displaystyle\hskip
11.38092pt\times\left(\hat{n}_{B}+\sqrt{\hat{n}_{B}^{2}+(2\pi
A_{f}(1-g^{2}/2\pi^{2}))^{2}/3^{5}}\right)\;\;.\;\;$

Inserting Eq. 59 in Eq. 54 we obtain the equation of state of strangeness-
neutral, isospin-symmetric matter in the high-temperature limit. The result in
the ideal gas limit ($g^{2}=0$) is shown in Fig. 17. Comparing with Fig. 4 one
sees that the density range covered at high temperature in current Taylor
series expansions with $\hat{\mu}_{B}\leq 2.5$, i.e. $n_{B}/T^{3}<0.6$
corresponds to the low density limit of an ideal gas, $\mu_{B}/T\sim
n_{B}/T^{3}$. The transition between the low and high density regions,
however, starts to set in for $n_{B}/T^{3}\gtrsim 0.6$.

##  Appendix C Speed of sound

The speed of sound is defined as

$c_{\vec{X}}^{2}=\left(\frac{\partial
p}{\partial\epsilon}\right)_{\vec{X}}\;,$ (61)

where $\vec{X}$ defines the external conditions, e.g. constant entropy per
particle number. We generalize this here to the case of QCD, where three
conserved currents need to be specified. We do so by choosing
$\vec{X}=(s/n_{B},n_{Q}/n_{B},n_{S})$ With this we obtain

$\displaystyle c^{2}_{\vec{X}}$ $\displaystyle=$
$\displaystyle\frac{N_{\vec{X}}}{(\epsilon+p-\mu_{S}n_{S})\ D_{\vec{X}}}\;,\;$
(62)

with

$\displaystyle N_{\vec{X}}$ $\displaystyle=$ $\displaystyle
s^{2}a_{s2}+s^{\prime}a_{s}+s\ (s_{B}a_{B}+s_{Q}a_{Q}+s_{S}a_{S})$ (63)
$\displaystyle+s_{B}^{2}a_{B2}+s_{Q}^{2}a_{Q2}+s_{S}^{2}a_{S2}+s_{B}s_{Q}a_{BQ}$
$\displaystyle+s_{B}s_{S}a_{BS}+s_{Q}s_{S}a_{QS}$ $\displaystyle D_{\vec{X}}$
$\displaystyle=$ $\displaystyle
s^{\prime}b_{s}+s_{B}^{2}b_{B2}+s_{Q}^{2}b_{Q2}+s_{S}^{2}b_{S2}$
$\displaystyle+s_{B}s_{Q}b_{BQ}+s_{B}s_{S}b_{BS}+s_{Q}s_{S}b_{QS}$

Here we use the abbreviation $s_{i}$ for the derivatives of the entropy with
respect to the chemical potential $\mu_{i}$, i.e. $s_{i}=\partial
s/\partial\mu_{i}$.

Using $n_{Q}/n_{B}=r$, the various coefficients appearing in Eq. 63 are given
by

$\displaystyle a_{s2}$ $\displaystyle=$
$\displaystyle-(\chi_{11}^{BS})^{2}\chi_{2}^{Q}+2\chi_{11}^{BQ}\chi_{11}^{BS}\chi_{11}^{QS}-\chi_{2}^{B}(\chi_{11}^{QS})^{2}$
$\displaystyle-(\chi_{11}^{BQ})^{2}\chi_{2}^{S}+\chi_{2}^{B}\chi_{2}^{Q}\chi_{2}^{S}$
$\displaystyle a_{s}$ $\displaystyle=$
$\displaystyle-(r^{2}b_{Q2}+rb_{BQ}+b_{B2})n_{B}^{2}$
$\displaystyle-\frac{1}{2}(rb_{QS}+b_{BS})n_{B}n_{S}$ $\displaystyle a_{B}$
$\displaystyle=$ $\displaystyle(2b_{B2}+rb_{BQ})n_{B}+\frac{1}{2}b_{BS}n_{S}$
$\displaystyle a_{Q}$ $\displaystyle=$
$\displaystyle(b_{BQ}+2rb_{Q2})n_{B}+\frac{1}{2}b_{QS}n_{S}$ $\displaystyle
a_{S}$ $\displaystyle=$ $\displaystyle(b_{BS}+rb_{QS})n_{B}+b_{S2}n_{S}$
$\displaystyle a_{B2}$ $\displaystyle=$
$\displaystyle-r^{2}\chi_{2}^{S}n_{B}^{2}+r\chi_{11}^{QS}n_{B}n_{S}$
$\displaystyle a_{Q2}$ $\displaystyle=$
$\displaystyle-\chi_{2}^{S}n_{B}^{2}+\chi_{11}^{BS}n_{B}n_{S}$ (64)
$\displaystyle a_{S2}$ $\displaystyle=$
$\displaystyle(2r\chi_{11}^{BQ}-r^{2}\chi_{2}^{B}-\chi_{2}^{Q})n_{B}^{2}$
$\displaystyle a_{BQ}$ $\displaystyle=$ $\displaystyle
2r\chi_{2}^{S}n_{B}^{2}-(r\chi_{11}^{BS}+\chi_{11}^{QS})n_{B}n_{S}$
$\displaystyle a_{BS}$ $\displaystyle=$
$\displaystyle-2r(\chi_{11}^{QS}-r\chi_{11}^{BS})n_{B}^{2}+(\chi_{2}^{Q}-r\chi_{11}^{BQ})n_{B}n_{S}$
$\displaystyle a_{QS}$ $\displaystyle=$ $\displaystyle
2(\chi_{11}^{QS}-r\chi_{11}^{BS})n_{B}^{2}+(r\chi_{2}^{B}-\chi_{11}^{BQ})n_{B}n_{S}$

$\displaystyle b_{s}$ $\displaystyle=$ $\displaystyle a_{s2}$ $\displaystyle
b_{B2}$ $\displaystyle=$
$\displaystyle(\chi_{11}^{QS})^{2}-\chi_{2}^{Q}\chi_{2}^{S}$ $\displaystyle
b_{Q2}$ $\displaystyle=$
$\displaystyle(\chi_{11}^{BS})^{2}-\chi_{2}^{B}\chi_{2}^{S}$ $\displaystyle
b_{S2}$ $\displaystyle=$
$\displaystyle(\chi_{11}^{BQ})^{2}-\chi_{2}^{B}\chi_{2}^{Q}$ $\displaystyle
b_{BQ}$ $\displaystyle=$ $\displaystyle
2(\chi_{11}^{BQ}\chi_{2}^{S}-\chi_{11}^{BS}\chi_{11}^{QS})$ $\displaystyle
b_{BS}$ $\displaystyle=$ $\displaystyle
2(\chi_{11}^{BS}\chi_{2}^{Q}-\chi_{11}^{BQ}\chi_{11}^{QS})$ $\displaystyle
b_{QS}$ $\displaystyle=$ $\displaystyle
2(\chi_{11}^{QS}\chi_{2}^{B}-\chi_{11}^{BQ}\chi_{11}^{BS})$ (65)

We note that all cumulants $\chi_{11}^{ab}$ and $\chi_{2}^{c}$, with $a,\
b=B,\ Q,\ S$, appearing in Eqs. 64 and 65 are defined as $\mu_{i}$, $i=a,b$,
derivatives of $P$ and as such are functions of the chemical potentials, i.e
$\chi_{11}^{ab}\equiv\chi_{11}^{ab}(\mu_{B},\mu_{Q},\mu_{S})$. They need to be
chosen such that the constraint $\vec{X}$ is satisfied.

For $\vec{\mu}=0$ Eq. 63 reduces to $N_{\vec{\mu}=0}=s^{2}a_{s2}$ and
$D_{\vec{\mu}=0}=s^{\prime}b_{s}$. All other coefficients appearing in Es. 65
and 64 are proportional to $\mu^{2}$.

As discussed in Appendix A some care has to be taken when rewriting the
derivatives of e.g. entropy density with respect to $T$,

$\displaystyle s^{\prime}\equiv\frac{\partial s}{\partial
T}\Bigg{|}_{\vec{\mu}}$ $\displaystyle=$ $\displaystyle
T^{2}\left(3\hat{s}+T\frac{\partial\hat{s}}{\partial
T}\Bigg{|}_{\vec{\hat{\mu}}}-\sum_{i}\hat{\mu}_{i}\frac{\partial\hat{s}}{\partial\hat{\mu}_{i}}\Bigg{|}_{T}\right).$
(66)

For $\vec{\mu}=0$ we obtain the well known expression

$\left.c^{2}_{s}\right|_{\vec{\mu}=0}=\left.\frac{s}{Ts^{\prime}}\right|_{\vec{\mu}=0}\;.$
(67)

Furthermore, for an ideal gas Eq. 62 gives, of course, $c_{s}^{2}=1/3$ for
arbitrary $(\mu_{B},\mu_{Q},\mu_{S})$.

We are interested here in the speed of sound in strangeness-neutral matter. In
that case all terms proportional to $n_{S}$ appearing in Eqs. 62 and 64 can be
set to zero.

## References

  * Braun-Munzinger  _et al._ [2016] P. Braun-Munzinger, V. Koch, T. Schäfer, and J. Stachel, Properties of hot and dense matter from relativistic heavy ion collisions, Phys. Rept. 621, 76 (2016), arXiv:1510.00442 [nucl-th] .
  * Middeldorf-Wygas  _et al._ [2022] M. M. Middeldorf-Wygas, I. M. Oldengott, D. Bödeker, and D. J. Schwarz, Cosmic QCD transition for large lepton flavor asymmetries, Phys. Rev. D 105, 123533 (2022), arXiv:2009.00036 [hep-ph] .
  * Karsch  _et al._ [2000] F. Karsch, E. Laermann, and A. Peikert, The Pressure in two flavor, (2+1)-flavor and three flavor QCD, Phys. Lett. B 478, 447 (2000), arXiv:hep-lat/0002003 .
  * Borsanyi  _et al._ [2014] S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, and K. K. Szabo, Full result for the QCD equation of state with 2+1 flavors, Phys. Lett. B 730, 99 (2014), arXiv:1309.5258 [hep-lat] .
  * Bazavov  _et al._ [2014] A. Bazavov _et al._ (HotQCD), Equation of state in (2+1)-flavor QCD, Phys. Rev. D 90, 094503 (2014), arXiv:1407.6387 [hep-lat] .
  * Gavai and Gupta [2001] R. V. Gavai and S. Gupta, Quark number susceptibilities, strangeness and dynamical confinement, Phys. Rev. D 64, 074506 (2001), arXiv:hep-lat/0103013 .
  * Allton  _et al._ [2002] C. R. Allton, S. Ejiri, S. J. Hands, O. Kaczmarek, F. Karsch, E. Laermann, C. Schmidt, and L. Scorzato, The QCD thermal phase transition in the presence of a small chemical potential, Phys. Rev. D 66, 074507 (2002), arXiv:hep-lat/0204010 .
  * D’Elia and Lombardo [2003] M. D’Elia and M.-P. Lombardo, Finite density QCD via imaginary chemical potential, Phys. Rev. D 67, 014505 (2003), arXiv:hep-lat/0209146 .
  * de Forcrand and Philipsen [2002] P. de Forcrand and O. Philipsen, The QCD phase diagram for small densities from imaginary chemical potential, Nucl. Phys. B 642, 290 (2002), arXiv:hep-lat/0205016 .
  * Borsanyi  _et al._ [2012] S. Borsanyi, G. Endrodi, Z. Fodor, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo, QCD equation of state at nonzero chemical potential: continuum results with physical quark masses at order $\mu^{2}$, JHEP 08, 053, arXiv:1204.6710 [hep-lat] .
  * Guenther  _et al._ [2017] J. N. Guenther, R. Bellwied, S. Borsanyi, Z. Fodor, S. D. Katz, A. Pasztor, C. Ratti, and K. K. Szabó, The QCD equation of state at finite density from analytical continuation, Nucl. Phys. A 967, 720 (2017), arXiv:1607.02493 [hep-lat] .
  * Ejiri  _et al._ [2006] S. Ejiri, F. Karsch, E. Laermann, and C. Schmidt, The Isentropic equation of state of 2-flavor QCD, Phys. Rev. D 73, 054506 (2006), arXiv:hep-lat/0512040 .
  * Bazavov  _et al._ [2017] A. Bazavov _et al._ , The QCD Equation of State to $\mathcal{O}(\mu_{B}^{6})$ from Lattice QCD, Phys. Rev. D 95, 054504 (2017), arXiv:1701.04325 [hep-lat] .
  * Bollweg  _et al._ [2022a] D. Bollweg, J. Goswami, O. Kaczmarek, F. Karsch, S. Mukherjee, P. Petreczky, C. Schmidt, and P. Scior (HotQCD), Taylor expansions and Padé approximants for cumulants of conserved charge fluctuations at nonvanishing chemical potentials, Phys. Rev. D 105, 074511 (2022a), arXiv:2202.09184 [hep-lat] .
  * Gavai  _et al._ [2005] R. V. Gavai, S. Gupta, and S. Mukherjee, The Speed of sound and specific heat in the QCD plasma: Hydrodynamics, fluctuations and conformal symmetry, Phys. Rev. D 71, 074013 (2005), arXiv:hep-lat/0412036 .
  * Floerchinger and Martinez [2015] S. Floerchinger and M. Martinez, Fluid dynamic propagation of initial baryon number perturbations on a Bjorken flow background, Phys. Rev. C 92, 064906 (2015), arXiv:1507.05569 [nucl-th] .
  * Bazavov  _et al._ [2020] A. Bazavov _et al._ , Skewness, kurtosis, and the fifth and sixth order cumulants of net baryon-number distributions from lattice QCD confront high-statistics STAR data, Phys. Rev. D 101, 074502 (2020), arXiv:2001.08530 [hep-lat] .
  * Dimopoulos  _et al._ [2022] P. Dimopoulos, L. Dini, F. Di Renzo, J. Goswami, G. Nicotra, C. Schmidt, S. Singh, K. Zambello, and F. Ziesché, Contribution to understanding the phase structure of strong interaction matter: Lee-Yang edge singularities from lattice QCD, Phys. Rev. D 105, 034513 (2022), arXiv:2110.15933 [hep-lat] .
  * Bazavov  _et al._ [2019] A. Bazavov _et al._ (HotQCD), Chiral crossover in QCD at zero and non-zero chemical potentials, Phys. Lett. B 795, 15 (2019), arXiv:1812.08235 [hep-lat] .
  * Mazur [2021] L. Mazur, Topological Aspects in Lattice QCD, Ph.D. thesis, Bielefeld University 10.4119/unibi/2956493 (2021).
  * Bollweg  _et al._ [2022b] D. Bollweg, L. Altenkort, D. A. Clarke, O. Kaczmarek, L. Mazur, C. Schmidt, P. Scior, and H.-T. Shu, HotQCD on multi-GPU Systems, PoS LATTICE2021, 196 (2022b), arXiv:2111.10354 [hep-lat] .
  * Bollweg  _et al._ [2021] D. Bollweg, J. Goswami, O. Kaczmarek, F. Karsch, S. Mukherjee, P. Petreczky, C. Schmidt, and P. Scior (HotQCD), Second order cumulants of conserved charge fluctuations revisited: Vanishing chemical potentials, Phys. Rev. D 104, 10.1103/PhysRevD.104.074512 (2021), arXiv:2107.10011 [hep-lat] .
  * [23] AnalysisToolbox: A set of Python tools for analyzing physics data, in particular targeting lattice QCD, https://github.com/LatticeQCD/AnalysisToolbox.
  * Virtanen  _et al._ [2020] P. Virtanen _et al._ , SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods 17, 261 (2020).
  * Haque  _et al._ [2014] N. Haque, A. Bandyopadhyay, J. O. Andersen, M. G. Mustafa, M. Strickland, and N. Su, Three-loop HTLpt thermodynamics at finite temperature and chemical potential, JHEP 05, 027, arXiv:1402.6907 [hep-ph] .
  * Mogliacci  _et al._ [2013] S. Mogliacci, J. O. Andersen, M. Strickland, N. Su, and A. Vuorinen, Equation of State of hot and dense QCD: Resummed perturbation theory confronts lattice data, JHEP 12, 055, arXiv:1307.8098 [hep-ph] .
  * Bazavov  _et al._ [2013] A. Bazavov, H. T. Ding, P. Hegde, F. Karsch, C. Miao, S. Mukherjee, P. Petreczky, C. Schmidt, and A. Velytsky, Quark number susceptibilities at high temperatures, Phys. Rev. D 88, 094021 (2013), arXiv:1309.2317 [hep-lat] .
  * Ding  _et al._ [2015] H. T. Ding, S. Mukherjee, H. Ohno, P. Petreczky, and H. P. Schadler, Diagonal and off-diagonal quark number susceptibilities at high temperatures, Phys. Rev. D 92, 074043 (2015), arXiv:1507.06637 [hep-lat] .
  * Bellwied  _et al._ [2015a] R. Bellwied, S. Borsanyi, Z. Fodor, S. D. Katz, A. Pasztor, C. Ratti, and K. K. Szabo, Fluctuations and correlations in high temperature QCD, Phys. Rev. D 92, 114505 (2015a), arXiv:1507.04627 [hep-lat] .
  * Prosperi  _et al._ [2007] G. M. Prosperi, M. Raciti, and C. Simolo, On the running coupling constant in QCD, Prog. Part. Nucl. Phys. 58, 387 (2007), arXiv:hep-ph/0607209 .
  * Clarke  _et al._ [2021] D. A. Clarke, O. Kaczmarek, F. Karsch, A. Lahiri, and M. Sarkar, Sensitivity of the Polyakov loop and related observables to chiral symmetry restoration, Phys. Rev. D 103, L011501 (2021), arXiv:2008.11678 [hep-lat] .
  * Pelissetto and Vicari [2013] A. Pelissetto and E. Vicari, Relevance of the axial anomaly at the finite-temperature chiral transition in QCD, Phys. Rev. D 88, 105018 (2013), arXiv:1309.5446 [hep-lat] .
  * Kaczmarek  _et al._ [2011] O. Kaczmarek, F. Karsch, E. Laermann, C. Miao, S. Mukherjee, P. Petreczky, C. Schmidt, W. Soeldner, and W. Unger, Phase boundary for the chiral transition in (2+1) -flavor QCD at small values of the chemical potential, Phys. Rev. D 83, 014504 (2011), arXiv:1011.3130 [hep-lat] .
  * Karsch [2019] F. Karsch, Critical behavior and net-charge fluctuations from lattice QCD, PoS CORFU2018, 163 (2019), arXiv:1905.03936 [hep-lat] .
  * Bellwied  _et al._ [2015b] R. Bellwied, S. Borsanyi, Z. Fodor, J. Günther, S. D. Katz, C. Ratti, and K. K. Szabo, The QCD phase diagram from analytic continuation, Phys. Lett. B 751, 559 (2015b), arXiv:1507.07510 [hep-lat] .
  * Bonati  _et al._ [2018] C. Bonati, M. D’Elia, F. Negro, F. Sanfilippo, and K. Zambello, Curvature of the pseudocritical line in QCD: Taylor expansion matches analytic continuation, Phys. Rev. D 98, 054510 (2018), arXiv:1805.02960 [hep-lat] .
  * Adamczyk  _et al._ [2017] L. Adamczyk _et al._ (STAR), Bulk Properties of the Medium Produced in Relativistic Heavy-Ion Collisions from the Beam Energy Scan Program, Phys. Rev. C 96, 044904 (2017), arXiv:1701.07065 [nucl-ex] .
  * D’Elia [2019] M. D’Elia, High-Temperature QCD: theory overview, Nucl. Phys. A 982, 99 (2019), arXiv:1809.10660 [hep-lat] .
  * Borsanyi  _et al._ [2020] S. Borsanyi, Z. Fodor, J. N. Guenther, R. Kara, S. D. Katz, P. Parotto, A. Pasztor, C. Ratti, and K. K. Szabo, QCD Crossover at Finite Chemical Potential from Lattice Simulations, Phys. Rev. Lett. 125, 052001 (2020), arXiv:2002.02821 [hep-lat] .
  * Shuryak [1978] E. V. Shuryak, Theory of Hadronic Plasma, Sov. Phys. JETP 47, 212 (1978).
  * Chin [1978] S. A. Chin, Transition to Hot Quark Matter in Relativistic Heavy Ion Collision, Phys. Lett. B 78, 552 (1978).
  * Aoki  _et al._ [2022] Y. Aoki _et al._ (Flavour Lattice Averaging Group (FLAG)), FLAG Review 2021, Eur. Phys. J. C 82, 869 (2022), arXiv:2111.09849 [hep-lat] .