Finite-temperature expansion of the dense-matter equation of state

Debora Mroczek Nanxi Yao Katherine Zine Jacquelyn Noronha-Hostler Illinois Center for Advanced Studies of the Universe, Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA Veronica Dexheimer Center for Nuclear Research, Department of Physics, Kent State University, Kent, OH 44243 USA Alexander Haber Department of Physics, Washington University in St. Louis, 63130 Saint Louis, MO, USA Elias R. Most TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA Walter Burke Institute for Theoretical Physics, California Institute of Technology, Pasadena, CA 91125, USA

(May 2, 2024)

Abstract

In this work we provide a new, well-controlled expansion of the equation of state of dense matter from zero to finite temperatures ( $T$ ), while covering a wide range of charge fractions ( $Y_{Q}$ ), from pure neutron to isospin symmetric nuclear matter. Our expansion can be used to describe neutron star mergers and core-collapse supernova explosions using as a starting point neutron star observations, while maintaining agreement with laboratory data, in a model independent way. We suggest new thermodynamic quantities of interest that can be calculated from theoretical models or directly inferred by experimental data that can help constrain the finite $T$ equation of state. With our new method, we can quantify the uncertainty in our finite $T$ and $Y_{Q}$ expansions in a well-controlled manner without making assumptions about the underlying degrees of freedom. We can reproduce results from a microscopic equation of state up to $T=100$ MeV for baryon chemical potential $\mu_{B}\gtrsim 1100$ MeV ( $\sim 1-2\ n_{\rm sat}$ ) within $5\%$ error, with even better results for larger $\mu_{B}$ and/or lower $T$ . We investigate the sources of numerical and theoretical uncertainty and discuss future directions of study.

I Introduction

The dynamical description of the formation of neutron stars and their mergers (specifically in the post-merger phase) requires a multidimensional equation of state (EOS) in the phase space of finite temperature, baryon number density, and electric charge fraction, $\left\{T,n_{B},Y_{Q}^{\rm QCD}\right\}$ . The charge fraction is defined as the quantum chromodynamics (QCD) – hadronic and quark – electric charge, $Q$ , per baryon, $B$ , $Y_{Q}^{\rm QCD}=\frac{Q^{\rm QCD}}{B}$ , or dividing be the volume, the hadronic and quark charge density per baryon density, $Y_{Q}^{\rm QCD}=\frac{n_{Q}^{\rm QCD}}{n_{B}}$ . Its value usually goes from $Y_{Q}^{\rm QCD}=0$ for pure neutron matter to $Y_{Q}^{\rm QCD}=0.5$ to isospin symmetric matter containing the same number of protons and neutrons.

In fully evolved (beyond the proto-neutron star stage) neutron stars, $Y_{Q}^{\rm QCD}=0.01-0.2$ , with charge neutrality being enforced by leptons, $Y_{Q}^{\rm net}=Y_{Q}^{\rm QCD}+Y_{Q}^{\rm lep}=0$ . In this case, $Y_{Q}^{\rm QCD}$ changes as a function of $n_{B}$ and the temperature is low enough when compared to the Fermi energy that it can be approximated by $T=0$ . Additionally, the leptons are considered to be in weak-interaction ( $\beta$ -)equilibrium with the baryons (protons, neutrons, hyperons, etc.). In this case, the EOS can be described by a simple one-dimensional function of $n_{B}$ . Thus, even the best information on the mass, radius, and tidal deformability of neutron stars from astrophysical observations (e.g., by NICER Miller et al. (2019); Riley et al. (2019); Miller et al. (2021); Riley et al. (2021) or LIGO Abbott et al. (2019)) only constrain the EOS in the limit of zero temperature and $\beta$ -equilibrium. These constraints are insufficient for extrapolating the EOS in the additional dimensions required for dynamical calculations, such as numerical simulations of supernova explosions and neutron star mergers. A multidimensional description of the EOS requires laboratory data from experiments capable of probing larger $Y_{Q}^{\rm QCD}$ and finite $T$ , see Ref. Kumar et al. (2023) for a complete review.

Numerical relativity calculations of merging binary neutron stars have incorporated finite-temperate effects in different ways Bauswein et al. (2010); Raithel and Paschalidis (2023a); Figura et al. (2020); Perego et al. (2019). One approach is to take a cold neutron star EOS and add a thermal contribution to the EOS based on a semi-degenerate, ideal neutron gas Bauswein et al. (2010); Figura et al. (2020); Villa-Ortega et al. (2023). The drawback is that a constant adiabatic index is used, which may not approximate full microscopic EOS calculations well. Another possibility is to incorporate finite $T$ effects via an effective particle mass motivated by the symmetry energy expansion Constantinou et al. (2015); Raithel et al. (2021); Most and Raithel (2021); Fields et al. (2023); Raithel and Paschalidis (2023a); Jacobi et al. (2023); Raithel and Paschalidis (2023b). In this approach, it is assumed that the underlying degrees-of-freedom consist of only protons, neutrons, and leptons, excluding hyperons and/or quarks (see also Blacker et al. (2024, 2023)). Most commonly, fully tabulated EOS are used that depend on the nuclear composition and temperature. However, these are usually coupled to a fixed set of nuclear physics parameters. Furthermore, only a few EOS available, and even fewer with phase transitions or exotic degrees of freedom. A number of studies exist that compare the consequences of these microscopic assumptions on finite $T$ on the neutron star EOS Figura et al. (2020); Wei et al. (2021); Raduta (2022). The drawback in this case is that the degrees of freedom and functional form of the $\beta$ -equilibrated, cold neutron star EOS are fixed.

Rather, what we would like to construct is a generic framework that allows for the systematic study of finite $T$ effects without assuming specific degrees-of-freedom and phases of matter within the $\beta$ -equilibrated, cold neutron star EOS. This framework should be thermodynamically consistent such that the resulting EOS is causal, stable, and respects thermodynamic identities. Two other conditions that we would like to fulfill are: 1) all information about the EOS can be found at $T=0$ , such that we do not require information from regimes of the QCD phase diagram that are hard to access and 2) stronger connections to heavy-ion collision physics, considering new data available from ongoing heavy-ion collision experiments Abdallah et al. (2022a) that are providing new insights into the cold, dense EOS Lovato et al. (2022); Oliinychenko et al. (2023); Omana Kuttan et al. (2023); Yao et al. (2023). Additionally, new heavy-ion experiments are currently being built, such as the CBM at FAIR Almaalol et al. (2022) and FRIB400 Sorensen et al. (2024) that will provide significantly more data in years to come.

In this work, we have developed an entirely new finite temperature expansion of the QCD EOS to be used in simulations of neutron star mergers. The premise is that given a $\beta-$ equilibrated EOS for cold nuclear matter, we can then expand across both in the $Y_{Q}$ direction using a modified symmetry energy expansion (see Yao et al. (2023)) and in the finite $T$ direction to obtain a three-dimensional (3D) EOS that can be used directly in numerical relativity simulations. Additionally, we match our new expanded EOS to a finite-temperature crust at low densities using a method that ensures thermodynamic consistency. The code developed for this work is called FiniteT and will be available online upon publication Mroczek et al. (2024).

I.1 Executive Summary

Methodology: in this paper we develop a new method for generating a 3D EOS based on a given cold, $\beta$ -equilibrated EOS (1D). Our approach relies on an existing expansion (the symmetry energy expansion) and two entirely new expansions (the finite $T$ expansion and the $s/n_{B}$ expansion). Here we provide an brief overview of the theoretical framework discussed in Sec. II.

•

Going from the baryon number density dependent zero temperature, $\beta$ -equilibrated EOS to arbitrary charge fraction. This step is based on previous work Yao et al. (2023) on the symmetry energy expansion and requires knowledge about the EOS at zero temperature, around saturation density for symmetric nuclear matter which can be extracted from low-energy nuclear physics experiments and effective theories.
•

Going from the zero temperature, arbitrary baryon number density and charge fraction EOS to finite temperature, assuming knowledge about the EOS across all charge fractions. We propose a new, well-controlled expansion of the pressure to finite temperatures that requires input from the cold EOS obtained in the previous step and a heat capacity term at zero temperature that is dependent only on the charge fraction and baryon number density.
•

Obtaining the charge fraction and baryon number density dependence of the heat capacity. We propose yet another expansion to account for the charge fraction dependence of the heat capacity. We expand around a fixed charged fraction that can be matched to that of nuclei used in heavy-ion collision experiments. This expansion provides a direct connection to heavy-ion data via the finite temperature EOS of (approximately) symmetric nuclear matter.

Main conclusions: The primary focus of this manuscript is to formalize the framework behind a new open-source code called FiniteT Mroczek et al. (2024) to be released upon the publication of this work. Because the finite $T$ expansion is performed using chemical potentials, its accuracy across a grid of chemical potentials is at the percent level even up to $T\sim 100$ MeV for the core of neutron stars. However, switching to a grid of $n_{B},Y_{Q}$ does introduce some numerical error from various sources that we outline in Sec. IV. We also discuss strategies to improve this error in future work. Our new framework also allows for us to quantify certain features of the QCD phase diagram. For instance, our finite $T$ expansion is driven by the heat capacity $\partial s/\partial T|_{T=0}$ that we calculate in different limits such as heavy-ion collisions and neutron star mergers. From this quantification of the heat capacity in Sec. III.4, we determined that in a relativistic mean-field EOS Alford et al. (2022) the EOS was more sensitive to finite $T$ effects for heavy-ion collisions than in neutron stars. In future work, we can compare different features of the QCD phase diagram using the thermodynamic variables necessary for our 3D framework in order to better quantify thermal effects.

II Theoretical Framework

In the following, we describe the equations that we use to take a zero-temperature, $\beta$ -equilibrated EOS and expand, first to slices along $Y_{Q}=$ const., then to finite temperature, in order to create a 3D EOS.

Before discussing the details of these expansions, we remind the reader of a few thermodynamic relations that are vital to our work. The first is the Gibbs-Duhem relation:

p+\varepsilon=sT+\sum_{X}n_{X}\mu_{X},

(1)

where $p$ is the pressure, $\varepsilon$ the energy density, $s$ the entropy density, $T$ the temperature, and $n$ and $\mu$ are the number densities and chemical potentials associated with the different $X$ conserved charges in the system. For instance, when discussing leptons, there is only one conserved charge, electric charge, (assuming that neutrinos can escape from the system) so there is only contribution of $n_{Q}\mu_{Q}$ . However, for a system that contains baryons that carry electric charge, one has two contributions: $n_{B}\mu_{B}+n_{Q}\mu_{Q}$ . These can be rewritten in terms of $Y_{Q}$ such that the last term in Eq. 1 becomes $n_{B}\left(\mu_{B}+Y_{Q}\mu_{Q}\right)$ . Furthermore, in heavy-ion collisions, where it is generally believed that strangeness ( $S$ ) is conserved¹¹1This is because of the extremely short time scales of $\tau_{\rm HIC}\sim 10^{-22}$ s that are significantly shorter than time scales for the quickest weak decays to occur, i.e., $\tau_{\rm weak}\gtrsim 10^{-10}$ s. In contrast, neutron star mergers occur on the order of milliseconds, i.e., $\tau_{\rm NSM}\sim 10^{-6}$ s such that weak decays play a role., one has a third contribution of $n_{B}\mu_{B}+n_{Q}\mu_{Q}+n_{S}\mu_{S}$ .

Because we work with two conserved charges in this paper ( baryon number and electric charge), we introduce the following notation for the chemical potentials:

\vec{\mu}=\left\{\mu_{B},\mu_{Q}\right\},

(2)

and the densities:

\vec{n}=\left\{n_{B},n_{Q}\right\},

(3)

that is regularly used in our equations. For much of this work, it is quite important to understand what are the “free variables” (sometimes also called “fixed variables”) vs the dependent variables in each case. This is equivalent to working in different ensembles. The relevant thermodynamical potential in high-energy physics is the grand canonical potential, $\Omega/V=-p$ , where $V$ is the volume. Thus, most microscopic models that describe high-energy systems are computed in terms of of $\left\{T,\vec{\mu}\right\}$ , the natural parameters in the grand canonical potential.

If we take derivatives of the pressure with respect to the chemical potentials, then we obtain susceptibilities of the system:

\chi_{ijk}^{BSQ}=\frac{\partial^{i+j+k}p}{\partial^{i}\mu_{B}\partial^{j}\mu_{% S}\partial^{k}\mu_{Q}}\Bigg{|}_{T},

(4)

where the indices $i,j,k$ can run from 0 to $\infty$ but, in practice, we are only aware of these derivatives being calculated from QCD up to $8^{th}$ order Borsanyi et al. (2018); Bazavov et al. (2020); Borsányi et al. (2023). Susceptibilities are important because they can be indicators of phase transitions. For instance, $\chi_{1}$ has a jump in a first-order phase transition and $\chi_{2}$ diverges at a critical point (second-order phase transition). Additionally, they can be calculated from lattice QCD and used to reconstruct the finite $T$ /low density equation of state that is used for heavy-ion collisions Monnai et al. (2019); Noronha-Hostler et al. (2019).

Here is is important to point out that the first susceptibilities relate to the respective densities of a conserved charge, i.e.,

n_{X}=\chi_{1}^{X}=\frac{\partial p}{\partial\mu_{X}}\Big{|}_{T,\mu_{Y\neq X}},

(5)

where $Y$ are any other conserved charges in the system that are not $X$ , such that if we know the pressure at a fixed point in $\left\{T,\vec{\mu}\right\}$ space, then we are able to recover the densities. Note that when taking a derivative of $p$ with respect to $\mu_{X}$ , the derivative will pick up the sign of the conserved charge $X$ . If we isolate the partial pressure $p_{i}$ (i.e., the pressure contribution just from a specific species $i$ ) then this derivative provides the charge number density for a conserved charge $X$ for species $i$

n_{X,i}=\frac{\partial p_{i}}{\partial\mu_{X}}\Big{|}_{T,\mu_{Y\neq X}},

(6)

that picks up the sign of $X_{i}$ . In QCD, relevant conserved charges are electric charge $Q$ , baryon number $B$ , and strangeness $S$ such that we can define electric charge number densities $n_{Q,i}$ , baryon number density $n_{B,i}$ (or simply $n_{i}$ ), and strangeness number density $n_{S,i}$ where each may be either positive or negative depending on the particle $i$ . As an example, if we have the partial pressure for the electrons, $p_{e^{-}}$ , then $Q_{e^{-}}=-1$ and $n_{Q,e^{-}}<0$ . Hyperons can have a mixture of positive and negative charge number densities. As an example, the $\Sigma^{-}$ has $Q_{\Sigma^{-}}=-1$ , $B_{\Sigma^{-}}=+1$ , $S_{\Sigma^{-}}=-1$ such that $n_{Q,\Sigma^{-}}<0$ , $n_{B,\Sigma^{-}}>0$ , $n_{S,\Sigma^{-}}<0$ .

The higher-order susceptibilities ( $i+j+k>2$ ) encode other important information about the EOS. For instance, the second baryon susceptibility must be positive $\chi_{2}^{B}>0$ for thermodynamically stable matter. The higher-order susceptibilities can be related to moments of the distribution of net-particle number Karsch (2012) and be used to search for first or second-order phase transitions Stephanov (2009, 2011).

Derivatives of the pressure with respect to temperature are also important in our calculations. Entropy is the first derivative of the pressure when holding the chemical potentials of the system fixed, i.e.,

s=\frac{\partial p}{\partial T}\Bigg{|}_{\mu_{X}},

(7)

and the temperature derivative of the entropy is the heat capacity

C_{T}=\frac{\partial s}{\partial T}\Bigg{|}_{\mu_{X}}=\frac{\partial^{2}p}{% \partial T^{2}}\Bigg{|}_{\mu_{X}},

(8)

where it also can be rewritten as the second derivative of the pressure.

The third law of thermodynamics states that the entropy of a closed system approaches a constant value at $T=0$ but, in practice, the entropy generally vanishes at $T=0$ , i.e., $s(T=0)=0$ . An exception to $s(T=0)\neq 0$ is if there is residual entropy resulting from the system being stuck in a configuration where the energy is not minimized at $T=0$ , or if the state of minimized energy is not unique. At the time of writing this paper, we are not aware of any microscopic model that produces residual entropy for neutron stars and, therefore, we argue that as $T\rightarrow 0$ we expect $s\rightarrow 0$ . However, it may be interesting to explore this possibility in microscopic models, especially ones that could potentially consider a glass or solid state within a neutron star.

Note that certain derivatives are more naturally taken on a grid of $\left\{s,\vec{n}\right\}$ (for instance, the speed of sound). The chemical potentials are the amount of energy required to add a new particle of that specific conserved quantity to the system. For example, when it comes to baryons it is the amount of energy required to add another baryon to the system. Thus, the equation for the chemical potential at $T=0$ can be written as

\mu_{X}=\frac{\partial\varepsilon}{\partial n_{X}}\Bigg{|}_{n_{Y\neq X,}}.

(9)

II.1 Removing lepton contributions

We use the symmetry energy expansion to extrapolate a $\beta$ -equilibrium EOS to a range of $Y^{\rm QCD}_{Q}$ . The first step when starting with a $\beta$ -equilibrated EOS is to isolate the QCD (hadronic and quark) contribution. To do so, we first discuss charge neutrality within neutron stars.

The net charge density of hadrons/quarks and leptons can be defined as:

	$\displaystyle n_{Q}^{QCD}$	$\displaystyle=$	$\displaystyle\sum_{i}^{\rm had}n_{Q,i}+\sum_{i}^{\rm quarks}n_{Q,i},$		(10)
	$\displaystyle n_{Q}^{\rm lep}$	$\displaystyle=$	$\displaystyle\sum_{i}^{\rm lep}n_{Q,i},$		(11)

where $n_{Q,i}$ is the electric charge number density for an individual hadron, quark, or lepton. In principle, both positive and negative charges are possible, e.g., the electric charge of a positron is $Q_{e^{+}}=+1$ and the electric charge of an electron is $Q_{e^{-}}=-1$ that lead to $n_{Q,e^{+}}>0$ and $n_{Q,e^{-}}<0$ , respectively. Overall, hadron and quark contributions to the net charge are positive (although down, strange quarks and negative hyperons contribute negatively), i.e., $n_{Q}^{\rm had},~{}n_{Q}^{\rm quark}>0$ , and lepton contributions to the net charge are negative, i.e., $n_{Q}^{\rm lep}<0$ .

For stability, a neutron star should be electrically neutral so the net charge density of hadrons and quarks $n_{Q}^{\rm QCD}$ is identically balanced by the net charge density of leptons $n_{Q}^{\rm lep}$ , implying

n_{Q}^{\rm net}=n_{Q}^{\rm QCD}+n_{Q}^{\rm lep}=0.

(12)

For the simplest possible description of neutron stars, containing only neutrons, protons, and electrons (known as npe matter), Eq. (12) simplifies to

\displaystyle n_{Q}^{\rm QCD}\rightarrow n_{Q,p},\quad n_{Q}^{\rm lep}% \rightarrow n_{Q,e},\quad n_{Q,p}=-n_{Q,{e^{-}}}.

(13)

In this case, $\beta$ -equilibrium is reached through a balance between electron capture and neutron decay,

	$\displaystyle p+e^{-}$	$\displaystyle\leftrightarrow n+\nu_{e},$		(14)
	$\displaystyle n$	$\displaystyle\leftrightarrow p+e^{-}+\nu_{e},$		(15)

where, in the neutrino free-streaming regime (e.g., at low T), it is typically assumed that the mean free-path of neutrinos is large enough that they can escape ( $n_{\nu_{e}}=0$ ) and thus the reverse rates for the equations above are zero. The, at zero-temperature, $\beta$ -equilibrium it is required that the chemical potentials follow the relation

\mu_{n}=\mu_{p}+\mu_{e},

(16)

and muonic contributions in weak $\beta$ -equilibrium are included via

\mu_{\mu}=\mu_{e},

(17)

due to electron-muon converting weak processes. Note that this is in reference to the zero-temperature, $\beta-$ equilibrated EOS only, as in differences can arise in, e.g., the post-merger phase of neutron star mergers Alford and Harris (2018); Alford et al. (2021); Loffredo et al. (2023); Alford et al. (2023).

Eq. (12) is generic but the EOS may be much more complex than simple npe matter, such that $n_{Q}^{\rm QCD}$ can include contributions from various hadrons, such as $\Lambda,\Sigma,\Xi,\Delta$ , etc. Additionally, other leptons can contribute to $n_{Q}^{\rm lep}$ , such as $\mu^{-}$ . Thus, in this work we use super/subscripts of “QCD” and “lep” in order to distinguish QCD vs QED (quantum electrodynamics) contributions to various thermodynamic observables. Within the total EOS, there are contributions from both QCD and leptons to the pressure, energy density, entropy density, and baryon density,

$\displaystyle p^{\rm tot}$	$\displaystyle=$	$\displaystyle p^{\rm QCD}+p^{\rm lep},$	(18)
$\displaystyle\varepsilon^{\rm tot}$	$\displaystyle=$	$\displaystyle\varepsilon^{\rm QCD}+\varepsilon^{\rm lep},$	(19)
$\displaystyle s^{\rm tot}$	$\displaystyle=$	$\displaystyle s^{\rm QCD}+s^{\rm lep},$	(20)
$\displaystyle n_{B}^{\rm tot}$	$\displaystyle=$	$\displaystyle n_{B}^{\rm QCD}.$	(21)

The lepton contributions can also be determined for any generic EOS from Eq. (12). As long as $Y_{Q}$ , $n_{B}$ , and which leptons contribute to the EOS are known, then $n_{Q}^{\rm QCD}$ and $n_{Q}^{\rm lep}$ can be determined, and it is possible to calculate all the leptonic contributions find the total EOS, as described in Eqs. (18-21). One final comment is that most EOS assume charge neutrality for the full 3D phase-space such that Eqs. (18-21) are also applied at finite $T$ .

The symmetry energy expansion that is detailed in the next section is only valid for the QCD sector of the EOS. The QED sector does not contribute to the EOS of nuclear experiments and, therefore, any contributions from the leptons must be removed from descriptions intended for neutron stars before applying the symmetry energy expansion.

To remove the lepton contribution for the zero-temperature, $\beta$ -equilibrated EOS, we use the following algorithm:

1.

Calculate the (free Fermi gas) leptonic EOS across a wide range of $\left\{n_{Q}^{\rm lep}\right\}$ at $T=0$ ;
2.

Determine $n_{Q}^{\rm lep}=-n_{Q}^{\rm QCD}=-n_{B}Y_{Q}^{\rm QCD}$ for each point in the original EOS (we discuss the determination of $Y_{Q}^{\rm QCD}$ in the following);
3.

Obtain $p^{\rm QCD}$ , $\varepsilon^{\rm QCD}$ from the neutron star EoS and Eqs. (18-19) for a wide range of $\left\{n_{B}^{\rm lep}\right\}$ .

Once we have the QCD contributions to the EOS, then the symmetry energy expansion can be applied.

II.1.1 EOS contributions at $T=0$

The $T=0$ limit for a free Fermi gas of spin $1/2$ leptons of one species is:

$\displaystyle\varepsilon^{\rm lep}$	$\displaystyle=$	$\displaystyle\sum_{i}\frac{1}{\pi^{2}}\Bigg{[}\left(\frac{1}{8}m_{i}^{2}k_{F,i% }+\frac{1}{4}k_{F,i}^{3}\right)\sqrt{k_{F,i}^{2}+m_{i}^{2}}$	(22)
	$\displaystyle-$	$\displaystyle\frac{1}{8}m_{i}^{4}\ln{\frac{k_{F,i}+\sqrt{k_{F,i}^{2}+m_{i}^{2}% }}{m_{i}}}\Bigg{]}\ ,$	(22)
$\displaystyle p^{\rm lep}$	$\displaystyle=$	$\displaystyle\sum_{i}\frac{1}{3}\frac{1}{\pi^{2}}\Bigg{[}\left(\frac{1}{4}k_{F% ,i}^{3}-\frac{3}{8}m_{i}^{2}k_{F,i}\right)\sqrt{k_{F,i}^{2}+m_{i}^{2}}$	(23)
	$\displaystyle+$	$\displaystyle\frac{3}{8}m_{i}^{4}\ln{\frac{k_{F,i}+\sqrt{k_{F,i}^{2}+m_{i}^{2}% }}{m_{i}}}\Bigg{]}\ ,$	(23)
$\displaystyle n^{\rm lep}$	$\displaystyle=$	$\displaystyle\sum_{i}\frac{k_{F,i}^{3}}{3\pi^{2}}\ ,$	(24)

where $i=e^{-},\mu^{-}$ , $k_{F,i}$ is the Fermi momentum of particle $i$ , and $m_{i}$ the mass of particle $i$ .

II.1.2 EOS contributions for $T>0$

After the entire 3D EOS is generated for the QCD contribution, we add back the lepton contribution, needed to describe all astrophysical scenarios. But now this means adding leptons to a 3D EOS (at finite $T$ and out of $\beta$ -equilibrium). At finite $T$ we can use an ideal gas of non-interacting fermions. For such a system, we can define a general normalized chemical potential

\hat{\mu}_{i}=B_{i}\frac{\mu_{B}}{T}+S_{i}\frac{\mu_{S}}{T}+Q_{i}\frac{\mu_{Q}% }{T},

(25)

for the particle $i$ where the quantum numbers of the particle such as the baryon number $B_{i}$ , strangeness $S_{i}$ , and electric charge $Q_{i}$ are taken into account. When we consider only leptons then this general chemical potential simplifies to

\hat{\mu}_{i}^{\rm lep}=Q_{i}\frac{\mu_{Q}}{T},

(26)

where it picks up the sign of the electric charge from $Q_{i}$ . At finite temperatures one incorporates the Fermi-Dirac distribution into the integrals

f_{i}^{\rm lep}=\frac{1}{e^{E_{i}/T-\hat{\mu}_{i}^{\rm lep}}+1},

(27)

where the energy of particle $i$ is defined as

E_{i}=\sqrt{k^{2}+m_{i}^{2}}.

(28)

The relevant thermodynamic quantities are then the number density

n^{\rm lep}=\frac{1}{\pi^{2}}\sum_{i}\int_{0}^{\infty}dk~{}{k^{2}}f_{i}^{\rm lep},

(29)

the energy density

\displaystyle\varepsilon^{\rm lep}

\displaystyle=

\displaystyle\sum_{i}\frac{1}{\pi^{2}}\int_{0}^{\infty}dk~{}{k^{2}}E_{i}\;f_{i% }^{\rm lep},

(30)

the pressure

\displaystyle p^{\rm lep}

\displaystyle=

\displaystyle\frac{1}{\pi^{2}}\sum_{i}\int_{0}^{\infty}dk~{}{k^{4}}\frac{1}{E_% {i}}f_{i}^{\rm lep},

(31)

and the entropy

s^{\rm lep}=\frac{1}{T}\left[\varepsilon^{\rm lep}+p^{\rm lep}-n^{\rm lep}\mu_% {Q}\right],

(32)

were the latter was obtained from the Gibbs-Duhem relation.

We now have all the thermodynamic variables needed to add the electron and muon contributions to the 3D EOS. We calculate $n_{\rm Q,QCD}$ directly from our 3D EOS at a point in $(T,n_{B},Y_{Q})$ and then can solve the equation

n_{\rm Q,lep}(T,\mu_{Q})=-n_{\rm Q,QCD}(T,n_{B},Y_{Q}),

(33)

at a fixed $T$ to find the corresponding $\mu_{Q}$ needed to ensure charge neutrality. With the knowledge of $T$ and $\mu_{Q}$ , we can solve Eqs. (30-32) to determine the lepton contribution to the EOS. Then, the lepton contributions are added to the 3D EOS by using Eqs. (18-20).

II.2 Going away from $\beta$ -equilibrium

It is also possible to probe the QCD EOS using flow data from low-energy heavy-ion experiments that can constrain the EOS up to $\sim 5$ $n_{\rm sat}$ at finite temperature Sorensen et al. (2024), and other low-energy nuclear experiments that measure properties of the zero-temperature EOS around $n_{\rm sat}$ Kumar et al. (2023). In heavy-ion collisions where identical ions are used, the proton to nucleon fraction, $Z/A$ , is constant, as electric charge is conserved. In fact, this fraction is equivalent to $Y_{Q}^{\rm QCD}$ , since leptons do not influence the EOS of nuclear systems. Heavy-ion experiments that collide stable nuclei probe $Y_{Q}^{\rm QCD}$ in the range $Y_{Q}^{\rm QCD}=\left\{0.38,0.5\right\}$ , where ${}^{238}U$ corresponds to the most neutron rich nucleus and in contrast ${}^{16}O$ has exactly the same number of neutrons and protons. In astrophysical systems, including neutron stars and supernovae, $Y_{Q}^{\rm QCD}$ is typically within $\left\{0,0.5\right\}$ Lattimer and Prakash (2000). To describe and connect these different environments, we make the following definitions:

•

Isospin symmetric nuclear matter (SNM). In this limit, $Y_{Q}^{\rm QCD}=0.5$ . Information in this regime comes either from nuclear experiments Kumar et al. (2023) or calculations such as chiral effective field theory Tews et al. (2013); Drischler et al. (2019). Because the system is isospin symmetric, $\mu_{Q}=0$ . Leptons are not included.
•

Isospin asymmetric nuclear matter (ANM). In this case pertinent to astrophysics, the charge fraction is a function of baryon density $Y_{Q}^{\rm QCD}(n_{B})$ . At $T=0$ , matter is in $\beta$ -equilibrium, which allow us to determine $\mu_{e}=-\mu_{Q}$ . The charge density from QCD is identically equal to that of leptons: $n_{Q}^{\rm QCD}=n_{Q}^{\rm lep}$ to ensure that the star is electrically neutral, $n_{Q}^{\rm net}=0$ .
•

Pure neutron matter (PNM). In this limit, only neutrons are present in the EOS such that $Y_{Q}^{\rm QCD}=0$ and no leptons are required $n_{Q}^{\rm lep}=0$ . PNM is only a theoretical limit but can be relevant to the expansion schemes described in this section. Information in this regime comes from chiral effective field theory.

We proceed by introducing a parametric form of the binding energy of nucleons valid in the intermediate regime between PNM and SNM. We do so by expanding in $n_{B}$ and $Y_{Q}^{\rm QCD}$ . The difference between the energy (per baryon) within these two limits (SNM and PNM) is known as the symmetry energy,

E_{\rm sym}(n_{B})=\frac{1}{A}\left(E_{\rm PNM}(n_{B})-E_{\rm SNM}(n_{B})% \right),

(34)

where the total energy $E$ can be related to the energy density $\varepsilon$ via

\frac{\varepsilon}{n_{B}}=\frac{E}{A}+m_{n},

(35)

where $m_{n}$ is the mass of the nucleons.

Neutron stars in $\beta$ -equilibrium (ANM), which present a finite value of $Y_{Q}^{\rm QCD}(n_{B})$ that depends on the baryon density, can be used as a starting point for an expansion. First, we can expand around SNM in terms of isospin asymmetry

\delta_{\rm iso}=1-2Y_{Q}^{\rm QCD},

(36)

where the expansion of the energy is Baldo and Burgio (2016)

E_{\rm sym}(n_{B})\delta_{\rm iso}^{2}+\mathcal{O}(\delta_{\rm iso}^{4})=\frac% {1}{n_{B}}\left(\varepsilon_{\rm ANM}-\varepsilon_{\rm SNM}\right),

(37)

where typically terms only up to $\delta_{\rm iso}^{2}$ are taken and any Taylor expansion coefficients are absorbed by $E_{\rm sym}(n_{B})$ .

It is then useful to take a second expansion of Eq. (37) around $n_{\rm sat}$ where nuclear physics properties are known best, such that

	$\displaystyle E_{\rm sym}(n_{B})=E_{\rm sym,sat}+\frac{L_{\rm sym}}{3}\left(% \frac{n_{B}-n_{\rm sat}}{n_{\rm sat}}\right)+$
	$\displaystyle\frac{K_{\rm sym}}{18}\left(\frac{n_{B}-n_{\rm sat}}{n_{\rm sat}}% \right)^{2}+\frac{J_{\rm sym}}{162}\left(\frac{n_{B}-n_{\rm sat}}{n_{\rm sat}}% \right)^{3}+$
	$\displaystyle\mathcal{O}\left(\frac{n_{B}-n_{\rm sat}}{n_{\rm sat}}\right)^{4},$		(38)

where

$\displaystyle E_{\rm sym,sat}$	$\displaystyle=$	$\displaystyle\frac{1}{A}\left(E_{\rm PNM}(n_{\rm sat})-E_{\rm SNM}(n_{\rm sat}% )\right),$	(39)
$\displaystyle L_{\rm sym}$	$\displaystyle=$	$\displaystyle 3n_{\rm sat}\left(\frac{dE_{\rm sym}}{dn_{B}}\right)_{n_{\rm sat% }},$	(40)
$\displaystyle K_{\rm sym}$	$\displaystyle=$	$\displaystyle 9n_{\rm sat}^{2}\left(\frac{d^{2}E_{\rm sym}}{dn_{B}^{2}}\right)% _{n_{\rm sat}},$	(41)
$\displaystyle J_{\rm sym}$	$\displaystyle=$	$\displaystyle 27n_{\rm sat}^{3}\left(\frac{d^{3}E_{\rm sym}}{dn_{B}^{3}}\right% )_{n_{\rm sat}}.$	(42)

It should be clear then that the higher-order coefficients play a more significant role at large $n_{B}$ . Currently, experimental constraints exist on $E_{\rm sym,sat}$ and $L_{\rm sym}$ and theoretical constraints exist for $K_{\rm sym}$ and $J_{\rm sym}$ . Recent work has study the constraints on these coefficientsTews et al. (2017); Zhang et al. (2018); Li et al. (2019); Xie and Li (2020); Adhikari et al. (2021); Reed et al. (2021); Sun et al. (2023) and an alternative expansion scheme Imam et al. (2022).

From this point, we drop terms above the order of $\delta_{\rm iso}^{2}$ and above the order of $\left(\frac{n_{B}-n_{\rm sat}}{n_{\rm sat}}\right)^{3}$ , obtaining:

		$\displaystyle\left[E_{\rm sym,sat}+\frac{L_{\rm sym}}{3}\left(\frac{n_{B}-n_{% \rm sat}}{n_{\rm sat}}\right)+\frac{K_{\rm sym}}{18}\left(\frac{n_{B}-n_{\rm sat% }}{n_{\rm sat}}\right)^{2}+\frac{J_{\rm sym}}{162}\left(\frac{n_{B}-n_{\rm sat% }}{n_{\rm sat}}\right)^{3}\right]\left(1-Y_{Q}^{\rm QCD}(n_{B})\right)^{2}$		(43)
		$\displaystyle=\frac{1}{n_{B}}\left(\varepsilon_{\rm ANM}-\varepsilon_{\rm SNM}% \right),$		(44)

where we can now expand from a $\beta$ -equilibrated EOS corresponding to a $Y_{Q}^{\rm QCD}(n_{B})$ into SNM where $Y_{Q}^{\rm QCD}=0.5$ . Note that linear terms in Eq. (37) always vanish. See Yao et al. (2023) for more details on the derivation of Eq. (43).

We can further generalize this equation to any $Y_{Q}^{\rm const}$ slice that can be used to reconstruct the EOS along a regular grid of $Y_{Q}$

	$\displaystyle\varepsilon_{Y_{Q}^{\rm const}}$	$\displaystyle=$	$\displaystyle\varepsilon_{\rm ANM}-4n_{B}\left[E_{\rm sym,sat}+\frac{L_{\rm sym% }}{3}\left(\frac{n_{B}}{n_{\rm sat}}-1\right)+\frac{K_{\rm sym}}{18}\left(% \frac{n_{B}}{n_{\rm sat}}-1\right)^{2}+\frac{J_{\rm sym}}{162}\left(\frac{n_{B% }}{n_{\rm sat}}-1\right)^{3}\right]$		(45)
		$\displaystyle\times$	$\displaystyle\left[\left(Y_{\rm Q,QCD}^{\rm const}-Y_{\rm Q,QCD}\right)+\left(% Y_{\rm Q,QCD}^{2}-\left(Y_{\rm Q,QCD}^{\rm const}\right)^{2}\right)\right]\ .$		(45)

We have yet to develop a model-independent framework for obtaining $Y_{Q}(n_{B})$ in a neutron star. However, based on the symmetry energy expansion, we can parameterize $Y_{Q}(n_{B})$ at $\beta$ -equilibrium by expanding the proton faction around $n_{\rm sat}$ :

\displaystyle Y_{\rm Q,QCD}^{\rm\beta,sym}(n_{B})

\displaystyle=

\displaystyle\frac{1}{16}\left[8-\frac{\pi^{4/3}n_{B}}{2^{1/3}\lambda}+\left(% \frac{\pi}{2}\right)^{2/3}\frac{\lambda}{E_{\rm sym}^{3}}\right]~{},

(46)

where $\lambda$ is given by

\lambda=\left(-24E_{\rm{\rm sym}}^{6}n_{B}+\sqrt{2}\sqrt{288E_{\rm{\rm sym}}^{% 12}n_{B}^{2}+\pi^{2}E_{\rm{\rm sym}}^{9}n_{B}^{3}}\right)^{1/3}.

(47)

This formula was derived in Yao et al. (2023) for the generic case of an unconstrained value of $Y_{Q}$ and four symmetry energy coefficients (previous work with three symmetry energy coefficients and assuming a small $Y_{Q}$ can be found in Mendes et al. (2021)). One should note that Eq. (46) is the only real solution that was obtained in Yao et al. (2023), although imaginary solutions also exist.

Once we have $\varepsilon_{Y_{Q}^{\rm const}}$ , we can use thermodynamic relations to obtain the remaining thermodynamic quantities. For numerical relativity simulations at $T=0$ , we require the following set of thermodynamic state variables:

\left\{\varepsilon,p,\mu_{B},\mu_{Q}\right\}(n_{B},Y_{Q}^{\rm const}),

(48)

taken on a grid of fixed $(n_{B},Y_{Q}^{\rm const})$ . In order to obtain the chemical potentials, we rewrite the Gibbs-Duhem equation for the case of two conserved quantities, baryon number and electric charge, in the $T=0$ limit,

\varepsilon=-p+n_{B}\left(\mu_{B}+Y_{Q}\mu_{Q}\right),

(49)

where inside the parenthesis we have the system’s Gibbs free energy per baryon Hempel et al. (2013). Along a slice of $Y_{Q}=\frac{n_{Q}}{n_{B}}=$ const., we can derive

\frac{\partial\varepsilon}{\partial n_{B}}\Big{|}_{Y_{Q}}=\mu_{B}+Y_{Q}\mu_{Q},

(50)

\displaystyle\mu_{B}

\displaystyle=

\displaystyle\frac{\partial\varepsilon}{\partial n_{B}}\Big{|}_{Y_{Q}}-Y_{Q}% \mu_{Q},

(51)

and rewriting Eq. 49 in terms of $n_{Q}$

\varepsilon=-p+n_{B}\mu_{B}+n_{Q}\mu_{Q},

(52)

\frac{\partial\varepsilon}{\partial n_{Q}}\Big{|}_{Y_{Q}}=\mu_{Q}+\mu_{B}\frac% {\partial n_{B}}{\partial n_{Q}}\Big{|}_{Y_{Q}}=\mu_{Q}-\frac{\mu_{B}}{Y_{Q}},

(53)

\displaystyle\mu_{Q}

\displaystyle=

\displaystyle\frac{\partial\varepsilon}{\partial n_{Q}}\Big{|}_{Y_{Q}}+\frac{% \mu_{B}}{Y_{Q}},

(54)

where we used that

\displaystyle\dfrac{\partial n_{B}}{\partial n_{Q}}\Big{|}_{Y_{Q}}

\displaystyle=\dfrac{\partial Y_{Q}}{\partial n_{Q}}\Big{|}_{n_{B}}\dfrac{% \partial n_{B}}{\partial Y_{Q}}\Big{|}_{n_{Q}}=\dfrac{1}{n_{B}}\left(-\dfrac{n% _{B}}{Y_{Q}}\right)=-\dfrac{1}{Y_{Q}}.

(55)

The speed of sound can then be calculated along slices of $Y_{Q}^{\rm const}$ using

c_{s}^{2}=\left(\frac{\partial p}{\partial\varepsilon}\right)_{T=0}\ .

(56)

Note that the $c_{s}^{2}$ relation becomes more complicated if the grid is in $\mu_{B},\mu_{Q}$ instead of $n_{B},Y_{Q}$ , see Parotto et al. (2020); Parotto (2019) for a discussion.

II.3 Finite $T$ expansion

We make use of the fact that we are describing matter at large $\mu_{B}$ (or in other words, large $n_{B}$ ) such that an expansion in $T$ is in the regime where $T\ll\mu_{B}$ . In fact, such an expansion is significantly more well-controlled than what is typically done in lattice QCD, where an expansion in $\mu_{B}/T$ to reach finite $n_{B}$ are currently used up to $\mathcal{O}(\mu_{B}/T)=3.5$ Parotto et al. (2022); Borsanyi et al. (2022). In contrast, for neutron star mergers, the maximum temperatures are up to $T<100$ MeV (and even lower for supernova explosions) whereas typical chemical potentials are $\mu_{B}\gtrsim 1000$ MeV, such that $T/\mu_{B}\lesssim 0.1$ Perego et al. (2019); Most et al. (2020); Hammond et al. (2021).

In the following, we derive an expansion in terms of $T$ , while keeping $\mu_{B}=$ const. Note that the calculations in this section are performed at a fixed $\mu_{Q}$ =const., and the $T=0$ dependence of $\mu_{Q}$ is determined by the symmetry energy expansion that discussed in Sec. II.2. Here we consider a dimensionful expansion in $T$ . We find that in this case we only require the second term in the series, i.e., the $T^{2}$ coefficient, that is equivalent to a heat capacity defined at constant $\mu_{B},\mu_{Q}$ . We also derive analytical expressions for most of the relevant thermodynamic observables, specifically s, $n_{B}$ , $Y_{Q}$ , and $\varepsilon$ . Of particular importance to the finite $T$ expansion is the heat capacity $C_{T}$ , defined in Eq. (8). Thus, we must also determine a way to obtain $C_{T}$ at a particular slice in $Y_{Q}$ for a particular $n_{B}$ , i.e., $C_{T}(n_{B},Y_{Q})$ . We derive a new expansion in entropy over baryon density, $s/n_{B}$ , that builds a pathway to make connections to future heavy-ion data. The expansion in $s/n_{B}$ yields $C_{T}$ at fixed $n_{B}$ and $Y_{Q}$ , whereas the finite $T$ expansion has been performed at fixed $\mu_{B}$ and $\mu_{Q}$ . Thus, we also study the consequences of switching between a thermodynamic ensembles of fixed $(n_{B},Y_{Q})$ and $(\mu_{B},\mu_{Q})$ .

Finally, we discuss the consequences of EOS that assume a linear contribution to the pressure in order to obtain a finite $T$ EOS. This linear pressure contribution leads to a finite entropy at vanishing temperatures, which has non-trivial consequences that we point out here.

II.3.1 Expansion in $T$ for fixed $\vec{\mu}$

As a first step, we compute a Taylor series of the pressure, expanding around $T=0$ at fixed $\vec{\mu}$ ,

p(T,\vec{\mu})=p_{T=0}+\frac{\partial p}{\partial T}\biggr{\rvert}_{T=0,\vec{% \mu}}T+\frac{1}{2}\frac{\partial^{2}p}{\partial T^{2}}\biggr{\rvert}_{T=0,\vec% {\mu}}T^{2}+\frac{1}{3!}\frac{\partial^{3}p}{\partial T^{3}}\biggr{\rvert}_{T=% 0,\vec{\mu}}T^{3}+\mathcal{O}(T^{4}),

(57)

where we have expanded up to $T^{3}$ . To be clear in terms of notation, when we write fixed $\vec{\mu}$ this implies both chemical potentials are fixed individually, not that the magnitude of the vector quantity itself is fixed.

Let us recall that the entropy density is defined as $s={\partial p}/{\partial T}|_{\vec{\mu}}$ , which implies that the first term in Eq. (57) is the entropy at $T=0$ . In fact, because the entropy is obtained from the derivative of the pressure at fixed $\vec{\mu}=$ const., it is advantageous for us to calculate the finite $T$ expansion at $\vec{\mu}=$ const. instead of fixed densities. Precisely because of this advantage, we choose to perform our finite $T$ expansion at fixed $\vec{\mu}=$ const. even though numerical relativity simulations will later require a grid in $n_{B},Y_{Q}$ . If we were to rewrite our finite $T$ expansion, holding $n_{B},Y_{Q}$ constant, then extra terms would appear due to the change in thermodynamic ensembles, such that the linear $T$ term would not vanish.

Since we assume the entropy to vanish at $T=0$ , the linear $T$ term should be zero. We can then rewrite the expansion in terms of the entropy density up to third order as

p(T,\vec{\mu})=p_{T=0}+\frac{1}{2}\frac{\partial s}{\partial T}\biggr{\rvert}_% {T=0,\vec{\mu}}T^{2}+\frac{1}{6}\frac{\partial^{2}s}{\partial T^{2}}\biggr{% \rvert}_{T=0,\vec{\mu}}T^{3},

(58)

where $\frac{\partial s}{\partial T}\bigr{\rvert}_{T=0,\vec{\mu}}$ is related to the heat capacity (see the introduction to Sec. II).

Next, we make use the relativistic mean field (RMF) framework developed in Alford et al. (2022) that consists of protons and neutrons coupled to $\sigma$ , $\omega$ , and $\rho$ mesons and is informed by chiral effective field theory. This model can be calculated at finite $T$ , $\mu_{B}$ , $\mu_{Q}$ or $T$ , $n_{B}$ , $Y_{Q}$ wherein we can calculate the coefficients in the series in Eq. (58), understand how they vary with $Y_{Q}$ , and check the accuracy of our expansion in a realistic microscopic model. Later, in Sec. III.4, we calculate the numerical accuracy of our approach, and in Sec. IV we have a discussion on the types of error that enter our approach and ways to improve it for the future.

II.3.2 Connecting $ds/dT$ at different $Y_{Q}$ slices

Up until this point, we have assumed SNM when calculating the coefficients of the finite $T$ expansion, but merger and supernovae simulations require the EOS across a range of $Y_{Q}$ . Therefore, we also need to estimate the dependence of the heat capacity term $ds/dT$ on $Y_{Q}$ . We can then calculate the finite $T$ expansion not just for SNM but across a range of $\mu_{B},\mu_{Q}$ .

As discussed in the previous section, it is more natural to calculate the finite temperature dependence of the EOS at fixed $\vec{\mu}$ . However, heavy-ion collisions are sensitive to the dependence of $s/n_{B}$ at a fixed $Y_{Q}$ . Thus, in order to explore this connection, it is advantageous to calculate $\frac{ds}{dT}$ at fixed $n_{B}$ and $Y_{Q}$ , and study the implications of switching between fixed $\vec{\mu}$ and fixed $n_{B}$ , $Y_{Q}$ .

We derive here the finite temperature expansion in entropy across different values of $Y_{Q}$ and find a method to calculate $\frac{ds}{dT}(Y_{Q})$ along $Y_{Q}=$ const. Note that in the following we assume that $\mu_{S}=0$ . In reality, strangeness neutrality can play a large role in the EOS of heavy-ion collisions such that one must enforce on average $\langle n_{S}\rangle=0$ , which then leads to a finite $\mu_{S}$ at finite $\mu_{B}$ or $\mu_{Q}$ . We leave an investigation of the influence of strangeness for a future work.

At finite $T$ , natural quantities of interest are isentropic hypersurfaces, i.e., regions in the 3-dimensional space of parameters that present the same $s/n_{B}$ . This follows from ideal fluid descriptions of dynamic systems, which remain in such regions as the systems cool or heat up, since both entropy and baryon number are conserved. From heavy-ion collisions, isentropes can be extracted at freeze-out using thermal models Alba et al. (2014, 2020) and compared to lattice QCD results Guenther et al. (2017). Isentropes also take on unique properties around phase transitions Karthein et al. (2021) such that they can concentrate near a critical point, an effect known as critical lensing Dore et al. (2022); Stephanov (2004); Nonaka and Asakawa (2005); Asakawa et al. (2008).

Heavy-ion collisions often collide nuclei that have typical values of $Y_{Q}^{\rm{HIC}}=\frac{Z}{A}\sim 0.4$ , but there is some choice in the exact $Y_{Q}^{\rm{HIC}}$ . The recent isobar run at RHIC Abdallah et al. (2022b) ran ${}^{44}_{96}Ru$ and ${}^{40}_{96}Zr$ beams that have $Y_{Q}^{Ru}=0.46$ and $Y_{Q}^{Zr}=0.42$ , respectively. The global value of $Y_{Q}^{\rm{HIC}}$ remains constant throughout the expansion (note local fluctuations are possible Greif et al. (2018); Fotakis et al. (2020); Carzon et al. (2022); Plumberg et al. (2023), but we ignore that discussion for this work) such that any data used from HIC from a specific species of colliding heavy-ions can be taken with a constant $Y_{Q}^{\rm{HIC}}$ .

Let us start by defining the entropy over baryon number difference, $\Delta s/n_{B}$ , between PNM and HIC:

	$\displaystyle\frac{\Delta s}{n_{B}}\left(T,n_{B}\right)=$	$\displaystyle\frac{s_{\rm PNM}}{n_{B}}\left(T,n_{B},Y_{Q}=0\right)\,$
	$\displaystyle-$	$\displaystyle\frac{s_{\rm{HIC}}}{n_{B}}\left(T,n_{B},Y_{Q}=Y_{Q}^{\rm{HIC}}% \right),$		(59)

where all quantities are inherently at finite $T$ here because we assume entropy vanishes at $T=0$ (although how entropy approaches $T\rightarrow 0$ is what we are seeking to understand). Then, we can expand $s/n_{B}~{}(T,n_{B},Y_{Q})$ around $Y_{Q}=Y_{Q}^{\rm{HIC}}=$ const. just as it is typically done for the symmetry energy expansion. To do so, we define a new quantity

\delta_{\rm HIC}=1-\frac{Y_{Q}}{Y_{Q}^{\rm{HIC}}},

(60)

to make the expansion around $Y_{Q}^{\rm{HIC}}$ possible. Our $\delta_{\rm HIC}$ is similar to the typical isospin asymmetry parameter $\delta_{\rm iso}$ that expands around SNM defined in Eq. (36), but we purposely change the expansion (that is typically done around $Y_{Q}=0.5$ ) to be around a generic constant value of $Y_{Q}^{\rm{HIC}}$ that can correspond to $Z/A$ from a specific species used in heavy-ion collisions. Assuming for this study that it is enough to account only for the even orders of $\delta_{\rm HIC}$ (as it is done for $\delta_{\rm iso}$ ), our expansion at arbitrary $\delta_{\rm HIC}$ becomes:

\tilde{S}(T,n_{B},\delta_{\rm HIC})=\frac{s_{\rm{HIC}}}{n_{B}}(T,n_{B},\delta_% {\rm HIC}=0)+\frac{1}{2}\tilde{S}_{\rm{HIC},2}(T,n_{B})\delta_{\rm HIC}^{2}+% \mathcal{O}\left(\delta_{\rm HIC}^{4}\right),

(61)

where $\tilde{S}(T,n_{B},\delta_{\rm HIC})$ is then at some value of $Y_{Q}(n_{B})$ and

\tilde{S}_{\rm{HIC},2}(T,n_{B})\equiv\frac{\partial^{2}\tilde{S}(T,n_{B},% \delta_{\rm HIC})}{\partial\delta_{\rm HIC}^{2}}\Bigg{|}_{\delta_{\rm HIC}=0},

(62)

where $T$ , $n_{B}$ are fixed, $\delta_{\rm HIC}=0$ implies heavy-ion collisions, and $\delta_{\rm HIC}=1$ implies PNM. From this point on we drop higher-order terms and rewrite Eq. (61) in terms of $Y_{Q}$ instead of $\delta_{\rm HIC}$ for ease of interpretation:

	$\displaystyle\tilde{S}(T,n_{B},Y_{Q})$	$\displaystyle=\frac{s_{\rm{HIC}}}{n_{B}}(T,n_{B},Y_{Q}=Y_{Q}^{\rm{HIC}})$
		$\displaystyle+\frac{1}{2}\tilde{S}_{\rm{HIC},2}(T,n_{B})\left(1-\frac{Y_{Q}}{Y% _{Q}^{\rm{HIC}}}\right)^{2}.$		(63)

Now, we can take the temperature derivative

\frac{\partial\tilde{S}(T,n_{B},Y_{Q})}{\partial T}\Bigg{|}_{T=0}=\frac{1}{n_{% B}}\frac{\partial s_{\rm{HIC}}(T,n_{B},Y_{Q})}{\partial T}\Bigg{|}_{T=\delta_{% \rm HIC}=0}+\frac{1}{2}\left(1-\frac{Y_{Q}}{Y_{Q}^{\rm{HIC}}}\right)^{2}\frac{% \partial^{3}\tilde{S}_{\rm{HIC},2}(T,n_{B},\delta_{\rm HIC})}{\partial T% \partial\delta_{\rm HIC}^{2}}\Bigg{|}_{T=\delta_{\rm HIC}=0},

(64)

where we can pull out $n_{B}$ and $Y_{Q}$ terms since we are taking this expression at constant $n_{B}$ and $Y_{Q}$ . Using Eq. (64), if we know the first term, $\frac{\partial s_{\rm{HIC}}(T,n_{B},Y_{Q})}{\partial T}\Big{|}_{T=\delta_{\rm HIC% }=0}$ , and the second term, $\frac{\partial^{3}\tilde{S}_{\rm{HIC},2}(T,n_{B})}{\partial T\partial\delta_{% \rm HIC}^{2}}\Big{|}_{T=\delta_{\rm HIC}=0}$ , we can calculate $\frac{\partial s}{\partial T}\big{|}_{T=0}$ along any slice of $Y_{Q}=$ const., which is what we need for the finite $T$ expansion of the pressure.

Since we require the information about $\frac{\partial s}{\partial T}\big{|}_{T=0,\vec{\mu}}$ along a grid of fixed $\vec{\mu}$ in Eq. (58), but the expansion that we have derived above is in terms of fixed $n_{B},Y_{Q}$ , i.e., $\frac{\partial s}{\partial T}\big{|}_{T=0,n_{B},Y_{Q}}$ , we must perform a mapping between the two quantities, i.e.,

\frac{\partial s}{\partial T}\Big{|}_{T=0,\vec{\mu}}\underset{?}{% \leftrightarrow}\frac{\partial s}{\partial T}\Big{|}_{T=0,n_{B},Y_{Q}},

(65)

. In Appendix A, we show explicitly the conversion between the these two thermodynamic ensembles and demonstrate that the correction term goes to zero at $T=0$ . Thus, in this work we argue that

\frac{\partial s}{\partial T}\Big{|}_{T=0,\vec{\mu}}\sim\frac{\partial s}{% \partial T}\Big{|}_{T=0,n_{B},Y_{Q}},

(66)

and ignore any correction terms that show up at finite temperatures.

II.3.3 Extracting $s/n_{B}~{}(T,n_{B},Y_{Q})$ from HIC data

It may be possible in future work to use HIC data to constrain these two terms: $\frac{\partial s_{\rm{HIC}}(T,n_{B},Y_{Q}=Y_{Q}^{\rm{HIC}})}{\partial T}$ and $\frac{\partial^{3}\tilde{S}_{\rm{HIC},2}(T,n_{B})}{\partial T\partial\delta_{% \rm HIC}^{2}}\big{|}_{\delta_{\rm HIC}=0}$ . For instance, at chemical freeze-out one can use a model for the equation of state of the hadron resonance gas (HRG) to extract the temperature and chemical potentials using identified particles yields, ratios, and fluctuations. Typically, an ideal HRG model is used with the entire particle list from the particle data group (PDG) Alba et al. (2020). However, other, more realistic models can be used instead, such as a van der Waals EOS Poberezhnyuk et al. (2019), lattice QCD directly (at least for fluctuations) Borsanyi et al. (2013, 2014), or other effective models. In fact, any framework that provides microscopic information about the densities of specific particle species could then be used to extract these temperatures and chemical potentials from experimental data given a specific beam energy $\sqrt{s_{NN}}$ and a colliding species that fixes $Y_{Q}^{\rm{HIC}}=Z/A$ . Of course, if the model does not have the correct degrees of freedom and interactions it may not fit the data well. Presumably various assumptions could be tested against the particle yields, ratios, and fluctuations to determine the best fit, although we are unaware of such a study at this time.

Regardless of the underlying model, the general procedure is the same. The EOS must be 4D in terms of $T,\mu_{B},\mu_{S},\mu_{Q}$ wherein $\mu_{S}$ and $\mu_{Q}$ are constrained by

	$\displaystyle\langle n_{S}\rangle$	$\displaystyle=$	$\displaystyle 0,$		(67)
	$\displaystyle\langle n_{Q}\rangle$	$\displaystyle=$	$\displaystyle Y_{Q}^{\rm{HIC}}\langle n_{B}\rangle,$		(68)

such that we can always determine

\mu_{S}\left(T,\mu_{B}\right),\quad\mu_{Q}\left(T,\mu_{B}\right),

(69)

from $T$ and $\mu_{B}$ . Next, a minimum $\chi^{2}$ fit is performed using data from specific ion-ion (A-A) collisions at a specific beam energy $\sqrt{s_{NN}}$ and centrality. The result of the minimum $\chi^{2}$ fit provides the extracted freeze-out values of $T^{\rm FO},\mu_{B}^{\rm FO}$ from a given $\sqrt{s_{NN}}$ . Next, with this pair of $T^{\rm FO},\mu_{B}^{\rm FO}$ one can calculate EOS properties for that specific $\sqrt{s_{NN}}$ and $Y_{Q}^{\rm{HIC}}$ . Central collisions are considered the best for such a study since they result in the largest system size and presumably the closest to the infinite volume approximation of the Grand Canonical Ensemble.

At very large $\sqrt{s_{NN}}$ , the nuclei are Lorentz contracted and are very thin along the beam direction, passing through each other nearly instantaneously. As a result, in this limit, there is no time for baryons (i.e. valence quarks) to be stopped within the collision, such that the global baryon number of the system is extremely small $\langle n_{B}\rangle\rightarrow 0$ . This is also the regime where the collision reaches the highest temperatures. However, at lower $\sqrt{s_{NN}}$ , the nuclei must be treated as 3D objects, as they are traveling more slowly and take longer to pass through each other, allowing for enough time to capture baryons. For these intermediate $\sqrt{s_{NN}}$ , $n_{B}$ is large and the initial temperature is lower. Finally, at very low $\sqrt{s_{NN}}$ , the nuclei may not pass through each other entirely, but rather stick together. In this regime, $n_{B}$ is smaller, and the initial temperature is very low (may be so low that the quark-gluon plasma phase is no longer reached). Thus, the systems created in heavy-ion collisions can vary significantly as a function of $\sqrt{s_{NN}}$ and extend across a wide range on the QCD phase diagram.

Ion	$Y_{Q}^{\rm{HIC}}$	Data available
${}^{8}_{16}O$	0.5	not yet
${}^{29}_{63}Cu$	0.46	Alver et al. (2011)
${}^{44}_{96}Ru$	0.458	preliminary
${}^{40}_{96}Zr$	0.42	preliminary
${}^{79}_{198}Au$	0.4	Abelev et al. (2009)
${}^{92}_{238}U$	0.38	Abdallah et al. (2023)

Table 1: List of ion species (with the format of

{}^{Z}_{A}ion

, where

Z

is the number of protons and

A

is the number of nucleons) ran at

\sqrt{s_{NN}}=200

GeV at RHIC and their corresponding electric charge fraction.

Within heavy-ion collisions, we can also consider other possibilities to further refine the EOS that is studied. As previously mentioned, different ion species correspond to different $Y_{Q}^{\rm{HIC}}$ . In fact, a wide range of $Y_{Q}^{\rm{HIC}}$ values have already been run at $\sqrt{s_{NN}}=200$ GeV, where the baryon chemical potential is approximately $\mu_{B}\sim 20$ MeV. The different ions and their corresponding $Y_{Q}^{\rm{HIC}}$ are shown in Table 1. Additionally, we can also compare different centrality classes or rapidity windows, which also affect the temperatures and chemical potentials reached at freeze-out. With the upcoming FAIR fixed target experiment at GSI, there is an opportunity to run a variety of species with different $Y_{Q}^{\rm{HIC}}$ , EOS information can be extracted across a range of densities, temperatures, and $\delta_{\rm HIC}$ that would then provide constraints for Eq. (64).

One key difference between heavy-ion collisions and astrophysical scenarios is that in heavy-ion collision net-strangeness is conserved and, since the initial state has no strangeness, strangeness neutrality is enforced. Note that does not imply zero strangeness but just that the number of strange particles and anti-particles are equal to each other $N_{s}=N_{\bar{s}}$ . In contrast, in neutron star mergers there is no such conservation, and the star may violate strangeness neutrality. One potential way to avoid this issue is to only study yields, ratios, and fluctuations of light particles such as pions, protons, or light nuclei such as deuterons. However, one would have to consider how to map isentropes extracted using the strangeness neutrality constraint in Eq. 67 into isentropes applicable for neutron star mergers that would not have this constraint. We also leave that study for a future work.

II.4 Obtaining other thermodynamic quantities from $p(T,\vec{\mu})$

Once the finite temperature pressure is obtained, the other thermodynamic quantities can easily be calculated using:

$\displaystyle s(T,\vec{\mu})$	$\displaystyle=$	$\displaystyle\frac{dp}{dT}\biggr{\rvert}_{\vec{\mu}},$	(70)
$\displaystyle n_{B}(T,\vec{\mu})$	$\displaystyle=$	$\displaystyle\frac{\partial p}{\partial\mu_{B}}\biggr{\rvert}_{T,\mu_{Q}},$	(71)
$\displaystyle n_{Q}(T,\vec{\mu})$	$\displaystyle=$	$\displaystyle\frac{\partial p}{\partial\mu_{Q}}\biggr{\rvert}_{T,\mu_{B}},$	(72)
$\displaystyle Y_{Q}(T,\vec{\mu})$	$\displaystyle=$	$\displaystyle\frac{n_{Q}}{n_{B}},$	(73)
$\displaystyle\varepsilon$	$\displaystyle=$	$\displaystyle-p+sT+n_{B}\mu_{B}+n_{Q}\mu_{Q},$	(74)

where the last equation comes from the Gibbs-Durhem relation in Eq. (52) at finite temperature. Eqs. (70-72) can also be rewritten in terms of our series expansion, as shown in Appendix B. There, we show that the finite $T$ expansion requires a total $4$ coefficients in $\mathcal{O}(T^{2})$ and a total of $7$ coefficients in $\mathcal{O}(T^{3})$ to obtain all the necessary thermodynamic observables $\left\{s,n_{B},n_{Q},Y_{Q},\varepsilon\right\}$ analytically. The last two Eqs. (73-74) are straightforward to calculate once $\left\{s,n_{B},n_{Q}\right\}$ are obtained.

In order to extract the correct $Y_{Q}$ -dependent behavior, other derivatives are required. Specifically, to obtain $n_{B}(T,\vec{\mu})$ and $Y_{Q}(T,\vec{\mu})$ . In Appendix C we work through the calculations such that $\left\{s,n_{B},n_{Q}\right\}$ can be calculated analytically. We find that obtaining the entropy dependence on $Y_{Q}$ is very straightforward and only requires higher-order derivatives of the entropy with respect to the temperature. However, $n_{B}$ and $Y_{Q}$ would require a number of new derivatives with respect to $\mu_{B}$ or $\mu_{Q}$ . Considering the potential error in determining these fourth-order (or even high-order) derivatives, we instead advocate for simply taking the first-order derivatives in Eq. (71-72) numerically across the finite $T$ version of the EOS, once it is generated across a fixed table and using Eq. (73) to calculate $Y_{Q}$ .

II.5 Crust

While we have motivated this work with neutron star mergers in mind, up until this point the expansion scheme is generic and could be used in supernova explosions and nuclear experiments (although the range of $Y_{Q}$ values would be closer to SNM in some regions/stages of supernovae and heavy-ion collisions). However, the limit of $n_{B}\rightarrow 0$ in stars has significant differences from that in heavy-ion collisions. Thus, the EOS in this limit must be specific to the type of system that is being modeled. In Ref. Yao et al. (2023), a discussion on the method to attach an appropriate EOS for heavy-ion collisions when $n_{B}\rightarrow 0$ is discussed. Here, we focus on neutron star mergers.

The crust of a neutron star contains multiple different layers that depend strongly on $n_{B}$ Lattimer and Prakash (2001); Chamel and Haensel (2008). The outermost layers of a neutron star have an EOS that is dominated by a degenerate Fermi gas of electrons, with a small contribution to the energy density from nuclei. In this regime, the number density of protons, $n_{p}$ , within the nuclei is exactly equal to the number density of electrons, $n_{e}$ , in the system, i.e., $n_{p}=n_{e}$ , to ensure that the star is stable. As $n_{B}$ increases, the nuclei that minimize the energy become more and more neutron rich until eventually the neutron drip line is reached. Beyond the neutron drip line, there is a superfluid state of free neutrons and nuclei, then as $n_{B}$ continues to increase the nuclei are squeezed and form pasta phases. Eventually, nuclei are no longer stable at very high $n_{B}\sim n_{\rm sat}$ and one anticipates that N-body interactions between nucleons dominate the EOS.

Numerical relativity simulations requires the microscopic information from all of these rich phases of matter that take place from the outer core up to nearly $n_{\rm sat}$ . Thus, in this work, we match a microscopic EOS for the crust up to $n_{B}=n_{\rm sw}$ , where $n_{\rm sw}$ is the switching density given in units of $n_{\rm sat}$ , and at higher $n_{B}$ we use a functional form of the EOS that can be changed. In our framework, $n_{\rm sw}<1n_{\rm sat}$ because we allow the user to input the symmetry energy coefficients themselves. If one were to take $n_{\rm sw}\geq 1n_{\rm sat}$ , then the symmetry energy coefficients would already be defined since they are calculated at $n_{\rm sat}$ . Later in Sec. III, we discuss the practical implementation of a microscopic crust and its matching to the functional form of the EOS.

We should note the exact details related to the crust can modify the mass-radius of a neutron star by about $~{}10\%$ Grill et al. (2014); Pais et al. (2016); Ferreira and Providência (2020a, b). However, the inclusion of a crust is necessary to obtain the correct shape of the mass-radius sequence, such that it bends toward large radii at low masses. In contrast, compact objects that do not include a crust reach low radii for low masses, mimicking quark stars.

III Numerical Implementation

In order to construct the EOS across the full range of $\left\{T,n_{B},Y_{Q}\right\}$ required for numerical relativity simulations, we need to apply different types of EOS in different regions. Additionally, the construction of the 3D EOS will require different types of expansions across $Y_{Q}$ , across $T$ , and across $Y_{Q}$ and $T$ simultaneously. In this section, we discuss how we obtain the EOS in these different regimes numerically and also how we match the EOS in the different regimes to each other. Furthermore, we check the numerical accuracy of our EOS as it is being expanded in both $T$ and $Y_{Q}$ .

In Fig. 1 we sketch out the EOS models used in different regimes of the phase diagram as well as show the different expansion schemes. The finite $T$ regime comes from the finite $T$ expansion and we estimate the errors from this expansion in Fig. 1. At $T=0$ , our EOS varies across the $n_{B}$ and $Y_{Q}$ plane. At low $n_{B}$ , we use a microscopic crust based on a table that is an input to the code. At intermediate densities we interpolate the crust and switch to an interpolated version of the crust instead of the table. At high $n_{B}$ , we make use of a parametric form of the EoS that presents structure in $c_{s}^{2}$ , in agreement with different high-energy data. The behavior of this structure at $\beta$ -equilibrium is propagated across different $Y_{Q}$ using the symmetry energy expansion.

Refer to caption — Figure 1: Schematic diagram of the EOS and expansion schemes used in the different regions of the phase diagram. The EOS must cover a wide range in temperatures, $T$ , baryon number densities, $n_{B}$ , and charge fractions, $Y_{Q}$ . Pure neutron matter (PNM) is equivalent to $Y_{Q}=0$ , symmetric nuclear matter (SNM) is equivalent to $Y_{Q}=0.5$ , and a $\beta$ -equilibrated, cold neutron star is asymmetric nuclear matter (ANM) that falls between these two limits. Typical heavy-ion collisions (HIC) have values of $Y_{Q}\sim 0.38-0.5$ . The predicted numerical error for our $T$ expansion is shown vs $n_{B}$ . The symmetry energy expansion is performed around $Y_{Q}^{\rm HIC}$ in the direction of $Y_{Q}\rightarrow 0$ and the finite $T$ expansion is performed around $T=0$ in the direction of finite $T$ . At low densities we use a 3D crust EOS. For $T\leq 50$ MeV (green region) and $50\leq T\leq 100$ MeV (yellow region) the difference between the expanded pressure and the reference RMF EoS on the $(n_{B},Y_{Q})$ plane is below 20% and 40%, respectively. Note that the error bounds reflect the compounded error in the expansion when switching from a thermodynamical basis in $(\mu_{B},\mu_{Q})$ to $(n_{B},Y_{Q})$ .

In order to ensure that the EOS is thermodynamically consistent, we carefully match the different regimes of the EOS. In the next section we lay out all the steps in our algorithm and also define the free parameters and inputs needed in our framework.

III.1 Numerical algorithm and flow chart

In this section, we outline our numerical algorithm in order to produce a 3D table that is thermodynamically consistent for use in numerical relativity codes. The summary of our algorithm is shown in Fig. 2 together with the input parameters. The following sections provide more thorough details of the steps described here.

Our numerical scheme begins with a $\beta$ -equilibrated, cold neutron star EOS. At $T=0$ along $\beta$ -equilibrium, a microscopic crust is taken from $n_{B}=0$ up to $n_{\rm sw}$ where $n_{\rm sw}<n_{\rm sat}$ .

At high $n_{B}$ , we have two options for the core EOS. A user could input a table of $c_{s}^{2}(n_{B})$ at $\beta$ -equilibrium, which the code would then match to the crust, or the user could build in structure by hand within the EOS (this is the default). For the option of built-in structures in $c_{s}^{2}$ , it is possible to incorporate multiple different structures at specific points in $n_{B}$ , i.e., $n_{i},n_{i+1},\dots$ up until the maximum baryon density $n_{B}^{\rm max}$ considered for numerical relativity. Examples of these structures in $c_{s}^{2}$ are bumps, jumps, kinks, plateaus, steep rises or falls, which follow from a variety of physical motivations (see e.g., Tan et al. (2020, 2022a, 2022b); Mroczek et al. (2023)). In between $n_{\rm sw}<n_{B}<n_{i}$ , we take the crust EOS and interpolate its speed of sound, $c_{s}^{2}$ , but do not use its corresponding $Y_{Q}^{\beta,\rm crust}(n_{B})$ in this regime (instead $Y_{Q}^{\beta,\rm sym}(n_{B})$ comes from the symmetry energy expansion, as shown in Eq. (46)). Alternatively, a user can also input a table of $Y_{Q}^{\beta,\rm tab}(n_{B})$ that could be used. However, we caution that a poor choice in $Y_{Q}^{\beta,\rm tab}(n_{B})$ combined with specific symmetry energy coefficients may not correctly obtain saturation properties. The default is to obtain the $Y_{Q}^{\beta,\rm sym}(n_{B})$ from the symmetry energy expansion, but if a user provides a table, it will be matched to the $Y_{Q}^{\beta,\rm crust}(n_{B})$ in the same way as discussed below.

If we were to take the $Y_{Q}^{\beta,\rm crust}(n_{B})$ from the crust up to $n_{\rm sat}$ , we would be restricted to specific symmetry energy coefficients. However, by matching $Y_{Q}(n_{B})$ below $n_{\rm sat}$ , the symmetry energy coefficients can be varied, although we always check if saturation properties are maintained for SNM (as was was done in Ref. Yao et al. (2023)). At $n_{\rm sw}$ the $Y_{Q}^{\beta,\rm crust}(n_{B})$ is matched to the $Y_{Q}^{\beta,\rm sym}(n_{B})$ using a hyperbolic tangent switching function. For the $\beta$ -equilibrated EOS, the crust is matched to the interpolated (or user provided) EOS at $n_{\rm sw}$ using a $\tanh$ in $c_{s}^{2}(n_{B})$ that allows for a smooth transition between the microscopic EOS and the interpolated version. Then we can integrate $c_{s}^{2}(n_{B})$ using a second-order Runge Kutta to obtain $p(\varepsilon)$ , and from there use thermodynamic properties to obtain all necessary thermodynamic state variables. Since the matching at $n_{\rm sw}$ is based on the same EOS (one is a table and the other is interpolated) this matching should be extremely smooth. We then output the $\beta$ -equilibrated EOS.

Once our new $\beta$ -equilibrated EOS and $Y_{Q}(n_{B})$ is obtained, we can then use the symmetry energy expansion in Eq. (45) to expand to different slices in $Y_{Q}^{\rm const}$ . Following the expansion across $Y_{Q}$ , we calculate $c_{s}^{2}(n_{B},Y_{Q})_{\rm expn}$ across the 2D $T=0$ plane. Then, we return to our original 2D crust and use that table directly up to $n_{B}\leq n_{\rm sw}$ and switch to $c_{s}^{2}(n_{B},Y_{Q})_{\rm expn}$ for $n_{B}>n_{\rm sw}$ . Along slices of $Y_{Q}=$ const., we can integrate $c_{s}^{2}(n_{B},Y_{Q})_{\rm expn}$ upwards from $n_{\rm sw}$ using the second-order Runge Kutta and the thermodynamic information at $n_{\rm sw}$ from the crust table as initial conditions for the Runge Kutta. Using this method, we obtain $p(n_{B},Y_{Q})$ and $\varepsilon(n_{B},Y_{Q})$ and from there all other necessary thermodynamic observables can be calculated. From this point we can check saturation properties, and we output if the EOS fits within these known constraints or not.

Our final step is the development of the finite $T$ EOS. To do so we must reformat our EOS along a fixed grid in $\vec{\mu}$ space using a 2D interpolation. Then, we can calculate $p(T,\vec{\mu})$ and $s(T,\vec{\mu})$ analytically, after which we can calculate $n_{B}(T,\vec{\mu})$ and $n_{Q}(T,\vec{\mu})$ numerically. Then, we add the contribution of leptons back into the EOS, ensuring local charge neutrality at each point in the grid. We can then interpolate and reformat the grid in terms of $T,n_{B},Y_{Q}$ and calculate $c_{s}^{2}$ across this grid. Finally, we are able to convert our units back to astrophysical units and covert to the format required for numerical relativity simulations, which is then recorded as an output for the user.

III.2 Crust

In this section, we detail the changes and calculations that need to be performed on the crust EOS in order to make it compatible with the core EOS. These changes primarily consist of unit conversions but also include the calculations of certain thermodynamic variables. Additionally, we need to find the $\beta$ -equilibrated EOS, since this is typically not automatically provided from the crust.

The crust we use here is SLy4 Schneider et al. (2017) that is available online at SLy on the StellarCollapse website. In principle, any other crust could be used, although it would need to be in the same format as StellarCollapse. In Appendix D, we discuss the format of SLy4 and unit conversion for the various thermodynamic variables relevant to this work.

III.2.1 Finding $\beta$ -equilibrium

Cold, isolated neutron stars reach a state of $\beta$ -equilibrium (see Sec. II.1) wherein a relationship between the chemical potentials (assuming neutrinos are not trapped) forms, shown in Eq. (16). The 3D table provided does not supply the EOS at $\beta$ -equilibrium. However, we can reconstruct it using the relation in Eq. (16). Our algorithm for finding $\beta$ -equilibrium calculates the chemical potential that is a shift out-of- $\beta$ -equilibrium, $\mu_{\Delta\beta}$ at each point in the table, i.e.,

\mu_{\Delta\beta}=\mu_{n}-(\mu_{p}+\mu_{e}),

(75)

where the $\mu_{\Delta\beta}=0$ at $\beta$ -equilibrium. Then we employ a 1D interpolation and root-finder at a fixed $n_{B}$ at $T=0$ to find the corresponding $Y_{Q}(n_{B})$ where $\mu_{\Delta\beta}=0$ . Repeating this method across the entire range of $n_{B}$ we reconstruct the full $\beta$ -equilibrium $Y_{Q}(n_{B})$ and from this point can interpolate the table to obtain the remaining thermodynamic quantities at $\beta$ -equilibrium.

III.3 Crust-core matching

In this work, we want to ensure that the matching between the crust and the core is thermodynamically consistent, such that it does not introduce phase transitions, nor leads to significant numerical error within numerical relativity simulations. Based on our algorithm in Fig. 2, we first perform a 1D matching of the crust at $\beta$ -equilibrium to the high $n_{B}$ region of the EOS. Then we smoothly switch from our 2D crust table to the 2D high $n_{B}$ (core) EOS that was reconstructed using the symmetry energy expansion. In order to perform the matching (for the symmetry energy expansion), we pick a point to switch from the crust to the high $n_{B}$ EOS. Here we match at a fixed point of $n_{B}^{sw}$ that is the same point for both the 1D matching at $\beta$ -equilibrium as well as the point in 2D where we switch between the crust and the expanded EOS along slices of $Y_{Q}=$ const.

In past work Most et al. (2019); Most and Raithel (2021), thermodynamic consistency from the crust-core transition was not fully ensured, but here we attempt to address that issue. As discussed previously in Sec. III.1, we have three different sections of our $T=0$ EOS: the crust (table), the interpolated crust (int) wherein $Y_{Q}(n_{B})$ at $\beta$ -equilibrium is replaced with Eq. (46)), and structure in $c_{s}^{2}$ that appears at high $n_{B}=n_{\rm strc}$ (strc). Then, both the 1D and 2D matching only occur at the boundary between the table EOS and the int EOS, since the connection between the int EOS and the strc EOS is made at the level of $c_{s}^{2}(n_{B})$ , such that one can simply integrate upward from a lower $n_{B}$ to obtain $p(\varepsilon)$ and no unintentional phase transitions are introduced in the process.

For the 1D problem along $\beta$ -equilibrium, a variety of approaches have been followed in the past. One of them is to use a $\tanh$ to avoid unintentionally introducing phase transitions (e.g., Tan et al. (2022a)), which is what we use here. Generically, we can achieve a smooth transition between two 1D functions that vary with $n_{B}$ by using:

f_{\rm tanh}(n_{B})\equiv S(n_{B})f_{\rm int}(n_{B})+\left[1-S(n_{B})\right]f_% {\rm crust}(n_{B})\,,

(76)

where $S(n_{B})$ can be any smoothing function but here we use a $\tanh$ :

S(n_{B}):=0.5+0.5\tanh\left[\frac{(n_{B}[n_{\rm sat}]-n_{\rm sw})}{\Delta n_{B% }}\right]\,,

(77)

where $\Delta n_{B}$ determines the width of the smoothing region, the baryon density is given in terms of $[n_{\rm sat}]$ , and $n_{\rm sw}$ is the point in $n_{B}[n_{\rm sat}]$ where the matching happens. The $c^{2}_{s}$ is matched at $n_{B}=0.5n_{\rm sat}$ and $\Delta n_{B}=0.1$ .

Then, for the 2D connection between the crust table and the expanded 2D EOS (expn) using the symmetry energy expansion, we use a different approach. We simply switch from the $c_{s}^{2}(n_{B},Y_{Q})_{\rm table}$ at $n_{\rm sw}$ directly to the $c_{s}^{2}(n_{B},Y_{Q})_{\rm expn}$ . Then, we can once again apply the second-order Runge Kutta and integrate starting from $n_{\rm sw}$ up to higher $n_{B}$ to obtain the thermodynamic variables needed for the full EOS. This approach provides a fully smooth EOS when switching from the table to the expanded EOS.

III.3.1 Adding structure into $c_{s}^{2}$

One of our goals here is to add in structure into $c_{s}^{2}(n_{B})$ at $\beta$ -equilibrium and then use the symmetry energy expansion to allow that structure to propagate to different slices of $Y_{Q}$ . There are many ways to add structure in $c_{s}^{2}(n_{B})$ such as using Gaussian Processes Landry and Essick (2019); Han et al. (2021); Legred et al. (2021); Landry et al. (2020); Essick et al. (2020); Mroczek et al. (2023), adding in phase transitions through the constant speed of sound method Alford et al. (2013), or inserting functional forms using specific types of equations (e.g., Gaussians, polynomials, etc) Tews et al. (2018). The code is flexible enough to change the details of the structure at high densities or it would be straightforward enough to input an external $c_{s}^{2}(n_{B})$ that could be used instead at some chosen $n_{i}$ .

III.3.2 EOS reconstruction from a second-order Runge Kutta

Now that the description of $c_{s}^{2}(n_{B})$ is established for $\beta$ -equilibrium, we can integrate upwards from $n_{\rm sw}$ to obtain the remaining thermodynamic variables. Previous work Tews et al. (2018); Tan et al. (2020, 2022a) used an Eulerian integration that is sufficient only for studying properties at $\beta$ -equilibrium. However, because we will need to expand to other $Y_{Q}$ slices and calculate a number of thermodynamic derivatives, we have decided to use a second-order Runge-Kutta instead to improve numerical accuracy.

We have devised the following scheme to obtain $p$ , $\varepsilon$ , and $\mu_{B}$ from a starting point $n_{0}$ where the step size is $\Delta n_{B}$ . A Runge Kutta of second-order requires the initial point at $n_{B}=n_{0}$ , the mid point at $n_{B}=n_{0}+\Delta n_{B}/2$ , and the slope at the mid point $n_{B}=n_{0}+\Delta n_{B}/2$ . Thus, for a generic function $f(n_{B})$ one calculates the derivative at $n_{0}$

\displaystyle f^{\prime}_{0}(n_{0})=\frac{df_{0}(n_{0})}{dn_{B}},

(78)

using that to calculate the function at the mid point

f_{1}\left(n_{0}+\frac{\Delta n_{B}}{2}\right)=f_{0}(n_{0})+f^{\prime}_{0}(n_{% 0})\frac{\Delta n_{B}}{2},

(79)

and the slope at the midpoint

f^{\prime}_{1}\left(n_{0}+\frac{\Delta n_{B}}{2}\right)=\frac{df_{1}\left(n_{0% }+\frac{\Delta n_{B}}{2}\right)}{d\left(n_{0}+\frac{\Delta n_{B}}{2}\right)},

(80)

that can be used to provide our final value for the full step:

f_{2}\left(n_{0}+\Delta n_{B}\right)=f_{0}(n_{0})+f^{\prime}_{1}\left(n_{0}+% \frac{\Delta n_{B}}{2}\right)\Delta n_{B}.

(81)

Since we assume $\beta$ -equilibrium and charge neutrality, Eqs. 49, 51, and 56 simplify to

\varepsilon=-p+n_{B}\mu_{B},

(82)

\mu_{B}=\frac{d\varepsilon}{dn_{B}},

(83)

c_{s}^{2}=\frac{dp}{d\varepsilon},

(84)

that help us to calculate the slopes.

For our second-order Runge Kutta we are integrating over baryon density such that:

	$\displaystyle n_{B,1}$	$\displaystyle=$	$\displaystyle n_{B,0}+\frac{\Delta n_{B}}{2},$		(85)
	$\displaystyle n_{B,2}$	$\displaystyle=$	$\displaystyle n_{B,0}+\Delta n_{B},$		(86)

such that $n_{B,1}$ is our mid point and $n_{B,2}$ is the full step. The energy density at the mid point is then

	$\displaystyle\varepsilon_{1}$	$\displaystyle=$	$\displaystyle\varepsilon_{0}+\left(\frac{d\varepsilon}{dn_{B}}\right)_{0}\frac% {\Delta n_{B}}{2}$		(87)
		$\displaystyle=$	$\displaystyle\varepsilon_{0}+\mu_{B,0}\frac{\Delta n_{B}}{2},$		(87)

where

\mu_{B,i}=\frac{\varepsilon_{i}+p_{i}}{n_{B,i}}.

(88)

Then we also require the information about the pressure where

$\displaystyle p_{1}$	$\displaystyle=$	$\displaystyle p_{0}+\left(\frac{dp}{dn_{B}}\right)_{0}\frac{\Delta n_{B}}{2}$	(89)
	$\displaystyle=$	$\displaystyle p_{0}+\left(\frac{dp}{d\varepsilon}\frac{d\varepsilon}{dn_{B}}% \right)_{0}\frac{\Delta n_{B}}{2}$
	$\displaystyle=$	$\displaystyle p_{0}+c^{2}_{s}(n_{B,0})\mu_{B,0}\frac{\Delta n_{B}}{2},$

where

c_{s}^{2}(n_{B,0})=\frac{dp_{0}}{d\varepsilon_{0}}.

(91)

Putting this all together, at the mid point we can calculate the slopes:

	$\displaystyle\mu_{B,1}$	$\displaystyle=$	$\displaystyle\frac{e_{1}+p_{1}}{n_{B,1}},$		(92)
	$\displaystyle c_{s}^{2}(n_{B,1})$	$\displaystyle=$	$\displaystyle\frac{dp_{1}}{d\varepsilon_{1}}.$		(93)

Then, this allows up to construct the full time step for the energy, the pressure, the chemical potential, and speed of sound squared:

$\displaystyle\varepsilon_{2}$	$\displaystyle=$	$\displaystyle\varepsilon_{0}+\mu_{B,1}\cdot\Delta n_{B},$	(94)
$\displaystyle p_{2}$	$\displaystyle=$	$\displaystyle p_{0}+c^{2}_{s}(n_{B,1})\mu_{B,1}\cdot\Delta n_{B},$	(95)
$\displaystyle\mu_{B,2}$	$\displaystyle=$	$\displaystyle\frac{e_{2}+p_{2}}{n_{B,2}},$	(96)
$\displaystyle c_{s}^{2}(n_{B,2})$	$\displaystyle=$	$\displaystyle\frac{dp_{2}}{d\varepsilon_{2}},$	(97)

such that we have recovered all necessary thermodynamic quantities at $\beta$ -equilibrium. To ensure numerical accuracy, we use steps in baryon density of $\Delta n_{B}=0.001n_{\rm sat}$ .

For the 2D EOS, a similar process must be performed. Once again, we integrate over baryon density along a slice of $Y_{Q}=$ const. such that

	$\displaystyle n_{B,1}$	$\displaystyle=$	$\displaystyle n_{B,0}+\frac{\Delta n_{B}}{2},$		(98)
	$\displaystyle n_{B,2}$	$\displaystyle=$	$\displaystyle n_{B,0}+\Delta n_{B},$		(99)

and we can write the charge density as

	$\displaystyle n_{Q,1}$	$\displaystyle=$	$\displaystyle n_{Q,0}+\frac{\Delta n_{Q}}{2}=n_{Q,0}+Y_{Q}\frac{\Delta n_{B}}{% 2},$		(100)
	$\displaystyle n_{Q,2}$	$\displaystyle=$	$\displaystyle n_{Q,0}+Y_{Q}\Delta n_{B}.$		(101)

Since $Y_{Q}$ and $\Delta n_{B}$ are both constants, then $\Delta n_{Q}=Y_{Q}\Delta n_{B}=$ const.

The energy density at the mid point is then

\displaystyle\varepsilon_{1}

\displaystyle=

\displaystyle\varepsilon_{0}+\left(\frac{\partial\varepsilon}{\partial n_{B}}% \right)\Big{|}_{Y_{Q},0}\frac{\Delta n_{B}}{2},

(102)

where we no longer have $\mu_{B}\neq\frac{d\varepsilon}{dn_{B}}|_{Y_{Q}}$ but rather $\frac{d\varepsilon}{dn_{B}}|_{Y_{Q}}$ has contributions both from $\mu_{B}$ and $\mu_{Q}$ since both $n_{Q}$ and $n_{B}$ vary along slices of $Y_{Q}=$ const. (see Eq. 51). Thus, we simply take the derivative $\frac{d\varepsilon}{dn_{B}}|_{Y_{Q}}$ for the slope and do not substitute in the chemical potential as we did previously for the 1D case.

Then we also require the information about the pressure where

$\displaystyle p_{1}$	$\displaystyle=$	$\displaystyle p_{0}+\left(\frac{dp}{dn_{B}}\right)_{0}\frac{\Delta n_{B}}{2}$	(103)
	$\displaystyle=$	$\displaystyle p_{0}+\left(\frac{\partial p}{\partial\varepsilon}\frac{\partial% \varepsilon}{\partial n_{B}}\right)_{0}\frac{\Delta n_{B}}{2}$
	$\displaystyle=$	$\displaystyle p_{0}+c^{2}_{s}(n_{B,0})\left(\frac{\partial\varepsilon}{% \partial n_{B}}\right)\Big{\|}_{Y_{Q},0}\frac{\Delta n_{B}}{2}.$

The speed of sound squared can be calculated directly for $Y_{Q}=$ const.

c_{s}^{2}(n_{B,0})=\frac{\partial p_{0}}{\partial\varepsilon_{0}}\Big{|}_{Y_{Q% }=\rm const},

(104)

since both $\Delta n_{B}=$ const. and $\Delta n_{Q}=$ const. Putting this all together, at the mid point we can calculate the slope

\displaystyle c_{s}^{2}(n_{B,1})

\displaystyle=

\displaystyle\frac{\partial p_{1}}{\partial\varepsilon_{1}}\Big{|}_{Y_{Q}=\rm const}.

(105)

Then, this allows up to construct the full time step for the energy, the pressure, and speed of sound squared

$\displaystyle\varepsilon_{2}$	$\displaystyle=$	$\displaystyle\varepsilon_{0}+\left(\frac{\partial\varepsilon}{\partial n_{B}}% \right)\Big{\|}_{Y_{Q},1}\cdot\Delta n_{B},$	(106)
$\displaystyle p_{2}$	$\displaystyle=$	$\displaystyle p_{0}+c^{2}_{s}(n_{B,1})\left(\frac{\partial\varepsilon}{% \partial n_{B}}\right)\Big{\|}_{Y_{Q},1}\cdot\Delta n_{B},$	(107)
$\displaystyle c_{s}^{2}(n_{B,2})$	$\displaystyle=$	$\displaystyle\frac{dp_{2}}{d\varepsilon_{2}}\Big{\|}_{Y_{Q}=\rm const},$	(108)

such that we have recovered most necessary thermodynamic quantities across the 2D EOS.

Once our full 2D EOS is calculated across the necessary $n_{B},Y_{Q}$ , we can then interpolate our 2D EOS such that we can take derivatives along trajectories of $n_{X}$ while holding $n_{Y\neq X}=$ const. and obtain $\mu_{B}$ and $\mu_{Q}$ . Note that we only require the Runge Kutta approach for $n_{B}>n_{\rm sw}$ since below that point we are using the table directly.

III.3.3 Interpolated EOS and saturation properties

After the symmetry energy expansion part of our algorithm is complete and we have reconstructed the entire 2D EOS from crust to core at $T=0$ , we can then determine if we have obtained reasonable saturation properties. From low-energy nuclear physics experiments, we know certain properties of nuclear matter at $n_{\rm sat}$ (see Kumar et al. (2023) for a discussion on the experimental and astrophysical constraints). Here we can check that our EOS for SNM ( $Y_{Q}=0.5$ ) reproduces saturation properties. The first property at $n_{\rm sat}$ that we check is the binding energy

B=\frac{\varepsilon_{\rm SNM}}{n_{\rm sat}}-m_{N},

(109)

that presents a minimum at $n_{B}=n_{\rm sat}$ , where one expects $n_{\rm sat}\sim 0.14-0.18\;fm^{-3}$ . Then, the binding energy at $n_{\rm sat}$ is typically around $B=-16~{}\rm{MeV}$ with a range of $-14$ to $-18$ MeV. One can also calculate the compressibility at $n_{\rm sat}$ , which is defined as

K=9\frac{\partial p}{\partial n_{B}}\Bigg{|}_{n_{B}=n_{\rm sat}},

(110)

wherein requiring $200<K<300$ MeV is a reasonable constraint. Our code checks if a given EOS passes these constraints and prints off a warning if one of these is out-of-bounds but then proceeds to the next steps.

III.4 Finite T expansion

Now that we have constructed the 2D EOS at T=0, we can begin our work on the finite $T$ expansion. There are a few questions that we wish to explore here:

•

How significant are the terms at $\mathcal{O}(T^{2})$ and $\mathcal{O}(T^{3})$ in the expansion?
•

Are the terms well-behaved across $Y_{Q}$ ? Does our expansion in $s/n_{B}$ accurately describe finite $T$ behavior along different $Y_{Q}$ slices?
•

How accurate is the finite $T$ expansion at $\mathcal{O}(T^{2})$ and $\mathcal{O}(T^{3})$ ?
•

Does our switch between a grid in fixed $\vec{\mu}$ to fixed $Y_{Q},n_{B}$ introduce significant numerical error?

Note that all these questions are numerical. Thus, we study them in the context of a single model – the relativistic mean field (RMF) model Alford et al. (2022), although our approach is generic and can work for any model not containing first-order phase transitions. This model describes protons and neutrons coupled to $\sigma$ , $\omega$ , and $\rho$ mesons. We choose the following values for the couplings: $g_{\sigma}=8.3965$ , $g_{\omega}=10.1845$ , and $g_{\rho}=10.9176$ and the remaining free parameters $b=0.00438046$ , $c=-0.00359399$ , and $b_{1}=5.18964$ (see Ref. Alford et al. (2022) for a description of the model Lagrangian and free parameters). These parameters are a representative sample of a large set of constrained parameter values, which reproduce the properties of uniform PNM obtained from chiral effective field theory and satisfy basic astrophysical constraints; namely, the maximum mass of a stable, isolated, slowly-rotating neutron star is at least two solar masses Demorest et al. (2010); Arzoumanian et al. (2018); Cromartie et al. (2020); Fonseca et al. (2021); Agazie et al. (2023) and the radii predicted for 1.4 and 2.1 solar-mass stars is within the one-sigma posterior obtained from NICER data Miller et al. (2019, 2021); Riley et al. (2019, 2021).

We start by picking a fixed point in $Y_{Q}=Y_{Q}^{\rm HIC}$ to study the finite $T$ expansion. In this limit, we demonstrate the accuracy of the finite $T$ expansion at $\mathcal{O}(T^{2})$ and $\mathcal{O}(T^{3})$ . Then, we calculate the $s/n_{B}$ expansion and demonstrate its accuracy across $Y_{Q}$ . Finally, we quantify the numerical error introduced by switching between a grid in fixed $\vec{\mu}$ to fixed $Y_{Q},n_{B}$ .

III.4.1 Terms at $\mathcal{O}(T^{2})$ and $\mathcal{O}(T^{3})$

Using the RMF model, we can calculate the $\mathcal{O}(T^{2})$ term $\partial s/\partial T|_{T=0}$ and the $\mathcal{O}(T^{3})$ term $\partial^{2}s/\partial T^{2}|_{T=0}$ . In order to calculate accurate derivatives we require step sizes in the temperature of $\Delta T=1$ MeV and must generate tables in ranges of $T=0$ to $100$ MeV. To do so, we calculate the terms on a fixed grid in $\mu_{B}$ and $\mu_{Q}$ for all points at $T=0$ . Then, we can perform a root-finding method to reconstruct $\partial s/\partial T|_{T=0}$ and $\partial^{2}s/\partial T^{2}|_{T=0}$ along lines of $Y_{Q}=$ const.

In Fig. 3 we plot the second-order term $\partial s/\partial T|_{T=0}$ on the top and the third-order term $\partial^{2}s/\partial T^{2}|_{T=0}$ on the bottom for a fixed $Y_{Q}=0.4$ (the value corresponding to Au-Au heavy-ion collisions). The second-order term is positive definite, which is what one should expect, because at finite temperature entropy is non-zero. Thus, the derivative of $\partial s/\partial T$ at the limit of $T=0$ must be positive to ensure that as one goes from $T=0$ to any $T$ , entropy is increased. We find that the $\partial s/\partial T|_{T=0}$ term increases monotonically with $\mu_{B}$ and has a relatively smooth behavior.

When we study the third-order term $\partial^{2}s/\partial T^{2}|_{T=0}$ , we find that it is at least 5-orders of magnitude smaller than the second-order term for this specific $Y_{Q}=0.4$ . We find that the third-order term is nearly zero at large $\mu_{B}$ and increases with decreasing $\mu_{B}$ , although it is not necessarily monotonic, having some small oscillations at very large $\mu_{B}$ that are due to numerical error from taking a third-order derivative or difficulties calculating the entropy at very low $T$ in the RMF EOS. Since it has the opposite qualitative behavior as the second-order term, one should expect that the second-order term would play the largest role at low $\mu_{B}$ . However, for this specific realization of the EOS at $Y_{Q}=0.4$ , one would need to go to extremely large $T$ in order to see an influence of the third-order term. Thus, any potential numerical noise in the third-order term is of little relevance except at very high $T$ . Later on we quantify more precisely the influence of the third-order term at specific temperatures in Sec. IV and discuss potential numerical noise in this term.

In Fig. 4 we plot the second-order term $\partial s/\partial T|_{T=0}$ on the top and the third-order term $\partial^{2}s/\partial T^{2}|_{T=0}$ on the bottom but this time we compare their behavior across $Y_{Q}$ . Qualitatively, we find a very similar behavior for $\partial s/\partial T|_{T=0}$ across $Y_{Q}$ , but the overall magnitude changes. When comparing across slices in $Y_{Q}$ , we find that the $Y_{Q}\rightarrow 0$ limit has the smallest overall values of $\partial s/\partial T|_{T=0}$ whereas the magnitude of $\partial s/\partial T|_{T=0}$ increases with increasing $Y_{Q}$ . Thus, from the second-order term alone, we anticipate that SNM is more sensitive to $T$ dependence compared to PNM. In other words, heavy-ion collisions should be more sensitive to temperature effects than neutron star mergers.

The third-order term in Fig. 4 (bottom) is overall not positive definite, nor is it monotonic, unlike the second-order term that has both properties. As we saw previously for the case of $Y_{Q}=0.4$ , at low $\mu_{B}$ generally the $\partial^{2}s/\partial T^{2}|_{T=0}$ term is the largest, especially for small $Y_{Q}$ . Comparing Fig. 3 to Fig. 4, we find that for $Y_{Q}=0.4$ the $\partial^{2}s/\partial T^{2}|_{T=0}$ result is at least 2-orders of magnitude smaller than for small $Y_{Q}\lesssim 0.1$ . In fact, when plotting in Fig. 4, it appears that for $Y_{Q}$ close to SNM that the results are nearly 0. Comparing the second and third-order terms we can conclude that SNM is more strongly $T$ dependent and likely has the most straight-forward $T$ dependence (the third-order term is probably mostly negligible), at least in this RMF model. On the other hand, for PNM and ANM the system does not have as strong of an overall $T$ dependence, but some non-trivial behavior may show up due to the contribution of the third-order term, especially at low $\mu_{B}$ , high $T$ . However, further investigation must be made into the numerical accuracy of the third-order term to better understand if the wiggles that appear are due to numerical error or physics.

III.4.2 Proof-of-principles for HIC with $Y_{Q}^{\rm HIC}=0.4$

Now that we have an idea of what the $\mathcal{O}(T^{2})$ and $\mathcal{O}(T^{3})$ terms look like, we can use these terms to reconstruct the pressure at finite $T$ . As a proof of principle, we use $Y_{Q}^{\rm HIC}=0.4$ because both $\partial s/\partial T|_{T=0}(\mu_{B})$ and $\partial^{2}s/\partial T^{2}|_{T=0}(\mu_{B})$ appeared well-behaved in this limit in Fig. 3. Another motivation for choosing $Y_{Q}^{\rm HIC}=0.4$ as a starting point is that we also need to expand around $Y_{Q}^{\rm HIC}$ in order to obtain a $Y_{Q}$ dependence for $\partial s/\partial T|_{T=0}(Y_{Q})$ . We wanted to pick a $Z/A$ from heavy-ions collisions that is as close to typical $Y_{Q}$ values reached in neutron stars as possible, such that one has to expand the smallest distance possible in $\delta_{\rm HIC}$ .

In Fig. 5 (top panel) we show a comparison between the pressure $p$ in MeV/fm³ vs $\mu_{B}$ at $T=0$ vs $T=100$ MeV for the RMF EOS. Then, our finite $T$ expansion in Eq. (58) is shown including either just the second-order term (black solid line) or both the second- and third-order terms from Fig. 3 (green dotted line). We find that the inclusion of just the second-order term already provides an extremely accurate representation of the $T=100$ MeV EOS from $\mu_{B}=1000$ to $2000$ MeV. Only a very small deviation can be seen around $\mu_{B}\sim 1000$ MeV. When we include the third-order term we produce no visible difference in the results.

We can quantify the deviation from the original RMF EOS at $T=100$ MeV by calculating the following relative $\%$ error of deviation from the pressure:

\frac{|\Delta p|}{p_{\rm RMF}}\times 100\%=\frac{p_{\rm expn}-p_{\rm RMF}}{p_{% \rm RMF}}\times 100\%,

(111)

where $p_{\rm RMF}$ is directly from the RMF model and $p_{\rm expn}$ is from our finite $T$ expansion. The results of this comparison are shown in the middle panel in Fig. 5, where we find that for $\mu_{B}\gtrsim 1200$ MeV our error is only at the $\mathcal{O}(1\%)$ or even significantly less. However, for smaller $\mu_{B}$ we do find larger deviations from the original RMF table because there may be possible contributions from the liquid-gas phase transition. In this figure we also find oscillations in the relative $\%$ error, which we suspect arise from issues in the original RMF model in accurately calculating the entropy at low $T$ , which would in turn affect our pressure expansion.

Finally, in the bottom panel of Fig. 5 we study the contribution of the second-order term:

p_{2}\equiv\frac{\partial s}{\partial T}\big{|}_{T=0}T^{2},

(112)

vs the third-order term:

p_{3}\equiv\frac{\partial^{2}s}{\partial T^{2}}\big{|}_{T=0}T^{3},

(113)

by plotting the ratio $p_{3}/p_{2}$ vs $\mu_{B}$ MeV. For $Y_{Q}=0.4$ the $p_{3}$ contribution is heavily suppressed and even at its largest is 4 orders of magnitude smaller than $p_{2}$ such that it plays no role in the expansion. Thus, the small deviation that we find at low $\mu_{B}$ may either be fixed in future work by improved numerical accuracy of the entropy in the RMF model or by possibility exploring fourth-order terms.

III.4.3 $s/n_{B}$ expansion in $Y_{Q}$

Now that we have shown that our finite $T$ expansion works well in a specific limit of $Y_{Q}=0.4$ , we test our $s/n_{B}$ expansion to see if it allows us to properly capture the $Y_{Q}$ behavior of $\partial s/\partial T|_{T=0}$ . Looking at the $Y_{Q}$ dependence of $\partial s/\partial T|_{T=0}$ in Fig. 4, it seems that it should be relatively straightforward to capture this behavior.

As a first step, we plot $s/n_{B}$ versus the parameter $\delta_{\rm HIC}$ expanded around $Y_{Q}^{\rm HIC}$ for a fixed temperature $T=1$ MeV along lines of $n_{B}\;[n_{\rm sat}]=$ const. The result of $s/n_{B}$ versus $\delta_{\rm HIC}$ is shown in Fig. 6 in the top panel. We find that $s/n_{B}$ is extremely flat versus $\delta_{\rm HIC}$ for different $n_{B}$ slices, although at low $n_{B}$ some small amount of $\delta_{\rm HIC}$ dependence begins to appear. The small amount of $\delta_{\rm HIC}$ dependence indicates a clear decrease in $s/n_{B}$ as one increases $\delta_{\rm HIC}$ (and that change increases more quickly for larger $\delta_{\rm HIC}$ ). Thus, from this behavior we anticipate a negative second derivative of $\delta_{\rm HIC}^{2}$ such that it is more pronounced also at small $n_{B}$ .

Next, in Fig. 6 in the bottom panel we show the second derivative, $\partial^{2}(s/n_{B})/\partial\delta_{\rm HIC}^{2}|_{\delta_{\rm HIC}=0}$ , versus temperature $T$ . We find that $\partial^{2}(s/n_{B})/\partial\delta_{\rm HIC}^{2}|_{\delta_{\rm HIC}=0}$ is always negative, which implies that $s/n_{B}$ decreases with increasing $\delta_{\rm HIC}$ or, in other words, $s/n_{B}$ is larger (at a fixed $T,n_{B}$ ) for SNM than for PNM. This qualitative effect was indeed confirmed in the top panel of Fig. 6. Using $\partial^{2}(s/n_{B})/\partial\delta_{\rm HIC}^{2}|_{\delta_{\rm HIC}=0}$ we can quantify how much and how quickly $s/n_{B}$ decreases as one approaches more asymmetric nuclear matter. The other effect that we see in Fig. 6 (bottom) is that the $\partial^{2}(s/n_{B})/\partial\delta_{\rm HIC}^{2}|_{\delta_{\rm HIC}=0}$ term has a temperature dependence that strongly depends on $n_{B}$ , whereas we find at high $n_{B}$ there is little temperature dependence of $\partial^{2}(s/n_{B})/\partial\delta_{\rm HIC}^{2}|_{\delta_{\rm HIC}=0}$ .

Putting this all together, we can now finally plot in Fig. 7 the third-derivative that is needed for the $s/n_{B}$ expansion in Eq. (64). The overall term is negative as we would expect from the combination of Fig. 4 and Eq. (64). Fig. 4 demonstrates that $\partial s/\partial T|_{T=0}$ must decrease with decreasing $Y_{Q}$ , which is only possible from Eq. (64) if the term $\partial^{3}(s/n_{B})/\partial\delta_{\rm HIC}^{2}\partial T|_{\delta_{\rm HIC% }=T=0}$ is negative. Additionally, we find that the overall magnitude of $\partial^{3}(s/n_{B})/\partial\delta_{\rm HIC}^{2}\partial T|_{\delta_{\rm HIC% }=T=0}$ is the largest for small $n_{B}$ , such that small $n_{B}$ is be most sensitive to changes in $Y_{Q}$ . The fact that low $n_{B}$ is more sensitive to $Y_{Q}$ changes is another source of potential error that may appear at low $n_{B}$ .

We can test how well our expansion works by studying how $\partial(s/n_{B})/\partial T|_{T=0}$ varies with $Y_{Q}$ in the RMF model and then try to reproduce its results for a small value of $Y_{Q}$ using Eq. (64). The result of the $s/n_{B}$ expansion at finite $T$ can be found in Fig. 8 where we used $Y_{Q}^{\rm HIC}=0.4$ as our starting point and expanded to $Y_{Q}=0.05$ . We show the result for the baryon density range of $n_{B}=0.5-6\;n_{\rm sat}$ . The difference between $\partial(s/n_{B})/\partial T|_{T=0}$ for $Y_{Q}=0.4$ and $Y_{Q}=0.05$ is already very small and only shows small deviations at low $n_{B}$ , whereas for high $n_{B}$ there is almost no $Y_{Q}$ dependence. The results of our $s/n_{B}$ expansion are shown compared to the $Y_{Q}=0.05$ result from RMF and the $s/n_{B}$ expanded results reproduce the behavior of $\partial(s/n_{B})/\partial T|_{T=0}$ vs $n_{B}$ very well. At low $n_{B}$ , we do see a small deviation, but this regime also corresponds to low $\mu_{B}$ , such that we already anticipate more numerical error there from the finite temperature expansion itself.

Finally, in Fig. 9 we summarize our findings by performing error quantification for the accuracy of our expansion from Eq. (64) across the $\mu_{B},\mu_{Q}$ plane at three different temperatures $T=30,50,100$ MeV in the top, middle, and bottom panels, respectively. We have a density plot across $\mu_{B},\mu_{Q}$ where the colors indicate the $\frac{|\Delta p|}{p_{\rm RMF}}\times 100\%$ as defined in Eq. (111). Our results for $T=30$ MeV show that most of the error is less than $1\%$ , although some higher deviations appear around $\mu_{B}\sim 1000$ MeV. We also do not see almost any $\mu_{Q}$ dependence of the error, rather it depends almost entirely on $\mu_{B}$ . A similar story is seen at $T=50$ MeV although larger deviations begin to appear around $\mu_{B}\sim 1100$ MeV but higher $\mu_{B}$ ’s also demonstrate error only at the $\sim 1\%$ level. Even up to $T\sim 100$ MeV our expansion still works quite well above $\mu_{B}\gtrsim 1150$ MeV where we are only seeing deviations at the percent level. However, it is clear that around $\mu_{B}\gtrsim 1000$ MeV and $T=100$ MeV that error can rise as high as $15-20\%$ , consistent with our findings in Fig. 5.

III.4.4 Conversion between a fixed grid in $\vec{\mu}$ vs a fixed grid in $n_{B},Y_{Q}$

It is significantly more natural to perform our finite $T$ expansion across a grid of fixed $T,\vec{\mu}$ than across a grid of $T,n_{B},Y_{Q}$ , even though the latter are more natural variables in hydrodynamical simulations. In this section, we discuss the consequences of performing our expansion in $T,\vec{\mu}$ and then the need to redo our EOS grid in terms of $T,n_{B},Y_{Q}$ for it to be usable in numerical relativity simulations.

We saw already in Fig. 5 that in the range of $\mu_{B}=1000-2000$ MeV the finite $T$ expansion can reproduce the pressure extremely well at $T=100$ MeV, although small deviations begin at $\mu_{B}\lesssim 1150$ MeV. One complication to this picture, however, is that a given range of $n_{B}$ has a non-trivial relationship with $\mu_{B}$ at finite $T$ . In other words, if we take the range of $\mu_{B}=1000-2000$ MeV and then calculate the baryon density at different temperatures, we find that

n_{B}(T=100,\mu_{B}=1000)\neq n_{B}(T=0,\mu_{B}=1000),

(114)

where the units are in MeV. While this may seem an obvious statement, it has important consequences for error quantification in our approach.

Since the error quantification in Fig. 5 was done at fixed $\mu_{B}$ , the corresponding $n_{B}$ for that range of $\mu_{B}$ is not the same at fixed $T$ . Rather, as we go to higher $T$ than we actually require a lower range in $\mu_{B}$ in order to produce the same range in $n_{B}$ . To see this effect clearly, we fix our range in $n_{B}=1-5n_{\rm sat}$ , which is a very reasonable range reached within neutron star mergers (see, e.g., Most et al. (2019, 2020); Perego et al. (2019) for example ranges). Then we plot the extent of the $n_{B}=1-5n_{\rm sat}$ range in $\mu_{B}$ at different temperatures and indicate integer units of $n_{B}[n_{\rm sat}]$ using tick marks. The results of the range of $n_{B}$ mapped onto $\mu_{B}$ ranges at different temperatures is shown in Fig. 10.

In Fig. 10 we find that a fixed range in $n_{B}=1-5n_{\rm sat}$ shifts to lower $\mu_{B}$ as one reaches higher $T$ . The implication of this effect is that for a fixed $n_{B}$ the error increases with increasing $T$ not only from the expansion itself but also because one shift to lower $\mu_{B}$ where the error is larger. However, this shift is a stronger effect for heavy-ion collisions (i.e., $Y_{Q}=0.4$ ) than for ANM in neutron stars (shown for $Y_{Q}=0.02$ ). One should also note that this effect is not that large at temperatures of $T\lesssim 50$ MeV so depending on the maximum $T$ reached within a merger it may or may not play a role.

We can also see this effect numerically using the expansion derived in Appendix B in Eqs. (129,134). If we just take the expansion of $n_{B}$ up to $\mathcal{O}(T^{2})$ shown in Eq. (134), we find that the change of $n_{B}$ with $T$ at a fixed $\vec{\mu}$ can be determined as:

n_{B}(T,\vec{\mu})-n_{B}(\vec{\mu})\biggr{\rvert}_{T=0}=\frac{1}{2}\frac{% \partial^{2}s(T,\vec{\mu})}{\partial T\partial\mu_{B}}\biggr{\rvert}_{T=0,\mu_% {Q}}T^{2},

(115)

where the sign of $\frac{\partial^{2}s(T,\vec{\mu})}{\partial T\partial\mu_{B}}\bigr{\rvert}_{T=0% ,\mu_{Q}}$ tells you if $n_{B}$ increases or decreases with $T$ . However, in Fig. 4 (top) we can already get a sense of sign of this derivative. In Fig. 4 (top) we see that $\frac{\partial s(T,\vec{\mu})}{\partial T}\bigr{\rvert}_{T=0,Y_{Q}}$ monotonically increases with $\mu_{B}$ such that the derivative is most likely always positive (note that there $Y_{Q}=$ const. rather than $\mu_{Q}$ needed for our expansion, however we expect that difference to be small). Thus, this is consistent with our finding that at a fixed $T,\vec{\mu}$ then $n_{B}$ should always increase with T.

With this in mind, we can then understand the relative error percentage of the pressure as a function of $n_{B},Y_{Q}$ shown in Fig. 11. Fig. 11 is made after calculating the thermodynamic derivatives in Eqs. (71-72) wherein numerical error at low $\mu_{B}$ can influence regions of higher $\mu_{B}$ when taking a derivative. Thus, we find that in Fig. 11 there is a clear increase in the relative error percentage of the pressure across $n_{B},Y_{Q}$ . Above $n_{B}\gtrsim 3n_{\rm sat}$ the error still remains small and is under $10\%$ even up to $T=100$ MeV. Similar to the behavior across $\mu_{Q}$ in Fig. 9, we do not see almost any $Y_{Q}$ dependence to the error in Fig. 11 as well.

In order to better understand the low $n_{B}$ error, let us compare $T=30$ MeV. In Fig. 9 for $\mu_{B}\gtrsim 1000$ MeV most of the error is at the percent level. Then, $T=30$ MeV and $\mu_{B}\gtrsim 1000$ corresponds to $n_{B}\sim 1.5n_{\rm sat}$ (there is little shift in the $n_{B}(\mu_{B})$ relationship at that low of a $T$ ). However, the numerical error in Fig. 11 for $T=30$ MeV and $n_{B}\sim 1.5n_{\rm sat}$ is quite large varying from $10-50\%$ . This increase arises primarily from the numerical derivatives and possibly even from boundary issues (taking numerical derivatives close to where our knowledge ends).

IV Error Quantification and Future Improvements

In this section, we break down the different sources of error in our approach and methods that can be used in future work to improve the error. To be clear, we will begin with the assumption that the $\beta$ -equilibrated $T=0$ EOS is known and only discuss error from the expansion itself. Additionally, we should clarify that there are three different types of error that arise in our approach:

•

Numerical Error. The first source of error is the most intuitive to understand. We require a number of derivatives in our approach, some that are even up to third-order. An especially significant challenge is that these derivatives are often temperature derivatives that must be taken in the limit of $T=0$ . However, a known issue in microscopic EOS codes is that the Fermi-Dirac distribution is difficult to solve in this regime, leading to numerical noise. Another potential source of numerical error is taking derivatives close to the edge of a grid.
•

Expansion Error. We use three different expansions in our approach and only consider terms up to $\mathcal{O}(T^{2})$ for the finite $T$ expansion, $\mathcal{O}(\delta_{\rm HIC}^{2})$ for the symmetry energy and isentrope expansions. Error can arise from higher-order contributions in these expansions.
•

Uncertainty in coefficients. The input for our 3D expansion requires knowledge of coefficients and functions (see Fig. 2 for a summary). The symmetry energy coefficients and the functions required for the finite $T$ expansion can come from experiments or be calculated from a given microscopic EOS (although some numerical error will be introduced). However, if we wish to use a generic $\beta$ -equilibrated EOS at $T=0$ then these inputs will not be known and introduce a new source of error.

IV.1 Numerical Error

The numerical error in our approach is a relatively small contribution but does likely introduce “wave-like” structure when comparing our expanded pressure at $T=100$ MeV to the original RMF EOS in the middle-panel of Fig. 5 and that leads to a similar effect at very high $T=100$ MeV in Fig. 9. The primary issue is that it is difficult to calculate the entropy accurately in the RMF model at low $T$ , which is required for our second and third-order terms in Fig. 3-4. That numerical error may be leading to the extra wiggles that appear in Fig. 4 (bottom) for $\partial^{2}s/\partial T^{2}|_{T=0}$ . There are ways to improve these calculations in microscopic models and we plan to explore them in a future work. However, considering how well our $\mathcal{O}(T^{2})$ results reproduce the pressure up to $T=100$ MeV in Fig. 5 and Fig. 9, we believe that the numerical error at least in $\partial s/\partial T|_{T=0}$ is likely quite small on a grid of $\vec{\mu}$ .

Another source of numerical error arises after the 3D pressure is calculated and we obtain the other thermodynamic observables in the EOS such as $s$ , $\varepsilon$ , $n_{B}$ , $Y_{Q}$ . We must take first-order numerical derivatives of the pressure to recover these thermodynamic variables (see Sec. II.4) or calculate these analytically with further information about $T=0$ (see Appendix B). In this work we calculated the thermodynamic variables numerically and find a significant increase in the error when ones switches to a grid of $n_{B},Y_{Q}$ . The numerical arises for a few different reasons: the shift in $n_{B}$ to lower $\mu_{B}$ at finite $T$ , the higher error from the expansion up to $\mathcal{O}(T^{2})$ in that regime, and potentially boundary effects when taking numerical derivatives since the error is most pronounced close to the edge of our grid. Thus, our error increases as one converts our grid from $T,\vec{\mu}$ into $T,n_{B},Y_{Q}$ .

IV.2 Expansion Error

We can first discuss error that can arise in keeping terms in the finite $T$ expansion only up to $\mathcal{O}(T^{2})$ . To do so we can consider order-of-magnitude contributions for the second and third-order terms in Eq. (58).

We should note that these are dimensionful terms such that they should be considered in conjunction with their $T$ expansion as well. Here we estimate the percentage contribution from the inclusion of the third-order term. To do so, we first define the order-of-magnitude of the maximum contribution to the second and third-order terms, i.e.,

	$\displaystyle\frac{\partial s}{\partial T}\bigg{\|}_{T=0}$	$\displaystyle\sim$	$\displaystyle 10^{5},$
	$\displaystyle\frac{\partial^{2}s}{\partial T^{2}}\bigg{\|}_{T=0}$	$\displaystyle\sim$	$\displaystyle 10^{2},$

Then we can take the ratios of the entire second and third-order terms (including their $T$ contributions), i.e.,

	$\displaystyle\frac{\frac{\partial s}{\partial T}\big{\|}_{T=0}T^{2}}{\frac{% \partial^{2}s}{\partial T^{2}}\big{\|}_{T=0}T^{3}}$	$\displaystyle\sim$	$\displaystyle\frac{10^{2}\;T^{3}}{10^{5}\;T^{2}}$		(116)
		$\displaystyle=$	$\displaystyle 10^{-3}\;T.$		(116)

Thus, for a contribution from the third-order term on the order of $1\%$ , one requires a temperature of at least $T\sim 10$ MeV. For a contribution on the order of $10\%$ , one requires a temperature of at least $T\sim 100$ MeV. Although, we remind the reader here that these are the maximum contributions from the third-order terms and for other regions of the EOS the contribution is sometimes orders of magnitude smaller. Typical neutron star mergers reach up to temperature of about $T\sim 50$ MeV but some simulations find temperatures up to $T\sim 100$ MeV Most et al. (2019, 2020).

However, we provide a word of caution here. In the RMF model we consider hadrons only and no phase transitions exist between these hadrons and quarks. Thus, the order-of-magnitude or the behavior of these contributions may change for models that have more complex phase structures. In fact, one could even consider using these terms to calculate a radius-of-convergence to obtain a phase transition, but we leave that for future studies.

The symmetry energy expansion includes both uncertainty in the coefficients and also the limitation of expanding only up to $\mathcal{O}(\delta^{2}_{\rm iso})$ and the fact that the expansion is around $n_{B}\sim n_{\rm sat}$ such that error increases at large $n_{B}$ . The $\delta_{\rm iso}$ expansion error is likely very small in our approach because the $Y_{Q}$ range needed for neutron star mergers is very close to the $\beta$ -equilibrated value such that one does not need to expand far away. The $n_{B}$ error is more of a concern because the expansion is around $n_{B}\sim n_{\rm sat}$ whereas neutron stars may reach densities up to $[4-10]n_{\rm sat}$ . In this work we use four coefficients but, as was shown in Ref. Yao et al. (2023), error can already appear at higher $n_{B}$ . Thus, we argue that the symmetry energy expansion warrants further investigation in future work.

We also have the isentropic expansion across $\delta_{\rm HIC}$ . We found in this work that works very well with terms up $\delta^{2}_{\rm HIC}$ , although the low $n_{B}$ regime appears to be most sensitive to this expansion. It may be possible to explore higher terms up to $\delta^{4}_{\rm HIC}$ in future work but this appears to be one of the smallest sources of error in this work. Lastly, we point out that we ignored linear contributions to the $s/n_{B}$ expansion. However, in principle that is only a valid assumption if one expands around SNM i.e. as was done in Eq. 36. We plan to relax that assumption in future work to test any influence that it may have had in our results.

IV.3 Uncertainty in coefficients

We have the following input coefficients and functions of $n_{B}$ :

•

coefficients: $E_{\rm sym},\ L_{\rm sym},\ J_{\rm sym},\ K_{\rm sym}$ ,
•

function: $\frac{\partial s_{\rm{HIC}}(T,n_{B},Y_{Q}=Y_{Q}^{\rm{HIC}})}{\partial T}\Big{|% }_{T=0}(n_{B})$ ,
•

function: $\frac{\partial^{3}\tilde{S}_{\rm{HIC},2}(T,n_{B})}{\partial T\partial\delta_{% \rm HIC}^{2}}\Big{|}_{T=\delta_{\rm HIC}=0}(n_{B})$ ,

where the four coefficients come from the symmetry energy expansion and at least the first two of them ( $E_{\rm sym},\ L_{\rm sym}$ ) have bounds set by various experiments Li et al. (2019); Adhikari et al. (2021); Reed et al. (2021). The remaining two ( $J_{\rm sym},\ K_{\rm sym}$ ) have theoretical bounds Tews et al. (2017) but these bounds are very broad still. However, as discussed in Ref. Yao et al. (2023) one can also impose causality and stability constraints on $c_{s}^{2}(n_{B})$ at SNM that further restricts the range of these coefficients. It is also possible to impose astrophysical constraints Xie and Li (2020).

The two functions were first introduced in this work and, therefore, have not yet been constrained from theoretical and experimental work. However, since it is possible to extract information about isentropes from heavy-ion collisions and one has the flexibility to vary $Y_{Q}^{\rm HIC}$ by changing $Z/A$ , there is hope for eventual experimental constraints. Furthermore, future work is warranted on better understanding these functions using a variety of microscopic models. In this work, we have assumed that these functions are known from a model and quantified the error under that assumption.

IV.4 Final error quantification

In Fig. 12 we summarize our findings for the error quantification across $n_{B}$ at different slices in the temperature $T$ . The error is averaged over the relevant range of $Y_{Q}$ for neutron star mergers i.e. $Y_{Q}=0.02-0.2$ and shown on a log scale to emphasize the accuracy across $n_{B}$ . We find an overall trend of the error decreasing with $n_{B}$ , consistent with previous plots. Using $10\%$ error as a gauge, we find that the lower temperatures of $T=10,30$ MeV reach $10\%$ error at $n_{B}\sim 2n_{\rm sat}$ whereas the highter temperatures of $T=50,100$ MeV reach $10\%$ error at $n_{B}\sim 3n_{\rm sat}$ .

Note that the error that appears in Fig. 5 is not the same as what appears in Fig. 12 for a number of reasons that we outline here. These are the different types of error in our analysis of the 3D table:

1.

The finite $T$ expansion up to $\mathcal{O}(T^{2})$ , shows largest deviations at low $\mu_{B}$ . This error can be improved once entropy at low $T$ is numerically more accurate within the RMF model, allowing us to check contributions from higher-order terms in the series.
2.

The symmetry energy expansion has uncertainty in its coefficients, uncertainty in the knowledge of $Y_{Q}(n_{B})$ at $\beta$ -equilibrium, and because it is taken up to $\mathcal{O}(\delta_{\rm iso}^{2})$ and expanded around $n\sim n_{\rm sat}$ . Because the range in $Y_{Q}$ required for numerical relativity simulations is small, the $\mathcal{O}(\delta_{\rm iso}^{2})$ does not produce significant amounts of error. The uncertainty in the symmetry energy coefficients is mitigated by using known constraints and requiring causality and stability. However, the expansion around $n\sim n_{\rm sat}$ leads to unavoidable error at large $n_{B}$ .
3.

The $s/n_{B}$ up to $\mathcal{O}(\delta_{\rm HIC}^{2})$ shows largest deviations at low $n_{B}$ . Since the error is nearly independent of $Y_{Q}$ at a fixed $n_{B}$ , this error is likely negligible.
4.

A fixed $n_{B}$ shifts to lower $\mu_{B}$ at finite $T$ . Because the expansion is worse at lower $\mu_{B}$ then this shift means that at finite $T$ there is a larger $n_{B}$ regime with high error. This is technically not an error that can be improved on, but rather a physical consequence of $n_{B}$ at finite $T$ . However, improvements in 1. will help to reduce the error that arise from this issue.
5.

Numerical calculation of the EOS from $c_{s}^{2}(n_{B})$ . We mitigate this error by using a second-order Runge Kutta.
6.

Calculation of numerical derivatives $s,\ n_{B},\ n_{Q}$ at finite $T$ . There are multiple challenges taking these derivatives as discussed Sec. IV.1. Future work will explore an analytical method to mitigate this issue.

Of these six sources of error, working to resolve 1. would most likely have the largest impact on the uncertainty quantification because it would also aid with 3. and might make it possible to calculate the derivatives in 6. analytically. Additionally, we argue that improving upon 1. should be prioritized because our largest uncertainties occur at low $n_{B}$ and they are in part driven by 1. The best starting point to fix this error is to improve the entropy calculations in the RMF model at low $T$ and studying the influence of higher-order contributions in $T$ .

We also argued that 6. is a significant source of error in this work. In future work, we could improve on this issue by creating a grid in $\vec{\mu}$ well-beyond the regime of validity of the model to avoid boundary issues and trying other numerical approaches to calculate these derivatives. There is also the possibility of analytically calculating these state variables, but, as previously stated, that requires knowledge of new, potentially challenging derivatives at $T=0$ that have not yet been calculated.

It may also be worth thinking about techniques other than an expansion around $n_{\rm sat}$ that is used for the symmetry energy expansion. We do not find large error at large $n_{B}$ in our model. However, if one would like to extrapolate these results to supernova or heavy-ion collisions, then the compounded expansion in both $\delta_{\rm iso}$ and $n_{B}$ could lead to larger issues with uncertainty at large $n_{B}$ .

In contrast, the error arising from expansions in $\delta_{\rm HIC}$ and $\delta_{iso}$ are likely some of the smallest source of error and do not warrant as much attention. Additionally, we estimate that the numerical error is very minimal from 5.

V Conclusions and Outlook

In this paper we present a new approach to incorporating finite $T$ effects for an arbitrary cold, $\beta$ -equilibrated neutron star EOS. Our approach results in a 3D EOS across $T,n_{B},Y_{Q}$ using the symmetry energy expansion that has already been well-established in nuclear physics but introduces two new expansions: the finite $T$ expansion and the entropy over baryon number expansion. We find that the finite $T$ expansion only requires one non-zero term beyond the $T=0$ EOS to reproduce the pressure at $T=100$ MeV at $\mu_{B}\gtrsim 1150$ MeV for a fixed $Y_{Q}$ , using an relativistic mean field model as a benchmark. We then break down all the different sources of error in our approach and discuss methods to continue to decrease the error in future work. This framework provides a clear path for quantifying uncertainty in the finite $T$ EOS relevant for, e.g., neutron star mergers and outlines a connection to heavy-ion collision experiments through the entropy over baryon number expansion.

We have developed an open-source code Mroczek et al. (2024) called FiniteT that will be released upon publication of this work that can take an input $\beta$ -equilibrium EOS or build a functional form of a $\beta$ -equilibrated EOS, which is then expanded into the phase space relevant for numerical relativity simulations of neutron star mergers. The code matches this 3D expanded EOS to a crust at low $n_{B}$ and produces a table in a format compatible with numerical relativity codes. Our code also outputs the $\beta$ -equilibrated EOS, verifies if known nuclear saturation properties are respected, and checks causality and stability of the EOS across the full 3D EOS. Next steps include testing within numerical relativity simulations, especially with a focus on varying finite $T$ effects, grid size convergence tests, and studying how the error in our expansion affects numerical relativity simulations (see also Raithel et al. (2022)).

While our focus in this work has been primarily on neutron star mergers, these techniques could be applied to both heavy-ion collisions and supernova studies. For both heavy-ion collisions and neutron star mergers, the uncertainty quantification would need to be studied with more focus on large $Y_{Q}$ values (for instance, supernovae may reach values of $Y_{Q}>0.5$ (e.g., Oertel et al. (2017); Kumar et al. (2020)), that we did not explore in this work. For heavy-ion collisions, we would not want to match to a crust, as we did in this work, but rather to the lattice QCD EOS and/or a hadron resonance gas (see, e.g., Noronha-Hostler et al. (2019); Monnai et al. (2019)). Additionally, in heavy-ion collisions strangeness plays a large role in the EOS and would need to be included in the expansion as well. An additional unexplored avenue is the impact of approximations in the symmetry energy and how this would impact nuclear reaction rates. This may affect especially supernovae simulations, where more consideration of neutrino transport effects would need to be considered Mezzacappa et al. (2020). Similar considerations also apply to neutrino-driven composition changes in neutron star mergers, including potential imprints of neutrino bulk viscosity, which may be of a comparable order Alford et al. (2018); Hammond et al. (2021, 2023); Most et al. (2021, 2022); Zappa et al. (2022); Espino et al. (2023); Loffredo et al. (2023).

At this point, we have only benchmarked our approach using a relativistic mean field EOS, but we plan to perform future studies on EOS that contain hyperons and/or quarks. We do not anticipate that our finite $T$ expansion can handle first-order phase transitions due to the divergence of thermodynamic state variables at the transition, but cross-over phase transitions into quark phases should be possible if the higher-order terms in the finite $T$ expansion remain small.

VI Acknowledgements

We acknowledge support from the support from the US-DOE Nuclear Science Grant No. DE-SC0023861 and the National Science Foundation under grants PHY1748621, MUSES OAC-2103680, NP3M PHY-2116686, and PHY2309210. D.M is supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE – 1746047 and the Illinois Center for Advanced Studies of the Universe Graduate Fellowship. We also acknowledge support from the Illinois Campus Cluster, a computing resource that is operated by the Illinois Campus Cluster Program (ICCP) in conjunction with the National Center for Supercomputing Applications (NCSA), which is supported by funds from the University of Illinois at Urbana-Champaign. A.H. is partly supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under Award No. #DE-FG02-05ER41375.

References

Miller et al. (2019) M. C. Miller et al., Astrophys. J. Lett. 887, L24 (2019), arXiv:1912.05705 [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. (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] .
Abbott et al. (2019) B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. X 9, 011001 (2019), arXiv:1805.11579 [gr-qc] .
Kumar et al. (2023) R. Kumar et al. (MUSES), (2023), arXiv:2303.17021 [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] .
Raithel and Paschalidis (2023a) C. A. Raithel and V. Paschalidis, (2023a), arXiv:2312.14046 [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] .
Perego et al. (2019) A. Perego, S. Bernuzzi, and D. Radice, Eur. Phys. J. A 55, 124 (2019), arXiv:1903.07898 [gr-qc] .
Villa-Ortega et al. (2023) V. Villa-Ortega, A. Lorenzo-Medina, J. Calderón Bustillo, M. Ruiz, D. Guerra, P. Cerdá-Duran, and J. A. Font, (2023), arXiv:2310.20378 [gr-qc] .
Constantinou et al. (2015) C. Constantinou, B. Muccioli, M. Prakash, and J. M. Lattimer, Annals Phys. 363, 533 (2015), arXiv:1507.07874 [nucl-th] .
Raithel et al. (2021) C. Raithel, V. Paschalidis, and F. Özel, Phys. Rev. D 104, 063016 (2021), arXiv:2104.07226 [astro-ph.HE] .
Most and Raithel (2021) E. R. Most and C. A. Raithel, Phys. Rev. D 104, 124012 (2021), arXiv:2107.06804 [astro-ph.HE] .
Fields et al. (2023) J. Fields, A. Prakash, M. Breschi, D. Radice, S. Bernuzzi, and A. da Silva Schneider, Astrophys. J. Lett. 952, L36 (2023), arXiv:2302.11359 [astro-ph.HE] .
Jacobi et al. (2023) M. Jacobi, F. M. Guercilena, S. Huth, G. Ricigliano, A. Arcones, and A. Schwenk, Mon. Not. Roy. Astron. Soc. 527, 8812 (2023), arXiv:2309.09233 [astro-ph.HE] .
Raithel and Paschalidis (2023b) C. A. Raithel and V. Paschalidis, Phys. Rev. D 108, 083029 (2023b), arXiv:2306.13144 [astro-ph.HE] .
Blacker et al. (2024) S. Blacker, H. Kochankovski, A. Bauswein, A. Ramos, and L. Tolos, Phys. Rev. D 109, 043015 (2024), arXiv:2307.03710 [astro-ph.HE] .
Blacker et al. (2023) S. Blacker, A. Bauswein, and S. Typel, Phys. Rev. D 108, 063032 (2023), arXiv:2304.01971 [astro-ph.HE] .
Wei et al. (2021) J.-B. Wei, G. F. Burgio, A. R. Raduta, and H. J. Schulze, Phys. Rev. C 104, 065806 (2021), arXiv:2112.05323 [nucl-th] .
Raduta (2022) A. R. Raduta, Eur. Phys. J. A 58, 115 (2022), arXiv:2205.03177 [nucl-th] .
Abdallah et al. (2022a) M. S. Abdallah et al. (STAR), Phys. Lett. B 827, 137003 (2022a), arXiv:2108.00908 [nucl-ex] .
Lovato et al. (2022) A. Lovato et al., (2022), arXiv:2211.02224 [nucl-th] .
Oliinychenko et al. (2023) D. Oliinychenko, A. Sorensen, V. Koch, and L. McLerran, Phys. Rev. C 108, 034908 (2023), arXiv:2208.11996 [nucl-th] .
Omana Kuttan et al. (2023) M. Omana Kuttan, J. Steinheimer, K. Zhou, and H. Stoecker, Phys. Rev. Lett. 131, 202303 (2023), arXiv:2211.11670 [hep-ph] .
Yao et al. (2023) N. Yao, A. Sorensen, V. Dexheimer, and J. Noronha-Hostler, (2023), arXiv:2311.18819 [nucl-th] .
Almaalol et al. (2022) D. Almaalol et al., (2022), arXiv:2209.05009 [nucl-ex] .
Sorensen et al. (2024) A. Sorensen et al., Prog. Part. Nucl. Phys. 134, 104080 (2024), arXiv:2301.13253 [nucl-th] .
Mroczek et al. (2024) D. Mroczek, N. Yao, and K. Zine, “Finite temperature expansion,” https://github.com/nucthUIUC/Finite_T (2024).
Alford et al. (2022) M. G. Alford, L. Brodie, A. Haber, and I. Tews, Phys. Rev. C 106, 055804 (2022), arXiv:2205.10283 [nucl-th] .
Borsanyi et al. (2018) S. Borsanyi, Z. Fodor, J. N. Guenther, S. K. Katz, K. K. Szabo, A. Pasztor, I. Portillo, and C. Ratti, JHEP 10, 205 (2018), arXiv:1805.04445 [hep-lat] .
Bazavov et al. (2020) A. Bazavov et al., Phys. Rev. D 101, 074502 (2020), arXiv:2001.08530 [hep-lat] .
Borsányi et al. (2023) S. Borsányi, Z. Fodor, J. N. Guenther, S. D. Katz, P. Parotto, A. Pásztor, D. Pesznyák, K. K. Szabó, and C. H. Wong, (2023), arXiv:2312.07528 [hep-lat] .
Monnai et al. (2019) A. Monnai, B. Schenke, and C. Shen, Phys. Rev. C 100, 024907 (2019), arXiv:1902.05095 [nucl-th] .
Noronha-Hostler et al. (2019) J. Noronha-Hostler, P. Parotto, C. Ratti, and J. M. Stafford, Phys. Rev. C 100, 064910 (2019), arXiv:1902.06723 [hep-ph] .
Karsch (2012) F. Karsch, Central Eur. J. Phys. 10, 1234 (2012), arXiv:1202.4173 [hep-lat] .
Stephanov (2009) M. A. Stephanov, Phys. Rev. Lett. 102, 032301 (2009), arXiv:0809.3450 [hep-ph] .
Stephanov (2011) M. A. Stephanov, Phys. Rev. Lett. 107, 052301 (2011), arXiv:1104.1627 [hep-ph] .
Alford and Harris (2018) M. G. Alford and S. P. Harris, Phys. Rev. C 98, 065806 (2018), arXiv:1803.00662 [nucl-th] .
Alford et al. (2021) M. G. Alford, A. Haber, S. P. Harris, and Z. Zhang, Universe 7, 399 (2021), arXiv:2108.03324 [nucl-th] .
Loffredo et al. (2023) E. Loffredo, A. Perego, D. Logoteta, and M. Branchesi, Astron. Astrophys. 672, A124 (2023), arXiv:2209.04458 [astro-ph.HE] .
Alford et al. (2023) M. G. Alford, A. Haber, and Z. Zhang, (2023), arXiv:2306.06180 [nucl-th] .
Lattimer and Prakash (2000) J. M. Lattimer and M. Prakash, Phys. Rept. 333, 121 (2000), arXiv:astro-ph/0002203 .
Tews et al. (2013) I. Tews, T. Krüger, K. Hebeler, and A. Schwenk, Phys. Rev. Lett. 110, 032504 (2013), arXiv:1206.0025 [nucl-th] .
Drischler et al. (2019) C. Drischler, K. Hebeler, and A. Schwenk, Phys. Rev. Lett. 122, 042501 (2019), arXiv:1710.08220 [nucl-th] .
Baldo and Burgio (2016) M. Baldo and G. F. Burgio, Prog. Part. Nucl. Phys. 91, 203 (2016), arXiv:1606.08838 [nucl-th] .
Tews et al. (2017) I. Tews, J. M. Lattimer, A. Ohnishi, and E. E. Kolomeitsev, Astrophys. J. 848, 105 (2017), arXiv:1611.07133 [nucl-th] .
Zhang et al. (2018) N.-B. Zhang, B.-A. Li, and J. Xu, Astrophys. J. 859, 90 (2018), arXiv:1801.06855 [nucl-th] .
Li et al. (2019) B.-A. Li, P. G. Krastev, D.-H. Wen, and N.-B. Zhang, Eur. Phys. J. A 55, 117 (2019), arXiv:1905.13175 [nucl-th] .
Xie and Li (2020) W.-J. Xie and B.-A. Li, Astrophys. J. 899, 4 (2020), arXiv:2005.07216 [astro-ph.HE] .
Adhikari et al. (2021) D. Adhikari et al. (PREX), Phys. Rev. Lett. 126, 172502 (2021), arXiv:2102.10767 [nucl-ex] .
Reed et al. (2021) B. T. Reed, F. J. Fattoyev, C. J. Horowitz, and J. Piekarewicz, Phys. Rev. Lett. 126, 172503 (2021), arXiv:2101.03193 [nucl-th] .
Sun et al. (2023) B. Sun, S. Bhattiprolu, and J. M. Lattimer, (2023), arXiv:2311.00843 [nucl-th] .
Imam et al. (2022) S. M. A. Imam, N. K. Patra, C. Mondal, T. Malik, and B. K. Agrawal, Phys. Rev. C 105, 015806 (2022), arXiv:2110.15776 [nucl-th] .
Mendes et al. (2021) M. Mendes, A. Cumming, C. Gale, and F. J. Fattoyev, in 16th Marcel Grossmann Meeting on Recent Developments in Theoretical and Experimental General Relativity, Astrophysics and Relativistic Field Theories (2021) arXiv:2110.11077 [nucl-th] .
Hempel et al. (2013) M. Hempel, V. Dexheimer, S. Schramm, and I. Iosilevskiy, Phys. Rev. C 88, 014906 (2013), arXiv:1302.2835 [nucl-th] .
Parotto et al. (2020) P. Parotto, M. Bluhm, D. Mroczek, M. Nahrgang, J. Noronha-Hostler, K. Rajagopal, C. Ratti, T. Schäfer, and M. Stephanov, Phys. Rev. C 101, 034901 (2020), arXiv:1805.05249 [hep-ph] .
Parotto (2019) P. Parotto, Characterization of the Transition Region in the QCD Phase Diagram, Ph.D. thesis, Houston U. (2019).
Parotto et al. (2022) P. Parotto, S. Borsanyi, Z. Fodor, J. N. Guenther, R. Kara, S. D. Katz, A. Pasztor, C. Ratti, and K. K. Szabo, PoS LATTICE2021, 460 (2022), arXiv:2112.00083 [hep-lat] .
Borsanyi et al. (2022) S. Borsanyi, J. N. Guenther, R. Kara, Z. Fodor, P. Parotto, A. Pasztor, C. Ratti, and K. K. Szabo, Phys. Rev. D 105, 114504 (2022), arXiv:2202.05574 [hep-lat] .
Most et al. (2020) E. R. Most, L. Jens Papenfort, V. Dexheimer, M. Hanauske, H. Stoecker, and L. Rezzolla, Eur. Phys. J. A 56, 59 (2020), arXiv:1910.13893 [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] .
Alba et al. (2014) P. Alba, W. Alberico, R. Bellwied, M. Bluhm, V. Mantovani Sarti, M. Nahrgang, and C. Ratti, Phys. Lett. B 738, 305 (2014), arXiv:1403.4903 [hep-ph] .
Alba et al. (2020) P. Alba, V. M. Sarti, J. Noronha-Hostler, P. Parotto, I. Portillo-Vazquez, C. Ratti, and J. M. Stafford, Phys. Rev. C 101, 054905 (2020), arXiv:2002.12395 [hep-ph] .
Guenther et al. (2017) J. N. Guenther, R. Bellwied, S. Borsanyi, Z. Fodor, S. D. Katz, A. Pasztor, C. Ratti, and K. K. Szabó, Nucl. Phys. A 967, 720 (2017), arXiv:1607.02493 [hep-lat] .
Karthein et al. (2021) J. M. Karthein, D. Mroczek, A. R. Nava Acuna, J. Noronha-Hostler, P. Parotto, D. R. P. Price, and C. Ratti, Eur. Phys. J. Plus 136, 621 (2021), arXiv:2103.08146 [hep-ph] .
Dore et al. (2022) T. Dore, J. M. Karthein, I. Long, D. Mroczek, J. Noronha-Hostler, P. Parotto, C. Ratti, and Y. Yamauchi, Phys. Rev. D 106, 094024 (2022), arXiv:2207.04086 [nucl-th] .
Stephanov (2004) M. A. Stephanov, Prog. Theor. Phys. Suppl. 153, 139 (2004), arXiv:hep-ph/0402115 .
Nonaka and Asakawa (2005) C. Nonaka and M. Asakawa, Phys. Rev. C 71, 044904 (2005), arXiv:nucl-th/0410078 .
Asakawa et al. (2008) M. Asakawa, S. A. Bass, B. Muller, and C. Nonaka, Phys. Rev. Lett. 101, 122302 (2008), arXiv:0803.2449 [nucl-th] .
Abdallah et al. (2022b) M. Abdallah et al. (STAR), Phys. Rev. C 105, 014901 (2022b), arXiv:2109.00131 [nucl-ex] .
Greif et al. (2018) M. Greif, J. A. Fotakis, G. S. Denicol, and C. Greiner, Phys. Rev. Lett. 120, 242301 (2018), arXiv:1711.08680 [hep-ph] .
Fotakis et al. (2020) J. A. Fotakis, M. Greif, C. Greiner, G. S. Denicol, and H. Niemi, Phys. Rev. D 101, 076007 (2020), arXiv:1912.09103 [hep-ph] .
Carzon et al. (2022) P. Carzon, M. Martinez, M. D. Sievert, D. E. Wertepny, and J. Noronha-Hostler, Phys. Rev. C 105, 034908 (2022), arXiv:1911.12454 [nucl-th] .
Plumberg et al. (2023) C. Plumberg et al., in 30th International Conference on Ultrarelativstic Nucleus-Nucleus Collisions (2023) arXiv:2312.07415 [hep-ph] .
Poberezhnyuk et al. (2019) R. Poberezhnyuk, V. Vovchenko, A. Motornenko, M. I. Gorenstein, and H. Stoecker, Phys. Rev. C 100, 054904 (2019), arXiv:1906.01954 [hep-ph] .
Borsanyi et al. (2013) S. Borsanyi, Z. Fodor, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo, Phys. Rev. Lett. 111, 062005 (2013), arXiv:1305.5161 [hep-lat] .
Borsanyi et al. (2014) S. Borsanyi, Z. Fodor, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo, Phys. Rev. Lett. 113, 052301 (2014), arXiv:1403.4576 [hep-lat] .
Alver et al. (2011) B. Alver et al. (PHOBOS), Phys. Rev. C 83, 024913 (2011), arXiv:1011.1940 [nucl-ex] .
Abelev et al. (2009) B. I. Abelev et al. (STAR), Phys. Rev. C 79, 034909 (2009), arXiv:0808.2041 [nucl-ex] .
Abdallah et al. (2023) M. S. Abdallah et al. (STAR), Phys. Rev. C 107, 024901 (2023), arXiv:2208.00653 [nucl-ex] .
Lattimer and Prakash (2001) J. M. Lattimer and M. Prakash, Astrophys. J. 550, 426 (2001), arXiv:astro-ph/0002232 .
Chamel and Haensel (2008) N. Chamel and P. Haensel, Living Rev. Rel. 11, 10 (2008), arXiv:0812.3955 [astro-ph] .
Grill et al. (2014) F. Grill, H. Pais, C. Providência, I. Vidaña, and S. S. Avancini, Phys. Rev. C 90, 045803 (2014), arXiv:1404.2753 [nucl-th] .
Pais et al. (2016) H. Pais, A. Sulaksono, B. K. Agrawal, and C. Providência, Phys. Rev. C 93, 045802 (2016), arXiv:1603.06899 [nucl-th] .
Ferreira and Providência (2020a) M. Ferreira and C. Providência, Universe 6, 220 (2020a).
Ferreira and Providência (2020b) M. Ferreira and C. Providência, Phys. Rev. D 102, 103003 (2020b), arXiv:2010.05588 [nucl-th] .
Tan et al. (2020) H. Tan, J. Noronha-Hostler, and N. Yunes, Phys. Rev. Lett. 125, 261104 (2020), arXiv:2006.16296 [astro-ph.HE] .
Tan et al. (2022a) H. Tan, T. Dore, V. Dexheimer, J. Noronha-Hostler, and N. Yunes, Phys. Rev. D 105, 023018 (2022a), arXiv:2106.03890 [astro-ph.HE] .
Tan et al. (2022b) H. Tan, V. Dexheimer, J. Noronha-Hostler, and N. Yunes, Phys. Rev. Lett. 128, 161101 (2022b), arXiv:2111.10260 [astro-ph.HE] .
Mroczek et al. (2023) D. Mroczek, M. C. Miller, J. Noronha-Hostler, and N. Yunes, (2023), arXiv:2309.02345 [astro-ph.HE] .
Schneider et al. (2017) A. S. Schneider, L. F. Roberts, and C. D. Ott, Phys. Rev. C 96, 065802 (2017), arXiv:1707.01527 [astro-ph.HE] .
(93) “SLy4 with NSE (3335 nuclides),” https://stellarcollapse.org/EOS/SLy4_3335_rho391_temp163_ye66.h5.bz2, 2023.
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] .
Landry and Essick (2019) P. Landry and R. Essick, Phys. Rev. D 99, 084049 (2019), arXiv:1811.12529 [gr-qc] .
Han et al. (2021) M.-Z. Han, J.-L. Jiang, S.-P. Tang, and Y.-Z. Fan, Astrophys. J. 919, 11 (2021), arXiv:2103.05408 [hep-ph] .
Legred et al. (2021) I. Legred, K. Chatziioannou, R. Essick, S. Han, and P. Landry, Phys. Rev. D 104, 063003 (2021), arXiv:2106.05313 [astro-ph.HE] .
Landry et al. (2020) P. Landry, R. Essick, and K. Chatziioannou, Phys. Rev. D 101, 123007 (2020), arXiv:2003.04880 [astro-ph.HE] .
Essick et al. (2020) R. Essick, P. Landry, and D. E. Holz, Phys. Rev. D 101, 063007 (2020), arXiv:1910.09740 [astro-ph.HE] .
Alford et al. (2013) M. G. Alford, S. Han, and M. Prakash, Phys. Rev. D 88, 083013 (2013), arXiv:1302.4732 [astro-ph.SR] .
Tews et al. (2018) I. Tews, J. Carlson, S. Gandolfi, and S. Reddy, Astrophys. J. 860, 149 (2018), arXiv:1801.01923 [nucl-th] .
Demorest et al. (2010) P. B. Demorest, T. Pennucci, S. Ransom, M. Roberts, and J. Hessels, nature 467, 1081 (2010).
Arzoumanian et al. (2018) Z. Arzoumanian, P. Baker, A. Brazier, S. Burke-Spolaor, S. Chamberlin, S. Chatterjee, B. Christy, J. M. Cordes, N. J. Cornish, F. Crawford, et al., The Astrophysical Journal 859, 47 (2018).
Cromartie et al. (2020) H. T. Cromartie, E. Fonseca, S. M. Ransom, P. B. Demorest, Z. Arzoumanian, H. Blumer, P. R. Brook, M. E. DeCesar, T. Dolch, J. A. Ellis, et al., Nature Astronomy 4, 72 (2020).
Fonseca et al. (2021) E. Fonseca et al., Astrophys. J. Lett. 915, L12 (2021), arXiv:2104.00880 [astro-ph.HE] .
Agazie et al. (2023) G. Agazie et al. (NANOGrav), Astrophys. J. Lett. 951, L9 (2023), arXiv:2306.16217 [astro-ph.HE] .
Raithel et al. (2022) C. A. Raithel, P. Espino, and V. Paschalidis, Mon. Not. Roy. Astron. Soc. 516, 4792 (2022), arXiv:2206.14838 [astro-ph.HE] .
Oertel et al. (2017) M. Oertel, M. Hempel, T. Klähn, and S. Typel, Rev. Mod. Phys. 89, 015007 (2017), arXiv:1610.03361 [astro-ph.HE] .
Kumar et al. (2020) A. Kumar, H. C. Das, S. K. Biswal, B. Kumar, and S. K. Patra, Eur. Phys. J. C 80, 775 (2020), arXiv:2005.08320 [nucl-th] .
Mezzacappa et al. (2020) A. Mezzacappa, E. Endeve, O. E. Bronson Messer, and S. W. Bruenn, Liv. Rev. Comput. Astrophys. 6, 4 (2020), arXiv:2010.09013 [astro-ph.HE] .
Alford et al. (2018) M. G. Alford, L. Bovard, M. Hanauske, L. Rezzolla, and K. Schwenzer, Phys. Rev. Lett. 120, 041101 (2018), arXiv:1707.09475 [gr-qc] .
Hammond et al. (2023) P. Hammond, I. Hawke, and N. Andersson, Phys. Rev. D 107, 043023 (2023), arXiv:2205.11377 [astro-ph.HE] .
Most et al. (2021) E. R. Most, S. P. Harris, C. Plumberg, M. G. Alford, J. Noronha, J. Noronha-Hostler, F. Pretorius, H. Witek, and N. Yunes, Mon. Not. Roy. Astron. Soc. 509, 1096 (2021), arXiv:2107.05094 [astro-ph.HE] .
Most et al. (2022) E. R. Most, A. Haber, S. P. Harris, Z. Zhang, M. G. Alford, and J. Noronha, (2022), arXiv:2207.00442 [astro-ph.HE] .
Zappa et al. (2022) F. Zappa, S. Bernuzzi, D. Radice, and A. Perego, (2022), 10.1093/mnras/stad107, arXiv:2210.11491 [astro-ph.HE] .
Espino et al. (2023) P. L. Espino, P. Hammond, D. Radice, S. Bernuzzi, R. Gamba, F. Zappa, L. F. L. Micchi, and A. Perego, (2023), arXiv:2311.00031 [astro-ph.HE] .

Appendix A Converting between different thermodynamical bases

We want to express $\dfrac{\partial s}{\partial T}\Big{|}_{\vec{\mu}}$ through the thermodynamical basis $(T,n_{B},Y_{Q})$ . In this basis, the differential of the entropy is

ds=\dfrac{\partial s}{\partial T}\Big{|}_{n_{B},Y_{Q}}dT+\dfrac{\partial s}{% \partial n_{B}}\Big{|}_{T,Y_{Q}}dn_{B}+\dfrac{\partial s}{\partial Y_{Q}}\Big{% |}_{n_{B},T}dY_{Q},

(117)

so we can calculate

\dfrac{\partial s}{\partial T}\Big{|}_{\vec{\mu}}=\dfrac{\partial s}{\partial T% }\Big{|}_{n_{B},Y_{Q}}+\dfrac{\partial s}{\partial n_{B}}\Big{|}_{T,Y_{Q}}% \dfrac{\partial n_{B}}{\partial T}\Big{|}_{\vec{\mu}}+\dfrac{\partial s}{% \partial Y_{Q}}\Big{|}_{n_{B},T}\dfrac{\partial Y_{Q}}{\partial T}\Big{|}_{% \vec{\mu}}.

(118)

We can then write Eq. (118) concisely in terms of a matrix product,

\dfrac{\partial s}{\partial T}\Big{|}_{\vec{\mu}}=\dfrac{\partial s}{\partial T% }\Big{|}_{n_{B},Y_{Q}}+\begin{pmatrix}\dfrac{\partial s}{\partial n_{B}}\Big{|% }_{T,Y_{Q}}&0\\ 0&\dfrac{\partial s}{\partial Y_{Q}}\Big{|}_{n_{B},T}\\ \end{pmatrix}\begin{pmatrix}\dfrac{\partial n_{B}}{\partial T}\Big{|}_{\vec{% \mu}}\\ \dfrac{\partial Y_{Q}}{\partial T}\Big{|}_{\vec{\mu}}\end{pmatrix}.

(119)

The elements of the column vector specify trajectories of constant $\vec{\mu}$ on the $(n_{B},Y_{Q})$ plane. These trajectories are given by the constraint equations,

	$\displaystyle 0$	$\displaystyle=d\mu_{B}=\dfrac{\partial\mu_{B}}{\partial T}\Big{\|}_{n_{B},Y_{Q}% }dT+\dfrac{\partial\mu_{B}}{\partial n_{B}}\Big{\|}_{T,Y_{Q}}dn_{B}+\dfrac{% \partial\mu_{B}}{\partial Y_{Q}}\Big{\|}_{T,n_{B}}dY_{Q},$		(120)
	$\displaystyle 0$	$\displaystyle=d\mu_{Q}=\dfrac{\partial\mu_{Q}}{\partial T}\Big{\|}_{n_{B},Y_{Q}% }dT+\dfrac{\partial\mu_{Q}}{\partial n_{B}}\Big{\|}_{T,Y_{Q}}dn_{B}+\dfrac{% \partial\mu_{Q}}{\partial Y_{Q}}\Big{\|}_{T,n_{B}}dY_{Q},$		(121)

implying that

	$\displaystyle 0$	$\displaystyle=\dfrac{\partial\mu_{B}}{\partial T}\Big{\|}_{n_{B},Y_{Q}}+\dfrac{% \partial\mu_{B}}{\partial n_{B}}\Big{\|}_{T,Y_{Q}}\dfrac{\partial n_{B}}{% \partial T}\Big{\|}_{\vec{\mu}}+\dfrac{\partial\mu_{B}}{\partial Y_{Q}}\Big{\|}_% {T,n_{B}}\dfrac{\partial Y_{Q}}{\partial T}\Big{\|}_{\vec{\mu}},$		(122)
	$\displaystyle 0$	$\displaystyle=\dfrac{\partial\mu_{Q}}{\partial T}\Big{\|}_{n_{B},Y_{Q}}+\dfrac{% \partial\mu_{Q}}{\partial n_{B}}\Big{\|}_{T,Y_{Q}}\dfrac{\partial n_{B}}{% \partial T}\Big{\|}_{\vec{\mu}}+\dfrac{\partial\mu_{Q}}{\partial Y_{Q}}\Big{\|}_% {T,n_{B}}\dfrac{\partial Y_{Q}}{\partial T}\Big{\|}_{\vec{\mu}}.$		(123)

Again, we write Eq. 122 as a matrix product,

-\begin{pmatrix}\dfrac{\partial\mu_{B}}{\partial T}\Big{|}_{n_{B},Y_{Q}}\\ \dfrac{\partial\mu_{Q}}{\partial T}\Big{|}_{n_{B},Y_{Q}}\end{pmatrix}=\begin{% pmatrix}\dfrac{\partial\mu_{B}}{\partial n_{B}}\Big{|}_{T,Y_{Q}}&\dfrac{% \partial\mu_{B}}{\partial Y_{Q}}\Big{|}_{T,n_{B}}\\ \dfrac{\partial\mu_{Q}}{\partial T}\Big{|}_{n_{B},Y_{Q}}&\dfrac{\partial\mu_{Q% }}{\partial Y_{Q}}\Big{|}_{T,n_{B}}\end{pmatrix}\begin{pmatrix}\dfrac{\partial n% _{B}}{\partial T}\Big{|}_{\vec{\mu}}\\ \dfrac{\partial Y_{Q}}{\partial T}\Big{|}_{\vec{\mu}}\end{pmatrix}.

(124)

We can now invert the matrix to obtain

\begin{pmatrix}\dfrac{\partial n_{B}}{\partial T}\Big{|}_{\vec{\mu}}\\ \dfrac{\partial Y_{Q}}{\partial T}\Big{|}_{\vec{\mu}}\end{pmatrix}=-\begin{% pmatrix}\dfrac{\partial\mu_{B}}{\partial n_{B}}\Big{|}_{T,Y_{Q}}&\dfrac{% \partial\mu_{B}}{\partial Y_{Q}}\Big{|}_{T,n_{B}}\\ \dfrac{\partial\mu_{Q}}{\partial T}\Big{|}_{n_{B},Y_{Q}}&\dfrac{\partial\mu_{Q% }}{\partial Y_{Q}}\Big{|}_{T,n_{B}}\end{pmatrix}^{-1}\begin{pmatrix}\dfrac{% \partial\mu_{B}}{\partial T}\Big{|}_{n_{B},Y_{Q}}\\ \dfrac{\partial\mu_{Q}}{\partial T}\Big{|}_{n_{B},Y_{Q}}\end{pmatrix}.

(125)

Finally, we substitute Eq. 125 into Eq. 119 to obtain

\dfrac{\partial s}{\partial T}\Big{|}_{\vec{\mu}}=\dfrac{\partial s}{\partial T% }\Big{|}_{n_{B},Y_{Q}}-\begin{pmatrix}\dfrac{\partial s}{\partial n_{B}}\Big{|% }_{T,Y_{Q}}&0\\ 0&\dfrac{\partial s}{\partial Y_{Q}}\Big{|}_{n_{B},T}\\ \end{pmatrix}\begin{pmatrix}\dfrac{\partial\mu_{B}}{\partial n_{B}}\Big{|}_{T,% Y_{Q}}&\dfrac{\partial\mu_{B}}{\partial Y_{Q}}\Big{|}_{T,n_{B}}\\ \dfrac{\partial\mu_{Q}}{\partial T}\Big{|}_{n_{B},Y_{Q}}&\dfrac{\partial\mu_{Q% }}{\partial Y_{Q}}\Big{|}_{T,n_{B}}\end{pmatrix}^{-1}\begin{pmatrix}\dfrac{% \partial\mu_{B}}{\partial T}\Big{|}_{n_{B},Y_{Q}}\\ \dfrac{\partial\mu_{Q}}{\partial T}\Big{|}_{n_{B},Y_{Q}}\end{pmatrix}.

(126)

If we now evaluate this expression at $T=0$ ,

\dfrac{\partial s}{\partial T}\Big{|}_{\mu_{X},T=0}=\dfrac{\partial s}{% \partial T}\Big{|}_{n_{B},Y_{Q},T=0}-\begin{pmatrix}\dfrac{\partial s}{% \partial n_{B}}\Big{|}_{T=0,Y_{Q}}&0\\ 0&\dfrac{\partial s}{\partial Y_{Q}}\Big{|}_{n_{B},T=0}\\ \end{pmatrix}\begin{pmatrix}\dfrac{\partial\mu_{B}}{\partial n_{B}}\Big{|}_{T=% 0,Y_{Q}}&\dfrac{\partial\mu_{B}}{\partial Y_{Q}}\Big{|}_{T=0,n_{B}}\\ \dfrac{\partial\mu_{Q}}{\partial T}\Big{|}_{n_{B},Y_{Q},T=0}&\dfrac{\partial% \mu_{Q}}{\partial Y_{Q}}\Big{|}_{T=0,n_{B}}\end{pmatrix}^{-1}\begin{pmatrix}% \dfrac{\partial\mu_{B}}{\partial T}\Big{|}_{n_{B},Y_{Q},T=0}\\ \dfrac{\partial\mu_{Q}}{\partial T}\Big{|}_{n_{B},Y_{Q},T=0}\end{pmatrix},

(127)

knowing that for a system with a non-degenerate ground state for all $\mu_{B}$ and $\mu_{Q}$ we expect $s=0$ at $T=0$ , the terms $\dfrac{\partial s}{\partial n_{B}}\Big{|}_{T=0,Y_{Q}}$ and $\dfrac{\partial s}{\partial Y_{Q}}\Big{|}_{n_{B},T=0}$ should vanish, and thus $\dfrac{\partial s}{\partial T}\Big{|}_{\mu_{X},T=0}=\dfrac{\partial s}{% \partial T}\Big{|}_{n_{B},Y_{Q},T=0}$ should hold in the $T=0$ limit. We verified this numerically for the RMF model shown in this work.

Appendix B Other thermodynamic observables for the finite $T$ expansion

In fact, it is possible to derive the other thermodynamic quantities directly from Eq. (58). In Eq. (58) we consider terms only up to $\mathcal{O}(T^{3})$ in the pressure, which is then reflected in our equations below such that we drop high-order terms. Beginning with the entropy:

	$\displaystyle s(T,\vec{\mu})$	$\displaystyle=$	$\displaystyle\frac{dp(T,\vec{\mu})}{dT}\biggr{\rvert}_{\vec{\mu}}$		(128)
		$\displaystyle=$	$\displaystyle 2\frac{\partial s(T,\vec{\mu})}{\partial T}\biggr{\rvert}_{T=0,% \vec{\mu}}T+\frac{\partial s^{2}(T,\vec{\mu})}{\partial T^{2}}\biggr{\rvert}_{% T=0,\vec{\mu}}T^{2}+\frac{1}{3}\frac{\partial^{3}s(T,\vec{\mu})}{\partial T^{3% }}\biggr{\rvert}_{T=0,\vec{\mu}}T^{3}+\mathcal{O}(T^{4}).$		(128)

The same can be done for $n_{B}$ :

$\displaystyle n_{B}(T,\vec{\mu})$	$\displaystyle=$	$\displaystyle\frac{\partial p(T,\vec{\mu})}{\partial\mu_{B}}\biggr{\rvert}_{T,% \mu_{Q}}$	(129)
	$\displaystyle=$	$\displaystyle\frac{\partial p(\vec{\mu})}{\partial\mu_{B}}\biggr{\rvert}_{T=0,% \mu_{Q}}+\underbrace{\frac{\partial s(T,\vec{\mu})}{\partial\mu_{B}}\biggr{% \rvert}_{T=0,\mu_{Q}}}_{=0}T+\frac{1}{2}\frac{\partial^{2}s(T,\vec{\mu})}{% \partial T\partial\mu_{B}}\biggr{\rvert}_{T=0,\mu_{Q}}T^{2}+\frac{1}{6}\frac{% \partial^{3}s(T,\vec{\mu})}{\partial T^{2}\partial\mu_{B}}\biggr{\rvert}_{T=0,% \mu_{Q}}T^{3}+\mathcal{O}(T^{4})$
	$\displaystyle=$	$\displaystyle n_{B}(\vec{\mu})\biggr{\rvert}_{T=0}+\frac{1}{2}\frac{\partial^{% 2}s(T,\vec{\mu})}{\partial T\partial\mu_{B}}\biggr{\rvert}_{T=0,\mu_{Q}}T^{2}+% \frac{1}{6}\frac{\partial^{3}s(T,\vec{\mu})}{\partial T^{2}\partial\mu_{B}}% \biggr{\rvert}_{T=0,\mu_{Q}}T^{3}+\mathcal{O}(T^{4}),$

where $\partial s/\partial\mu_{B}$ is zero in the limit of $T=0$ since the entropy is zero.

And also for $n_{Q}$ :

	$\displaystyle n_{Q}(T,\vec{\mu})$	$\displaystyle=$	$\displaystyle\frac{\partial p(T,\vec{\mu})}{\partial\mu_{Q}}\biggr{\rvert}_{T,% \mu_{B}}$		(130)
		$\displaystyle=$	$\displaystyle n_{Q}(\vec{\mu})\biggr{\rvert}_{T=0}+\frac{1}{2}\frac{\partial^{% 2}s(T,\vec{\mu})}{\partial T\partial\mu_{Q}}\biggr{\rvert}_{T=0,\mu_{B}}T^{2}+% \frac{1}{6}\frac{\partial^{3}s(T,\vec{\mu})}{\partial T^{2}\partial\mu_{Q}}% \biggr{\rvert}_{T=0,\mu_{B}}T^{3}+\mathcal{O}(T^{4}).$		(130)

We can also calculate higher-order susceptibilities:

	$\displaystyle\chi_{X}^{n}(T,\vec{\mu})$	$\displaystyle=$	$\displaystyle\frac{\partial^{n}p(T,\vec{\mu})}{(\partial\mu_{X})^{n}}\biggr{% \rvert}_{T,\mu_{Y\neq X}}$		(131)
		$\displaystyle=$	$\displaystyle\chi_{X}^{n}(\vec{\mu})\bigr{\rvert}_{T=0}+\frac{1}{2}\frac{% \partial^{n+1}s(T,\vec{\mu})}{\partial T(\partial\mu_{X})^{n}}\biggr{\rvert}_{% T=0,\mu_{Y\neq X}}T^{2}+\frac{1}{6}\frac{\partial^{n+2}s(T,\vec{\mu})}{% \partial T^{2}(\partial\mu_{X})^{n}}\biggr{\rvert}_{T=0,\mu_{Y\neq X}}T^{3}+% \mathcal{O}(T^{4}),$		(131)

where $X=B,Q$ .

For a number of EOS it may be sufficient to only expand the pressure up to $\mathcal{O}(T^{2})$ and in that case we summarize the corresponding equations:

$\displaystyle p(T,\vec{\mu})$	$\displaystyle=$	$\displaystyle p_{T=0}(\vec{\mu})+\frac{1}{2}\frac{\partial s(T,\vec{\mu})}{% \partial T}\biggr{\rvert}_{T=0,\vec{\mu}}T^{2}+\mathcal{O}(T^{3}),$	(132)
$\displaystyle s(T,\vec{\mu})$	$\displaystyle=$	$\displaystyle 2\frac{\partial s(T,\vec{\mu})}{\partial T}\biggr{\rvert}_{T=0,% \vec{\mu}}T+\frac{\partial^{2}s(T,\vec{\mu})}{\partial^{2}T}\biggr{\rvert}_{T=% 0,\vec{\mu}}T^{2}+\mathcal{O}(T^{3}),$	(133)
$\displaystyle n_{B}(T,\vec{\mu})$	$\displaystyle=$	$\displaystyle n_{B}(\vec{\mu})\biggr{\rvert}_{T=0}+\frac{1}{2}\frac{\partial^{% 2}s(T,\vec{\mu})}{\partial T\partial\mu_{B}}\biggr{\rvert}_{T=0,\mu_{Q}}T^{2}+% \mathcal{O}(T^{3}),$	(134)
$\displaystyle n_{Q}(T,\vec{\mu})$	$\displaystyle=$	$\displaystyle n_{Q}(\vec{\mu})\biggr{\rvert}_{T=0}+\frac{1}{2}\frac{\partial^{% 2}s(T,\vec{\mu})}{\partial T\partial\mu_{Q}}\biggr{\rvert}_{T=0,\mu_{B}}T^{2}+% \mathcal{O}(T^{3}).$	(135)

The number of required coefficients are listed for $\mathcal{O}(T^{2})$ and $\mathcal{O}(T^{3})$ in Table. 2. We find that for the finite $T$ expansion, we require only $4$ coefficients for up to $\mathcal{O}(T^{2})$ and $7$ coefficients up to $\mathcal{O}(T^{3})$ .

Thermodynamics	$\mathcal{O}(T^{2})$	$\mathcal{O}(T^{3})$
Pressure	$\frac{\partial s(T,\vec{\mu})}{\partial T}\bigr{\rvert}_{T=0,\vec{\mu}}$	$\frac{\partial s(T,\vec{\mu})}{\partial T}\bigr{\rvert}_{T=0,\vec{\mu}}$
		$\frac{\partial^{2}s(T,\vec{\mu})}{\partial^{2}T}\bigr{\rvert}_{T=0,\vec{\mu}}$
Entropy	$\frac{\partial s(T,\vec{\mu})}{\partial T}\bigr{\rvert}_{T=0,\vec{\mu}}$	$\frac{\partial s(T,\vec{\mu})}{\partial T}\bigr{\rvert}_{T=0,\vec{\mu}}$
	$\frac{\partial^{2}s(T,\vec{\mu})}{\partial^{2}T}\bigr{\rvert}_{T=0,\vec{\mu}}$	$\frac{\partial^{2}s(T,\vec{\mu})}{\partial^{2}T}\bigr{\rvert}_{T=0,\vec{\mu}}$
		$\frac{\partial^{3}s(T,\vec{\mu})}{\partial^{3}T}\bigr{\rvert}_{T=0,\vec{\mu}}$
Baryon Density	$\frac{\partial^{2}s(T,\vec{\mu})}{\partial T\partial\mu_{B}}\bigr{\rvert}_{T=0% ,\mu_{Q}}$	$\frac{\partial^{2}s(T,\vec{\mu})}{\partial T\partial\mu_{B}}\bigr{\rvert}_{T=0% ,\mu_{Q}}$
		$\frac{\partial^{3}s(T,\vec{\mu})}{\partial T^{2}\partial\mu_{B}}\bigr{\rvert}_% {T=0,\mu_{Q}}$
Charge Density	$\frac{\partial^{2}s(T,\vec{\mu})}{\partial T\partial\mu_{Q}}\bigr{\rvert}_{T=0% ,\mu_{B}}$	$\frac{\partial^{2}s(T,\vec{\mu})}{\partial T\partial\mu_{Q}}\bigr{\rvert}_{T=0% ,\mu_{B}}$
		$\frac{\partial^{3}s(T,\vec{\mu})}{\partial T^{2}\partial\mu_{Q}}\bigr{\rvert}_% {T=0,\mu_{B}}$
Total Unique Terms	4	7

Table 2: List of all needed thermodynamic coefficients up to orders

T^{2}

and

T^{3}

for the analytical finite

T

expansion.

Appendix C $Y_{Q}$ expansion of entropy derivatives

Following the $\delta_{\rm HIC}$ expansion of the isentropes in Sec. II.3.2 we can obtain the $Y_{Q}$ dependence of the derivatives shown in Table 2. We have already derived the $\partial s/\partial T$ term in Eq. (64). Thus, we begin with the $\partial^{2}s/\partial T^{2}$ by applying another $T$ derivative on Eq. (64)

\displaystyle\frac{\partial^{2}\tilde{S}(T,n_{B},Y_{Q})}{\partial T^{2}}\Bigg{% |}_{T=0,n_{B}}

\displaystyle=\frac{1}{n_{B}}\frac{\partial^{2}s_{\rm{HIC}}(T,n_{B},Y_{Q})}{% \partial T^{2}}\Bigg{|}_{T=\delta_{\rm HIC}=0,n_{B}}+\frac{1}{2}\left(1-\frac{% Y_{Q}}{Y_{Q}^{\rm{HIC}}}\right)^{2}\frac{\partial^{4}\tilde{S}_{\rm{HIC},2}(T,% n_{B})}{\partial T^{2}\partial\delta_{\rm HIC}^{2}}\Bigg{|}_{T=\delta_{\rm HIC% }=0,n_{B}}.

(136)

In fact, we can generically write:

\displaystyle\frac{\partial^{n}\tilde{S}(T,n_{B},Y_{Q})}{\partial T^{n}}\Bigg{% |}_{T=0,n_{B}}

\displaystyle=\frac{1}{n_{B}}\frac{\partial^{n}s_{\rm{HIC}}(T,n_{B},Y_{Q})}{% \partial T^{n}}\Bigg{|}_{T=\delta_{\rm HIC}=0,n_{B}}+\frac{1}{2}\left(1-\frac{% Y_{Q}}{Y_{Q}^{\rm{HIC}}}\right)^{2}\frac{\partial^{n+2}\tilde{S}_{\rm{HIC},2}(% T,n_{B})}{\partial T^{n}\partial\delta_{\rm HIC}^{2}}\Bigg{|}_{T=\delta_{\rm HIC% }=0,n_{B}},

(137)

where $n$ is the order of the derivative. However, the numerical stability of higher-order derivatives may become quite challenging, especially if one would like to connect to heavy-ion data.

Appendix D Unit conversion

The SLy4 data table is provided along a fixed grid of

\log_{10}T,\log_{10}\rho_{B},Y_{Q},

(138)

where they are in the following units:

T\;MeV,\rho_{B}\;\left[\frac{g}{cm^{3}}\right],Y_{Q},

(139)

where $\rho_{B}$ is the baryon mass density and $Y_{Q}$ is dimensionless. However, in this work we assume natural units $\hbar=c=1$ , which are more natural for calculations of a relativistic EOS. Thus, our first step is to convert all the quantities in the EOS into natural units.

Even for natural units there is a choice of writing variables in terms of some power of MeV, fm or a mixture of the two. All EOS calculations within our code are performed in powers of MeV. However, for plots, it is normally advantageous to write in a mixture of variables. We have compiled a summary table of all the needed thermodynamic variables in Tab. 3 and compared their units and format between the table, the code, and the plots. For our final EOS that we output to be used within numerical relativity simulations, we convert all variables back to the original table format. We also provide the user the option to output the table in $\left[MeV\right]$ as well.

Variable	Table	Code	Plots
$T$	$\log_{10}$ and $\left[MeV\right]$	$\left[MeV\right]$	$\left[MeV\right]$
$n_{B}$	$\rho_{B}$ $\left[\frac{g}{cm^{3}}\right]$	$n_{B}$ $\left[MeV^{3}\right]$	$\left[1/fm^{3}\right]$
$Y_{Q}$	$Y_{Q}$	$Y_{Q}$	$Y_{Q}$
$\mu_{B}$	$\mu_{n}+m_{N}$ and $\left[MeV\right]$	$\left[MeV\right]$	$\left[MeV\right]$
$\mu_{Q}$	$\mu_{p}$ - $\mu_{n}$ $\left[MeV\right]$	$\left[MeV\right]$	$\left[MeV\right]$
$p$	$\log_{10}$ and $\left[\frac{dyn}{cm^{2}}\right]$	$\left[MeV^{4}\right]$	$\left[MeV/fm^{3}\right]$
$\varepsilon$	$\log_{10}$ and $\left[\frac{erg}{g}\right]$	$\left[MeV^{4}\right]$	$\left[MeV/fm^{3}\right]$
$c_{s}^{2}$	$\left[\frac{cm^{2}}{s^{2}}\right]$	$c$	$c$
$\mu_{n}$	$\left[MeV\right]$	$\left[MeV\right]$	$\left[MeV\right]$
$\mu_{p}$	$\left[MeV\right]$	$\left[MeV\right]$	$\left[MeV\right]$
$\mu_{e}$	$\left[MeV\right]$	$\left[MeV\right]$	$\left[MeV\right]$

Table 3: List of thermodynamic variables required within the code, their units and format obtained from the crust table of SLy4, the units used in the code, and the units used for plots.

From the table, we are provided with T, $\rho_{B}$ , $Y_{Q}$ , $p$ , $\varepsilon$ , $c_{s}^{2}$ , and the chemical potentials for protons $\mu_{p}$ , electrons $\mu_{e}$ , and neutrons $\mu_{n}$ . However, the variables required for our calculations are $T$ , $n_{B}$ , $Y_{Q}$ , $p$ , $\varepsilon$ , $c_{s}^{2}$ , $\mu_{B}$ , $\mu_{Q}$ . Here we walk through each variable, any required unit conversion, and other changes requires:

Baryon density The table provides the baryon mass density, $\rho_{B}$ , in units of $\left[g/cm^{3}\right]$ along a grid in $log_{10}$ , whereas we require the baryon number density in the code in $\left[MeV^{3}\right]$ and for our plots in $\left[fm^{-3}\right]$ . The first step is to remove the $\log_{10}$ such that one takes the input to the $10^{th}$ power: $\rho_{B}=10^{\rho_{B,log10}}$ . To convert from the mass density to number density, we do the following:

n_{B}[fm^{-3}]=\frac{\rho_{B}\left[g\;cm^{-3}\right]}{m_{N}[g]}\frac{1\;[cm^{3% }]}{(10^{13})^{3}\;[fm^{3}]},

(140)

where $m_{N}$ is the mass of the nucleons, which is $m_{N}[g]=940.6MeV\cdot 1.783\cdot 10^{-30}\left[\frac{g}{MeV}\right]$ in Sly4. Thus, the full conversion ends up being

n_{B}[fm^{-3}]=5.97\cdot 10^{16}\cdot 10^{\rho_{B,table}},

(141)

to obtain $n_{B}[fm^{-3}]$ . For the actual EOS calculations, we use only units of $\left[MeV\right]$ such that the conversion is

n_{B}[MeV^{3}]=197.3^{3}\cdot n_{B}[fm^{-3}].

(142)

Speed of sound For the speed of sound squared, $c_{s}^{2}$ , the table had units of $[cm^{2}/s^{2}]$ and we use natural units such that the speed of light is $c=1$ and $c_{s}^{2}$ is in terms of $c$ :

c_{s}^{2}=\frac{c^{2}_{s,table}\left[cm^{2}\;s^{-2}\right]}{c_{cgs}^{2}\left[% cm^{2}\;s^{-2}\right]},

(143)

where $c_{s}^{2}$ is the speed of sound squared in terms of $c$ , $c^{2}_{s,table}$ are the values from the table in cgs units of $\left[cm^{2}\;s^{-2}\right]$ , and $c_{cgs}^{2}$ is the speed of light in cgs units of $\left[cm^{2}\;s^{-2}\right]$ .

Pressure The pressure is provided in units of $[dyn/cm^{2}]$ and also in log scale of $log_{10}$ . To convert into our needed units of $\left[MeV^{4}\right]$

p\left[MeV^{4}\right]=F*10^{p_{table}}\;\left[\frac{dyn}{cm^{2}}\right],

(144)

where $p$ is the pressure in $\left[MeV^{4}\right]$ , $10^{p_{table}}$ is the pressure from the table in units of $\left[\frac{dyn}{cm^{2}}\right]$ , and $F$ is a conversion factor to convert the units, in this case $F=4.79363*10^{-27}\;\left[\frac{MeV^{4}cm^{2}}{dyn}\right]$ .

Energy The table provides the internal energy $E$ in cgs units of $\left[erg/g\right]$ wherein

\left[\frac{erg}{g}\right]=\left[\frac{dyn\cdot cm}{g}\right]=\left[\frac{cm^{% 2}}{s^{2}}\right].

(145)

In cgs units, the variables are not normalized by the speed of light $c$ such that we need to divide by the correct dimensions of $c$ in cgs units. Additionally, from the table the energy density provided is in log scale of $log_{10}$ . The internal energy from the table does not include the rest mass, $m_{N}$ , whereas for our calculations we require the total energy density $\varepsilon$ that does include the rest mass. Another change is that the values stored in the original table are shifted by a value in units of $\left[erg/g\right]$ . This energy shift is used to ensure that $log_{10}\varepsilon\;[erg/g]$ is positive and is a constant value of $1.91312955\cdot 10^{19}\;[erg/g]$ . To summarize, we can put together all of these effects to get the desired energy density in terms of $\left[MeV^{4}\right]$ :

\varepsilon\left[MeV^{4}\right]=m_{N}\left[MeV\right]*n_{B}\left[MeV^{3}\right% ]\\ \left(1+\frac{10^{E}\left[\frac{erg}{g}\right]-y\left[\frac{erg}{g}\right]}{c^% {2}\left[\frac{cm^{2}}{s^{2}}\right]}\right).

(146)

Here $\varepsilon$ is the energy density in units of $\left[MeV^{4}\right]$ and $m_{N}\left[MeV\right]$ is the mass of the nucleons that we defined previously.

Chemical potentials The table provides information about the microscopic chemical potentials, i.e., for a specific species (either electrons, protons, or neutrons) whereas we require the chemical potentials of conserved charges. Additionally, the table uses chemical potentials that are shifted by their mass that we denote as $\tilde{\mu}$ :

\mu\left[MeV\right]=\tilde{\mu}\left[MeV\right]+m_{N}\left[MeV\right],

(147)

where both $\mu$ and $\tilde{\mu}$ are in MeV. In order to get the electric charge chemical potential, we can simply use

\mu_{Q}=\mu_{p}-\mu_{n},

(148)

and the baryon chemical potential is

\mu_{B}=\mu_{n}.

(149)

From this point forward, we assume that all thermodynamic quantities are in units of $\left[MeV\right]$ within the code itself. Then only at the very end when the entire 3D is finalized, we convert back to the units required for numerical relativity simulations using the inverse of the formulas discussed here.

Finite-temperature expansion of the dense-matter equation of state

Abstract

I Introduction

I.1 Executive Summary

II Theoretical Framework

II.1 Removing lepton contributions

II.1.1 EOS contributions at T=0𝑇0T=0italic_T = 0

II.1.2 EOS contributions for T>0𝑇0T>0italic_T > 0

II.2 Going away from β𝛽\betaitalic_β-equilibrium

II.3 Finite T𝑇Titalic_T expansion

II.3.1 Expansion in T𝑇Titalic_T for fixed μ→→𝜇\vec{\mu}over→ start_ARG italic_μ end_ARG

II.3.2 Connecting d⁢s/d⁢T𝑑𝑠𝑑𝑇ds/dTitalic_d italic_s / italic_d italic_T at different YQsubscript𝑌𝑄Y_{Q}italic_Y start_POSTSUBSCRIPT italic_Q end_POSTSUBSCRIPT slices

II.4 Obtaining other thermodynamic quantities from p⁢(T,μ→)𝑝𝑇→𝜇p(T,\vec{\mu})italic_p ( italic_T , over→ start_ARG italic_μ end_ARG )

II.5 Crust

III Numerical Implementation

III.1 Numerical algorithm and flow chart

III.2 Crust

III.2.1 Finding β𝛽\betaitalic_β-equilibrium

III.3 Crust-core matching

III.3.1 Adding structure into cs2superscriptsubscript𝑐𝑠2c_{s}^{2}italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT

III.3.2 EOS reconstruction from a second-order Runge Kutta

III.3.3 Interpolated EOS and saturation properties

III.4 Finite T expansion

III.4.1 Terms at 𝒪⁢(T2)𝒪superscript𝑇2\mathcal{O}(T^{2})caligraphic_O ( italic_T start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) and 𝒪⁢(T3)𝒪superscript𝑇3\mathcal{O}(T^{3})caligraphic_O ( italic_T start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT )

III.4.2 Proof-of-principles for HIC with YQHIC=0.4superscriptsubscript𝑌𝑄HIC0.4Y_{Q}^{\rm HIC}=0.4italic_Y start_POSTSUBSCRIPT italic_Q end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_HIC end_POSTSUPERSCRIPT = 0.4

III.4.3 s/nB𝑠subscript𝑛𝐵s/n_{B}italic_s / italic_n start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT expansion in YQsubscript𝑌𝑄Y_{Q}italic_Y start_POSTSUBSCRIPT italic_Q end_POSTSUBSCRIPT

III.4.4 Conversion between a fixed grid in μ→→𝜇\vec{\mu}over→ start_ARG italic_μ end_ARG vs a fixed grid in nB,YQsubscript𝑛𝐵subscript𝑌𝑄n_{B},Y_{Q}italic_n start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT , italic_Y start_POSTSUBSCRIPT italic_Q end_POSTSUBSCRIPT

IV Error Quantification and Future Improvements

IV.1 Numerical Error

IV.2 Expansion Error

IV.3 Uncertainty in coefficients

IV.4 Final error quantification

V Conclusions and Outlook

VI Acknowledgements

References

Appendix A Converting between different thermodynamical bases

Appendix B Other thermodynamic observables for the finite T𝑇Titalic_T expansion

Appendix C YQsubscript𝑌𝑄Y_{Q}italic_Y start_POSTSUBSCRIPT italic_Q end_POSTSUBSCRIPT expansion of entropy derivatives

Appendix D Unit conversion

II.1.1 EOS contributions at $T=0$

II.1.2 EOS contributions for $T>0$

II.2 Going away from $\beta$ -equilibrium

II.3 Finite $T$ expansion

II.3.1 Expansion in $T$ for fixed $\vec{\mu}$

II.3.2 Connecting $ds/dT$ at different $Y_{Q}$ slices

II.4 Obtaining other thermodynamic quantities from $p(T,\vec{\mu})$

III.2.1 Finding $\beta$ -equilibrium

III.3.1 Adding structure into $c_{s}^{2}$

III.4.1 Terms at $\mathcal{O}(T^{2})$ and $\mathcal{O}(T^{3})$

III.4.2 Proof-of-principles for HIC with $Y_{Q}^{\rm HIC}=0.4$

III.4.3 $s/n_{B}$ expansion in $Y_{Q}$

III.4.4 Conversion between a fixed grid in $\vec{\mu}$ vs a fixed grid in $n_{B},Y_{Q}$

Appendix B Other thermodynamic observables for the finite $T$ expansion

Appendix C $Y_{Q}$ expansion of entropy derivatives