# Dense and Hot QCD at Strong Coupling

Tuna Demircik<EMAIL_ADDRESS>Asia Pacific Center for Theoretical
Physics, Pohang, 37673, Korea  Christian Ecker<EMAIL_ADDRESS>Institut für Theoretische Physik, Goethe Universität, Max-von-Laue-Str. 1,
60438 Frankfurt am Main, Germany  Matti Järvinen<EMAIL_ADDRESS>Asia Pacific Center for Theoretical Physics, Pohang, 37673, Korea

###### Abstract

We present a novel framework for the equation of state of dense and hot
Quantum Chromodynamics (QCD), which focuses on the region of the phase diagram
relevant for neutron star mergers and core-collapse supernovae. The model
combines predictions from the gauge/gravity duality with input from lattice
field theory, QCD perturbation theory, chiral effective theory and statistical
modeling. It is therefore, by construction, in good agreement with theoretical
constraints both at low and high densities and temperatures. The main
ingredients of our setup are the non-perturbative V-QCD model based on the
gauge/gravity duality, a van der Waals model for nucleon liquid, and the DD2
version of the Hempel-Schaffner-Bielich statistical model of nuclear matter.
By consistently combining these models, we also obtain a description for the
nuclear to quark matter phase transition and its critical endpoint. The
parameter dependence of the model is represented by three (soft, intermediate
and stiff) variants of the equation of state, all of which agree with
observational constraints from neutron stars and their mergers. We discuss
resulting constraints for the equation of state, predictions for neutron stars
and the location of the critical point.

quantum chromodynamics, gauge/gravity duality, equation of state

††preprint: APCTP Pre2021 - 034

##  I Introduction

Solving QCD at intermediate density and temperature is a long-standing open
problem. Recent and upcoming developments in relativistic heavy ion collision
experiments and astrophysical observations of compact stars as well as their
mergers urgently demand progress in the theoretic modelling of QCD at few
times nuclear saturation density $n_{s}=0.16\,\rm{fm}^{-3}$ and temperatures
up to around $100$ MeV. First principle approaches such as lattice QCD and
perturbation theory are not applicable in the relevant regime, and effective
field theory only works up to densities around the nuclear saturation density.
Despite recent progress in finite temperature chiral effective theory
computations at high loop-order Keller  _et al._ (2021), typical densities and
temperatures such as estimated by realistic binary neutron star merger
simulations Most  _et al._ (2019); Perego  _et al._ (2019); Endrizzi  _et al._
(2020); Hammond  _et al._ (2021); Figura  _et al._ (2021); Raithel  _et al._
(2021) remain currently out of reach.

A central quantity that is absolutely crucial in the modelling of compact
stars is the equation of state (EoS). Due to the aforementioned lack of first-
principles results, the QCD EoS at intermediate densities has currently large
uncertainties at zero temperature, and even less is known about the
temperature dependence. These uncertainties motivate us to formulate a novel
framework, which combines predictions from various different approaches in
different temperature and density regions of the QCD phase diagram where they
are expected to work best. The main idea is to use gauge/gravity duality to
model the physics at large and intermediate densities, and combine this with
effective theory at low densities. This combination allows us to establish a
unified description of QCD matter for a wide range of densities and
temperatures including, but not limited to, the ranges relevant for core-
collapse supernovae and neutron star mergers.

The QCD phase diagram is conjectured to include a critical point where the
nuclear to quark matter transition ends. Ongoing experiments at RHIC (the beam
energy scan) already probe the region of the phase diagram where the critical
point may lie Adam  _et al._ (2021). Future experiments at FAIR Ablyazimov
_et al._ (2017); Durante  _et al._ (2019) and NICA Sissakian and Sorin (2009)
will provide more detailed information about this region, and will reach
substantially higher densities, i.e., densities well above the nuclear
saturation density. Consequently, it is timely to improve theoretical
predictions for the location of the critical point and the EoS in its
vicinity. Our approach leads to EoSs which are in good agreement, among other
things, with lattice data at small density as well as with ab initio
calculations and observational data from neutron star measurements at finite
density and small temperature. Therefore we are able to obtain controlled
interpolations of the EoS to the theoretically challenging region of
intermediate densities and temperatures, and sound estimates for the location
of the critical point.

Figure 1:  A schematic diagram that shows the construction for the temperature
dependence in the model. See text for details.

The model, which we construct in this article, is a thermodynamically
consistent combination of three main approaches: the holographic V-QCD model
Järvinen and Kiritsis (2012); Ishii  _et al._ (2019), an adjusted van der
Waals model of nuclear matter Vovchenko (2020), and the nuclear theory model
Hempel-Schaffner-Bielich (HS) EoS Hempel and Schaffner-Bielich (2010) with DD2
relativistic mean field theory interactions Typel  _et al._ (2010). We also
use the Akmal-Pandharipande-Ravenhall (APR) model Akmal  _et al._ (1998a) of
cold nuclear matter 111We use APR rather than HS(DD2) for cold nuclear matter
near the saturation density because if we used HS(DD2) only, the resulting
EoSs would be so stiff that they would conflict with the LIGO/Virgo
measurement of the tidal deformability for GW170817. and a meson gas model
near the QCD crossover region (see Fig. 1), but they play a smaller role in
the final EoS.

We use V-QCD Järvinen and Kiritsis (2012); Ishii  _et al._ (2019), which is a
non-perturbative model for QCD that is based on the gauge/gravity duality, at
intermediate and high densities, i.e., densities above the nuclear saturation
density. This model allows us to describe both the nuclear and quark matter
phases and therefore the phase transition in a single framework. It has been
shown Ecker  _et al._ (2020); Jokela  _et al._ (2021a) to lead to a feasible
EoS for cold QCD matter which is in good agreement with all available data.
Its extension to finite temperature is however problematic due to a generic
limitation in gauge/gravity duality that the EoS in confined phases, including
the nuclear matter phase, is independent of temperature and therefore not
fully realistic. In the absence of other reliable ways to estimate the
temperature dependence of the EoS for dense nuclear matter we use essentially
the simplest approach: a van der Waals type model of nuclear matter, i.e., a
gas of nucleons and electrons with excluded volume corrections and an
effective potential. The effective potential is tailored for the model to
agree with V-QCD at zero temperature, and therefore the van der Waals model
gives an extrapolation of the V-QCD nuclear matter result to finite
temperatures. The third and final main constituent of the model is the nuclear
theory model at low densities (below and around nuclear saturation density).
In this article, we use the HS(DD2) model.

Fig. 1 shows a sketch of the building blocks and the resulting phase structure
of our model and their dependence on temperature and density, which is the
main topic of this article. This dependence is completed by the dependence on
the charge fraction $Y_{q}$, i.e., allowing deviation from beta-equilibrium,
which is required for realistic simulations of neutron star mergers and core-
collapse supernovae. For this dependence we use the prediction of the HS(DD2)
model in the NM phase and a simple model, arising from the pressure of free
electrons, in the QM phase. We will provide three tabulated variants of the
density, temperature and charge fraction dependent EoSs of this article in the
CompOSE database Typel  _et al._ (2015); com (a).

Finally, let us comment on how the model of this article is related to other
recent research on the EoS of QCD. The cold versions of the model, i.e.,
hybrid EoSs using V-QCD for dense nuclear and quark matter at zero temperature
and at beta equilibrium, and various models for the low density nuclear
matter, were established and studied in Ecker  _et al._ (2020); Jokela  _et
al._ (2021a). They have been successfully applied in simulations of binary
neutron star mergers Ecker  _et al._ (2020) and the study of rapidly rotating
isolated stars Demircik  _et al._ (2021). A related work Chesler  _et al._
(2019) combined the HS(DD2) model, as well as two other general purpose models
for the EoS, directly with the quark matter EoS from V-QCD, producing models
for the EoS with dependence on temperature and charge fraction outside beta
equilibrium. In the current article, we extend this work, among other things,
by including the description of nuclear matter from gauge/gravity duality and
its temperature dependence, which is our main focus, and also allows a
consistent description of the mixed nuclear-quark matter phase. A recent
article Fujimoto  _et al._ (2021) also considered models in the vdW class for
hot and dense nuclear matter. Their focus was matching the predictions of
chiral effective theory with the vdW models at low density, and use the vdW
models to extend the results to higher densities and temperatures, whereas in
the current work we use the vdW model to extrapolate the holographic cold
hybrid EoS (which is feasible for all densities) to finite temperatures.

The article is organized as follows. In Sec. II, we review the various
building blocks of our approach. In Sec. III, we show how the results for the
EoS from these building blocks can be combined into a single thermodynamically
consistent EoS. In Sec. IV we analyze the EoS, predictions for the critical
point, and predictions for nonrotating neutron stars. Our conclusions and
future directions are given in Sec. V. In three appendices we provide
technical details on the construction of our EoSs (App. A), compare the vdW
model to the nuclear matter part of V-QCD (App. B) and explain how we
determine our predictions for the QCD critical point (App. C).

Throughout this article we use Planck units where $c=\hbar=k_{B}=1$.

##  II Building blocks of the model

We start the discussion of our model by briefly reviewing the various building
blocks of the model.

###  II.1 Holographic V-QCD

A central ingredient in our setup is the use of the gauge/gravity duality. It
maps, in general, strongly coupled four dimensional field theory to classical
five dimensional gravity. Therefore challenging questions on the field theory
side can be solved by carrying out a simple classical analysis on the dual
gravity side. In this work we use the V-QCD model in regions which do not
admit a weakly coupled description in terms of quarks, gluons or hadrons.
V-QCD is an effective gauge/gravity model in the sense that it contains a
relatively large number of parameters that are tuned to match with QCD data
from experiments, lattice analysis, and perturbation theory. This is
particularly useful in the context of hot and dense QCD as there is plenty of
lattice data available at low density, and gauge/gravity models can be used to
extrapolate this data to higher densities at which first-principles methods
are not available DeWolfe  _et al._ (2011); Knaute  _et al._ (2018); Critelli
_et al._ (2017); Grefa  _et al._ (2021). Applying these ideas to V-QCD has
been shown to lead to feasible and well constrained EoSs for both dense quark
Jokela  _et al._ (2019) and nuclear matter Ishii  _et al._ (2019); Ecker  _et
al._ (2020); Jokela  _et al._ (2021a). This means the model is also able to
describe the nuclear to quark matter phase transition in a single framework.

Let us then explain briefly how V-QCD is constructed. We discuss here only the
main features of the model and refer the reader to Järvinen (2021) for a
complete review with precise definitions. The model contains both a gluon
sector, given by the improved holographic QCD model (five-dimensional dilaton
gravity) Gursoy and Kiritsis (2008); Gursoy  _et al._ (2008a), and a flavor
sector arising through a pair of space filling flavor branes Bigazzi  _et al._
(2005); Casero  _et al._ (2007). The flavors are dynamical: full backreaction
of the branes to the geometry is included formally by working in the Veneziano
limit Veneziano (1979) where both the number of colors $N_{c}$ and flavors
$N_{f}$ is large but their ratio is $\mathcal{O}(1)$ Järvinen and Kiritsis
(2012).

A basic feature of the gauge/gravity duality is that various phases in the
field theory map to different geometries in the five dimensional gravity
theory. In case of V-QCD, there are two possible geometries. The first is a
horizonless geometry ending at a “good” kind of IR singularity Gubser (2000)
and the second is “planar” black hole solution. The black holes can be
charged, which is interpreted as a nonzero baryon number arising from
deconfined quark matter Alho  _et al._ (2013, 2014). The gravitational
solution includes a scalar condensate in the confined phase, which is dual to
the chiral condensate of QCD therefore implementing chiral symmetry breaking
Bigazzi  _et al._ (2005); Casero  _et al._ (2007); Järvinen and Kiritsis
(2012). These geometries are therefore dual to a chirally broken confined
hadron gas phase and a chirally symmetric deconfined quark-gluon plasma phase,
respectively. For the plain confined geometry, the thermodynamics is trivial
in the sense that for example the pressure is independent of temperature and
chemical potential; we will use other methods in this region as we will
discuss below. For the deconfined phase, the temperature and entropy density
are calculated through black hole thermodynamics Gursoy  _et al._ (2008b,
2009); Alho  _et al._ (2013, 2014). Apart from hadron gas and quark matter, we
consider nuclear matter by employing the approach of Ishii  _et al._ (2019).
The nuclear matter phase is obtained in an approximation that is based on a
homogeneous five-dimensional bulk field in the confined horizonless geometry.
The five-dimensional action of the model is then obtained schematically as a
sum of three terms:

$\displaystyle
S_{\mathrm{V-QCD}}=S_{\mathrm{g}}+S_{\mathrm{f}}+S_{\mathrm{nm}}\,,$ (1)

where the first one ($S_{\mathrm{g}}$) is the action for improved holographic
QCD, i.e., the gluon sector of the model, the second one ($S_{\mathrm{f}}$) is
the flavor brane action important in the quark gluon plasma and quark matter
phases, and ($S_{\mathrm{nm}}$) is the action for homogeneous nuclear matter
derived in Ishii  _et al._ (2019).

We do not present the details of the actions here, but as we pointed out
above, they contain a relatively high number of parameters that need to be
tuned to match the model with QCD data. First, the model must agree with known
features of QCD such as asymptotic freedom, confinement, linear glueball and
meson trajectories, and chiral symmetry breaking. Second, those parameters
that are left free after considering such qualitative constraints, are fitted
to lattice data for thermodynamics of QCD Alho  _et al._ (2015); Jokela  _et
al._ (2019). Specifically, we use lattice data for the equation of state of
large-$N_{c}$ pure Yang-Mills Panero (2009) and data for $N_{c}=3$ QCD with
$N_{f}=2+1$ flavors at physical quark masses Borsanyi  _et al._ (2014, 2012)
at small baryon number density. Interestingly, the fit is stiff in the sense
that the dependence of the various model parameters is mild, but despite this
a fit of high quality is possible with all the parameter values remaining in a
natural range. At the moment, the fit uses flavor independent quarks with zero
mass. See the review Järvinen (2021) for a detailed discussion of this fit and
on the comparison of the model with QCD in general.

After determining the model parameters, V-QCD has a natural phase diagram that
includes both nuclear and quark matter Ishii  _et al._ (2019); Järvinen
(2021). The EoS for quark matter (QM) agrees with QCD lattice data at zero
density and finite temperature, with perturbative QCD for high values of the
baryon number chemical potential and/or temperature, and with constraints for
the EoS of cold quark matter EoS at low density Jokela  _et al._ (2019);
Chesler  _et al._ (2019). Notice however that we have not included electrons
or photons, and as the quark flavors are identical, there is no dependence on
charge fraction. For the full model EoS in the quark matter phase, we also
need to model these features, i.e., the electron pressure and the dependence
on charge fraction. We will discuss this in Sec. II.

In the nuclear matter phase, the model is feasible at zero temperature and can
be used to construct phenomenologically viable EoSs for cold QCD matter as was
done in Ecker  _et al._ (2020); Jokela  _et al._ (2021a, b). However, there is
an issue with the extension to finite temperature: the thermodynamics is
temperature independent in the confined phases, as we pointed out above. This
is a rather generic feature of gauge/gravity duality and arises due to taking
the limit of large $N_{c}$: the pressure of confined color singlet hadron
states, which are the constituents of the confined matter, is suppressed by
$1/N_{c}^{2}$ with respect to the pressure of the deconfined matter 222In
principle working in the Veneziano limit, where $N_{f}\sim N_{c}$ should
enhance the pressure of the confined matter. However it turns out that such
contributions are not captured only by backreacting the flavor action to the
geometry, but one should include nontrivial string loop diagrams Alho  _et
al._ (2015).. Due to this issue, the thermodynamics in the confined “hadron
gas” phase without nuclear matter is trivial, i.e., that of empty space. But
also in the nuclear matter phase, the temperature dependence is absent. While
this result may be a rather good zeroth-order approximation, it prevents us
from building a fully realistic EoS directly based on V-QCD. In order to cure
this issue, we will use a vdW model instead for the temperature dependence,
which we shall discuss next.

###  II.2 van der Waals Hadron Gas model (vdW-HG)

The vdW model consists of a bosonic and a fermionic sector. The ideal gas
pressure of the whole system is:

$\displaystyle p_{\mathrm{id}}(T,\\{\mu_{k}\\})$ $\displaystyle=$
$\displaystyle\sum_{i}p_{\mathrm{FD}}^{(i)}(T,\mu_{i},m_{i})$ (2)
$\displaystyle+\sum_{j}p_{\mathrm{BE}}^{(j)}(T,m_{j})+p_{\gamma}(T).\qquad$

The first term, $p_{\mathrm{FD}}^{(i)}$ is the fermionic contribution that
describes nucleons, antinucleons, electrons and positrons denoted by the index
$i\in\\{n,\bar{n},p,\bar{p},e,\bar{e}\\}$. The bosonic sector is given by the
last two terms, where $p_{\mathrm{BE}}^{(j)}$ and $p_{\gamma}$ are the
contributions of mesons and photons, respectively. We include all mesons with
$m_{j}\leq 1$ GeV from the particle data group listings Zyla  _et al._ (2020).
The photon contribution $p_{\gamma}(T)=\pi^{2}T^{4}/45$ is simply that of a
black-body photon gas. The pressure of a relativistic ideal Fermi (Bose) gas
is given by

$\displaystyle
p_{\mathrm{FD/BE}}^{(k)}(T,\mu_{k},m_{k})=\frac{g_{k}}{6\pi^{2}}\int_{0}^{\infty}\frac{p^{4}}{E_{k}}\frac{dp}{e^{(E_{k}-\mu_{k})/T}\pm
1},$ (3)

wherein the index $k$ denotes the fermion (boson) species with relativistic
dispersion relation $E_{k}=\sqrt{p^{2}+m_{k}^{2}}$, $+$ $(-)$ is for fermions
(bosons) and $g_{k}$ is the spin degeneracy factor. We also note that
different fermion and antifermion species have different chemical potentials
while $\mu_{k}=0$ for the bosonic sector.

It is well known that the ideal hadron gas (IHG) picture has vital
shortcomings. Perhaps most importantly, it fails in describing the ground
state of NM, viz., saturation of NM at $T=0$ and $n=n_{s}$. The reason for
this shortcoming is the fact that the IHG does not capture repulsive and
attractive van der Waals (vdW) interactions of nucleons. Repulsive vdW
interactions are often implemented via excluded volume (EV) corrections
describing a hard-core repulsion of nucleons. The thermodynamically consistent
formulation of an IHG with EV correction was developed in Rischke  _et al._
(1991) (see also Hagedorn and Rafelski (1980)). A more natural picture of vdW
interactions that realizes both the short range repulsive and the intermediate
range attractive interactions was formulated in Vovchenko  _et al._ (2017a, b)
(see Vovchenko (2020) for a more complete list of references). The parameters
of the vdW interactions are typically fixed by requiring consistency with the
saturation density $n_{s}$ and binding energy per nucleon of the ground state
$\epsilon_{B}/n_{b}=-16$MeV Vovchenko  _et al._ (2015); Redlich and Zalewski
(2016).

In our construction, we follow a different strategy and use holography as
guidance to model vdW interactions. In the rest of this subsection, we present
our implementation of repulsive vdW interactions via EV corrections by
employing the formulation in Rischke  _et al._ (1991). Attractive vdW
interactions are incorporated by a direct matching to the cold V-QCD hybrid
EoSs at beta-equilibrium described in section III.2.

The excluded volume corrected pressure is defined implicitly by introducing a
shifted chemical potential $\tilde{\mu}_{i}$:

$\displaystyle
p_{\mathrm{ex}}(T,\\{\mu_{k}\\})=p_{\mathrm{id}}(T,\\{\tilde{\mu}_{k}\\})\,;$
(4)
$\displaystyle\tilde{\mu}_{i}=\mu_{i}-v_{i}p_{\mathrm{ex}}(T,\\{\mu_{k}\\})\,.$
(5)

We choose $v_{p}=v_{\bar{p}}=v_{n}=v_{\bar{n}}=v_{0}$, $v_{e}=v_{\bar{e}}=0$
and set $v_{0}=0.56\,\text{fm}^{3}$. The choice of $v_{0}$ is motivated by
comparison with V-QCD EoSs at $T=0$. More details about our choice and a
comparison of different values for $v_{0}$ are given in Appendix B.

The number density is found by differentiating (4):

$\displaystyle n_{\mathrm{ex}}^{(i)}(T,\\{\mu_{k}\\})$ $\displaystyle=$
$\displaystyle\frac{\partial
p_{\mathrm{ex}}(T,\\{\mu_{k}\\})}{\partial\mu_{i}}$ (6) $\displaystyle=$
$\displaystyle\frac{n_{\mathrm{id}}^{(i)}(T,\\{\tilde{\mu}_{k}\\})}{1+v_{0}\sum_{l}n_{\mathrm{id}}^{(l)}(T,\\{\tilde{\mu}_{k}\\})},\;\;$

where $n_{\mathrm{id}}^{(i)}=\partial
p_{\mathrm{id}}(T,\\{\tilde{\mu}_{k}\\})/\partial\mu_{i}$ and $i$ denotes the
fermion species, while the index $l$ runs over only nucleons and antinucleons.
At this point, we fix the chemical potentials of antifermions as
$\mu_{\bar{i}}=-\mu_{i}$. The total number density for each fermion species is
then given by the difference between the corresponding particle and
antiparticle number densities:

$\tilde{n}_{\mathrm{ex}}^{(i)}(T,\\{\tilde{\mu}_{k}\\})=n_{\mathrm{ex}}^{(i)}(T,\\{\tilde{\mu}_{k}\\})-n_{\mathrm{ex}}^{(\bar{i})}(T,\\{\tilde{\mu}_{k}\\})\,,$
(7)

where $i\in\\{p,n,e\\}$. Requiring charge neutrality in addition,

$\tilde{n}_{\mathrm{ex}}^{(e)}(T,\\{\mu_{k}\\})=\tilde{n}_{\mathrm{ex}}^{(p)}(T,\\{\tilde{\mu}_{k}\\})\,,$
(8)

leaves only two free chemical potentials which we choose to be $\mu_{p}$ and
$\mu_{n}$. By using (6) and (7), the definitions for the baryon number density
$n_{b}$, the charge fractions $Y_{q}$ and charge neutrality condition are
obtained as

$\displaystyle n_{b}(T,\mu_{p},\mu_{n})$ $\displaystyle=$
$\displaystyle\tilde{n}_{\mathrm{ex}}^{(p)}(T,\mu_{p},\mu_{n})+\tilde{n}_{\mathrm{ex}}^{(n)}(T,\mu_{p},\mu_{n})\,,$
(9) $\displaystyle Y_{q}(T,\mu_{p},\mu_{n})$ $\displaystyle=$
$\displaystyle\tilde{n}_{\mathrm{ex}}^{(p)}(T,\mu_{p},\mu_{n})/n_{b}(T,\mu_{p},\mu_{n})\,.$
(10)

For future convenience, we conclude this section by giving the definition of
the free energy which is the natural thermodynamic potential in the canonical
ensemble:

$f_{\mathrm{ex}}(T,n_{b},Y_{q})=\sum_{i}n_{\mathrm{ex}}^{(i)}\mu_{i}(T,n_{b},Y_{q})-p_{\mathrm{ex}}(T,n_{b},Y_{q}),$
(11)

where the thermodynamics is expressed in terms of variables $n_{b}$ and
$Y_{q}$ instead of $\mu_{p}$ and $\mu_{n}$.

###  II.3 HS(DD2) model

HS(DD2) is a commonly used EoS that was originally developed to simulate core-
collapse of supernovae Hempel and Schaffner-Bielich (2010); Typel  _et al._
(2010). It provides reliable modelling of NM below and around the saturation
density, with a consistent description for the transition from nonuniform to
uniform nuclear matter.

HS(DD2) consists of two different sectors. The first sector models light and
heavy nuclear clusters in nuclear statistical equilibrium with EV correction.
The second sector describes unbound nucleons in relativistic mean-field (RMF)
theory. HS(DD2) includes all possible light nuclei (e.g. deuterium, tritium,
etc.) in addition to $\alpha$-particles and heavy nuclei up to mass number
$A\sim 330$.

The RMF approach is used to model interactions of unbound nucleons with the
exchange of $\sigma$, $\omega$ and $\rho$ mesons. At low densities, the
interactions become negligible and the system reduces to an ideal Fermi-Dirac
gas of nucleons. Photons are added separately as a free Bose gas and the
Wigner-Seitz approximation is employed for the Coulomb interaction between
electrons and nuclei. EV effects are implemented in a way to treat unbound
nucleons and nuclei in a different manner: while the volume of all baryons is
excluded for nuclei, unbound nucleons only feel the volume of nuclei since the
interaction among them is already modelled by the RMF model.

The parameters of the RMF model are the masses of the nucleons and the mesons,
and the coupling constants. In HS(DD2), the so called TMA parameter set is
used Toki  _et al._ (1995). For the masses of nuclei, the experimental data
Audi  _et al._ (2003) is used when it is available, otherwise the data is
taken from nuclear calculations Geng  _et al._ (2005).

By construction, HS(DD2) is a thermodynamically consistent model for both
nonuniform matter of light and heavy clusters within the gas of unbound
nucleons at low density, uniform matter described by the RMF model above
densities higher than $n_{s}$, and the transition between them. The resulting
EoS are provided in three dimensional ($T$,$n_{b}$,$Y_{q}$) tabular form
covering a wide range. The EoS tables are publicly available on the CompOSE
database Typel  _et al._ (2015); com (a). Different variants of the model were
constructed by using different sets of RMF parameters Hempel  _et al._ (2012);
Steiner  _et al._ (2013).

##  III Combining the building blocks

In this section we discuss in detail how the various building blocks of the
model are combined into a unified model (see Fig. 1).

###  III.1 Construction of the cold EoS

We start from the construction of the nuclear matter EoS at zero temperature,
which mostly follows Refs. Ecker  _et al._ (2020) and Jokela  _et al._
(2021a). That is, we first construct cold hybrid EoSs by combining the
predictions of V-QCD with nuclear theory models at beta equilibrium.

Our approach for nuclear matter in V-QCD, which is based on a homogeneous bulk
field, is natural and realistic for densities well above $n_{s}$. At such high
densities the average distance between neighboring nucleons are comparable or
smaller than their diameters, so that their wave functions overlap, and
approximating the system as homogeneous matter is expected to work. However,
this is not case at densities below $n_{s}$, wherein homogeneous
approximations will break down. This is not a problem, since in this region
traditional NM models have proven to be reliable and feasible. In this sense
traditional NM models and homogeneous NM in V-QCD complement each other. In
Ecker  _et al._ (2020); Jokela  _et al._ (2021a), this idea was implemented to
construct hybrid EoSs: low density EoSs from various traditional nuclear
theory models were combined with the high density V-QCD NM EoSs, matching them
continuously at an intermediate density around $1.5n_{s}$ to $2n_{s}$.

A potential weakness of the construction of Ecker  _et al._ (2020); Jokela
_et al._ (2021a) is that the matching point between the low and high density
models introduces, in effect, a second order phase transition, where for
example the speed of sound is discontinuous. While such discontinuities are
ubiquitous in constructions of EoSs in the literature and appear in particular
in the commonly used polytropic EoSs, we find that for the extension to finite
temperature, it is better if the cold EoS has continuous speed of sound. In
order to smooth it out, we consider an improved matching setup with two
separate matching densities in the same region of densities slightly above the
saturation density. The EoS in the intermediate region between the two
matching densities is chosen to have a speed of sound which is linear in the
baryon number density, whereas the low and high density regions are treated as
before. The slope of the speed of sound in the intermediate region is
determined by requiring continuity of the speed of sound, i.e., third order
phase transitions at both the transition densities. See Appendix A.1 for
details. The resulting EoS (free energy as a function of baryon number
density) is denoted below as $f_{\mathrm{cold}}(n_{b})$.

In this work, we will use three variants of hybrid EoSs for nuclear matter at
zero temperature. We select the Akmal-Pandharipande-Ravenhall (APR) model
Akmal  _et al._ (1998b) for the low density regime and three variants of V-QCD
at high density: these variants are defined by the data fits 5b, 7a, and 8b of
Jokela  _et al._ (2019, 2021a). These choices represent the leftover parameter
dependence of V-QCD after the comparison with lattice data and taking into
account observational constraints. If we used the approach of Ecker  _et al._
(2020); Jokela  _et al._ (2021a) with a single transition density $1.6n_{s}$
the hybrids with these three choices would be exactly the soft, intermediate
and stiff variants of the V-QCD(APR) EoS published in the CompOSE database com
(b, c, d). Here the stiffness refers to a property of dense nuclear matter:
the stiff EoS reaches a noticeably higher speed of sound than the soft one. In
the improved approach of this article, we choose the two transition densities
to be $1.4n_{s}$ and $1.8n_{s}$. That is, the choices of the cold nuclear
matter EoS in this article are practically the same as the published
V-QCD(APR) variants, but the kink (say, in pressure as a function of density)
at the transition density has been smoothed out.

###  III.2 Matching the vdW model with holography

We then discuss how the EoS for cold nuclear matter is extrapolated to finite
temperatures, and also outside beta equilibrium. As pointed out above, we use
the vdW model for the temperature dependence in the dense nuclear matter
phase. This approach is motivated by the fact that, for an appropriate choice
of the excluded volume $v_{0}$, the EoSs of the vdW model and the nuclear
matter of the V-QCD model are relatively close (see Appendix B). Therefore
adding only a small, mostly attractive potential is required to match the vdW
model with the cold hybrid EoS. Moreover, we add the dependence on the charge
fraction outside beta equilibrium by using the HS(DD2) model. In principle,
one could use the vdW model also for the charge fraction. The vdW model is
however too simple to satisfy the experimental constraints for nuclear matter
below and around saturation density. In particular, the symmetry energy is too
low – this is known to happen for a free gas, and excluded volume effects
alone are not sufficient to improve the result Baldo and Burgio (2016). We
could in principle use the vdW model at higher densities, were experimental
constraints do not apply, but for simplicity we adopt the $Y_{q}$ dependence
from HS(DD2) everywhere in the nuclear matter regime.

The above adjustments are taken into account by redefining the free energy of
the dense NM phase as

$\displaystyle
f_{\mathrm{vdW}}(T,n_{b},Y_{q})=f_{\mathrm{ex}}(T,n_{b},Y_{q})+\Delta
f(n_{b},Y_{q}),$ (12)

where $\Delta f$ models the mean contribution from an attractive potential. It
is simpler to specify the free energy difference $\Delta f$ directly rather
than start from the definition of the potential. Because the contribution from
the potential separates (see, e.g., Vovchenko (2020)), the two ways of
formulating this contribution are practically equivalent. That is, we take
333We use $T_{\mathrm{min}}=0.1$ MeV instead of $T=0$ when evaluating these
definitions numerically because this value is the lowest available in the data
grid of the HS(DD2) EoS.

$\displaystyle\Delta f(n_{b},Y_{q})$ $\displaystyle=$ $\displaystyle
f_{\mathrm{cold}}(n_{b})$ (13) $\displaystyle-
f_{\mathrm{ex}}(T=0,n_{b},Y_{q})$
$\displaystyle+f_{\mathrm{HS(DD2)}}(T=0,n_{b},Y_{q})$ $\displaystyle-
f_{\mathrm{HS(DD2)}}(T=0,n_{b},Y_{q}^{\mathrm{eq}}(n_{b})),$

where $f_{\mathrm{cold}}$ is the free energy of the one dimensional cold
hybrid EoS constructed as discussed above, $f_{\mathrm{HS(DD2)}}$ is the free
energy of the HS(DD2) model, and $Y_{q}^{\mathrm{eq}}$ is the value at beta-
equilibrium for the HS(DD2) EoS.

The first two lines on the right hand side of (13) adjust the dependence of
the EoS on $n_{b}$ such that it matches with the cold hybrid EoS at zero
temperature. The last two lines in (13) adjust the dependence on $Y_{q}$ such
that it agrees with that of the HS(DD2) EoS at low temperatures, without
changing the EoS at beta-equilibrium.

###  III.3 Transition between the vdW model and HS(DD2)

At low densities, such as relevant for example in the crust of the neutron
stars, the construction based to the holographic model and vdW gas is too
simple to be reliable, so we use instead directly the HS(DD2) EoS. We
implement this by switching from the vdW model of (12), which already borrows
the $Y_{q}$ dependence from HS(DD2), smoothly to the exact HS(DD2) EoS at a
well chosen transition density. However, before this is possible it is
necessary to adjust the HS(DD2) EoS by adding the contribution from the mesons
of QCD as indicated in Fig. 1. This contribution is important only in the
region of low density and high temperatures (i.e., close to the transition
temperature in QCD), which is far from the regime relevant for neutron stars
and core-collapse supernovae. We however add this contribution since it
affects the study of the critical point, which we carry out below. We write
the “improved” HS(DD2) EoS as

$\displaystyle\hat{f}_{\mathrm{HS(DD2)}}(T,n_{b},Y_{q})$ $\displaystyle=$
$\displaystyle f_{\mathrm{HS(DD2)}}(T,n_{b},Y_{q})$ (14)
$\displaystyle+\sum_{j}p_{\mathrm{BE}}^{(j)}(T,m_{j})$

where the sum goes over all mesons from the particle data group Zyla  _et al._
(2020) with masses below 1 GeV.

After this modification, we define the final nuclear matter EoS as

$\displaystyle f_{\mathrm{NM}}(T,n_{b},Y_{q})$ $\displaystyle=$
$\displaystyle[1-w(n_{b})]\hat{f}_{\mathrm{HS(DD2)}}(T,n_{b},Y_{q})$ (15)
$\displaystyle+w(n_{b})f_{\mathrm{vdW}}(T,n_{b},Y_{q})$

where the weight function is

$w(n_{b})=\frac{1}{2}\left[1+\tanh\left(\frac{\log(n_{b}/n_{0})}{1.75}\right)\right]=\frac{\left(n_{b}/n_{0}\right)^{8/7}}{1+\left(n_{b}/n_{0}\right)^{8/7}}$
(16)

with $n_{0}\approx 0.0694n_{s}$. The numerical coefficients were chosen such
that the transition from HS(DD2) to vdW is smooth for all temperatures and
charge fractions.

###  III.4 The mixed phase and the critical point

The final step in our construction is to combine the NM and QM components into
a single EoS. In order to do this, we first need to adjust the V-QCD QM
result: the EoSs constructed in Jokela  _et al._ (2019) does neither include
dependence on the charge fraction nor electron pressure. It would be possible
to compute the charge fraction dependence from the model directly, but this
would require a significant extension of the model, which is beyond the scope
of this article in which the focus is on temperature dependence. Therefore we
resort to approximations.

There are two simple approximation schemes: First is to assume that the free
energy of strongly interacting matter only depends on the total baryon number,
given as the sum over the quark number densities as
$n_{b}=(n_{u}+n_{d}+n_{s})/3$. In this case the free energy of QM arises as
the sum over the electromagnetic contribution and the V-QCD pressure:

$\displaystyle f_{\mathrm{QM}}(T,n_{b},Y_{q})$ $\displaystyle=$ $\displaystyle
f_{e\bar{e}\gamma}(T,Y_{q}n_{b})$ (17)
$\displaystyle+f_{\mathrm{V-QCD}}(T,n_{b})\ ,$

where the electron density is $n_{e}=Y_{q}n_{b}$ by charge neutrality. The
electromagnetic term $f_{e\bar{e}\gamma}$ is estimated as the sum of the ideal
gas free energies of electrons, positrons, and photons at the given electron
density.

The second scheme (which was used in Chesler  _et al._ (2019)) assumes that
the free energy arises as a direct sum of the free energies of different quark
flavors with equal amount of down and strange quarks. This gives

$\displaystyle f_{\mathrm{QM}}(T,n_{b},Y_{q})$ $\displaystyle=$ $\displaystyle
f_{e\bar{e}\gamma}(T,Y_{q}n_{b})$ (18)
$\displaystyle+\frac{1}{3}f_{\mathrm{V-QCD}}\left(T,(1+Y_{q})n_{b}\right)$
$\displaystyle+\frac{2}{3}f_{\mathrm{V-QCD}}\left(T,(1-Y_{q}/2)n_{b}\right)$

where the second (third) line is the contribution from up (down) type quarks.
Notice that this approach includes a simple approximation for the symmetry
energy of quark matter, which assumes no interactions between the different
quark flavors, and therefore corresponds to the probe limit where the
backreaction of the flavors to the gluon dynamics is neglected. However, there
is no reason to expect that the effect of the backreaction is small. In this
article, we use the simplest approximation of (17).

The final EoS and phase diagram, including the mixed phase between NM and QM,
is then found by carrying out a Gibbs construction (see Appendix A.2 for
details). As it turns out, two different regimes can be clearly identified
from the result. At low temperatures, there is a very strong first order phase
transition, which becomes weaker with increasing temperature. At higher
temperatures there is a weak first order transition. We interpret this weak
transition as the signal of crossover: continuity over the phases is not
possible because this would require a precise match between the EoSs of the
meson gas (in the NM phase) with the V-QCD QM EoS. We have not tried to carry
out such matching here; this is left for future work. The transition between
the two regimes is therefore interpreted as the critical endpoint of the
nuclear to quark matter transition line. We will illustrate the mixed phase
and the critical point in detail in Sec. IV and in Appendix C.

Figure 2: Contours of the speed of sound squared in beta-equilibrium for soft
(left), intermediate (middle), stiff (right) EoSs. Solid black lines separate
the mixed phase from NM and QM phases; yellow stars mark the locations of
critical points whose numerical values are listed in Table 1.

##  IV Results

We then analyze the thermodynamic properties of the constructed EoSs. We have
carried out the steps outlined in Sec. III for all three versions of the cold
hybrid EoSs, leading likewise to three versions of the final equation of state
which depends on baryon number density, temperature, and charge fraction. The
naming of the models is inherited from the cold EoSs, so that we refer to the
three EoSs as “soft”, “intermediate” and “stiff”, according to the stiffness
(i.e., basically the values of the speed of sound) of the EoS in the region of
dense nuclear matter. The stiffness is directly related to the maximal radii
and masses of neutron stars described by the corresponding EoS.

Before going to the analysis of our EoSs, let us summarize how they agree with
known constraints from various sources. Since the V-QCD QM model was fitted to
lattice data, good agreement with lattice results at low density, including
the first nontrivial Taylor expansion coefficient in the chemical potential,
is guaranteed above the crossover temperature, $T\gtrsim 150$ MeV. Because the
meson gas contribution is added to HS(DD2) at lower temperatures, agreement
with QCD in this opposite region is obtained as well: it is known that even
with our simple approximation, good agreement with lattice data, including
higher order Taylor coefficients in chemical potential, is found Karsch  _et
al._ (2003a, b); Vovchenko  _et al._ (2017a). Moreover, the QM EoS agrees by
construction Alho  _et al._ (2013, 2014) with leading perturbative QCD results
both at high temperatures and high chemical potentials (see Jokela  _et al._
(2019); Chesler  _et al._ (2019) for explicit comparison); since V-QCD is a
strongly coupled model, more detailed agreement with higher order perturbative
results is not possible. The hybrid EoSs, which we use for the cold EoSs,
agree by construction with models at low density such as chiral effective
theory computations. Actually, the cold hybrid EoSs are in excellent agreement
(see Jokela  _et al._ (2021a, b)) at all densities with model independent
constructions of the EoSs such as polytropic interpolations between the known
low and high density limits Annala  _et al._ (2018a); Most  _et al._ (2018);
Annala  _et al._ (2020, 2021) (see also Komoltsev and Kurkela (2021)), and
consequently also agree with constraints from neutron star measurements and
from the GW170817 merger event Abbott  _et al._ (2017a, b); see, e.g.,
Rezzolla  _et al._ (2018). We will also discuss these constraints explicitly
below. Lastly, our model is in good agreement with finite temperature
calculations in chiral effective theory as we shall discuss in more detail
below.

In Fig. 2 we plot the adiabatic speed of sound squared for our three models in
beta equilibrium as function of baryon number density and temperature. Solid
black lines at intermediate densities represent phase boundaries between
baryonic, mixed and quark phase, while dashed lines correspond to the
artificial phase boundary of our construction at low densities and
temperatures close to the crossover between baryonic and quark phase in QCD
(see Appendix C). The speed of sound squared can be expressed in terms of
first derivatives as

$c_{s}^{2}=\frac{n^{2}\left[\left(\frac{s}{n_{b}}-\frac{\partial s}{\partial
n_{b}}\right)^{2}+\frac{\partial s}{\partial
T}\frac{\partial\hat{\mu}}{\partial n_{b}}\right]}{\frac{\partial s}{\partial
T}\left(\epsilon+p\right)},$ (19)

where we assumed charge neutrality, and all quantities are assumed to be
functions of $T$, $n_{b}$, and $Y_{q}$. The chemical potential is defined as
$\hat{\mu}=\mu_{b}+Y_{q}\mu_{le}=\partial f/\partial n_{b}$, where $\mu_{b}$
($\mu_{le}$) is the baryon number (electron lepton number) chemical potential.
In all three models $c_{s}^{2}$ exceeds the value of $c_{s,\rm{CFT}}^{2}=1/3$
in conformal field theory close to the onset of the mixed phase. The maximal
values for $c_{s}^{2}$ are $0.42$, $0.47$, and $0.59$ for soft, intermediate
and stiff variants, respectively. The contours in the mixed phase are also
determined by Eq. (19), but notice that this expression is not the physical
speed of sound in this phase.

Yellow stars in Fig. 2 mark the location ($T_{c}$,$n_{bc}$) of the critical
point in the respective models whose precise values are listed in Table 1. The
critical point has been analyzed in various models in the literature, and the
results for the location vary in a wide range depending on the model Stephanov
(2004); DeWolfe  _et al._ (2011); Ayala  _et al._ (2021); Aryal  _et al._
(2021). Recent results in a simpler holographic approach DeWolfe  _et al._
(2011), which extrapolates results for thermodynamics of QCD from lattice QCD
to higher values of baryon chemical potential by using a bottom-up setup, are
given by $\\{T_{c},\mu_{bc}\\}=\\{112,612\\}$ MeV Knaute  _et al._ (2018) and
$\\{T_{c},\mu_{bc}\\}=\\{89,724\\}$ MeV Critelli  _et al._ (2017); Grefa  _et
al._ (2021). These numbers are in the same ballpark with ours: we obtain in
average slightly higher critical temperatures and lower critical chemical
potentials. Notice also that the critical point in all three variants is
outside the regime probed by the second phase of the beam energy scan at RHIC
Yang (2017), but will be probed in future experiments at FAIR and at NICA.
Finally, our numbers lie close to the chemical freeze-out curve extracted from
heavy-ion experiments (see, e.g., Adamczyk  _et al._ (2017); An  _et al._
(2022)). Our numbers are actually slightly below the experimental data, but
consistent with the curve if the precision of the data and our approach are
taken into account.

Figure 3: Latent heat as function of $T$ are shown for the three variants of
EoSs. The location of critical endpoint is determined via the condition
$\Delta\epsilon(T_{c},n_{bc})=0$.

In Fig. 3 we plot the latent heat
$\Delta\epsilon=\epsilon_{\text{QM}}-\epsilon_{\text{NM}}$, i.e., the
difference between the energy density in the quark phase
$\epsilon_{\text{QM}}$ and the nuclear matter phase $\epsilon_{\text{NM}}$, as
function of the temperature. For the three models we have analysed the soft
(stiff) model leads to the smallest (largest) latent heat at small
temperatures. Curiously, at $T\approx 100$ MeV all three models lead to
approximately the same value of $\Delta\epsilon$. The vanishing of the latent
heat $\Delta\epsilon=0$ determines the location of the critical point. Notice
however that as we pointed out above, the phase transition in our model is
always of first order $\Delta\epsilon$ does not vanish exactly but becomes
small above certain temperatures, and the point where $\Delta\epsilon=0$ is
obtained via extrapolation (see Appendix C for details).

The temperature $T_{c}$ of the critical point is correlated with the stiffness
of the respective model: larger stiffness results in lower values of $T_{c}$.
However, there is no clear relation between the critical density and
stiffness, as the highest value of $n_{bc}$ is found for the intermediate
model. In this case we interpret the variation of $n_{bc}$ as a rough measure
of the precision for the value of the critical density.

Figure 4: Pressure a function of the baryon number density in beta-equilibrium
for different values of the temperatures. The lower (upper) bounds of the
coloured bands represent the soft (stiff) model, while central curves
correspond to the intermediate version.

In Fig. 4 we plot the pressure as function of the baryon number density at
different values of the temperature in beta equilibrium. The lower (upper)
bounds of the coloured bands represent the soft (stiff) model, while the
central curves correspond to the intermediate version. The results in the
nuclear matter phase for $n_{b}/n_{s}<1.4$ are the same for all the three
models as there is not input from holography in this region. The uppermost
curve ($T=150$ MeV) shows the EoS in the quark matter phase entirely described
by V-QCD, which is therefore slightly different for all three models even at
low $n_{b}$. The plateau at intermediate densities is a manifestation of the
strong first-order phase transition of the V-QCD model. The transition density
decreases with increasing temperature, and the transition becomes weaker, as
also can be seen from Fig. 3.

Figure 5:  Pressure for the intermediate model for different values of the
temperatures. Solid curves are beta equilibrium, while upper and lower bounds
of the hatched coloured bands represent maximum and minimum values of the
pressure at given temperature and density.

To illustrate impact of $Y_{q}$ outside beta-equilibrium we plot in Fig. 5 the
range in pressure that is covered by the intermediate model at various values
of the temperature. Solid curves show the pressure in beta equilibrium, while
upper and lower bounds of the hatched coloured bands are maximum and minimum
values, respectively at given density and temperature. Overall the pressure is
lowest close to and increases away from beta-equilibrium, except close to the
onset of the mixed phase, where the pressure can take smaller values also away
from beta-equilibrium. The variation of the pressure as function of $Y_{q}$ is
largest in the mixed phase.

Figure 6: Thermal index in beta-equilibrium. Notation for the bands as in Fig.
4.

A quantity that is useful to explore the thermal contributions to the energy
density and the pressure is the thermal index defined as

$\Gamma_{\rm{th}}(T,n_{b},Y_{q})=1+\frac{p(T,n_{b},Y_{q})-p(0,n_{b},Y_{q})}{\epsilon(T,n_{b},Y_{q})-\epsilon(0,n_{b},Y_{q})}\,.$
(20)

In Fig. 6 we show the thermal index on the beta equilibrium slice for various
values of the temperature. As in Fig. 4, the upper (lower) bound of the
colored bands represent the stiff (soft) model, while the central curves
represent the intermediate case. The dashed grey line is the thermal index of
a free Fermi gas $\Gamma^{(FD)}_{\rm{th}}=5/3$ which is independent of density
and temperature. Any deviations of $\Gamma_{\rm{th}}(T,n_{b},Y_{q})$ from
$\Gamma^{(FD)}_{\rm{th}}$ are due to thermal interaction effects. Results for
the thermal index were computed at low temperatures and up to densities
between $1n_{s}$ and $2n_{s}$ by using chiral effective theory (CET) in
Carbone and Schwenk (2019); Keller  _et al._ (2021). It is important to
compare our EoS to these results in particular because the temperature
dependence in our setup in this region is based on the vdW setup, which is not
guaranteed to be realistic enough to agree with the CET calculations.
Interestingly, our result show good overall agreement with the CET
predictions, i.e., deviations are mostly below 10% in the relevant region.
Below $\sim 0.5n_{s}$ the CET predicts values above 1.6 which are slightly
higher than our $T=20$ MeV curve. Also the CET analysis predicts mild decrease
of the thermal index with increasing density, whereas in our model the index
mildly increases with increasing density. In neutron star simulations a
constant thermal index $\Gamma_{\rm{th}}=1-2$ is often assumed to mimic such
finite temperature effects Bauswein  _et al._ (2010); De Pietri  _et al._
(2020); Xie  _et al._ (2020), where $\Gamma_{\rm{th}}\approx 1.7$ has been
argued Figura  _et al._ (2020) to best approximate the dynamical and
thermodynamical behavior of neutron star merger simulations with microscopic
prescription of finite-temperature effects. The thermal index in our
construction remains also well within these bounds at temperatures relevant in
such simulations, except in and close to the mixed phase, where
$\Gamma_{\rm{th}}$ can take values smaller smaller than one. As expected,
thermal interaction effects become more important at higher temperature, where
the deviations of the thermal index from $\Gamma^{(FD)}_{\rm{th}}$ are
largest.

Finally, in Fig. 7 we show the mass-radius relation of cold isolated non-
rotating neutron stars for the soft (red), intermediate (green) and stiff
(blue) model, together with the relevant observational constraints. Solid
parts of the curves represent purely baryonic stars, while dashed parts belong
to stars with quark matter cores, which in our model turn out to be unstable.
Circles mark the maximum mass of stable non-rotating stars
($M_{\mathrm{TOV}}$) of the respective model. For the soft and intermediate
model $M_{\mathrm{TOV}}$ is determined by the onset of the phase transition at
which the star becomes unstable to black hole collapse. For the stiff model
the maximum mass is already reached in the baryonic phase and the phase
transition is realised only in the unstable branch of the mass-radius
sequence.

Figure 7: Mass-radius relation of non-rotating stars.

The green band shows the result $M_{\mathrm{TOV}}=2.08\pm 0.07M_{\odot}$ from
direct mass measurements of the pulsar J0740+6620 Cromartie  _et al._ (2019);
Fonseca  _et al._ (2021), which sets a lower bound on the maximum mass of
nonrotating neutron stars. Pink ellipses are radius measurements from the
NICER experiment for the pulsars J0030+0451 Riley  _et al._ (2019); Miller
_et al._ (2019) and J0740+6620 Miller  _et al._ (2021); Riley  _et al._
(2021), while cyan area is from the measurement of the X-ray binary 4U
1702-429 Nättilä  _et al._ (2017). All three examples pass the observational
constraints shown in the plot in addition to the constraint on the tidal
deformability $\Lambda_{1.4}<580$ of a $M=1.4M_{\odot}$ star deduced from the
analysis of GW170817 by LIGO/Virgo (low-spin prior at 90% confidence level)
Abbott  _et al._ (2018) (see also Table 1 and the detailed analysis in Jokela
_et al._ (2021b)).

Table 1:  EoS and neutron star properties in beta equilibrium.  Model | $\frac{n_{\text{bc}}}{n_{\text{s}}}$ | $\frac{\mu_{bc}}{\text{MeV}}$ | $\frac{T_{\text{c}}}{\text{MeV}}$ | $\frac{M_{\mathrm{TOV}}}{M_{\odot}}$ | $\frac{R_{\mathrm{e},1.4}}{\mathrm{km}}$ | $\Lambda_{1.4}$
---|---|---|---|---|---|---
soft | $0.46$ | $485$ | $128$ | $2.02$ | $12.41$ | $483$
interm. | $0.62$ | $575$ | $118$ | $2.14$ | $12.50$ | $511$
stiff | $0.32$ | $565$ | $112$ | $2.34$ | $12.64$ | $560$
HS(DD2) | – | – | – | $2.45$ | $13.2$ | $686$

##  V Conclusion and Outlook

In this article, we presented a novel framework for the EoS of hot and dense
QCD which combines ingredients from various approaches in different regions of
the phase diagram including gauge/gravity duality, van der Waals model of
nuclear matter, statistical models, and relativistic mean field theory. The
aim was to establish an EoS which uses the best available modeling in each of
the regimes. An essential new input here was the holographic V-QCD model,
which we used to cover the region of intermediate densities where computing
theoretical predictions is particularly hard. We presented three versions of
the EoS which are in good agreement with QCD data and constraints from
measurements of neutron stars and neutron star mergers. Using these models we
derived, apart from the properties of the EoS itself, predictions for the
location of the critical end point of the nuclear to quark matter transition.

Some details in our approach merit further study. It is expected (see, e.g.,
Alford  _et al._ (2008)) that at low temperatures in the quark matter phase,
pairing of quarks takes place, leading to potentially complex phase diagram
with various paired, color superconducting phases. Recently, there has been a
lot of interest in the analysis of such phases in gauge/gravity duality Chen
_et al._ (2010); Basu  _et al._ (2011); Bitaghsir Fadafan  _et al._ (2018);
Faedo  _et al._ (2019); Ghoroku  _et al._ (2019); Henriksson  _et al._ (2019).
Future work will study the inclusion of such phases in the V-QCD model and
their effect on the EoS.

While the main focus in this article was the temperature dependence, we also
included the dependence on the charge fraction following the HS(DD2) model in
the nuclear matter phase, and by using a simple approximation assuming a
strongly interacting component of free energy only depending on the total
baryon number (and with free electron gas) in the quark matter phase. This
latter approach can be improved by including proper flavor dependence in the
holographic model, so that one can consider states with unequal amount of
different quark flavors. This extension of the model is the topic of ongoing
research.

The inclusion of flavor dependence will also help to generalize the model to
analyze transport properties of QCD matter. The strongly interacting
components of conductivities and viscosities of QM at high density were solved
in Hoyos  _et al._ (2020, 2021a). Including flavor dependence will help to
properly analyze the correlators of electric and weak currents, which are
necessary to estimate the effects of electron and neutrino transport, known to
be important, among other things, for neutron star cooling, in core-collapse
supernovae, and for the behavior of ejected matter in neutron star mergers
(see, e.g., Schmitt and Shternin (2018)).

Furthermore in this work the treatment of the cross-over at low densities and
intermediate temperature was rather simple minded: the EoSs actually have
(very weak) first order phase transition instead of a crossover. It will be
interesting to study if the EoS can be improved in this region by matching the
meson gas pressure with the high-temperature holographic quark-gluon plasma
pressure more carefully, e.g., by extending the methods of Alho  _et al._
(2015) to finite density.

We should stress that the framework presented here admits immediate natural
generalizations where one replaces some of the building blocks of the model by
other approaches. As for the strongly coupled gauge/gravity duality model
there is (to our knowledge) currently no alternative to V-QCD that would allow
to repeat the analysis as done in this article. This field is however evolving
rapidly Hoyos  _et al._ (2016); Annala  _et al._ (2018b); Bitaghsir Fadafan
_et al._ (2020); Mamani  _et al._ (2020); Bitaghsir Fadafan  _et al._ (2021);
Ghoroku  _et al._ (2021); Hoyos  _et al._ (2021b). For example, very recently
it was demonstrated that a setup with nuclear matter in the Witten-Sakai-
Sugimoto model may lead to realistic neutron stars Kovensky  _et al._ (2021).
For the temperature dependence of dense nuclear matter we used the vdW model
for two reasons: it is simple, and agrees reasonably well with the holographic
model at zero temperature. A simple model may be the best guess for the
temperature dependence at densities around and above the saturation density,
where none of the known approaches is reliable. However, also other choices
are possible. We have checked that using (for example) the temperature
dependence of the HS(DD2) model instead only leads to rather mildly modified
EoSs. It would be also interesting to study the Carnahan-Starling
generalization of the excluded volume effect in the vdW model, recently
considered in Fujimoto  _et al._ (2021), since it may improve the temperature
dependence of the EoS at low but nonzero temperatures near the saturation
density, where ab initio results are available. Another simple generalization
of our approach is to use some other general purpose EoSs for the low density
region than the HS(DD2) model.

Apart from varying the building blocks of the model, a future study could
further explore the parameter space of the model presented in this article.
This would mean, at least, including additional variants of the V-QCD model
Jokela  _et al._ (2019) an studying the effect of varying the excluded volume
parameter $v_{0}$ in the vdW model.

The EoSs constructed in this article will be published in the standard format
in the CompOSE database, and can therefore immediately to be used in state-of-
the-art simulations of neutron star mergers and core-collapse supernovae.
Indeed, there is growing interest in effects arising due to the temperature
dependence and the phase transition in merger simulations Most  _et al._
(2019); Bauswein  _et al._ (2019); Raithel  _et al._ (2021); Prakash  _et al._
(2021). Work on applying the EoSs of this article in such simulations is
already in progress.

###### Acknowledgements.

We thank U. Gürsoy, N. Jokela, E. Kiritsis, L. Rezzolla, W. van der Schee, and
A. Vuorinen for useful discussions and comments. T. D. and M. J. have been
supported by an appointment to the JRG Program at the APCTP through the
Science and Technology Promotion Fund and Lottery Fund of the Korean
Government. T. D. and M. J. have also been supported by the Korean Local
Governments – Gyeongsangbuk-do Province and Pohang City – and by the National
Research Foundation of Korea (NRF) funded by the Korean government (MSIT)
(grant number 2021R1A2C1010834). T. D. would like to thank the Department of
Physics of Boğaziçi University for their hospitality during his visit. C. E.
acknowledges 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.

##  Appendix A Technical details on the EoS

In this Appendix, we discuss minor technical details which are not important
for our conclusions, but need to be discussed for all the results to be fully
reproducible.

###  A.1 Details of the construction of the cold EoS

We first give more details on the construction of the cold hybrid EoS. Recall
that the EoS of cold nuclear matter, as explained in Sec. III.1, requires
matching between EoSs from a nuclear theory model (in this article, the APR
model) at low density with V-QCD nuclear matter at higher density. As we
explained in main text, in order for the cold EoS to be smoother, we choose
two matching densities, instead of the single density approach of Ecker  _et
al._ (2020); Jokela  _et al._ (2021a), and connect the speeds of sound from
the two approaches though a linear interpolation of the squared speed of sound
$c_{s}^{2}$ as a function of the baryon number density $n_{b}$, therefore
avoiding a discontinuity in the speed of sound.

To be precise, the matching is carried out as follows. For
$n_{b}<1.4n_{s}\equiv n_{\mathrm{tr}}^{(1)}$ we use the APR EoS as such. We
then take the speed of sound for $n_{\mathrm{tr}}^{(1)}<n_{b}<1.8n_{s}\equiv
n_{\mathrm{tr}}^{(2)}$ to be

$\left(c_{s}(n_{b})\right)^{2}=\left(c_{s}(n_{\mathrm{tr}}^{(1)})\right)^{2}+\kappa\,(n-n_{\mathrm{tr}}^{(1)})$
(21)

where the slope $\kappa$ is a free parameter. The rest of the thermodynamic
functions are then obtained by integration for
$n_{\mathrm{tr}}^{(1)}<n_{b}<n_{\mathrm{tr}}^{(2)}$ so that $\kappa$ is the
only free parameter. We then require that the pressure, the baryon number
chemical potential, and the speed of sound are continuous at
$n_{b}=n_{\mathrm{tr}}^{(2)}$. These three conditions then determine $\kappa$
as well as the two parameters of V-QCD nuclear matter $c_{b}$ and $b$ (these
parameters appear in the nuclear matter action $S_{nm}$ of (1), see Jokela
_et al._ (2021a) for their definitions). These latter two parameters were also
determined by matching in the simpler approach of Ecker  _et al._ (2020);
Jokela  _et al._ (2021a).

We also modify the cold hybrid EoS in the crust region, $n_{b}\ll n_{s}$,
before using it to construct the adjusted vdW EoS through (13), in order to
remove some noise from the vdW EoS. This noise would appear because the fine
details of the beta-equilibrium, zero temperature APR EoS which turns out to
be inconsistent with the $Y_{q}$-dependence of the HS(DD2) EoS at low density.
One might wonder why we need to do this, since we anyhow replace the vdW EoS
with the exact HS(DD2) through (15) at low densities. However, since we use a
smooth weight function instead of an abrupt cutoff, noise in the vdW EoS at
low densities would remain in the final matched EoS even if it would be
heavily suppressed. We therefore implement the correction by using the HS(DD2)
at beta equilibrium for the cold EoS $f_{\mathrm{cold}}(n_{b})$ when
$n_{b}<0.008n_{s}$. At this value the pressures of the APR and HS(DD2) EoSs
cross, and we further introduce a small shift in the baryon chemical potential
for $n_{b}<0.008n_{s}$ in order to make it continuous at $n_{b}=0.008n_{s}$.
That is, we in effect introduce a second order phase transition. We stress
however that these modifications only remove noise from heavily suppressed
terms at low density in our final results.

###  A.2 Details of the construction of the three dimensional EoS

We then discuss some technical details on the determination of the final three
dimensional EoS and its components.

We carried out some minor but nontrivial modifications in the three
dimensional nuclear matter EoS $f_{\mathrm{NM}}(T,n_{b},Y_{q})$. The
transition from non-uniform to uniform nuclear matter in HS(DD2) proceeds via
a first order phase transition at around $0.3n_{s}\lesssim n_{b}\lesssim
0.5n_{s}$, at small temperatures $T\leq 5$ MeV, and values of $Y_{q}$ (mostly)
far from beta equilibrium Hempel and Schaffner-Bielich (2010). Due to the
smooth matching of (15), the mixed phase from this transition causes the
resulting nuclear matter to be inconsistent in this narrow range of
parameters, so that a thermodynamically unstable region appears. In order to
fix this, we carried out a simple one-dimensional Maxwell construction for
each value of $T$ and $Y_{q}$ which removed the inconsistency.

Let us also specify how the EoS for the mixed phase between the NM and QM
phases was calculated. It arises from requiring full chemical equilibrium, so
that the coexisting phases have the same pressure as well as baryon number and
electron lepton number chemical potentials

$\displaystyle\mu_{le}$ $\displaystyle=$
$\displaystyle\frac{1}{n_{b}}\\!\left.\frac{\partial f}{\partial
Y_{q}}\right|_{n_{b},T}\ ,$ (22) $\displaystyle\quad\mu_{b}$ $\displaystyle=$
$\displaystyle\left.\frac{\partial f}{\partial
n_{b}}\right|_{n_{b}Y_{q},T}=\left.\frac{\partial f}{\partial
n_{b}}\right|_{Y_{q},T}-Y_{q}\mu_{le}\ .$ (23)

In practice, the mixed phase is found as follows. Demanding equilibrium
between the phases, we need to solve the following set of equations:

$\displaystyle p_{\mathrm{NM}}(T,n_{b}^{(1)},Y_{q}^{(1)})$ $\displaystyle=$
$\displaystyle p_{\mathrm{QM}}(T,n_{b}^{(2)},Y_{q}^{(2)})\ ,$ (24)
$\displaystyle\mu_{b}^{(\mathrm{NM})}(T,n_{b}^{(1)},Y_{q}^{(1)})$
$\displaystyle=$
$\displaystyle\mu_{b}^{(\mathrm{QM})}(T,n_{b}^{(2)},Y_{q}^{(2)})\ ,$ (25)
$\displaystyle\mu_{le}^{(\mathrm{NM})}(T,n_{b}^{(1)},Y_{q}^{(1)})$
$\displaystyle=$
$\displaystyle\mu_{le}^{(\mathrm{QM})}(T,n_{b}^{(2)},Y_{q}^{(2)})\ .$ (26)

We have three conditions and four variables $n_{b}^{(i)}$, $Y_{q}^{(i)}$ so
the solution will involve one parameter which we call $\gamma$. The solution
defines two curves on the $(n_{b},Y_{q})$-plane, parameterized in terms of
$\gamma$, and a mapping between the curves. The mixed phase is found between
the curves. The temperature is a “trivial” parameter in these equations and we
will not denote the dependence on it explicitly below. The construction can be
carried out independently for each value of the temperature.

The mixed phase is then a mixture of NM and QM matter in the equilibrium
defined by Eqs. (24)–(26). The thermodynamic functions are most easily written
in a parametric representation using the $\gamma$ parameter and the volume
fraction $\alpha$ of the NM phase. That is, we may write

$\displaystyle n_{b}(\alpha,\gamma)$ $\displaystyle=$ $\displaystyle\alpha
n_{b}^{(1)}(\gamma)+(1-\alpha)n_{b}^{(2)}(\gamma)$ (27) $\displaystyle
Y_{q}(\alpha,\gamma)n_{b}(\alpha,\gamma)$ $\displaystyle=$
$\displaystyle\alpha Y_{q}^{(1)}(\gamma)n_{b}^{(1)}(\gamma)$ (28)
$\displaystyle\\!\\!+(1-\alpha)Y_{q}^{(2)}(\gamma)n_{b}^{(2)}(\gamma)$
$\displaystyle f(\alpha,\gamma)$ $\displaystyle=$ $\displaystyle\alpha
f_{\mathrm{NM}}(n_{b}^{(1)}(\gamma),Y_{q}^{(1)}(\gamma))$ (29)
$\displaystyle\\!\\!+(1-\alpha)f_{\mathrm{QM}}(n_{b}^{(2)}(\gamma),Y_{q}^{(2)}(\gamma))\
.$

separately for each temperature slice. Note that curves of constant $\gamma$
are straight lines on the $(n_{b},Y_{q}n_{b})$-plane. By construction, $p$,
$\mu_{b}$, and $\mu_{le}$ take constant values on these lines.

Notice that in the main text we focused on defining the free energy, which is
the natural thermodynamic potential in the canonical ensemble with the
parameters $T$, $n_{b}$, and $Y_{q}$. Determining first numerically the free
energy, and computing the other thermodynamic functions by using it as input,
is indeed enough to determine all thermodynamics. However, this procedure
requires taking numerical derivatives, which tend to increase numerical noise.
When computing the final results for other observables (such as the entropy
and the chemical potentials) we have therefore first computed analytically the
consequences of the various matching formulas, and computed these quantities
directly from the corresponding quantities of the EoSs being matched, avoiding
the use of numerical derivatives as much as possible.

Figure 8: Comparison of cold V-QCD EoSs and EoSs of hadron gas with excluded
volume correction. Left: normalized pressure ($p/\mu_{b}^{4}$) in terms of
baryon chemical potential $\mu_{b}$. Right: pressure in terms of baryon number
density $n_{b}$ in units of $n_{s}$. In both figures, stiff, intermediate,
soft variants are denoted by blue, green, red curves respectively and hadron
gas with excluded volume of $v_{0}=0.56\mathrm{fm}^{-3}$,
$v_{0}=1\mathrm{fm}^{-3}$and $v_{0}=1.5\mathrm{fm}^{-3}$ are exhibited by
solid, dotted, dashed black curves. In addition, ideal hadron gas result is
shown with the gray solid curve in the figure on the right.

##  Appendix B Comparison of the vdW EoS to the V-QCD NM EoS

A central motivation for the use of the vdW EoS in the article is its
agreement with the predicted EoS of cold NM by V-QCD. In this Appendix, we
study this by comparing the three different versions of the cold hybrid
V-QCD(APR) EoSs to simple vdW EoSs, i.e., those only with electrons, protons,
and neutrons with a constant excluded volume correction for the nucleons.

The EoSs are compared in Fig. 8. In both figures, stiff, intermediate, soft
variants of V-QCD EoSs up to the onset of the phase transition are showed via
blue, green, red curves respectively and simple vdW EoSs with EV parameter
$v_{0}=0.56\mathrm{fm}^{-3}$ (our choice), $v_{0}=1\mathrm{fm}^{-3}$,
$v_{0}=1.5\mathrm{fm}^{-3}$ are denoted by solid black, dotted, dashed curves.
The figure on the left shows the baryon chemical potential dependence of the
dimensionless pressure $p/\mu_{b}^{4}$, and the figure on the right shows the
pressure as a function of the baryon number density. The gray curve in the
right-hand side plot shows the ideal hadron gas result. This curve is not
shown in the left plot since it mostly stays outside the plotted range. It is
transparent from the figures that the chemical potential dependence is more
sensitive to $v_{0}$ at small values of the parameter.

We notice that the EV corrected EoS is relatively close to the cold hybrid
V-QCD EoSs, in particular when the pressure is plotted as a function of the
chemical potential. This happens in part because the left hand plot of Fig. 8
“zooms” into the region of higher densities where agreement is better. Notice
that the range of values of $v_{0}$ is typical for EV corrections in nuclear
matter, see, e.g., Vovchenko  _et al._ (2017b). The best fit between the EV
corrected pressure and V-QCD is found around $v_{0}\approx
1.5\mathrm{fm}^{-3}$. We however choose a smaller value $v_{0}\approx
0.56\mathrm{fm}^{-3}$, because for this value, the potential term of (13) is
more natural: this value corresponds to a potential which is attractive
(repulsive) at long (short) distances, whereas for larger $v_{0}$ the
potential would be attractive at both short and long distances with an
intermediate repulsive range in between. This unnatural behavior is reflected
in the nonmonotonic dependence of the pressure difference between the EV
corrected and V-QCD pressures as a function of $n_{b}$ in the left plot of
Fig. 8.

##  Appendix C Determination of the critical point

As discussed in Sec.III.4, the EoSs have a strong phase transition at low
temperature that becomes weaker at higher temperature. Continuity over the
phases is not possible due to limitation of our construction. Hence, we
interpret the weakening in the transition as signal for a crossover at high
temperature. Using this interpretation, it is possible to calculate the latent
heat and obtain an estimate for the location of the critical end point. The
latent heat is

$\Delta\epsilon(T)=\epsilon_{\text{QM}}(T,n_{b}^{(2)}\\!,Y_{q}^{\mathrm{eq}}(n_{b}^{(2)}))-\epsilon_{\text{NM}}(T,n_{b}^{(1)}\\!,Y_{q}^{\mathrm{eq}}(n_{b}^{(1)})),$
(30)

where $Y_{q}^{eq}$ is the value at beta-equilibrium and $n_{b}^{(2)}$,
$n_{b}^{(1)}$ are respectively the values for the baryon number density at the
onset and the end of the phase transition for a given value of $T$. Due to the
absence of a crossover, the transition lines get close, but do not intersect.
However, the latent heat as function of temperature shows a clean trend. By
using it, we perform a $5^{th}$-order polynomial fit to the data in the range
of strong first order transition, i.e., $T\lesssim 115$ MeV. The fit functions
that we obtained are

$\displaystyle\Delta\epsilon_{\mathrm{stiff}}$ $\displaystyle=$ $\displaystyle
1.37\times 10^{3}-1.41\times 10T+1.53\times 10^{-1}T^{2}$
$\displaystyle-4.36\times 10^{-3}T^{3}+6.26\times 10^{-5}T^{4}$
$\displaystyle-3.10\times 10^{-7}T^{5},$
$\displaystyle\Delta\epsilon_{\mathrm{interm}}$ $\displaystyle=$
$\displaystyle 1.05\times 10^{3}-8.85T+5.01\times 10^{-2}T^{2}$
$\displaystyle-9.04\times 10^{-4}T^{3}+1.28\times 10^{-5}T^{4}$
$\displaystyle-7.53\times 10^{-8}T^{5},$
$\displaystyle\Delta\epsilon_{\mathrm{soft}}$ $\displaystyle=$ $\displaystyle
9.23\times 10^{2}-7.07T+1.40\times 10^{-2}T^{2}$ $\displaystyle-2.90\times
10^{-4}T^{3}-3.40\times 10^{-6}T^{4}$ $\displaystyle+1.53\times 10^{-9}T^{5}.$

The critical temperature is then found by extrapolating to the point where the
latent heat vanishes: $\Delta\epsilon(T_{c})=0$. To determine value of the
baryon number density at critical end point $n_{bc}$, we use the geometric
mean of the density values on the two transition lines at $T_{c}$. We estimate
the uncertainty of this procedure to be below $0.05n_{s}$ for the soft EoS and
$0.1n_{s}$ for the intermediate and stiff EoSs by studying the variation
caused by varying the choice of the mean. While this uncertainty is arguably
sizeable, it is smaller than the differences between the numbers for the
different EoSs in Table 1. Then for the phase diagram, crossover (dashed)
lines are also computed by calculating the geometric mean of $n_{b}$ on the
transition lines at the temperature values above $T_{c}$.

## References

  * Keller  _et al._ (2021) J. Keller, C. Wellenhofer, K. Hebeler, and A. Schwenk, Phys. Rev. C 103, 055806 (2021), arXiv:2011.05855 [nucl-th] .
  * Most  _et al._ (2019) E. R. Most, L. J. Papenfort, V. Dexheimer, M. Hanauske, S. Schramm, H. Stöcker, and L. Rezzolla, Phys. Rev. Lett. 122, 061101 (2019), arXiv:1807.03684 [astro-ph.HE] .
  * Perego  _et al._ (2019) A. Perego, S. Bernuzzi, and D. Radice, Eur. Phys. J. A 55, 124 (2019), arXiv:1903.07898 [gr-qc] .
  * Endrizzi  _et al._ (2020) A. Endrizzi, A. Perego, F. M. Fabbri, L. Branca, D. Radice, S. Bernuzzi, B. Giacomazzo, F. Pederiva, and A. Lovato, Eur. Phys. J. A 56, 15 (2020), arXiv:1908.04952 [astro-ph.HE] .
  * Hammond  _et al._ (2021) P. Hammond, I. Hawke, and N. Andersson, Phys. Rev. D 104, 103006 (2021), arXiv:2108.08649 [astro-ph.HE] .
  * Figura  _et al._ (2021) A. Figura, F. Li, J.-J. Lu, G. F. Burgio, Z.-H. Li, and H. J. Schulze, Phys. Rev. D 103, 083012 (2021), arXiv:2103.02365 [gr-qc] .
  * Raithel  _et al._ (2021) C. Raithel, V. Paschalidis, and F. Özel, Phys. Rev. D 104, 063016 (2021), arXiv:2104.07226 [astro-ph.HE] .
  * Adam  _et al._ (2021) J. Adam _et al._ (STAR), Phys. Rev. Lett. 126, 092301 (2021), arXiv:2001.02852 [nucl-ex] .
  * Ablyazimov  _et al._ (2017) T. Ablyazimov _et al._ (CBM), Eur. Phys. J. A 53, 60 (2017), arXiv:1607.01487 [nucl-ex] .
  * Durante  _et al._ (2019) M. Durante _et al._ , Phys. Scripta 94, 033001 (2019), arXiv:1903.05693 [nucl-th] .
  * Sissakian and Sorin (2009) A. N. Sissakian and A. S. Sorin (NICA), J. Phys. G 36, 064069 (2009).
  * Järvinen and Kiritsis (2012) M. Järvinen and E. Kiritsis, JHEP 03, 002 (2012), arXiv:1112.1261 [hep-ph] .
  * Ishii  _et al._ (2019) T. Ishii, M. Järvinen, and G. Nijs, JHEP 07, 003 (2019), arXiv:1903.06169 [hep-ph] .
  * Vovchenko (2020) V. Vovchenko, Int. J. Mod. Phys. E 29, 2040002 (2020), arXiv:2004.06331 [nucl-th] .
  * Hempel and Schaffner-Bielich (2010) M. Hempel and J. Schaffner-Bielich, Nucl. Phys. A 837, 210 (2010), arXiv:0911.4073 [nucl-th] .
  * Typel  _et al._ (2010) S. Typel, G. Ropke, T. Klahn, D. Blaschke, and H. H. Wolter, Phys. Rev. C 81, 015803 (2010), arXiv:0908.2344 [nucl-th] .
  * Akmal  _et al._ (1998a) A. Akmal, V. R. Pandharipande, and D. G. Ravenhall, Phys. Rev. C 58, 1804 (1998a), arXiv:nucl-th/9804027 .
  * Note (1) We use APR rather than HS(DD2) for cold nuclear matter near the saturation density because if we used HS(DD2) only, the resulting EoSs would be so stiff that they would conflict with the LIGO/Virgo measurement of the tidal deformability for GW170817.
  * Ecker  _et al._ (2020) C. Ecker, M. Järvinen, G. Nijs, and W. van der Schee, Phys. Rev. D 101, 103006 (2020), arXiv:1908.03213 [astro-ph.HE] .
  * Jokela  _et al._ (2021a) N. Jokela, M. Järvinen, G. Nijs, and J. Remes, Phys. Rev. D 103, 086004 (2021a), arXiv:2006.01141 [hep-ph] .
  * Typel  _et al._ (2015) S. Typel, M. Oertel, and T. Klähn, Phys. Part. Nucl. 46, 633 (2015), arXiv:1307.5715 [astro-ph.SR] .
  * com (a) https://compose.obspm.fr (a).
  * Demircik  _et al._ (2021) T. Demircik, C. Ecker, and M. Järvinen, Astrophys. J. Lett. 907, L37 (2021), arXiv:2009.10731 [astro-ph.HE] .
  * Chesler  _et al._ (2019) P. M. Chesler, N. Jokela, A. Loeb, and A. Vuorinen, Phys. Rev. D 100, 066027 (2019), arXiv:1906.08440 [astro-ph.HE] .
  * Fujimoto  _et al._ (2021) Y. Fujimoto, K. Fukushima, Y. Hidaka, A. Hiraguchi, and K. Iida, (2021), arXiv:2109.06799 [nucl-th] .
  * DeWolfe  _et al._ (2011) O. DeWolfe, S. S. Gubser, and C. Rosen, Phys. Rev. D 83, 086005 (2011), arXiv:1012.1864 [hep-th] .
  * Knaute  _et al._ (2018) J. Knaute, R. Yaresko, and B. Kämpfer, Phys. Lett. B 778, 419 (2018), arXiv:1702.06731 [hep-ph] .
  * Critelli  _et al._ (2017) R. Critelli, J. Noronha, J. Noronha-Hostler, I. Portillo, C. Ratti, and R. Rougemont, Phys. Rev. D 96, 096026 (2017), arXiv:1706.00455 [nucl-th] .
  * Grefa  _et al._ (2021) J. Grefa, J. Noronha, J. Noronha-Hostler, I. Portillo, C. Ratti, and R. Rougemont, Phys. Rev. D 104, 034002 (2021), arXiv:2102.12042 [nucl-th] .
  * Jokela  _et al._ (2019) N. Jokela, M. Järvinen, and J. Remes, JHEP 03, 041 (2019), arXiv:1809.07770 [hep-ph] .
  * Järvinen (2021) M. Järvinen, (2021), arXiv:2110.08281 [hep-ph] .
  * Gursoy and Kiritsis (2008) U. Gursoy and E. Kiritsis, JHEP 02, 032 (2008), arXiv:0707.1324 [hep-th] .
  * Gursoy  _et al._ (2008a) U. Gursoy, E. Kiritsis, and F. Nitti, JHEP 02, 019 (2008a), arXiv:0707.1349 [hep-th] .
  * Bigazzi  _et al._ (2005) F. Bigazzi, R. Casero, A. L. Cotrone, E. Kiritsis, and A. Paredes, JHEP 10, 012 (2005), arXiv:hep-th/0505140 .
  * Casero  _et al._ (2007) R. Casero, E. Kiritsis, and A. Paredes, Nucl. Phys. B 787, 98 (2007), arXiv:hep-th/0702155 .
  * Veneziano (1979) G. Veneziano, Nucl. Phys. B 159, 213 (1979).
  * Gubser (2000) S. S. Gubser, Adv. Theor. Math. Phys. 4, 679 (2000), arXiv:hep-th/0002160 .
  * Alho  _et al._ (2013) T. Alho, M. Järvinen, K. Kajantie, E. Kiritsis, and K. Tuominen, JHEP 01, 093 (2013), arXiv:1210.4516 [hep-ph] .
  * Alho  _et al._ (2014) T. Alho, M. Järvinen, K. Kajantie, E. Kiritsis, C. Rosen, and K. Tuominen, JHEP 04, 124 (2014), [Erratum: JHEP 02, 033 (2015)], arXiv:1312.5199 [hep-ph] .
  * Gursoy  _et al._ (2008b) U. Gursoy, E. Kiritsis, L. Mazzanti, and F. Nitti, Phys. Rev. Lett. 101, 181601 (2008b), arXiv:0804.0899 [hep-th] .
  * Gursoy  _et al._ (2009) U. Gursoy, E. Kiritsis, L. Mazzanti, and F. Nitti, JHEP 05, 033 (2009), arXiv:0812.0792 [hep-th] .
  * Alho  _et al._ (2015) T. Alho, M. Järvinen, K. Kajantie, E. Kiritsis, and K. Tuominen, Phys. Rev. D 91, 055017 (2015), arXiv:1501.06379 [hep-ph] .
  * Panero (2009) M. Panero, Phys. Rev. Lett. 103, 232001 (2009), arXiv:0907.3719 [hep-lat] .
  * Borsanyi  _et al._ (2014) S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, and K. K. Szabo, Phys. Lett. B 730, 99 (2014), arXiv:1309.5258 [hep-lat] .
  * Borsanyi  _et al._ (2012) S. Borsanyi, Z. Fodor, S. D. Katz, S. Krieg, C. Ratti, and K. Szabo, JHEP 01, 138 (2012), arXiv:1112.4416 [hep-lat] .
  * Jokela  _et al._ (2021b) N. Jokela, M. Järvinen, and J. Remes, (2021b), arXiv:2111.12101 [hep-ph] .
  * Note (2) In principle working in the Veneziano limit, where $N_{f}\sim N_{c}$ should enhance the pressure of the confined matter. However it turns out that such contributions are not captured only by backreacting the flavor action to the geometry, but one should include nontrivial string loop diagrams Alho  _et al._ (2015).
  * Zyla  _et al._ (2020) P. A. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Rischke  _et al._ (1991) D. H. Rischke, M. I. Gorenstein, H. Stoecker, and W. Greiner, Z. Phys. C 51, 485 (1991).
  * Hagedorn and Rafelski (1980) R. Hagedorn and J. Rafelski, Phys. Lett. B 97, 136 (1980).
  * Vovchenko  _et al._ (2017a) V. Vovchenko, M. I. Gorenstein, and H. Stoecker, Phys. Rev. Lett. 118, 182301 (2017a), arXiv:1609.03975 [hep-ph] .
  * Vovchenko  _et al._ (2017b) V. Vovchenko, A. Motornenko, P. Alba, M. I. Gorenstein, L. M. Satarov, and H. Stoecker, Phys. Rev. C 96, 045202 (2017b), arXiv:1707.09215 [nucl-th] .
  * Vovchenko  _et al._ (2015) V. Vovchenko, D. V. Anchishkin, M. I. Gorenstein, and R. V. Poberezhnyuk, Phys. Rev. C 92, 054901 (2015), arXiv:1506.05763 [nucl-th] .
  * Redlich and Zalewski (2016) K. Redlich and K. Zalewski, Acta Phys. Polon. B 47, 1943 (2016), arXiv:1605.09686 [cond-mat.quant-gas] .
  * Toki  _et al._ (1995) H. Toki, D. Hirata, Y. Sugahara, K. Sumiyoshi, and I. Tanihata, Nuclear Physics A 588, c357 (1995), proceedings of the Fifth International Symposium on Physics of Unstable Nuclei.
  * Audi  _et al._ (2003) G. Audi, A. Wapstra, and C. Thibault, Nuclear Physics A 729, 337 (2003), the 2003 NUBASE and Atomic Mass Evaluations.
  * Geng  _et al._ (2005) L.-S. Geng, H. Toki, and J. Meng, Prog. Theor. Phys. 113, 785 (2005), arXiv:nucl-th/0503086 .
  * Hempel  _et al._ (2012) M. Hempel, T. Fischer, J. Schaffner-Bielich, and M. Liebendorfer, Astrophys. J. 748, 70 (2012), arXiv:1108.0848 [astro-ph.HE] .
  * Steiner  _et al._ (2013) A. W. Steiner, M. Hempel, and T. Fischer, Astrophys. J. 774, 17 (2013), arXiv:1207.2184 [astro-ph.SR] .
  * Akmal  _et al._ (1998b) A. Akmal, V. R. Pandharipande, and D. G. Ravenhall, Phys. Rev. C 58, 1804 (1998b).
  * com (b) https://compose.obspm.fr/eos/197 (b).
  * com (c) https://compose.obspm.fr/eos/198 (c).
  * com (d) https://compose.obspm.fr/eos/199 (d).
  * Baldo and Burgio (2016) M. Baldo and G. F. Burgio, Prog. Part. Nucl. Phys. 91, 203 (2016), arXiv:1606.08838 [nucl-th] .
  * Note (3) We use $T_{\mathrm{min}}=0.1$ MeV instead of $T=0$ when evaluating these definitions numerically because this value is the lowest available in the data grid of the HS(DD2) EoS.
  * Karsch  _et al._ (2003a) F. Karsch, K. Redlich, and A. Tawfik, Eur. Phys. J. C 29, 549 (2003a), arXiv:hep-ph/0303108 .
  * Karsch  _et al._ (2003b) F. Karsch, K. Redlich, and A. Tawfik, Phys. Lett. B 571, 67 (2003b), arXiv:hep-ph/0306208 .
  * Annala  _et al._ (2018a) E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Phys. Rev. Lett. 120, 172703 (2018a), arXiv:1711.02644 [astro-ph.HE] .
  * Most  _et al._ (2018) E. R. Most, L. R. Weih, L. Rezzolla, and J. Schaffner-Bielich, Phys. Rev. Lett. 120, 261103 (2018), arXiv:1803.00549 [gr-qc] .
  * Annala  _et al._ (2020) E. Annala, T. Gorda, A. Kurkela, J. Nättilä, and A. Vuorinen, Nature Phys. 16, 907 (2020), arXiv:1903.09121 [astro-ph.HE] .
  * Annala  _et al._ (2021) E. Annala, T. Gorda, E. Katerini, A. Kurkela, J. Nättilä, V. Paschalidis, and A. Vuorinen, (2021), arXiv:2105.05132 [astro-ph.HE] .
  * Komoltsev and Kurkela (2021) O. Komoltsev and A. Kurkela, (2021), arXiv:2111.05350 [nucl-th] .
  * Abbott  _et al._ (2017a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017a), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ (2017b) B. P. Abbott _et al._ (LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, Insight-Hxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GW-EM, DES, DLT40, GRAWITA, Fermi-LAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, J-GEM, GROWTH, JAGWAR, CaltechNRAO, TTU-NRAO, NuSTAR, Pan-STARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKI-GW Follow-up, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR, SKA South Africa/MeerKAT), Astrophys. J. Lett. 848, L12 (2017b), arXiv:1710.05833 [astro-ph.HE] .
  * Rezzolla  _et al._ (2018) L. Rezzolla, E. R. Most, and L. R. Weih, Astrophys. J. Lett. 852, L25 (2018), arXiv:1711.00314 [astro-ph.HE] .
  * Stephanov (2004) M. A. Stephanov, Prog. Theor. Phys. Suppl. 153, 139 (2004), arXiv:hep-ph/0402115 .
  * Ayala  _et al._ (2021) A. Ayala, B. A. Zamora, J. J. Cobos-Martínez, S. Hernández-Ortiz, L. A. Hernández, A. Raya, and M. E. Tejeda-Yeomans, (2021), arXiv:2108.02362 [hep-ph] .
  * Aryal  _et al._ (2021) K. Aryal, C. Constantinou, R. L. S. Farias, and V. Dexheimer, Universe 7, 454 (2021), arXiv:2109.14787 [nucl-th] .
  * Yang (2017) C. Yang (STAR), Nucl. Phys. A 967, 800 (2017).
  * Adamczyk  _et al._ (2017) L. Adamczyk _et al._ (STAR), Phys. Rev. C 96, 044904 (2017), arXiv:1701.07065 [nucl-ex] .
  * An  _et al._ (2022) X. An _et al._ , Nucl. Phys. A 1017, 122343 (2022), arXiv:2108.13867 [nucl-th] .
  * Carbone and Schwenk (2019) A. Carbone and A. Schwenk, Phys. Rev. C 100, 025805 (2019), arXiv:1904.00924 [nucl-th] .
  * Bauswein  _et al._ (2010) A. Bauswein, H. T. Janka, and R. Oechslin, Phys. Rev. D 82, 084043 (2010), arXiv:1006.3315 [astro-ph.SR] .
  * De Pietri  _et al._ (2020) R. De Pietri, A. Feo, J. A. Font, F. Löffler, M. Pasquali, and N. Stergioulas, Phys. Rev. D 101, 064052 (2020), arXiv:1910.04036 [gr-qc] .
  * Xie  _et al._ (2020) X. Xie, I. Hawke, A. Passamonti, and N. Andersson, Phys. Rev. D 102, 044040 (2020), arXiv:2005.13696 [astro-ph.HE] .
  * Figura  _et al._ (2020) A. Figura, J. J. Lu, G. F. Burgio, Z. H. Li, and H. J. Schulze, Phys. Rev. D 102, 043006 (2020), arXiv:2005.08691 [gr-qc] .
  * Cromartie  _et al._ (2019) H. T. Cromartie _et al._ (NANOGrav), Nature Astron. 4, 72 (2019), arXiv:1904.06759 [astro-ph.HE] .
  * Fonseca  _et al._ (2021) E. Fonseca _et al._ , Astrophys. J. Lett. 915, L12 (2021), arXiv:2104.00880 [astro-ph.HE] .
  * Riley  _et al._ (2019) T. E. Riley _et al._ , Astrophys. J. Lett. 887, L21 (2019), arXiv:1912.05702 [astro-ph.HE] .
  * Miller  _et al._ (2019) M. Miller _et al._ , Astrophys. J. Lett. 887, L24 (2019), arXiv:1912.05705 [astro-ph.HE] .
  * Miller  _et al._ (2021) M. C. Miller _et al._ , Astrophys. J. Lett. 918, L28 (2021), arXiv:2105.06979 [astro-ph.HE] .
  * Riley  _et al._ (2021) T. E. Riley _et al._ , Astrophys. J. Lett. 918, L27 (2021), arXiv:2105.06980 [astro-ph.HE] .
  * Nättilä  _et al._ (2017) J. Nättilä, M. C. Miller, A. W. Steiner, J. J. E. Kajava, V. F. Suleimanov, and J. Poutanen, Astron. Astrophys. 608, A31 (2017), arXiv:1709.09120 [astro-ph.HE] .
  * Abbott  _et al._ (2018) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 121, 161101 (2018), arXiv:1805.11581 [gr-qc] .
  * Alford  _et al._ (2008) M. G. Alford, A. Schmitt, K. Rajagopal, and T. Schäfer, Rev. Mod. Phys. 80, 1455 (2008), arXiv:0709.4635 [hep-ph] .
  * Chen  _et al._ (2010) H.-Y. Chen, K. Hashimoto, and S. Matsuura, JHEP 02, 104 (2010), arXiv:0909.1296 [hep-th] .
  * Basu  _et al._ (2011) P. Basu, F. Nogueira, M. Rozali, J. B. Stang, and M. Van Raamsdonk, New J. Phys. 13, 055001 (2011), arXiv:1101.4042 [hep-th] .
  * Bitaghsir Fadafan  _et al._ (2018) K. Bitaghsir Fadafan, J. Cruz Rojas, and N. Evans, Phys. Rev. D 98, 066010 (2018), arXiv:1803.03107 [hep-ph] .
  * Faedo  _et al._ (2019) A. F. Faedo, D. Mateos, C. Pantelidou, and J. Tarrío, JHEP 05, 106 (2019), arXiv:1807.09712 [hep-th] .
  * Ghoroku  _et al._ (2019) K. Ghoroku, K. Kashiwa, Y. Nakano, M. Tachibana, and F. Toyoda, Phys. Rev. D 99, 106011 (2019), arXiv:1902.01093 [hep-th] .
  * Henriksson  _et al._ (2019) O. Henriksson, C. Hoyos, and N. Jokela, JHEP 09, 088 (2019), arXiv:1907.01562 [hep-th] .
  * Hoyos  _et al._ (2020) C. Hoyos, N. Jokela, M. Järvinen, J. G. Subils, J. Tarrio, and A. Vuorinen, Phys. Rev. Lett. 125, 241601 (2020), arXiv:2005.14205 [hep-th] .
  * Hoyos  _et al._ (2021a) C. Hoyos, N. Jokela, M. Järvinen, J. G. Subils, J. Tarrio, and A. Vuorinen, (2021a), arXiv:2109.12122 [hep-th] .
  * Schmitt and Shternin (2018) A. Schmitt and P. Shternin, Astrophys. Space Sci. Libr. 457, 455 (2018), arXiv:1711.06520 [astro-ph.HE] .
  * Hoyos  _et al._ (2016) C. Hoyos, D. Rodríguez Fernández, N. Jokela, and A. Vuorinen, Phys. Rev. Lett. 117, 032501 (2016), arXiv:1603.02943 [hep-ph] .
  * Annala  _et al._ (2018b) E. Annala, C. Ecker, C. Hoyos, N. Jokela, D. Rodríguez Fernández, and A. Vuorinen, JHEP 12, 078 (2018b), arXiv:1711.06244 [astro-ph.HE] .
  * Bitaghsir Fadafan  _et al._ (2020) K. Bitaghsir Fadafan, J. Cruz Rojas, and N. Evans, Phys. Rev. D 101, 126005 (2020), arXiv:1911.12705 [hep-ph] .
  * Mamani  _et al._ (2020) L. A. H. Mamani, C. V. Flores, and V. T. Zanchin, Phys. Rev. D 102, 066006 (2020), arXiv:2006.09401 [hep-th] .
  * Bitaghsir Fadafan  _et al._ (2021) K. Bitaghsir Fadafan, J. Cruz Rojas, and N. Evans, Phys. Rev. D 103, 026012 (2021), arXiv:2009.14079 [hep-ph] .
  * Ghoroku  _et al._ (2021) K. Ghoroku, K. Kashiwa, Y. Nakano, M. Tachibana, and F. Toyoda, Phys. Rev. D 104, 126002 (2021), arXiv:2107.14450 [hep-th] .
  * Hoyos  _et al._ (2021b) C. Hoyos, N. Jokela, and A. Vuorinen, (2021b), arXiv:2112.08422 [hep-th] .
  * Kovensky  _et al._ (2021) N. Kovensky, A. Poole, and A. Schmitt, (2021), arXiv:2111.03374 [hep-ph] .
  * Bauswein  _et al._ (2019) A. Bauswein, N.-U. F. Bastian, D. B. Blaschke, K. Chatziioannou, J. A. Clark, T. Fischer, and M. Oertel, Phys. Rev. Lett. 122, 061102 (2019), arXiv:1809.01116 [astro-ph.HE] .
  * Prakash  _et al._ (2021) A. Prakash, D. Radice, D. Logoteta, A. Perego, V. Nedora, I. Bombaci, R. Kashyap, S. Bernuzzi, and A. Endrizzi, Phys. Rev. D 104, 083029 (2021), arXiv:2106.07885 [astro-ph.HE] .