aainstitutetext: Fakultät für Physik, Universität Bielefeld, D-33615
Bielefeld, Germanybbinstitutetext: Institut für Theoretische Physik, Goethe
Universität Frankfurt, D-60438 Frankfurt am Main, Germany

# Equation of state and speed of sound of isospin-asymmetric QCD on the
lattice

B. B. Brandt b F. Cuteri a G. Endrődi<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We determine the QCD equation of state at nonzero temperature in the presence
of an isospin asymmetry between the light quark chemical potentials on the
lattice. Our simulations employ $N_{f}=2+1$ flavors of dynamical staggered
quarks at physical masses, using three different lattice spacings. The main
results, obtained at the individual lattice spacings, are based on a two-
dimensional spline interpolation of the isospin density, from which all
relevant quantities can be obtained analytically. In particular, we present
results for the pressure, the interaction measure, the energy and entropy
densities, as well as the speed of sound. Remarkably, the latter is found to
exceed its ideal gas limit deep in the pion condensed phase, the first account
of the violation of this limit in first principles QCD. Finally, we also
compute the phase diagram in the temperature – isospin density plane for the
first time. Even though the results are not continuum extrapolated and thus
not final, the data for all observables will be useful for the benchmarking of
effective theories and low-energy models of QCD and are provided in ancillary
files for simple reuse.

##  1 Introduction

The theory of the strong interactions is Quantum Chromodynamics (QCD),
featuring confinement of quarks and gluons at low energies, as well as
asymptotic freedom at high scales. Albeit radically different in their
properties, these two phases of strongly interacting matter are connected by a
smooth crossover transition at zero net quark density according to lattice QCD
simulations Aoki:2006we ; Bhattacharya:2014ara . How the dominant degrees of
freedom transform from composite objects (hadrons) to colored quarks and
gluons through this transition is described by the equation of state (EoS) of
the system. In particular, the EoS gives a complete description of equilibrium
QCD in terms of a relationship between thermodynamic observables including the
pressure, the energy density or the entropy density. The phenomenological
relevance of the EoS is manifold and ranges from heavy-ion physics to
astrophysics and cosmology. The above observables control the evolution of the
quark-gluon plasma in hydrodynamic models of heavy-ion collisions
Teaney:2001av ; Kolb:2003dz , the expansion of the early universe via the
Friedmann equations Boyanovsky:2006bf and also the mass and radius of stable
neutron stars through the Tolman-Oppenheimer-Volkoff equation Lattimer:2000nx
. In the latter context, a particularly relevant feature of the EoS is the
speed of sound $c_{s}$ of QCD matter and the related polytropic index
$\gamma$, which might serve as a proxy to distinguish stars with and without
deconfined quark matter cores Tews:2018kmu ; Annala:2019puf .

The above physical systems contain QCD matter in very different environments
and therefore require the knowledge of the EoS as a function of different
control parameters. These parameters include the temperature $T$, the chemical
potentials $\mu$ conjugate to the conserved charges, as well as further
variables like external electromagnetic fields. The relevant chemical
potentials are charge $\mu_{Q}$, baryon $\mu_{B}$ and strangeness $\mu_{S}$
chemical potentials. While the baryon chemical potential is in most cases
assumed to carry the dominant effect, in some cases $\mu_{Q}$ can play the
major role. This occurs for example for an early Universe featuring large
lepton flavour asymmetries Oldengott:2017tzj ; Wygas:2018otj ; Middeldorf-
Wygas:2020glx . Here the isentropic cosmological expansion leads to
substantial charge chemical potentials, triggering the onset of pion
condensation and producing characteristic signals for primordial black holes
and gravitational wave spectra Vovchenko:2020crk . Significant negative charge
chemical potentials also arise for QCD matter in neutron stars due to an
excess of down quarks over up quarks.

The charge chemical potential can be rewritten in terms of a nonzero isospin
chemical potential $\mu_{I}$. While lattice simulations with a generic
combination of quark chemical potentials suffer from the infamous complex
action (or sign) problem, QCD at pure isospin chemical potential, i.e., at
vanishing other chemical potential components, has a real action and is
amenable to direct Monte-Carlo simulations Son:2000xc . At low temperature,
this setting exhibits a second-order phase transition to a phase with a Bose-
Einstein condensate (BEC) of charged pions according to chiral perturbation
theory Son:2000xc . First lattice simulations with higher-than-physical quark
masses have qualitatively confirmed this expectation Kogut:2002tm ;
Kogut:2002zg ; Kogut:2004zg ; Endrodi:2014lja and also gave important insight
to the structure of the phase diagram in the $T$-$\mu_{I}$ plane
deForcrand:2007uz ; Cea:2012ev as well as pion dynamics at low $T$
Detmold:2012wc . In Refs. Brandt:2017oyy ; Brandt:2018omg we carried out a
systematic investigation of this system with physical quark masses and
determined the phase diagram in the continuum limit, revealing an interesting
interplay of chiral symmetry breaking, deconfinement and Bose-Einstein
condensation. The continuum phase diagram for the parameter space relevant for
this study is shown in Fig. 1.

Figure 1:  Phase diagram of isospin asymmetric QCD determined in Ref.
Brandt:2017oyy . The green band is the phase boundary to the phase with Bose-
Einstein condensation (BEC), the blue band is the prolongation of the chiral
crossover in the $T$-$\mu_{I}$ plane and the red data point marks the pseudo-
triple point, the meeting point of the crossover and the BEC phase. In yellow
we also show a sketch of the region, where the speed of sound exceeds its
conformal limit.

In this paper we determine the EoS throughout the phase diagram for a broad
range of temperatures and isospin chemical potentials. Generalizing our
approach at (approximately) zero temperature Brandt:2018bwq , we construct the
pressure, the energy and entropy densities, the interaction measure and the
speed of sound from the isospin density as primary observable. At low
temperatures and high $\mu_{I}$, we find that the speed of sound increases
above its conformal limit $1/\sqrt{3}$ (we use natural units, with the speed
of light set to unity). This is the first evidence for the explicit violation
in first principles QCD of this general bound expected from holography
Cherman:2009tw and, together with the polytropic index, which we compute as
well, might provide relevant information for the modeling of the EoS based on
neutron star radii and masses Tews:2018kmu ; Annala:2019puf . We also included
a sketch of the region where $c_{s}$ exceeds this conformal bound in Fig. 1.
Besides the EoS, we also use our results to draw the QCD phase diagram in the
temperature – isospin density plane. This result, together with the complete
tabulated EoS as shown in the plots, is available in the ancillary files
submitted to the arXiv along with the preprint of this paper. To facilitate
the use of the EoS in phenomenological models, we provide an accompanying data
publication datapub , including the physical observables and the uncertainties
for all temperatures and chemical potentials where results are available. The
results can be used as benchmarks for low-energy models and effective theories
of QCD as well as for comparing to functional approaches. First accounts of
our findings have been given in Refs. Vovchenko:2020crk ; Brandt:2017zck ;
Brandt:2018wkp ; Brandt:2021yhc .

##  2 Determination of the EoS

###  2.1 Simulation setup and main observables

In our simulations we use $N_{f}=2+1$ flavors of rooted staggered quarks with
two steps of stout smearing at physical quark masses and the tree-level
Symanzik improved gluon action. The line of constant physics for the bare
quark masses is taken from Ref. Borsanyi:2010cj . For the approach to the
continuum limit we use lattices with temporal extents $N_{t}=8,\,10$ and $12$
and aspect ratios of $N_{s}/N_{t}\approx 3$ (these ensembles have already been
used for the phase diagram Brandt:2017oyy ), together with a number of
additional ensembles at $T=0$. More details concerning the run parameters are
collected in appendix D. The isospin chemical potential $\mu_{I}$ enters the
light quark Dirac operator in an exponential form and is normalized such that
pion condensation sets in at zero temperature at $\mu_{I}=m_{\pi}/2$. As in
our previous studies Brandt:2017oyy ; Brandt:2018omg ; Brandt:2018bwq ;
Vovchenko:2020crk the simulations are performed including a pionic source
parameter $\lambda$, in the light quark mass matrix (see also Refs.
Kogut:2002tm ; Kogut:2002zg ; Endrodi:2014lja ), which serves as an infrared
regulator and triggers pion condensation in a finite volume. Physical results
are obtained by means of an extrapolation $\lambda\to 0$, which is facilitated
by improving the observables and reweighting the configurations. For details
on this improvement, see Refs. Brandt:2017oyy ; Brandt:2018omg ;
Brandt:2018bwq . For computing uncertainties we use the bootstrap procedure
with 1000 samples.

As we will see below, the main observable is the isospin density,

$n_{I}=\frac{T}{V}\frac{\partial\log\mathcal{Z}}{\partial\mu_{I}}\,,$ (1)

from which the full $\mu_{I}$-dependence of the EoS can be extracted. The
strategy for the EoS computation will be outlined below. $n_{I}$ can be
computed directly from the simulations as described in Ref. Brandt:2018bwq .
To perform the $\lambda$-extrapolations we use the improvement program
introduced in Ref. Brandt:2017oyy with the application to $n_{I}$ as explained
in Ref. Brandt:2018bwq ; Brandt:2018omg . This improvement program results in
fully controlled extrapolations and from now on we only discuss results which
have already been extrapolated to $\lambda=0$. We note that a well controlled
$\lambda$-extrapolation is of particular importance to facilitate the
following spline interpolations of $n_{I}$ in $T$ and in $\mu_{I}$.

###  2.2 The EoS from an interpolation of the isospin density

Apart from the isospin density $n_{I}$, the main task for the determination of
the EoS is the computation of the pressure $p$ and the interaction measure
$I$. All other relevant quantities, apart from the speed of sound (to be
discussed in Sec. 4) follow from these three quantities. In particular, the
energy and entropy densities are given by

$\epsilon=I+3p,\quad\quad s=\frac{\epsilon+p-\mu_{I}n_{I}}{T}\,.$ (2)

The EoS at vanishing chemical potential has been computed in the continuum
limit in various setups and by different collaborations, see e.g. Refs.
Borsanyi:2013bia ; HotQCD:2014kol . It is thus convenient to separate the
effects due to nonzero temperature from the modifications due to the presence
of a nonzero isospin chemical potential. This is possible for the quantities
of this section, for which the two types of contributions are added, but not
for the speed of sound, defined via directional derivatives in the
$T$-$\mu_{I}$ plane. The pressure and the interaction measure can be written
as

$\begin{array}[]{rcl}\displaystyle p(T,\mu_{I})&=&\displaystyle p(T,0)+\Delta
p(T,\mu_{I})\,,\vspace*{2mm}\\\ \displaystyle I(T,\mu_{I})&=&\displaystyle
I(T,0)+\Delta I(T,\mu_{I})\,,\end{array}$ (3)

where the modifications of the EoS due to the isospin chemical potential,
$\Delta p$ and $\Delta I$, are the objects of interest in our study. Whenever
we need to use results for $p(T,0)$ and $I(T,0)$, we use the results obtained
from a reanalysis of the data of Ref. Borsanyi:2013bia with the
parameterisation of $I$ discussed in section 3 of Ref. Borsanyi:2010cj . The
correct inclusion of the correlations of the associated parameters are of
particular relevance for the computation of the speed of sound. The details of
the $\mu_{I}=0$ data are discussed in appendix A.

A possible starting point for the computation of the EoS is the relation

$\frac{I(T,\mu)}{T^{4}}=T\frac{\partial}{\partial
T}\left(\frac{p(T,\mu_{I})}{T^{4}}\right)+\frac{\mu_{I}n_{I}(T,\mu_{I})}{T^{4}}\,.$
(4)

At vanishing chemical potential this equation is used to rewrite $I$ as a
derivative of the partition function with respect to the lattice scale
Engels:1990vr ; Blum:1994zf ; Engels:1996ag ; Aoki:2005vt . At nonzero
chemical potential, one can follow a similar strategy to calculate the
modifications of pressure and interaction measure using $\mu=0$ subtraction
(see, e.g., Ref. Allton:2003vx and Refs. Iida:2022hyy ; Itou:2022ebw ). The
direct application of this method at nonzero $\mu_{I}$ is discussed in
appendix B, where we will see that it leads to large uncertainties for the
interaction measure, due to the $\mu_{I}=0$ subtractions. An alternative,
which leads to more accurate results that we will present in the following, is
to use a two-dimensional smooth interpolation of the results for $n_{I}$ to
obtain the function $n_{I}(T,\mu_{I})$.111We note that in the infinite volume,
$n_{I}(T,\mu_{I})$ has a kink at the second order phase transition to the BEC
phase, which cannot be described by a smooth function. Here we are working in
a finite volume, where this kink is absent; it only reappears in the
thermodynamic limit. Therefore a description of this region using smooth
splines is possible. Since the isospin density is the derivative of the
pressure with respect to $\mu_{I}$, see Eq. (1), the modification of the
pressure can be computed from such an interpolation as Vovchenko:2020crk ;
Brandt:2018wkp ; Brandt:2017zck

$\Delta
p(T,\mu_{I})=\int_{0}^{\mu_{I}}d\mu^{\prime}_{I}\,n_{I}(T,\mu^{\prime}_{I})\,.$
(5)

Inserting this into Eq. (4) gives for the interaction measure (see also Ref.
Vovchenko:2020crk )

$\Delta
I(T,\mu_{I})=\mu_{I}n_{I}(T,\mu_{I})+\int_{0}^{\mu_{I}}d\mu_{I}^{\prime}\,\Big{[}T\frac{\partial}{\partial
T}-4\Big{]}n_{I}(T,\mu_{I}^{\prime})\,.$ (6)

The remaining task is to perform the interpolation and to obtain
$n_{I}(T,\mu_{I})$ as a two-dimensional function. Since the interpolation of
an unknown function based on a discrete set of points subject to statistical
uncertainties is an ill-posed inverse problem, the final interpolation will
not be unique. This is already true for the interpolations used for the
computation of the EoS at $\mu=0$ and will extend to the computation of the
modifications of the EoS due to nonzero $\mu_{I}$ discussed in the next
section. The task at hand is to obtain an interpolation which is close to the
actual physical solution while remaining as model-independent as possible. For
the purpose of model-independence we average over all possible two-dimensional
cubic spline interpolations with variable spline nodepoints (spline fits),
weighted with the goodness of the description of the data. Here the goodness
of the description is determined via the Akaike information criterion
Akaike1973InformationTA and we have included a term to suppress unwanted (and
unphysical) oscillatory solutions. Furthermore, the spline boundary conditions
are chosen carefully to include the mandatory physical information on the
interpolated function. The individual spline configurations are generated by a
spline Monte-Carlo already introduced in Ref. Brandt:2016zdy and discussed
further in appendix C, where we also show a set of representative examples for
this interpolation.

##  3 Results for the EoS

###  3.1 Thermodynamic observables at $N_{t}=8$

Figure 2:  Results for the modifications of the pressure (top left), the
interaction measure (top right), the energy density (bottom left) and the
entropy density (bottom right) due to nonzero isospin chemical potential on
the $8\times 24^{3}$ lattices. Uncertainties are not shown for better
visibility.

Figure 3:  Results for the modifications of the pressure (left) and the
interaction measure (right) versus $\mu_{I}$ for different temperatures
obtained on the $8\times 24^{3}$ lattices from the two-dimensional spline
interpolation for $n_{I}(T,\mu_{I})$ described in the text.

Given the two-dimensional interpolation for the isospin density
$n_{I}(T,\mu_{I})$, we can now proceed with the computation of the
thermodynamic quantities we are interested in, namely the pressure $p$, the
interaction measure $I$, as well as energy $\epsilon$ and entropy $s$
densities using Eqs. (5), (6), (3) and (2). The modifications of the
individual observables due to $\mu_{I}\neq 0$ are shown in Fig. 2. The
pressure shows the strongest changes due to $\mu_{I}$ for small
temperatures.222Note, that this is partly also due to the normalization by
$T^{4}$. This normalization will become singular when we approach the $T=0$
limit. Here we have chosen this normalization to consent to the one typically
used in the literature and to allow for easier comparisons.

To show the magnitude of uncertainties and to have a more quantitative
picture, we plot the pressure versus $\mu_{I}$ for different temperatures in
Fig. 3. Within the BEC phase the modification of the interaction measure shows
an initial increase with $\mu_{I}$ before it decreases towards larger
$\mu_{I}$ values. This phenomenon has already been observed at vanishing
temperatures Brandt:2018bwq ; Brandt:2021yhc , where it leads to a negative
interaction measure starting at around $\mu_{I}/m_{\pi}\approx 0.84$
Brandt:2021yhc , and is a clear signature of the presence of the BEC in the
EoS, see also Ref. Vovchenko:2020crk . We will discuss the full interaction
measure at $T\neq 0$ below. The decrease of $\Delta I$ is shifted to larger
$\mu_{I}$ values with increasing temperature, as can be seen from the right
panel of Fig. 3, where we show the interaction measure versus $\mu_{I}$ for a
few different temperatures. The decrease is no longer visible at temperatures
above around 150 MeV – note that this is just around the edge of the BEC phase
for this particular lattice spacing. For these temperatures $\Delta I$ shows a
peak in temperature direction if $\mu_{I}/m_{\pi}\gtrsim 0.5$, indicating the
strong influence of the phase transition on the interaction measure. The
strong change of $I$ in this region also translates to energy and entropy
densities. The former generically shows a strong increase with $\mu_{I}$,
which becomes less pronounced at larger temperatures outside of the BEC phase.
The modifications of the entropy density mainly follow those of the
interaction measure.

Figure 4:  Results for the pressure (top left), the interaction measure (top
right), the energy density (bottom left) and the entropy density (bottom
right) in units of the temperature obtained on the $8\times 24^{3}$ lattices
from the two-dimensional spline interpolation for $n_{I}(T,\mu_{I})$ described
in the text. Uncertainties are not shown for better visibility. Figure 5:
Results for the interaction measure versus the temperature for different
values of $\mu_{I}/m_{\pi}$ obtained on the $8\times 24^{3}$ lattices from the
two-dimensional spline interpolation for $n_{I}(T,\mu_{I})$ described in the
text. We also compare to the $\mu_{I}=0$ results from the parameterization of
Ref. Borsanyi:2013bia described in appendix A.

To get a feeling for the overall magnitudes, we plot the full results for $p$,
$I$, $\epsilon$ and $s$ in Fig. 4. To obtain the full results, we use the
parameterization from Ref. Borsanyi:2013bia with the coefficients introduced
in appendix A for the $\mu_{I}=0$ quantities from Eq. (3).333Note that this
procedure uses continuum results at $\mu_{I}=0$, but our results for the
modifications by $\mu_{I}$ obtained at non-vanishing lattice spacing. We
observe that the pressure generically rises with $T$ and $\mu_{I}$. The
exception is again the small $T$ and large $\mu_{I}$ region, where the
increase seen in the plot, however, stems from the normalization with the
temperature as explained above. A similar monotonic rise with $T$ and
$\mu_{I}$ is generically observed for the energy density, even though the
increase tends to become less pronounced for larger temperatures outside of
the BEC phase. Note that the entropy density vanishes in the zero temperature
limit for all values of $\mu_{I}$.

For the interaction measure we can clearly observe the decrease with $\mu_{I}$
deep in the BEC phase for small temperatures. Furthermore, we also see a
flattening and the onset of a peak-like structure for temperatures around the
BEC phase boundary. This is an interesting observation. A peak of the
interaction measure is present at $\mu_{I}=0$ Borsanyi:2013bia , see the gray
curve in Fig. 5, at a temperature of around 190 MeV, i.e., above the phase
transition. The observed flattening of $I$ might thus be the onset of this
peak structure, shifted towards smaller temperatures for increasing $\mu_{I}$.
If this is the case, the peak approaches the thermal crossover at smaller
values of $\mu_{I}$, until it becomes mostly consistent with the boundary of
the BEC phase, see Fig. 1, around $\mu_{I}/m_{\pi}\approx 0.8$. This is
further visualized in Fig. 5, where we show the interaction measure versus the
temperature for different $\mu_{I}/m_{\pi}$. Together with the shift of the
onset of the plateau, we also observe the development of a narrower peak
structure for larger $\mu_{I}$ values. This might be an effect of the second
order phase transition at the BEC phase boundary. A similar decrease of the
temperature of the maximum with the chemical potential is also observed at
nonzero baryon chemical potential, see Borsanyi:2012cr ; Guenther:2017hnx (as
well as the large modifications of the interaction measure around $T\approx
155$ MeV observed in Ref. Bollweg:2022fqq ), for instance.

###  3.2 Lattice artifacts at $N_{t}=8$, 10 and 12

So far we have discussed results obtained for one value of $N_{t}$,
corresponding to one particular lattice spacing for each value of $T$. To
obtain continuum results we have to increase $N_{t}$ while keeping the
physical temperature and the aspect ratio $N_{s}/N_{t}$ fixed. Unfortunately,
for the isospin density a well controlled continuum extrapolation is not
possible with the $N_{t}=8,\,10$ and 12 lattices currently at our disposal. To
show this, we plot the results for $n_{I}$ obtained from $8\times 24^{3}$,
$10\times 28^{3}$ and $12\times 36^{3}$ lattices for two different
temperatures in Fig. 6. For both temperatures, we observe that at small
$\mu_{I}$ the results of the $N_{t}=12$ lattices lie between those of the
$N_{t}=8$ and 10 lattices. This indicates that higher order lattice artifacts
are still present in this region and impede a proper continuum extrapolation.
This behavior is in agreement with the large lattice artifacts observed for
the leading order Taylor expansion coefficient in the direction of the isospin
chemical potential, see, e.g. Ref. Borsanyi:2011sw . The ordering of the
results from the different $N_{t}$ values remains even for larger values of
$\mu_{I}$ at $T=149$ MeV. At smaller temperatures the data rearranges in the
vicinity of the phase transition to the BEC phase, such that the magnitude of
$n_{I}$ increases from $N_{t}=12$ to $N_{t}=10$ and 8.

Figure 6:  Results for the isospin density obtained on $8\times 24^{3}$,
$10\times 28^{3}$ and $12\times 36^{3}$ lattices, together with the spline
interpolation ($T=122$ and $149$ MeV slices of the global two-dimensional
spline fit).

Figure 7:  Results for the pressure (left) and the interaction measure (right)
for temperatures $T=120$, 145 and 165 MeV (top to bottom) from $N_{t}=8$, 10
and 12 lattices.

While a direct continuum extrapolation for the isospin density is not possible
with the current dataset, we can still look at the magnitude of lattice
artifacts for the other observables related to the EoS. We show the results
for the pressure and the interaction measure from the lattices at different
values of $N_{t}$ for three different temperatures in Fig. 7. The plot shows
that the pressure, as the direct integral over isospin density, suffers from
similar lattice artifacts as $n_{I}$ and again the results from the $N_{t}=12$
lattice are located between those of the $N_{t}=8$ and 10 lattices. For the
interaction measure the situation is a bit different in the sense that all
results typically overlap within the (comparably large) uncertainties.

###  3.3 The phase diagram in the $T$-$n_{I}$ plane

Figure 8:  Results for the chiral crossover temperature (left) and the BEC
phase boundary (right) in the $T$-$n_{I}$ plane for different values of
$N_{t}$, based on our results from Ref. Brandt:2017oyy .

From the interpolation of $n_{I}$ we can also extract the phase diagram in the
$T$-$n_{I}$ plane, of which a preliminary version has been presented in Ref.
Brandt:2021yhc . To this end we use the phase boundaries in the $T$-$\mu_{I}$
plane from Ref. Brandt:2017oyy for the individual temporal extents, $N_{t}=8$
to 12, and determine the value of $n_{I}$ on these phase boundaries. The
results for the chiral crossover and the BEC phase boundary for all $N_{t}$
are shown in Fig. 8. While the chiral crossover does not show any lattice
spacing dependence, the BEC phase boundary exhibits slight lattice artifacts.

Figure 9:  Results for the phase diagram in the $T$-$n_{I}$ plane for
$N_{t}=12$ together with the results from next-to-leading order chiral
perturbation theory of Ref. Adhikari:2020kdn (red curve). The continuation of
the $N_{t}=12$ pion condensation boundary to the chiral perturbation theory
line merely serves to guide the eye.

Due to the Silver blaze property at $T=0$, the BEC phase transition starts at
$n_{I}=T=0$. The BEC phase boundary at small temperatures and isospin
denstities has been calculated via next-to-leading order chiral perturbation
theory Adhikari:2020kdn . We show those results up to $T=30$ MeV (where the
phase boundary in chiral perturbation theory starts to deviate significantly
from the lattice result) together with the phase boundaries from the
$N_{t}=12$ lattices in Fig. 9.

##  4 The speed of sound

###  4.1 Computation from the interpolation of the isospin density

Another interesting observable related to the EoS is the speed of sound
$c_{s}$. The isentropic speed of sound, i.e., the speed of a sound wave
travelling isentropically through the plasma, is defined as

$c_{s}^{2}=\left.\frac{\partial p}{\partial\epsilon}\right|_{\rm
iso}=\frac{\partial_{\xi}\,p}{\partial_{\xi}\,\epsilon}\,,$ (7)

where the subscript ${\rm iso}$ refers to the derivative taken in the
direction of isentropic trajectories in phase space, i.e., for QCD at generic
nonzero quark chemical potentials in the direction where

$\frac{s}{n_{q}}={\rm
const}\quad\forall\>\text{flavours}\>q\qquad(q=u,\,d,\,s\>\>\text{for}\>\>N_{f}=2+1)\,,$
(8)

and we have introduced the directional derivative $\partial_{\xi}$ in this
direction in parameter space. For a pure isospin chemical potential the only
relevant density is the isospin density, so that the condition (8) reduces to

$\frac{s}{n_{I}}={\rm
const}\,,\quad\text{or}\quad\partial_{\xi}\left(\frac{s}{n_{I}}\right)=0\,,$
(9)

and we can write the directional derivative as

$\partial_{\xi}=\vec{\xi}\cdot\left(\begin{array}[]{c}\partial_{T}\\\
\partial_{\mu_{I}}\end{array}\right)=\cos\alpha\,\partial_{T}+\sin\alpha\,\partial_{\mu_{I}}\,.$
(10)

Combining Eqs. (9) and (10) we can compute the angle $\alpha$ via

$\tan\alpha=-\left(n_{I}\frac{\partial s}{\partial T}-s\frac{\partial
n_{I}}{\partial T}\right)\>\left(n_{I}\frac{\partial
s}{\partial\mu_{I}}-s\frac{\partial n_{I}}{\partial\mu_{I}}\right)^{-1}\,,$
(11)

where all quantities and derivatives can be obtained analytically from the
spline interpolation for $n_{I}$ and the analytic form for the interaction
measure at $\mu_{I}=0$ from appendix A. Once $\alpha$ has been obtained, one
can similarly analytically compute the directional derivatives of $p$ and
$\epsilon$ in Eq. (7).

Another quantity of interest related to the speed of sound, in particular for
astrophysical and cosmological applications, is the polytropic index (see
Annala:2019puf , as well as the lectures Silbar:2003wm ; Sagert:2005fw )

$\gamma=\left.\frac{d\log p}{d\log\epsilon}\right|_{\rm
iso}=\frac{\epsilon}{p}\,c_{s}^{2}\,.$ (12)

In the conformal limit, approached by QCD at asymptotically large densities or
temperatures, it takes a value of $\gamma=1$, while in the hadronic regime
conformal symmetry is broken due to spontaneous chiral symmetry breaking,
leading to large values of $\gamma$ in the range of $\gamma\gtrsim 2$ (see the
discussion in Ref. Annala:2019puf ). Consequently, $\gamma$ can be seen as a
measure in the EoS to distinguish between regions of hadron dominated matter
(confined) or matter dominated by free quarks (quarkyonic/deconfined). In the
study of Ref. Annala:2019puf a value of $\gamma=1.75$ has been introduced to
distinguish between these two types of matter in neutron star cores. Finally,
we also look at the normalized trace anomaly Fujimoto:2022ohj ,

$\Delta=\frac{1}{3}-\frac{p}{\epsilon}=\frac{I}{3\,\epsilon}\,,$ (13)

which should be a number between $-2/3$ and $1/3$ due to causality and
thermodynamic stability. Furthermore, in Ref. Fujimoto:2022ohj it has been
argued that $\Delta\geq 0$.

###  4.2 Speed of sound at vanishing temperature

Before discussing the results for the isentropic speed of sound in the
parameter space of nonzero $(T,\,\mu_{I})$, it is instructive to look at the
limiting case of vanishing temperature. An initial study of the EoS at $T=0$
on a coarse lattice with $a\approx 0.29$ fm has already been presented in
Refs. Brandt:2018bwq ; Brandt:2021yhc . Here we will present new results for
the speed of sound at $T=0$, obtained on $24^{3}\times 32$ and $32^{3}\times
48$ lattices at lattice spacings of $a\approx 0.22$ fm and $a\approx 0.15$ fm,
respectively, including data up to $\mu_{I}/m_{\pi}\approx 1$. The results for
these lattice spacings have already been presented partly in Ref.
Brandt:2022fij where they also have been compared to the $a\approx 0.29$ fm
data.

Figure 10:  Left: Results for the isospin density together with the spline
interpolation at $T=0$ obtained on $24^{3}\times 32$ and $32^{3}\times 48$
lattices and lattice spacings of $a\approx 0.22$ and 0.15 fm. The yellow part
of the curve is obtained directly from the chiral perturbation theory
expression for $n_{I}$. Right: Results for the isentropic speed of sound at
$T=0$, obtained from the spline interpolation of the left panel. Also shown
are the chiral perturbation theory result (dashed yellow line) with the pion
decay constant obtained from the fit discussed in the text as well as the
conformal bound Cherman:2009tw (dashed gray line).

The starting point for the extraction of the EoS at zero temperature is again
the isospin density, from which one can obtain the pressure and, consequently,
all other thermodynamic quantities, using Eq. (5). Due to the Silver Blaze
property, the isospin density vanishes outside of the BEC phase at $T=0$. In
practice, the simulations are performed at a small but non-vanishing
temperature, so that residual temperature effects on $n_{I}$ need to be
corrected in the vicinity of the transition. As already done in Ref.
Brandt:2018bwq we use chiral perturbation theory Son:2000xc to correct for
these $T\neq 0$ effects. In particular, we fit the results for $n_{I}$ for the
two smallest values of $\mu_{I}$ within the BEC phase, i.e., we include the
data points up to $\mu_{I}/m_{\pi}\approx 0.65$, to the chiral perturbation
theory expression for $n_{I}$ (see Eq. (B1) in appendix B of Ref.
Brandt:2018bwq ). For our present ensembles this fit yields values of
$f_{\pi}=136(2)$ and $130(3)$ MeV, respectively, in good agreement with the
physical value and the result obtained from the fit in Ref. Brandt:2018bwq for
the lattice with $a\approx 0.29$ fm. To obtain a smooth interpolation for the
isospin density, we again perform a model independent spline interpolation of
the remaining data points via a spline Monte-Carlo as discussed above, where
all splines are matched to chiral perturbation theory.444For a smooth matching
up to the second $\mu_{I}$-derivative, relevant for $c_{s}$, the matching of
the spline is done at half the distance between the second data point in the
BEC phase and the BEC phase boundary and we have included 20 additional data
points generated from chiral perturbation theory, equally separated in the
remaining interval up to the second data point, in the fit. The latter is
relevant for a smooth matching of the second derivative. The resulting
interpolation of the isospin density is shown in the left panel of Fig. 10.

At $T=0$ the condition of Eq. (9) is trivially fulfilled since $s$ vanishes.
Thus, the directional derivative is equivalent to the $\mu_{I}$-derviative,
$\partial_{\xi}=\partial_{\mu_{I}}$. The resulting derivatives of $p$ and
$\epsilon$ in Eq. (7) can again be computed analytically. The results for the
square of the isentropic speed of sound are shown in the right panel of Fig.
10, together with the conformal bound Cherman:2009tw as a gray dashed line. We
observe that the squared speed of sound crosses the conformal bound at
$\mu_{I}/m_{\pi}\approx 0.643(4)$ and $\mu_{I}/m_{\pi}\approx 0.646(5)$, for
$a\approx 0.22$ and 0.15 fm, respectively, and reaches a peak at

$\left.\frac{\mu_{I}}{m_{\pi}}\right|_{{\rm
max}(c_{s})}=\left\\{\begin{array}[]{ll}0.76(4)&\quad\text{with}\quad{\rm
max}(c_{s}^{2})=0.55(3)\quad\text{for}\quad a\approx
0.22\>\text{fm}\vspace*{2mm}\\\ 0.78(5)&\quad\text{with}\quad{\rm
max}(c_{s}^{2})=0.57(3)\quad\text{for}\quad a\approx
0.15\>\text{fm}\,.\end{array}\right.$ (14)

We note that our finding of $c_{s}^{2}>1/3$, as well as the development of a
peak is in good agreement with recent results obtained in two-color QCD
Iida:2022hyy ; Itou:2022ebw (see also Kojo:2021hqh ). Furthermore, similar
peaks in the speed of sound appear in quarkyonic models McLerran:2018hbz ;
Jeong:2019lhv ; Kovensky:2020xif ; Kojo:2021ugu . At larger $\mu_{I}$, the
speed of sound decreases and, on general grounds and according to perturbation
theory Annala:2019puf , is expected to approach the conformal bound
asymptotically from below. For our values of $\mu_{I}$, we currently do not
see the decrease below the conformal bound. This would require simulations at
yet higher isospin chemical potentials.

Figure 11:  Left: Results for the isentropic speed of sound at $T=0$ versus
the energy density in the BEC phase. Also shown is the conformal bound
Cherman:2009tw (gray dashed line). Right: Results for the polytropic index
$\gamma$ at $T=0$. Also shown is the result from chiral perturbation theory
(yellow dashed line) and the “quark matter bound” introduced in Ref.
Annala:2019puf (gray dashed line). Bottom: Results for the normalized
conformal anomaly $\Delta$ at $T=0$, together with the result from chiral
perturbation theory (yellow dashed line).

To allow for contact with recent studies on the EoS in neutron stars
Annala:2019puf ; Somasundaram:2021clp ; Annala:2021gom ; Altiparmak:2022bke ;
Ecker:2022xxj ; Marczenko:2022jhl , we plot $c_{s}^{2}$ versus the energy
density in the upper left panel of Fig. 11. Comparing to typical energy
densities reached in the most massive neutron star cores (which are of the
order of $10^{3}\text{ MeV}/\text{fm}^{3}$, see e.g. Refs.
Somasundaram:2021clp ; Marczenko:2022jhl ), we see that the speed of sound
reaches values which are larger than the conformal bound already for around
one to two orders of magnitude smaller energy densities. In the upper right
panel of Fig. 11, we show the polytropic index versus $\mu_{I}$. Just at the
onset of pion condensation, $\gamma$ assumes a value of $\gamma\geq 2$,
greater than the “quark matter bound” introduced in Ref. Annala:2019puf , in
agreement with the prediction from chiral perturbation theory (yellow dashed
curve). It then increases with $\mu_{I}$ until it reaches its maximal value
$\gamma\approx 2.4$ around $\mu_{I}/m_{\pi}$ between 0.67 and 0.72, depending
on the lattice spacing. Further increasing $\mu_{I}$, the polytropic index
decreases below 1.75 and is seen to approach its conformal value of $\gamma=1$
asymptotically. We note that the crossing of the “quark matter bound” might
provide an alternative definition for the BEC-BCS crossover Brandt:2019hel ;
Cuteri:2021hiq , where effective degrees of freedom change from pions to
Cooper pairs of $u$ and $\bar{d}$ quarks.

Finally, the normalized trace anomaly $\Delta$ is plotted in the bottom panel
of Fig. 11. It starts at $1/3$ at the onset of the BEC phase, in good
agreement with chiral perturbation theory, and decreases towards larger
$\mu_{I}$. Eventually it becomes negative between $\mu_{I}/m_{\pi}$ of 0.85 to
0.9 on the border of our parameter interval. Confirming the prediction of
chiral perturbation theory, this shows a specific counter-example to the claim
that $\Delta$ would be strictly positive in QCD.

###  4.3 Speed of sound at nonzero temperatures

Figure 12:  Results for the isentropic speed of sound at $T\neq 0$, obtained
on the $N_{t}=8,\,10$ and 12 lattices (from top left to bottom left) in the
$T$-$\mu_{I}$ plane, excluding uncertainties.

Figure 13:  Results for the isentropic speed of sound at $T\neq 0$ at
different temperatures on the $8\times 24^{3}$ lattices (left panel) and at a
temperature of 120 MeV for $N_{t}=8$, 10 and 12 (right panel).

The results for the isentropic speed of sound at $T\neq 0$ for the different
lattices are shown in Fig. 12. For small temperatures, the speed of sound
initially decreases slightly in the vicinity of the BEC phase boundary, before
it starts to rise within the BEC phase. For $N_{t}=8$ and 10 it crosses the
conformal bound around $\mu_{I}/m_{\pi}\approx 0.75$ to 0.8, i.e., at a
somewhat larger $\mu_{I}$ value than at $T=0$. This is also visible in the
left panel of Fig. 13, where we show the speed of sound obtained on the
$N_{t}=8$ lattice for different temperatures including the uncertainties.
Contrary to the $T=0$ case, we do not observe a maximum for $c_{s}^{2}$,
likely due to the fact that it appears on the border or outside of our
$\mu_{I}$ range.

At this point we note that the speed of sound depends on the derivatives of
$n_{I}$, which are not well determined at the borders of our interpolation
region at large $\mu_{I}$, as well as at our largest and smallest temperatures
(at $\mu_{I}=0$ it is determined by the parameterization from Ref.
Borsanyi:2013bia , see the lower left panel of Fig. 15). Consequently, the
results for the speed of sound have to be taken with care beyond
$\mu_{I}/m_{\pi}\gtrsim 0.85$. To allow the reader to scrutinize the
uncertainties at low $T$ but large $\mu_{I}$, we show in the right panel of
Fig. 13 $c_{s}^{2}$ including its uncertainties for the different $N_{t}$ at a
temperature of $T=120$ MeV, close to the lower border of the temperature
range. For $N_{t}=12$, the speed of sound does not reach the conformal bound,
but we can still observe an increase of $c_{s}$ towards larger values of
$\mu_{I}$. We interpret that as a shift of the maximum towards larger values
of $\mu_{I}$ when we approach the continuum. The question of the presence of
the peak in the continuum limit can be answered once the results on finer
lattices and at larger $\mu_{I}$ become available.

Figure 14:  Left: Results for the isentropic speed of sound versus the energy
density at $T=122$ MeV obtained on the $8\times 24^{3}$ lattices. Also shown
is the conformal bound Cherman:2009tw (gray dashed line). Right: Results for
the polytropic index $\gamma$ at $T=122$ MeV. Also shown is the “quark matter
bound” introduced in Ref. Annala:2019puf (gray dashed line).

Taking a closer look at the results from $N_{t}=8$, our most accurate results
concerning the extraction of $c_{s}$, albeit being furthest from the
continuum, we show the squared speed of sound versus the energy density and
the polytropic index $\gamma$ versus $\mu_{I}$ for a small temperature of
$T=122$ MeV in Fig. 14. Comparing the results for $c_{s}$ versus the energy
density against those obtained at $T=0$, shown in Fig. 11, we observe that the
speed of sound crosses the conformal bound at about five times larger energy
densities. Contrary to what is seen at $T=0$, the polytropic index starts from
a comparably small value $\gamma=0.82(4)$, and does not increase directly at
the BEC phase boundary. The drastic increase towards larger values happens at
$\mu_{I}/m_{\pi}\approx 0.7$ and it crosses the “quark matter bound” at around
$\mu_{I}/m_{\pi}\approx 0.85$. Within our range of chemical potentials we see
no sign of a maximum or the onset of a plateau.

##  5 Discussion and conclusions

In this paper we studied the QCD equation of state at nonzero temperatures and
isospin chemical potentials using first-principles lattice simulations at
physical quark masses. The results are shown in Figs. 2 and 4. We observe a
strong rise of the pressure within the phase of condensing charged pions (BEC
phase) at small temperature, which becomes milder when approaching the
boundary of the BEC phase with $T$. We mention that our results for the
pressure might also be useful to constrain the EoS at other chemical
potentials. In our simulations the quark chemical potentials are set as
$\mu_{u}=-\mu_{d}=\mu_{I}$ and $\mu_{s}=0$. Since flipping the sign of
$\mu_{d}$ merely amounts to a phase change in the fermion determinant, it is
simple to see that for the partition functions
$\mathcal{Z}(\mu_{u}=-\mu_{d})>\mathcal{Z}(\mu_{u}=\mu_{d})$ holds. This
implies that the pressure $\propto\log\mathcal{Z}$ at nonzero $\mu_{I}$
provides an upper bound for the pressure in a setup with equal light quark
chemical potentials.

The interaction measure was found to initially rise in the BEC phase, before
it reaches a maximum and decreases. For small temperatures it eventually
becomes negative deep in the BEC phase, providing an explicit counter-example
to general positivity arguments Fujimoto:2022ohj . This effect diminishes as
$T$ grows, and around the high-$T$ boundary of the BEC phase the impact of
$\mu_{I}$ is to shift the behavior of $I(T)$ to lower temperatures. We also
determine the QCD phase diagram in the $T$-$\mu_{I}$ plane, which for
$N_{t}=12$ is shown in Fig. 9.

We have put a particular focus on the determination of the isentropic speed of
sound $c_{s}$, for which we show results both at zero and non-vanishing
temperatures. At $T=0$, $c_{s}$ initially increases, crosses the conformal
bound of $1/\sqrt{3}$ around $\mu_{I}/m_{\pi}\approx 0.64$ and reaches a peak
at $\mu_{I}/m_{\pi}\approx 0.77$ with a maximum of $c_{s}^{2}\approx 0.56$,
see Eq. (14) and Fig. 10. It then decreases again and is expected to approach
the conformal bound asymptotically from below. To our knowledge, this is the
first evidence for the explicit violation of this bound in first principles
QCD. We note that similarly large $c_{s}$ values have also been observed at
large $\mu_{B}$ with the functional renormalization group approach in the
point-like approximation for four-quark interactions Leonhardt:2019fua . A
remnant of the peak in $c_{s}$ remains visible at low temperatures (see Fig.
12), where it is shifted towards larger values of $\mu_{I}$. In Fig. 1 we also
show a sketch of the region where the speed of sound exceeds the conformal
bound in the phase diagram.

Besides $c_{s}$, we also computed the polytropic index $\gamma$, which has
been discussed as an indicator for the state of matter in the core of neutron
stars. The results are shown in Figs. 11 and 14 (right panels). At $T=0$, the
polytropic index starts from $\gamma=2$, the value predicted by chiral
perturbation theory Son:2000xc , and then drops below $\gamma=1.75$, the
“quark matter bound” introduced in Annala:2019puf , around
$\mu_{I}/m_{\pi}\approx 0.8$, eventually approaching the conformal value of
$\gamma=1$ asymptotically. At $T\neq 0$ the behavior is quite different. The
polytropic index starts from the $\mu_{I}=0$ value, around $\gamma=0.8$ to 1.0
– see Fig. 15 bottom right – from where it increases and approaches larger
values of $\gamma\gtrsim 1.75$ only around $\mu_{I}/m_{\pi}\gtrsim 0.85$.

Our tabulated results for the EoS will be useful for comparison to low-energy
models and effective theories of QCD. To this end, our data for the full EoS,
including a code to compute the observables, are available with the published
version of this paper. Note that our calculations rely on three different
lattice spacings, but owing to enhanced lattice artifacts at low temperatures
for certain observables, we did not carry out a full continuum extrapolation
here.

We finally comment on the consequences of our finding of the excess of the
speed of sound over the conformal bound for the modeling of the EoS of neutron
stars. With increasing amount of data on the masses and radii of the observed
neutron stars in the Universe, several groups started to extract information
on the QCD EoS at nonzero baryon density using this experimental data. In
these approaches, the EoS is typically constrained at small and large
densities from effective hadronic models (e.g. Gandolfi:2009fj ; Tews:2012fj )
and perturbative QCD (e.g. Gorda:2018gpy ), respectively and then interpolated
using a set of basis functions (e.g. Hebeler:2013nza ; Kurkela:2014vha ;
Annala:2017llu ; Most:2018hfd – see also Han:2021kjx ; Han:2022rug for non-
parametric interpolations using neural networks). Newer studies use a large
set of different types of basic functions and millions of different EoS
interpolations Annala:2019puf ; Somasundaram:2021clp ; Annala:2021gom ;
Altiparmak:2022bke ; Ecker:2022xxj ; Marczenko:2022jhl . While most recent
studies indicate that experimental constraints favour a stiff EoS with a speed
of sound that exceeds the conformal limit Somasundaram:2021clp ;
Altiparmak:2022bke ; Ecker:2022xxj ; Marczenko:2022jhl ; Han:2022rug , it has
often been considered as extreme for the EoS to develop large speeds of sound
of $c_{s}^{2}>0.5$ to 0.6 or even to have an EoS which exceeds the conformal
bound. In our study we provide direct evidence that an EoS with a speed of
sound of this magnitude exists in QCD at small temperatures. Thus, such
conditions for the EoS of cold dense QCD matter are certainly not unrealistic.

Acknowledgments:
The authors are grateful to Szabolcs Borsányi, Gergely Markó, Guy Moore and
Aleksi Vuorinen for useful discussions and to Kálmán Szabó for providing the
parameterization for the lattice QCD EoS at $\mu_{I}=0$. We also thank Prabal
Adhikari, Jens Oluf Andersen and Martin Mojahed for discussions and for
providing the chiral perturbation theory data from Ref. Adhikari:2020kdn . The
authors acknowledge support by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) through the CRC-TR 211 “Strong- interaction matter
under extreme conditions” – project number 315477589 – TRR 211. F.C.
acknowledges the support by the State of Hesse within the Research Cluster
ELEMENTS (Project ID 500/10.006). The authors also gratefully acknowledge the
Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this
project by providing computing time on the GCS Supercomputer SuperMUC-NG at
Leibniz Supercomputing Centre (www.lrz.de). Parts of the computations in this
work were performed on the GPU cluster at Bielefeld University and at Goethe-
HLR at Goethe-University Frankfurt. We thank the computing staff of both
institutions for their support.

##  Appendix A Zero density input to the EoS

$h_{0}$ | $h_{1}$ | $h_{2}$ | $f_{0}$ | $f_{1}$ | $f_{2}$ | $g_{1}$ | $g_{2}$
---|---|---|---|---|---|---|---
0.1396(26) | -0.179(9) | 0.035(1) | 2.76(65) | 6.79(24) | -5.29(17) | -0.47(19) | 1.04(17)
Table 1:  Parameters for the parameterization of the $\mu_{I}=0$ EoS of Eq.
(15), obtained from the reanalysis of the data of Ref. Borsanyi:2013bia as
discussed in the text.

To compute the full EoS from the decompositions of Eq. (3), we need input at
$\mu_{I}=0$. Here we use the parameterization

$\frac{I(T,0)}{T^{4}}=e^{-h_{1}/t-h_{2}/t^{2}}\,\Big{(}h_{0}+\frac{f_{0}\,\big{[}\tanh(f_{1}\,t+f_{2})+1\big{]}}{1+g_{1}t+g_{2}t^{2}}\Big{)},$
(15)

with $t=T/(200\>\text{MeV})$, which has been employed in Refs. Borsanyi:2010cj
; Borsanyi:2013bia . Since we need to take the full correlations between the
parameters into account for the correct computation of the uncertainties of
derivatives and integrals, we use the parameters obtained from a reanalysis of
the data from Ref. Borsanyi:2013bia . The resulting parameters are listed in
Tab. 1. We note the slight differences in the parameters compared to the ones
obtained in Ref. Borsanyi:2013bia . These differences in the reanalysis can be
attributed to flat directions in parameter space, as already mentioned in Ref.
Borsanyi:2013bia , and do not lead to significant changes in the description
of the data and the curves for the thermodynamic observables. For
completeness, we show the results for the pressure, the interaction measure,
the squared speed of sound and the polytropic index for this parameterization
in the parameter region relevant for this study in Fig. 15.

Figure 15:  Pressure, interaction measure, speed of sound squared and
polytropic index in the temperature range relevant for this study, obtained
from the parameterization at $\mu_{I}=0$ of Ref. Borsanyi:2013bia discussed in
the main text.

##  Appendix B EoS from direct interaction measure calculations

An alternative to computing the EoS from the interpolation of the isospin
density, is to rewrite the interaction measure via the derivative of the
partition function with respect to the lattice scale, along similar lines as
in Refs. Allton:2003vx ; Iida:2022hyy ; Itou:2022ebw . Using Eq. (6) and the
derivatives of the lattice parameters with respect to the lattice scale, one
obtains

$I(T,\mu_{I})=\frac{\partial\beta}{\partial\log a}\left\langle
S_{g}\right\rangle_{T,\mu_{I}}-\sum_{q}\frac{\partial(am_{q})}{a\,\partial\log
a}\left\langle\bar{\psi}\psi_{q}\right\rangle_{T,\mu_{I}}\,,$ (16)

where $\beta$ is the lattice coupling, $S_{g}$ the gauge action and $am_{q}$
the bare quark mass of flavour $q$ in lattice units. We note that all the
quantities appearing in Eq. (16) need to be renormalized properly, demanding,
for instance, the knowledge of the quantities at $T=0$, but nonzero $\mu_{I}$.
Another way to ensure a proper renormalization is to make use of the
decomposition of Eq. (3) and to compute only $\Delta I$ instead of $I$.
Defining generically $\Delta O=\left\langle
O\right\rangle_{T,\mu_{I}}-\left\langle O\right\rangle_{T,0}$, we obtain

$\Delta I(T,\mu_{I})=\frac{\partial\beta}{\partial\log a}\Delta
S_{g}-\sum_{q}\frac{\partial(am_{q})}{a\,\partial\log
a}\Delta\bar{\psi}\psi_{q}\,.$ (17)

The remaining task is the computation of $\Delta S_{g}$ and
$\Delta\bar{\psi}\psi_{q}$.

Figure 16:  Left: Results for the interaction measure obtained from a direct
computation explained in the text (blue points) and obtained from the two-
dimensional interpolation of the isospin density (red band) on $24\times 6$
lattices at a temperature of 136 MeV. Right: $\lambda$-extrapolation (blue
band) of the Symanzik improved gauge action (blue open circles) on a
$24^{3}\times 6$ lattice with $T=136$ MeV and $\mu_{I}/m_{\pi}=0.25$. We also
show the result at $\mu_{I}=0$ for the same temperature (black filled box and
black band) which needs to be subtracted from the $\lambda=0$ result to obtain
the contribution to $\Delta I$.

We show the results for $\Delta I(T,\mu_{I})$ versus $\mu_{I}$ obtained for
one particular temperature on a set of $24^{3}\times 6$ lattices in the left
panel of Fig. 16. For comparison we also show the result for $\Delta I$
obtained from the interpolation of the isospin density. We note that the
qualitative behavior of the two data sets is similar, at least as far as
comparison is possible due to the uncertainties for the direct computation,
which are at least an order of magnitude larger compared to the uncertainties
of the results obtained from the interpolation. The question is, where these
large uncertainties orginate from. The problem is the combination of the
$\lambda$-extrapolations at $\mu_{I}\neq 0$ in combination with the
subtraction of the $\mu_{I}=0$ value (despite the fact that such an
extrapolation is not necessary there). We show the $\lambda$-extrapolation of
the Symanzik improved gauge action in the right panel of Fig. 16 in comparison
to the result at $\mu_{I}=0$ at the same temperature, which when subtracted
give $\Delta S_{g}$. As is evident from the plot, the large uncertainties come
from the subtraction of two quantities of similar magnitude, so that relative
uncertainties are enhanced by orders of magnitude compared to the
uncertainties of the individual quantities. This is particularly pronounced
for the gauge action, but a similar behavior is also seen for the quark
condensates. We note, that such a subtraction is absent for the isospin
density, leading to way more accurate results for the quantity to interpolate
and, consequently, for the EoS.

##  Appendix C Isospin density at $\lambda=0$ and model-independent spline
interpolations

The basic data for the isospin density, Eq. (1), is obtained at non-vanishing
pion source, $\lambda\neq 0$. For the extrapolations to $\lambda=0$ we use the
improvement program from Refs. Brandt:2017oyy ; Brandt:2018bwq ;
Brandt:2018omg . The resulting $\lambda$-extrapolations are basically flat and
can be done using either a linear function in $\lambda^{2}$ or a constant. As
an estimate for the systematic uncertainty associated with the extrapolation,
we use the maximal deviation of the final result with either the extrapolation
using the alternative (linear or constant), functional form or any of the two
data-points at the lowest $\lambda$ values. Note, that this is particularly
important for the extraction of the equation of state, since an
underestimation of uncertainties might lead to unphysical fluctuations which
significantly affect the interpolation using spline fits.

For the extraction of the equation of state we use an average over all
possible cubic spline interpolations of the isospin density in the two-
dimensional parameter space $(T,\,\mu_{I})$ with less grid- than data points
(i.e., spline fits), weighted with an estimator for the ‘‘goodness’’ of the
spline fit. Note, that for our spline fits the positions of the spline
nodepoints generically do not coincide with the positions of the data points.
As already discussed in Ref. Brandt:2016zdy , this average for an observable
$A$ (for instance the isospin density $n_{I}(T,\mu_{I})$ for given values $T$
and $\mu_{I}$) can be written as

$\big{\llangle}A\big{\rrangle}=\sum_{N_{G}}\int
d^{N_{x}(N_{G})}x\,\,A(\vec{x})\,\exp\big{(}-S_{\rm
spl}(\vec{x},N_{G})\big{)}\,.$ (18)

Here $N_{G}$ is the total number of spline nodepoints and $N_{x}(N_{G})$ the
number of nodepoints which can be varied in the particular spline setup. Note
that $N_{x}$ and $N_{G}$ do not need to be equivalent (but always $N_{x}\leq
N_{G}$), since some of the nodepoint positions can be fixed. This is the case
for the nodepoints at the lower end of the $n_{I}$-splines in the $\mu_{I}$
direction, for instance, which are kept at $\mu_{I}=0$. In Eq. (18), $\vec{x}$
is the vector of (two-dimensional) nodepoint positions for the $N_{x}$
variable nodepoints. Note that typically the allowed range for the nodepoint
values is restricted, as outlined below. The action $S_{\rm
spl}(\vec{x},N_{G})$ represents the estimate for the ‘‘goodness’’ of the
spline fit. Possible choices have already been discussed in Ref.
Brandt:2017oyy . As the basic action we use the Akaike information criterion
Akaike1973InformationTA (see also Ref. Akaike:1998zah ),

$S_{\rm AIC}=2N_{P}+\chi^{2}\,,$ (19)

where $N_{P}$ is the number of parameters of the fit.

One of the major problems for any spline interpolation or spline fit is the
possible appearance of oscillatory solutions, i.e., solutions with additional
minima and maxima as the spline attempts to capture all of the datapoints.
These solutions can in particular be triggered by statistical fluctuations of
data points and are particularly problematic for the equation of state, since
additional unphysical minima and maxima might have strong effects on
quantities like the speed of sound. To suppress those solutions we include
another term in the action, following the spirit of Ref. Endrodi:2010ai . The
term signifies the stability of the spline solution under small variations of
the nodepoints. The parameters of the spline $f_{n}$ (with $n=1,\ldots,N_{P}$)
are given either by the value of the spline on one of the spline nodepoints or
by the derivatives on the nodepoints of the spline boundaries, depending on
the particular spline setup. If we vary the nodepoints slightly and have a
stable, non-oscillatory spline solution, we expect those values to not change
significantly. I.e., given a variation $\alpha$ of one of the nodepoints, here
with nodepoint index $k$ and the variation in direction $i$,

$(x^{\alpha}_{k})_{i}=(x_{k})_{i}+\epsilon\qquad\text{and}\qquad(x^{\alpha}_{l})_{j}^{\prime}=(x_{l})_{j}\>\>\forall
l\neq k,\,j\neq i,$ (20)

where $\epsilon$ is a small (not necessarily positive) number compared to the
typical distance between two datapoints, we expect the spline parameters
$f^{\alpha}_{n}$ to differ only slightly from the previous parameters $f_{n}$.
The parameter variation with respect to the typical statistical uncertainties
for the parameters can be estimated by

$\mathcal{D}^{\alpha}=\frac{1}{N_{P}}\sum_{n}\frac{|f^{\alpha}_{n}-f_{n}|}{\sigma(f_{n})}\,.$
(21)

Here $\sigma(f_{n})$ is the statistical uncertainty of parameter $f_{n}$,
obtained from applying the same spline fit to the individual bootstrap samples
for the data points. For stable fits we expect $\mathcal{D}^{\alpha}$ to be a
number not much larger than one, of course depending on the typical order of
magnitude of the statistical uncertainties and the typical change of the value
of two consecutive data points. To suppress unwanted oscillatory solutions in
the sum of Eq. (18), we add to the action the average of
$\mathcal{D}^{\alpha}$ over all possible spline variations $\alpha$,

$S_{\rm STAB}=\frac{1}{N_{G}}\sum_{\alpha}\mathcal{D}^{\alpha}\,,$ (22)

so that the total action is given by

$S_{\rm spl}=S_{\rm AIC}+\gamma S_{\rm STAB}$ (23)

with the tunable parameter $\gamma$.

The tunable parameters of the spline average outlined above concern the
possible numbers of nodepoints in each direction, possible constraints on the
nodepoint locations, the boundary conditions of the spline, i.e., for each
spline boundary one derivative for a cubic spline, as used here, and the
parameter $\gamma$ as well as the size of $\epsilon$. For the interpolation of
the isospin density, we use three to five nodepoints in each direction and
demand that always two data points reside between the outermost and the
consecutive nodepoint on each border of the grid and at least one data point
lies between two consecutive nodepoints in each direction. The application of
these constraints to the spline grids is straightforward if the data points
themselves form a rectangular grid. The outer nodepoints in $\mu_{I}$
direction are fixed at $\mu_{I}=0$ and on this whole boundary we impose
$n_{I}=0$. To account for the isospin density being an uneven function in
$\mu_{I}$, we also impose
$\partial^{2}n_{I}/(\partial\mu_{I})^{2}|_{\mu_{I}=0}=0$. The second
derivatives in all other directions have been kept as free parameters for the
spline fit and the positions of the outer gridpoints are allowed to vary.

To efficiently perform the sum from Eq. (18), we use Monte-Carlo methods as
proposed in Ref. Brandt:2017oyy . In particular, we employ a Metropolis
algorithm with a symmetric proposal probability for changes in the spline
nodepoints. For $\epsilon$ we have chosen a random number between a tenth of
the distance between the two nearest data points in positive or negative
direction in such a way, that the nodepoint remains between the two data
points after variation. To tune the parameter $\gamma$, we have performed
several runs starting from small values of $\gamma$ and monitored the
resulting splines in the Markov chain. We stopped increasing $\gamma$ when we
found that no significant oscillations leading to local minima/maxima
structures between two nodepoints showed up in the final average over spline
configurations. As the final value we choose $\gamma=10$. For our final
results we have first tuned the number of nodepoints to the optimal value by
performing 100 independent thermalisations allowing changes in the number of
nodepoints and the nodepoint locations with respect to the data points. We
then search for the spline with the lowest action and restrict ourselves to
these number of nodepoints and the intervals in which the spline nodepoints
reside with respect to the grid of data points. The final results are then
obtained from 100 splines obtained by 20 independent runs where we vary the
nodepoints in this constrained setup with 20000 thermalisation updates and
with 5 spline configurations each, separated by 10000 spline updates.

$N_{t}$ | $N_{G}^{(T)}$ | $N_{G}^{(\mu_{I})}$ |  $T$-intervals [MeV] |  $\mu_{I}$-intervals [MeV]
---|---|---|---|---
8 | 3 | 4 |  $<114$, $[142,147]$, $>174$ |  $[34,51]$, $[51,68]$, $>119$
10 | 3 | 4 |  $<114$, $[158,168]$, $>179$ |  $[51,68]$, $[68,85]$, $>119$
12 | 3 | 3 |  $<114$, $[139,149]$, $>179$ |  $[34,51]$, $>119$
Table 2: Number of spline nodepoints in the $T$ and $\mu_{I}$ directions and
the associated intervals in which the Monte-Carlo-generated nodepoints reside.

Figure 17:  Summary of the location of spline nodepoints in the final run.
Note that the spline nodepoints can only be generated in certain regions of
the parameter space, indicated by the accumulation of dots.

Figure 18:  Comparison of the results for the isospin density $n_{I}$ versus
$\mu_{I}$ for different temperatures (left) and $T$ for different $\mu_{I}$
(right) obtained from the simulations on $8\times 24^{3}$ lattices and the
associated spline interpolation described in the text.

The final numbers of nodepoints in the different directions
$N_{G}=N_{G}^{(T)}\times N_{G}^{(\mu_{I})}$, as well as the intervals in which
the nodepoint positions have been varied are given in Tab. 2. We show the
location of the 100 final nodepoint sets used in the analysis in Fig. 17. Note
that the outer spline nodepoints can lie anywhere outside of the data point
interval and that the outer nodepoint of the lower boundary in $\mu_{I}$
direction has been held fixed at $\mu_{I}=0$. The resulting interpolation for
$N_{t}=8$ lattices is shown for a set of temperatures in Fig. 18. The
uncertainties include the uncertainty due to the individual data points
(computed using the bootstrap procedure with 1000 samples) and from the Monte-
Carlo over spline interpolations. As a comment concerning the $T=0$ splines,
in this case we are speaking about a one-dimensional spline interpolation,
which is rather well behaved and uncritical concerning the spline Monte-Carlo.
In this case, the number of nodepoints and their location has been allowed to
vary freely up to the maximal possible number of nodepoints given the number
of available data points with the matching to chiral perturbation theory as
discussed in section 4.2.

##  Appendix D Simulation Details

We provide the run parameters for the simulations at $T\neq 0$ and $T=0$ in
Tab. 3. For each of the mentioned parameter values we have simulated at up to
five different values of the pion source parameter $\lambda$ with
$1.4>\lambda/m_{ud}\gtrsim 0.05$. At $T=0$ we have used a fixed number of
three different $\lambda$ values. For each of these parameter sets we
generated between 500 and 2000 trajectories, measuring observables (here the
isospin density $n_{I}$) on every fifth configuration. Typically we use two to
three independent chains to acquire the full statistics. We estimate
autocorrelations using the integrated autocorrelation time of the plaquette
expectation value. The values we obtain on the different ensembles are mostly
between 5 and 10 in molecular dynamics units, with the tendency to large
values for smaller temperatures and $\lambda$-values, as well as for larger
values of $\mu_{I}$ and $N_{t}$. While this indicates that typically two
consecutive configurations are correlated, the autocorrelation times obtained
for $n_{I}$, the main observable of our study, are much smaller. In rare cases
we observed autocorrelation times of the order of 20 in molecular dynamics
units, so that 4 consecutive measurements are correlated. Generically, we have
checked the error analysis by using binning prior to the bootstrap procedure,
but did not observe a significant dependence of the uncertainties on the
binsize. Furthermore, in the improved $\lambda$-extrapolation, values at
different $\lambda$ are combined using a linear function to extract the
$\lambda=0$ result, averaging out the fluctuations of the individual
simulation points.

$N_{t}$ | $N_{s}$ | $T[\text{MeV}]$ | $\beta$ | $a\mu_{I}$
---|---|---|---|---
8 | 24 | 114 | 3.5500 | 0.019 0.037 0.056 0.075 0.094 0.113 0.131
8 | 24 | 123 | 3.5750 | 0.018 0.034 0.052 0.070 0.075 0.087 0.105 0.122
8 | 24 | 132 | 3.6000 | 0.016 0.032 0.048 0.065 0.075 0.081 0.098 0.113
8 | 24 | 137 | 3.6120 | 0.016 0.031 0.047 0.063 0.079 0.094 0.109
8 | 24 | 142 | 3.6250 | 0.015 0.030 0.045 0.060 0.075 0.076 0.091 0.106
8 | 24 | 147 | 3.6370 | 0.015 0.029 0.044 0.058 0.073 0.088 0.102
8 | 24 | 152 | 3.6500 | 0.014 0.028 0.042 0.056 0.071 0.075 0.085 0.098
8 | 24 | 163 | 3.6750 | 0.013 0.026 0.039 0.053 0.066 0.075 0.079 0.092
8 | 24 | 174 | 3.7000 | 0.012 0.024 0.037 0.049 0.062 0.074 0.075 0.086
10 | 28 | 114 | 3.6250 | 0.015 0.030 0.045 0.060 0.075 0.090 0.105
10 | 28 | 122 | 3.6500 | 0.014 0.028 0.042 0.056 0.070 0.084 0.098
10 | 28 | 130 | 3.6750 | 0.013 0.026 0.039 0.052 0.065 0.078 0.091
10 | 28 | 139 | 3.7000 | 0.012 0.024 0.037 0.049 0.061 0.073 0.086
10 | 28 | 149 | 3.7250 | 0.011 0.023 0.034 0.046 0.057 0.069 0.080
10 | 28 | 158 | 3.7500 | 0.011 0.022 0.032 0.043 0.054 0.065 0.075
10 | 28 | 168 | 3.7750 | 0.010 0.020 0.030 0.040 0.051 0.061 0.071
10 | 28 | 179 | 3.8000 | 0.010 0.019 0.029 0.038 0.048 0.057 0.067
12 | 36 | 114 | 3.6900 | 0.012 0.025 0.037 0.050 0.063 0.075 0.088
12 | 36 | 122 | 3.7200 | 0.011 0.023 0.034 0.046 0.058 0.069 0.081
12 | 36 | 132 | 3.7500 | 0.010 0.021 0.032 0.043 0.054 0.064 0.075
12 | 36 | 140 | 3.7750 | 0.010 0.020 0.030 0.040 0.051 0.060 0.071
12 | 36 | 149 | 3.8000 | 0.009 0.019 0.028 0.038 0.048 0.057 0.067
12 | 36 | 158 | 3.8250 | 0.009 0.018 0.027 0.036 0.045 0.054 0.063
12 | 36 | 167 | 3.8500 | 0.008 0.017 0.025 0.034 0.043 0.051 0.060
12 | 36 | 179 | 3.8800 | 0.008 0.016 0.023 0.032 0.040 0.048 0.056
32 | 24 | 29 | 3.5500 | 0.067 0.074 0.082 0.093 0.112 0.130
48 | 32 | 27 | 3.6700 | 0.053 0.059 0.066 0.080 0.092 0.106
Table 3: Summary of simulation parameters used for the determination of the
EoS.

## References

  * (1) Y. Aoki, G. Endrődi, Z. Fodor, S.D. Katz and K.K. Szabó, _The Order of the quantum chromodynamics transition predicted by the standard model of particle physics_ , _Nature_ 443 (2006) 675 [hep-lat/0611014].
  * (2) T. Bhattacharya et al., _QCD Phase Transition with Chiral Quarks and Physical Quark Masses_ , _Phys. Rev. Lett._ 113 (2014) 082001 [1402.5175].
  * (3) D. Teaney, J. Lauret and E.V. Shuryak, _A Hydrodynamic Description of Heavy Ion Collisions at the SPS and RHIC_ , nucl-th/0110037.
  * (4) P.F. Kolb and U.W. Heinz, _Hydrodynamic description of ultrarelativistic heavy ion collisions_ , nucl-th/0305084.
  * (5) D. Boyanovsky, H.J. de Vega and D.J. Schwarz, _Phase transitions in the early and the present universe_ , _Ann. Rev. Nucl. Part. Sci._ 56 (2006) 441 [hep-ph/0602002].
  * (6) J.M. Lattimer and M. Prakash, _Neutron star structure and the equation of state_ , _Astrophys. J._ 550 (2001) 426 [astro-ph/0002232].
  * (7) I. Tews, J. Carlson, S. Gandolfi and S. Reddy, _Constraining the speed of sound inside neutron stars with chiral effective field theory interactions and observations_ , _Astrophys. J._ 860 (2018) 149 [1801.01923].
  * (8) E. Annala, T. Gorda, A. Kurkela, J. Nättilä and A. Vuorinen, _Evidence for quark-matter cores in massive neutron stars_ , _Nature Phys._ 16 (2020) 907 [1903.09121].
  * (9) I.M. Oldengott and D.J. Schwarz, _Improved constraints on lepton asymmetry from the cosmic microwave background_ , _EPL_ 119 (2017) 29001 [1706.01705].
  * (10) M.M. Wygas, I.M. Oldengott, D. Bödeker and D.J. Schwarz, _Cosmic QCD Epoch at Nonvanishing Lepton Asymmetry_ , _Phys. Rev. Lett._ 121 (2018) 201302 [1807.10815].
  * (11) 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 (2022) 123533 [2009.00036].
  * (12) V. Vovchenko, B.B. Brandt, F. Cuteri, G. Endrődi, F. Hajkarim and J. Schaffner-Bielich, _Pion Condensation in the Early Universe at Nonvanishing Lepton Flavor Asymmetry and Its Gravitational Wave Signatures_ , _Phys. Rev. Lett._ 126 (2021) 012701 [2009.02309].
  * (13) D.T. Son and M.A. Stephanov, _QCD at finite isospin density_ , _Phys. Rev. Lett._ 86 (2001) 592 [hep-ph/0005225].
  * (14) J.B. Kogut and D.K. Sinclair, _Quenched lattice QCD at finite isospin density and related theories_ , _Phys. Rev._ D66 (2002) 014508 [hep-lat/0201017].
  * (15) J.B. Kogut and D.K. Sinclair, _Lattice QCD at finite isospin density at zero and finite temperature_ , _Phys. Rev._ D66 (2002) 034505 [hep-lat/0202028].
  * (16) J.B. Kogut and D.K. Sinclair, _The Finite temperature transition for 2-flavor lattice QCD at finite isospin density_ , _Phys. Rev._ D70 (2004) 094501 [hep-lat/0407027].
  * (17) G. Endrődi, _Magnetic structure of isospin-asymmetric QCD matter in neutron stars_ , _Phys. Rev._ D90 (2014) 094501 [1407.1216].
  * (18) P. de Forcrand, M.A. Stephanov and U. Wenger, _On the phase diagram of QCD at finite isospin density_ , _PoS_ LAT2007 (2007) 237 [0711.0023].
  * (19) P. Cea, L. Cosmai, M. D’Elia, A. Papa and F. Sanfilippo, _The critical line of two-flavor QCD at finite isospin or baryon densities from imaginary chemical potentials_ , _Phys. Rev._ D85 (2012) 094512 [1202.5700].
  * (20) W. Detmold, K. Orginos and Z. Shi, _Lattice QCD at non-zero isospin chemical potential_ , _Phys. Rev._ D86 (2012) 054507 [1205.4224].
  * (21) B.B. Brandt, G. Endrődi and S. Schmalzbauer, _QCD phase diagram for nonzero isospin-asymmetry_ , _Phys. Rev. D_ 97 (2018) 054514 [1712.08190].
  * (22) B.B. Brandt and G. Endrődi, _Reliability of Taylor expansions in QCD_ , _Phys. Rev. D_ 99 (2019) 014518 [1810.11045].
  * (23) B.B. Brandt, G. Endrődi, E.S. Fraga, M. Hippert, J. Schaffner-Bielich and S. Schmalzbauer, _New class of compact stars: Pion stars_ , _Phys. Rev. D_ 98 (2018) 094510 [1802.06685].
  * (24) A. Cherman, T.D. Cohen and A. Nellore, _A Bound on the speed of sound from holography_ , _Phys. Rev. D_ 80 (2009) 066003 [0905.0903].
  * (25) B.B. Brandt, F. Cuteri and G. Endrődi, Dataset for “Equation of state and speed of sound of isospin-asymmetric QCD on the lattice", Bielefeld University (2023), https://pub.uni-bielefeld.de/record/2980217 [DOI:10.4119/UNIBI/2980217].
  * (26) B.B. Brandt, G. Endrődi and S. Schmalzbauer, _QCD at finite isospin chemical potential_ , _EPJ Web Conf._ 175 (2018) 07020 [1709.10487].
  * (27) B.B. Brandt, G. Endrődi and S. Schmalzbauer, _QCD at nonzero isospin asymmetry_ , _PoS_ Confinement2018 (2018) 260 [1811.06004].
  * (28) B.B. Brandt, F. Cuteri and G. Endrődi, _QCD thermodynamics at non-zero isospin asymmetry_ , _PoS_ LATTICE2021 (2022) 132 [2110.14750].
  * (29) S. Borsányi, G. Endrődi, Z. Fodor, A. Jakovác, S.D. Katz, S. Krieg et al., _The QCD equation of state with dynamical quarks_ , _JHEP_ 11 (2010) 077 [1007.2580].
  * (30) S. Borsányi, Z. Fodor, C. Hoelbling, S.D. Katz, S. Krieg and K.K. Szabó, _Full result for the QCD equation of state with 2+1 flavors_ , _Phys. Lett. B_ 730 (2014) 99 [1309.5258].
  * (31) HotQCD collaboration, _Equation of state in ( 2+1 )-flavor QCD_ , _Phys. Rev. D_ 90 (2014) 094503 [1407.6387].
  * (32) J. Engels, J. Fingberg, F. Karsch, D. Miller and M. Weber, _Nonperturbative thermodynamics of SU(N) gauge theories_ , _Phys. Lett. B_ 252 (1990) 625.
  * (33) T. Blum, L. Karkkainen, D. Toussaint and S.A. Gottlieb, _The beta function and equation of state for QCD with two flavors of quarks_ , _Phys. Rev. D_ 51 (1995) 5153 [hep-lat/9410014].
  * (34) J. Engels, R. Joswig, F. Karsch, E. Laermann, M. Lutgemeier and B. Petersson, _Thermodynamics of four flavor QCD with improved staggered fermions_ , _Phys. Lett. B_ 396 (1997) 210 [hep-lat/9612018].
  * (35) Y. Aoki, Z. Fodor, S.D. Katz and K.K. Szabó, _The Equation of state in lattice QCD: With physical quark masses towards the continuum limit_ , _JHEP_ 01 (2006) 089 [hep-lat/0510084].
  * (36) C.R. Allton, S. Ejiri, S.J. Hands, O. Kaczmarek, F. Karsch, E. Laermann et al., _The Equation of state for two flavor QCD at nonzero chemical potential_ , _Phys. Rev. D_ 68 (2003) 014507 [hep-lat/0305007].
  * (37) K. Iida and E. Itou, _Velocity of Sound beyond the High-Density Relativistic Limit from Lattice Simulation of Dense Two-Color QCD_ , 2207.01253.
  * (38) E. Itou and K. Iida, _Bump of sound velocity in dense 2-color QCD_ , in _39th International Symposium on Lattice Field Theory_ , 10, 2022 [2210.14385].
  * (39) H. Akaike, _Information theory and an extension of the maximum likelihood principle_ , 1973.
  * (40) B.B. Brandt and G. Endrődi, _QCD phase diagram with isospin chemical potential_ , _PoS_ LATTICE2016 (2016) 039 [1611.06758].
  * (41) S. Borsányi, G. Endrődi, Z. Fodor, S.D. Katz, S. Krieg, C. Ratti et al., _QCD equation of state at nonzero chemical potential: continuum results with physical quark masses at order $\mu^{2}$_, _JHEP_ 08 (2012) 053 [1204.6710].
  * (42) J.N. Günther, R. Bellwied, S. Borsányi, Z. Fodor, S.D. Katz, A. Pásztor et al., _The QCD equation of state at finite density from analytical continuation_ , _Nucl. Phys. A_ 967 (2017) 720 [1607.02493].
  * (43) D. Bollweg, D.A. Clarke, J. Goswami, O. Kaczmarek, F. Karsch, S. Mukherjee et al., _Equation of state and speed of sound of (2+1)-flavor QCD in strangeness-neutral matter at non-vanishing net baryon-number density_ , 2212.09043.
  * (44) S. Borsányi, Z. Fodor, S.D. Katz, S. Krieg, C. Ratti and K. Szabó, _Fluctuations of conserved charges at finite temperature from lattice QCD_ , _JHEP_ 01 (2012) 138 [1112.4416].
  * (45) P. Adhikari, J.O. Andersen and M.A. Mojahed, _Condensates and pressure of two-flavor chiral perturbation theory at nonzero isospin and temperature_ , _Eur. Phys. J. C_ 81 (2021) 173 [2010.13655].
  * (46) R.R. Silbar and S. Reddy, _Neutron stars for undergraduates_ , _Am. J. Phys._ 72 (2004) 892 [nucl-th/0309041].
  * (47) I. Sagert, M. Hempel, C. Greiner and J. Schaffner-Bielich, _Compact stars for undergraduates_ , _Eur. J. Phys._ 27 (2006) 577 [astro-ph/0506417].
  * (48) Y. Fujimoto, K. Fukushima, L.D. McLerran and M. Praszalowicz, _Trace anomaly as signature of conformality in neutron stars_ , 2207.06753.
  * (49) B.B. Brandt, F. Cuteri and G. Endrődi, _Equation of state and Taylor expansions at nonzero isospin chemical potential_ , in _39th International Symposium on Lattice Field Theory_ , 12, 2022 [2212.01431].
  * (50) T. Kojo and D. Suenaga, _Peaks of sound velocity in two color dense QCD: Quark saturation effects and semishort range correlations_ , _Phys. Rev. D_ 105 (2022) 076001 [2110.02100].
  * (51) L. McLerran and S. Reddy, _Quarkyonic Matter and Neutron Stars_ , _Phys. Rev. Lett._ 122 (2019) 122701 [1811.12503].
  * (52) K.S. Jeong, L. McLerran and S. Sen, _Dynamically generated momentum space shell structure of quarkyonic matter via an excluded volume model_ , _Phys. Rev. C_ 101 (2020) 035201 [1908.04799].
  * (53) N. Kovensky and A. Schmitt, _Holographic quarkyonic matter_ , _JHEP_ 09 (2020) 112 [2006.13739].
  * (54) T. Kojo, _Stiffening of matter in quark-hadron continuity_ , _Phys. Rev. D_ 104 (2021) 074005 [2106.06687].
  * (55) R. Somasundaram, I. Tews and J. Margueron, _Investigating Signatures of Phase Transitions in Neutron-Star Cores_ , 2112.08157.
  * (56) E. Annala, T. Gorda, E. Katerini, A. Kurkela, J. Nättilä, V. Paschalidis et al., _Multimessenger Constraints for Ultradense Matter_ , _Phys. Rev. X_ 12 (2022) 011058 [2105.05132].
  * (57) S. Altiparmak, C. Ecker and L. Rezzolla, _On the Sound Speed in Neutron Stars_ , _Astrophys. J. Lett._ 939 (2022) L34 [2203.14974].
  * (58) C. Ecker and L. Rezzolla, _A General, Scale-independent Description of the Sound Speed in Neutron Stars_ , _Astrophys. J. Lett._ 939 (2022) L35 [2207.04417].
  * (59) M. Marczenko, L. McLerran, K. Redlich and C. Sasaki, _Reaching percolation and conformal limits in neutron stars_ , 2207.13059.
  * (60) B.B. Brandt, F. Cuteri, G. Endrődi and S. Schmalzbauer, _The Dirac spectrum and the BEC-BCS crossover in QCD at nonzero isospin asymmetry_ , _Particles_ 3 (2020) 80 [1912.07451].
  * (61) F. Cuteri, B.B. Brandt and G. Endrődi, _Searching for the BCS phase at nonzero isospin asymmetry_ , _PoS_ LATTICE2021 (2022) 232 [2112.11113].
  * (62) M. Leonhardt, M. Pospiech, B. Schallmo, J. Braun, C. Drischler, K. Hebeler et al., _Symmetric nuclear matter from the strong interaction_ , _Phys. Rev. Lett._ 125 (2020) 142502 [1907.05814].
  * (63) S. Gandolfi, A.Y. Illarionov, K.E. Schmidt, F. Pederiva and S. Fantoni, _Quantum Monte Carlo calculation of the equation of state of neutron matter_ , _Phys. Rev. C_ 79 (2009) 054005 [0903.2610].
  * (64) I. Tews, T. Krüger, K. Hebeler and A. Schwenk, _Neutron matter at next-to-next-to-next-to-leading order in chiral effective field theory_ , _Phys. Rev. Lett._ 110 (2013) 032504 [1206.0025].
  * (65) T. Gorda, A. Kurkela, P. Romatschke, M. Säppi and A. Vuorinen, _Next-to-Next-to-Next-to-Leading Order Pressure of Cold Quark Matter: Leading Logarithm_ , _Phys. Rev. Lett._ 121 (2018) 202701 [1807.04120].
  * (66) K. Hebeler, J.M. Lattimer, C.J. Pethick and A. Schwenk, _Equation of state and neutron star properties constrained by nuclear physics and observation_ , _Astrophys. J._ 773 (2013) 11 [1303.4662].
  * (67) A. Kurkela, E.S. Fraga, J. Schaffner-Bielich and A. Vuorinen, _Constraining neutron star matter with Quantum Chromodynamics_ , _Astrophys. J._ 789 (2014) 127 [1402.6618].
  * (68) E. Annala, T. Gorda, A. Kurkela and A. Vuorinen, _Gravitational-wave constraints on the neutron-star-matter Equation of State_ , _Phys. Rev. Lett._ 120 (2018) 172703 [1711.02644].
  * (69) E.R. Most, L.R. Weih, L. Rezzolla and J. Schaffner-Bielich, _New constraints on radii and tidal deformabilities of neutron stars from GW170817_ , _Phys. Rev. Lett._ 120 (2018) 261103 [1803.00549].
  * (70) M.-Z. Han, J.-L. Jiang, S.-P. Tang and Y.-Z. Fan, _Bayesian Nonparametric Inference of the Neutron Star Equation of State via a Neural Network_ , _Astrophys. J._ 919 (2021) 11 [2103.05408].
  * (71) M.-Z. Han, Y.-J. Huang, S.-P. Tang and Y.-Z. Fan, _Plausible presence of new state in neutron stars with masses above 0.98MTOV_ , _Sci. Bull._ 68 (2023) 913 [2207.13613].
  * (72) H. Akaike, _Information Theory and an Extension of the Maximum Likelihood Principle_ , (New York), Springer Science+Business Media (1998), DOI.
  * (73) G. Endrődi, _Multidimensional spline integration of scattered data_ , _Comput. Phys. Commun._ 182 (2011) 1307 [1010.2952].