#  Simulating neutron stars with a flexible enthalpy-based equation of state
parametrization in SpECTRE

Isaac Legred  Department of Physics, California Institute of Technology,
Pasadena, California 91125, USA  LIGO Laboratory, California Institute of
Technology, Pasadena, California 91125, USA  Yoonsoo Kim  Department of
Physics, California Institute of Technology, Pasadena, California 91125, USA
Theoretical Astrophysics 350-17, California Institute of Technology, Pasadena,
CA 91125, USA  Nils Deppe  Theoretical Astrophysics 350-17, California
Institute of Technology, Pasadena, CA 91125, USA  Katerina Chatziioannou
Department of Physics, California Institute of Technology, Pasadena,
California 91125, USA  LIGO Laboratory, California Institute of Technology,
Pasadena, California 91125, USA  Francois Foucart  Department of Physics &
Astronomy, University of New Hampshire, 9 Library Way, Durham NH 03824, USA
François Hébert  Theoretical Astrophysics 350-17, California Institute of
Technology, Pasadena, CA 91125, USA  Lawrence E. Kidder  Cornell Center for
Astrophysics and Planetary Science, Cornell University, Ithaca, New York
14853, USA

###### Abstract

Numerical simulations of neutron star mergers represent an essential step
toward interpreting the full complexity of multimessenger observations and
constraining the properties of supranuclear matter. Currently, simulations are
limited by an array of factors, including computational performance and input
physics uncertainties, such as the neutron star equation of state. In this
work, we expand the range of nuclear phenomenology efficiently available to
simulations by introducing a new analytic parametrization of cold, beta-
equilibrated matter that is based on the relativistic enthalpy. We show that
the new _enthalpy parametrization_ can capture a range of nuclear behavior,
including strong phase transitions. We implement the enthalpy parametrization
in the SpECTRE code, simulate isolated neutron stars, and compare performance
to the commonly used spectral and polytropic parametrizations. We find
comparable computational performance for nuclear models that are well
represented by either parametrization, such as simple hadronic EoSs. We show
that the enthalpy parametrization further allows us to simulate more
complicated hadronic models or models with phase transitions that are
inaccessible to current parametrizations.

##  I Introduction

Multimessenger observations of the gravitational wave event GW170817 Abbott
_et al._ (2017a, b) have highlighted the role of neutron star binaries (BNS)
in probing the physics of dense matter, e.g., Dietrich  _et al._ (2021);
Chatziioannou (2020); Piekarewicz (2022). In addition, further astronomical
observations Abbott  _et al._ (2017a, 2020); Cromartie  _et al._ (2019);
Fonseca  _et al._ (2021); Miller  _et al._ (2019, 2021); Riley  _et al._
(2019, 2021); Antoniadis  _et al._ (2013) and terrestrial nuclear experiments
Adhikari  _et al._ (2021, 2022) have facilitated new insights into the
equation of state (EoS) of NS matter Abbott  _et al._ (2018); Essick (2021);
Pang  _et al._ (2021, 2021); Raaijmakers  _et al._ (2020, 2021); Landry  _et
al._ (2020); Legred  _et al._ (2021); Miller  _et al._ (2019, 2021).
Nonetheless significant uncertainty exists about the properties of dense
matter above nuclear saturation density111The saturation density of atomic
nuclei is determined via theory and experiments Drischler  _et al._ (2016);
here we fix a value for convenience., $\rho_{\mathrm{nuc}}\equiv 2.8\times
10^{14}\mathrm{g}/\mathrm{cm}^{3}$, which translates to uncertainty in the
properties of astrophysical NSs whose densities can reach $\sim
7\rho_{\mathrm{nuc}}\leavevmode\nobreak\ $Pang  _et al._ (2021); Legred  _et
al._ (2021).

The merger phase of a BNS coalescence carries the largest imprint of nuclear
matter and strong gravity and it can only be studied numerically. Numerical
relativity (NR) simulations of BNS coalescences through merger require solving
the equations of general relativistic magnetohydrodynamics (GRMHD)
simultaneously with the Einstein field equations and, possibly, the Boltzmann
equations for neutrino radiation transport Antón  _et al._ (2006); Font
(2008); Baumgarte and Shapiro (2010). The system of equations is closed with a
nuclear EoS. See e.g. Baiotti and Rezzolla (2017); Radice  _et al._ (2020);
Foucart (2020); Kyutoku  _et al._ (2021); Rezzolla and Zanotti (2013) for
reviews of the field. Such simulations have been used to interpret existing
signals, e.g. Margalit and Metzger (2017); Radice and Dai (2019); Shibata  _et
al._ (2019); Köppel  _et al._ (2019); Annala  _et al._ (2022); Camilletti  _et
al._ (2022), and targeted simulations will likely be an essential tool for
understanding future observations.

The most generic strategy for representing the nuclear EoS numerically is
piecewise, i.e., using independent expressions in different density or
pressure intervals. For example, interpolated tables of thermodynamic
quantities such as pressure and internal energy at every value of the density
and composition offer access to the widest range of nuclear behavior. However,
the temperature- and composition-dependent tables currently used, e.g. Typel
_et al._ (2022), have a significant memory footprint and evaluation requires
computationally expensive operations such as constant access to the table and
interpolation Siegel  _et al._ (2018). The latter may also be inaccurate (at
low order) or prone to unphysical oscillations for EoSs with discontinuities
or underresolved features (at high order). A related approach makes use of
piecewise parametrizations such as a _piecewise-polytrope_ Read  _et al._
(2009a), which is effectively a sparsely sampled table for the polytropic
exponent. Though it can capture a range of high-density behavior,
discontinuities in derivatives of thermodynamic quantities can degrade
simulation accuracy Foucart  _et al._ (2019); Raithel and Paschalidis (2022).

A different strategy is based on functional representations of the EoS that
stay smooth across density scales, such as a _single-polytropic_ or _spectral_
parametrization Lindblom (2010); Greif  _et al._ (2019); Foucart  _et al._
(2019); Raithel and Paschalidis (2022). Such parametrizations typically cannot
fully represent nuclear EoS models, as they are restricted to a finite number
of parameters in the density range of interest Greif  _et al._ (2019); Foucart
_et al._ (2019). On top of this, smoothness across density scales fails to
capture nuclear models that contain nuclear transitions to exotic degrees of
freedom.

In this study, we propose a new parametrization of the nuclear EoS that
bridges smooth and discontinuous models while balancing accuracy and
computational efficiency.222We use the term “model” to refer to a nuclear-
theoretic prediction and “parametrization” for a functional form for the EoS.
We parametrize the relativistic enthalpy Lindblom (1992) via a combination of
analytic polynomials and trigonometric functions. Unlike pressure, the
difference in enthalpy at densities
$[\rho_{\mathrm{nuc}},3\rho_{\mathrm{nuc}}]$ for two EoSs is typically small
compared to the enthalpy of either. The enthalpy can thus be effectively
written as a “baseline” part plus small corrections. We capitalize on this in
order to write the enthalpy as a polynomial, typically capturing $\sim 99\%$
of the EoS, plus small trigonometric corrections, bringing the fit accuracy to
1 in $10^{5}$. Such a decomposition can capture a wide range of phenomenology
with modest changes to the relevant parameters. In addition, further
thermodynamic quantities such as the pressure can be evaluated efficiently and
analytically.

We implement this parametrization in SpECTRE Deppe  _et al._ (2022a); Kidder
_et al._ (2017), a scalable next-generation multiphysics computational
astrophysics code that uses task-based parallelism Kale  _et al._ (2020). A
primary science target for SpECTRE is fast and accurate GRMHD simulations of
BNS coalescences. We use SpECTRE to test the _enthalpy parametrization_ on
isolated NSs in the Cowling approximation, i.e. we do not evolve the spacetime
Cowling (1941), while evolving the ideal GRMHD equations Baumgarte and Shapiro
(2010) with a discontinuous Galerkin-finite difference (DG-FD) hybrid scheme
Deppe  _et al._ (2022b, c). Though these simulations assume a static
spacetime, they still allow us to evaluate the role of the enthalpy
parametrization in questions of convergence, efficiency, and resolvability of
nuclear physics in simulations.

We show that the enthalpy parametrization is able to effectively represent a
wide range of nuclear behavior, while incurring small additional computational
costs relative to simpler parametrizations. After reviewing the general
requirements a parametrization must meet in Sec. II, we introduce the enthalpy
parametrization in Sec. III. We demonstrate that it can faithfully fit various
nuclear models ranging from smooth EoSs to phase transitions in Sec. IV. We
perform numerical simulations with SpECTRE and find that for resolutions of at
least $130$ m, the EoS evaluation cost is subdominant to other simulation
components. We also simulate hybrid stars with quark cores and find that such
simulations can be carried out stably with better-than-expected runtime
scaling properties under increasing resolution. We conclude with discussions
in Sec. V.

##  II EoS Parametrizations for Relativistic Simulations

###  II.1 General requirements

We begin with a general discussion of the requirements phenomenological
parametrizations of the nuclear EoS must meet for efficient use in numerical
simulations. These include (i) faithful representation of target nuclear
models, (ii) parametric extensibility, and (iii) computational performance
related to smoothness (to the extent allowed by the underlying nuclear
physics) and/or a fully analytic formalism.

The first requirement is that the parametrization is generic enough that it
can faithfully represent the target nuclear physics. While no standard
faithfulness metrics exist, a common test is the $L^{2}$ difference of
quantities of interest Lindblom (2010); Read  _et al._ (2009a). Nonetheless it
is unclear how different metrics relate, for example the $L^{2}$ difference of
the local polytropic indices and that of the mass-radius curve Lindblom and
Indik (2014); Foucart  _et al._ (2019). One particular challenge to smooth
parametrizations is modeling strong phase transitions Han and Steiner (2019);
Pang  _et al._ (2020). In general we would like a parametrization where,
whatever the metric, we can improve the fit via iterative approximation. In
principle this is available to any parametrization by adding more parameters
and smoothly changing parameter values. In practice, however, the functional
form of the parametrization may limit the accessible parameter space, as shown
in Wysocki  _et al._ (2020) for the spectral parametrization.

The second requirement is that the parametrization allows us to parametrically
explore a wide range of possible high-density behavior. This entails
continuously, and without significant fine-tuning, extending the
parametrization to produce EoSs that might differ from existing nuclear
models. An example of such an extension would be a parameter which controls
the pressure at a particular density and thus allows us to isolate the effect
of this density scale on macroscopic observables. Another benefit of such
continuous extensibility is that it allows us to construct a map from the EoS
to observables, e.g. Özel and Psaltis (2009). This approach has already been
successfully used in the case of binary black hole mergers to produce accurate
_surrogates_ of the map of binary configurations to gravitational waves
Blackman  _et al._ (2015); Varma  _et al._ (2019). A similar methodology could
be used to construct a surrogate for the post-merger gravitational-wave
signature of BNS mergers, whose EoS dependence is not well captured by a small
number of parameters Wijngaarden  _et al._ (2022); Breschi  _et al._ (2022).

At the same time, we consider practical requirements in terms of computational
performance: speed and accuracy of the relevant evaluations, and smoothness of
the thermodynamic quantities where possible. A fully analytic form for the EoS
and all the relevant thermodynamic quantities is a sufficient (but perhaps not
necessary) condition. Tabulated EoSs, while guaranteeing maximal flexibility,
fail in this regard. Consider, for example, primitive variable recovery.
Numerical simulations evolve the components of the stress-energy tensor which
are nonlinear functions of primitive variables such the rest-mass baryon
density $\rho$, pressure $p$, and specific internal energy $\epsilon$. This
process involves inverting the relation between the stress-energy tensor and
the primitive variables with root-finding routines during which the EoS, for
example $p(\epsilon)$, is evaluated repeatedly. For tabulated EoSs this
includes computing the temperature $T$ from $\epsilon$ via another root-find
and then computing $p(T)$ via a table lookup and interpolation. The EoS tables
are typically too large to store in the CPU caches and so the nested root-
finding routines require repeated loading of data from main memory, causing
significant overhead that dominates simulation cost Siegel  _et al._ (2018).

Another advantage of fully analytical parametrizations is that they enable
efficient computation of all necessary thermodynamic quantities in a
consistent way. Besides tabulated EoSs, this also applies to certain
parametrizations that require interpolation or numerical integration. For
example, the spectral parametrization allows for analytic evaluation but not
integration of $d\epsilon/d\rho$. Then, $\epsilon(\rho)$ is computed via a
computationally expensive numerical integral as high accuracy is required to
avoid thermodynamic inconsistency during primitive variable recovery. Even if
tables are used in simulations, ensuring smoothness and consistency requires
building higher-order interpolants (or sampling very densely). This
effectively amounts to constructing local parametrizations of the EoS which
satisfy some stitching constraints. Therefore, even the use of tables in NR
simulations stands to gain from understanding fully analytic representations
of the nuclear EoS.

###  II.2 Existing parametrizations of the EoS

The simplest parametrization of cold, beta-equilibrated, dense matter is a
single polytrope that prescribes a relationship between the rest-mass baryon
density $\rho$ and the pressure $p$

$p(\rho)=K\rho^{\Gamma}\,,$ (1)

where $\Gamma$ is the the polytropic exponent and $K$ is the polytropic
constant; both are independent of $\rho$. For example, a degenerate neutron
gas would obey a polytropic relation with $\Gamma=5/3$. Polytropes have a long
history in NS simulations, e.g., Shibata and Uryu (2000); Etienne  _et al._
(2008); Baiotti  _et al._ (2008); Duez  _et al._ (2008), and more recent code
tests, e.g., Radice  _et al._ (2014); DeBuhr  _et al._ (2018); Deppe  _et al._
(2022c), due to their simplicity, low computational cost, and the fact that
they allow for analytic evaluation of pressure, internal energy, specific
enthalpy, and rest-mass density. Nonetheless, their simplicity makes
polytropes incompatible with realistic EoS nuclear models, either hadronic
(for example, polytropes do not satisfy the same universal relations as
hadronic models Yagi and Yunes (2013)) or hybrid ones that include multiple
degrees of freedom.

Piecewise-polytropes Read  _et al._ (2009b) extend single-polytropes to
multiple polytropic segments at different densities, thereby decoupling low-
and high-density behavior. With enough piecewise segments, piecewise-
polytropes can also fit EoSs with strong phase transitions Ujevic  _et al._
(2022). While piecewise-polytropes retain some of the computational simplicity
of the single-polytrope and have been employed in BNS mergers Hotokezaka  _et
al._ (2011); Lackey  _et al._ (2014); Dietrich  _et al._ (2018, 2017), the
lack of smoothness across stitching boundaries tends to increase the
computational cost and decrease the accuracy Foucart  _et al._ (2019); Raithel
and Paschalidis (2022). Extensions to continuous polytropic indices O’Boyle
_et al._ (2020); Raithel and Paschalidis (2022) guarantee differentiability of
the pressure; however, it is unclear how to extend the parametrization to
guarantee further derivatives of the pressure exist at the stitching point.
Generically stitching two $\mathcal{C}^{n}$ functions to form a globally
$\mathcal{C}^{n}$ function requires matching $n+1$ derivatives, which may
require the introduction of functions to the parametrization of $p(\rho)$ for
example, which make it difficult to solve for $e(\rho)$ analytically.

Finally, the spectral parametrization Lindblom (2010) accurately reflects a
broad range of nuclear models while maintaining smoothness across density
scales. The parametrization has a similar form to a polytrope

$p(\rho)=K\rho^{\Gamma(\rho)}\,,$ (2)

but now $\Gamma(\rho)$ is expanded in a basis of smooth functions, typically a
polynomial. The spectral parametrization can successfully reproduce hadronic
nuclear models with a comparable number of parameters as polytropes, though it
cannot capture sharp changes in the speed of sound that are associated with
phase transitions Lindblom (2010); Foucart  _et al._ (2019). Compared to
piecewise polytropes and other EoS with discontinuities, the spectral
parametrization can lead to reduced computational cost in simulations Foucart
_et al._ (2019) for a given accuracy requirement, while remaining more
computationally intensive than pure polytropes. Our current implementation of
the spectral EoS balances faithfulness to nuclear models and computational
efficiency by expressing $\Gamma(\rho)$ as a a polynomial in $\log\rho$
Foucart  _et al._ (2019). More complex basis functions could improve
faithfulness, but they would come at the cost of computational efficiency
since computation of the internal energy requires a numeric integral whose
accuracy depends on how rapidly $\Gamma(\rho)$ varies.

The above discussion highlights the role of balancing faithfulness and
computational efficiency in selecting EoS parametrizations for numerical
simulations. While the single-polytrope is computationally efficient, it is
too restrictive in terms of nuclear physics. Piecewise-polytropes expand the
range of nuclear models accessible, but at the cost of longer runtimes and
loss of accuracy due to non-smoothness at the stitching boundaries. The
spectral parametrization strikes some balance, but performs optimally when few
parameters are used; it is therefore restricted to simple nuclear models.
Ultimately, we would prefer an EoS parametrization which is able to fit to a
problem-specific precision, matching the level of other errors in simulations
at the lowest possible cost. This motivates the introduction of a new
parametrization with increased flexibility to model a wider range of nuclear
EoSs without considerable performance losses.

##  III enthalpy parametrization of the EoS

In this section we introduce a new _enthalpy parametrization_ with a flexible
number of degrees of freedom that expands the range of microscopic physics we
are able to represent in numerical simulations. In the following, we work in
geometric units: $c=1$, $G=1$.

###  III.1 Parametrizing the enthalpy

The specific enthalpy of a system $h$ is defined as the enthalpy per unit
mass. In relativistic contexts it represents the energy required to inject a
unit of rest mass into the system while remaining in thermodynamic
equilibrium. The first law of thermodynamics requires that at zero temperature
$T$ and in $\beta-$equilibrium,

$h(\rho)\equiv\left(\frac{\partial
e}{\partial\rho}\right)_{T,\beta}=\frac{de}{d\rho}=\frac{p(\rho)+e(\rho)}{\rho}\,,$
(3)

where $e$ and $p$ are the energy density and pressure, while $\rho$ is the
rest-mass energy density of baryons.

We choose to directly parameterize the enthalpy for three primary reasons.
First, the enthalpy is a monotonically-increasing and slowly-varying function
of the baryon density, which is numerically beneficial. Second, the enthalpy
can be intuitively interpreted as a measure of the stiffness of the EoS: a
larger enthalpy corresponds to higher pressure and energy density. Third, and
importantly for hydrodynamic simulations, the enthalpy in cold, beta-
equibrilated matter is related to other thermodynamic quantities by linear
operations, which facilitates analytic calculations and avoids interpolation
or numerical integration.

From the first law, we have

$\displaystyle\frac{dh}{d\log\rho}$
$\displaystyle=\frac{de}{d\rho}+\frac{dp}{d\rho}-h$ (4)
$\displaystyle=\frac{dp}{d\rho}=\frac{dp}{de}\frac{de}{d\rho}$ (5)
$\displaystyle=hc_{s}^{2}\,,$ (6)

Equation (5) suggests that $dp/d\rho$ is zero if and only if $dh/d\rho$ is
zero. Equation (6) provides the motivation for our parametrization choices.
Consider, for example, a constant speed of sound $c_{s}=c_{s,0}$. Then

$c_{s,0}^{2}=c_{s}^{2}=\frac{dp}{de}\implies p=p_{0}+c_{s,0}^{2}\Delta e\,,$
(7)

with $p_{0}=p(e_{0})$ and $\Delta e\equiv e-e_{0}$. In this special case Eq.
(6) becomes

$\displaystyle h(\log\rho)$
$\displaystyle\propto\exp\left(c_{s,0}^{2}\log\rho\right)$
$\displaystyle\approx\rho_{0}\left[1+c_{s,0}^{2}\log\left(\rho/\rho_{0}\right)+\dots\right]\,,$
(8)

where $\rho_{0}$ is some fiducial density. Equation (III.1) suggests that if
$c_{s}^{2}$ is slowly varying,333In general, causality and stability bound
$0\leq c_{s}^{2}\leq 1$. the enthalpy can be approximated as exponential in
$\log\rho$. Moreover, a smaller $c_{s}^{2}$ accelerates the convergence of the
series of Eq. (III.1), though this also depends on the choice of density scale
$\rho_{0}$. We therefore choose the Taylor expansion in Eq. (III.1) as the
starting point of the enthalpy parametrization.

We further select $\log\rho/\rho_{0},$ as the independent variable of the
parametrization. This choice enables us to better resolve the low-density EoS.
Equation (III.1) further suggests that
$h(\log\rho)\propto\exp(c_{s,0}^{2}\log\rho)$ is analytically and
computationally simpler than $h(\rho)\propto\rho^{c_{s,0}^{2}}$ as the Taylor
expansion of $\rho^{c_{s,0}^{2}}$ converges more slowly than the expansion of
$\exp(c_{s,0}^{2}\log\rho)$ for non-integer $c_{s,0}$.

Lastly, a desirable property of the specific enthalpy is that it is continuous
across first-order phase transitions. This can be seen from Eq. (5) where
maintaining a constant pressure across the transition guarantees that the
enthalpy will be constant as well. This indicates that across certain weak
transitions the enthalpy can be expanded in a basis of continuous functions,
unlike, for example, a local polytropic exponent.

###  III.2 Decomposition

Motivated by Eq. (III.1), we introduce a parametrization of $h(\log\rho)$.
Given an EoS in some density region $\rho_{\min}\leq\rho\leq\rho_{\max}$ we
select a density scaling parameter $\rho_{0}\leq\rho_{\min}$ such that
$z\equiv\log(\rho/\rho_{0})$ is positive in the relevant density range.
Importantly, $\rho_{0}$ is not necessarily equal to $\rho_{\min}$ thus
introducing an additional parameter. We then write

$h(z)\approx h_{p}(z)\equiv\sum_{i=0}^{i_{\max}}\gamma_{i}z^{i}\,,$ (9)

where $h(z)$ is the target enthalpy and $h_{p}(z)$ is its approximation. This
polynomial decomposition is motivated by the previous observation that $h(z)$
is approximately exponential in $z$ for nearly constant speeds of sound,
corresponding to $\gamma_{i}\sim c_{s,0}^{2i}/i!$. The rapid convergence of
the $\gamma_{i}$ sequence indicates that the $i>i_{\max}$ terms will be small
provided that the speed of sound is slowly varying.

Given that $h(z)$ is positive and increasing and $z_{i}>0$, catastrophic
floating point cancellation in numerical calculations can be avoided by
restricting to $\gamma_{i}\geq 0$. This guarantees that each term
$\gamma_{i}z^{i}$ is a small and positive correction to previous terms.
Furthermore, the polynomial expansion of Eq. (9) can be efficiently and stably
evaluated with Horner’s method Press  _et al._ (2007). Allowing for more
general $\gamma_{i}$ is possible, but this comes at a risk of oscillatory
behavior and cancellation of large terms which make convergence predictions
difficult. The implications and rationale behind the choice to set
$\gamma_{i}\geq 0$ are further discussed in App. C.

A consequence of setting $\gamma_{i}\geq 0$ is that Eq. (9) is unable to model
certain EoSs, for example the case where $dh/dz=hc_{s}^{2}$ is not strictly
increasing, even with $i_{\max}\rightarrow\infty$. Such a non-monotonic speed
of sound could be encountered for complicated hadronic models or more
generically if non-hadronic degrees of freedom are introduced McLerran and
Reddy (2019); Tews  _et al._ (2018a); Kapusta and Welle (2021); Han and
Steiner (2019). We therefore augment Eq. (9) by decomposing $h_{t}(z)\approx
h(z)-h_{p}(z)$ as a Fourier series

$h_{t}(z)\equiv\sum_{j=1}^{j_{\max}}a_{j}\sin(jkz)+b_{j}\cos(jkz)\,,$ (10)

where $k$ sets the “wavelength scale” of the fit. In a Fourier series $k$ is
typically fixed to

$k=k_{F}\equiv\frac{2\pi}{z_{\max}-z_{\min}}=\frac{2\pi}{\log(\rho_{\max}/\rho_{\min})}\,,$
(11)

but here we vary it and find that $k\gtrapprox k_{F}$ leads to good fits. The
effect of perturbing $k$ around $k_{F}$ is small, as we explore in App. A. The
trigonometric expansion of Eq. (10) can also serve as a low-pass filter to
remove high-frequency oscillations from the tabulated EoS data that may not be
physical or computationally resolvable. In summary, the enthalpy
parametrization is

$h_{*}(z)\equiv h_{t}(z)+h_{p}(z)\approx h(z)\,.$ (12)

Figure 1: Results of a fit to an EoS drawn from a Gaussian process with the
enthalpy parametrization. We plot various thermodynamic quantities as a
function of $z$. The fit parameters are
$\rho_{\min}=\rho_{\operatorname{nuc}}$, $\rho_{\max}=7\rho_{\mathrm{nuc}}$,
$\rho_{0}=0.5\rho_{\operatorname{nuc}}$, $k=\pi/(\log(7))$, and
$i_{\max}=j_{\max}=10$. _Top Panel:_ The tabulated EoS $h$ (solid, orange) and
the total fit $h_{*}$ (solid, light blue). We also plot the polynomial fit to
the EoS $h_{p}$ (dashed, indigo). Both the total and the polynomial fit are
indistinguishable from the tabulated EoS by eye. _Second Panel:_ The residuals
of the total fit $h-h_{*}$. In this metric, the fit demonstrates excellent
agreement relative to $h-1=p/\rho+\epsilon\gtrsim 1\times 10^{-2}$ _Third
Panel:_ The trigonometric fit $h_{r}=h-h_{p}$. _Fourth Panel:_
$(1/h)dh/dz=c_{s}^{2}$, for both the tabulated EoS and the total fit.
Heuristically, the speed of sound has a comparable number of plateaus to the
number of obvious peaks in $h_{t}$.

In Fig. 1 we demonstrate the enthalpy parametrization fit of Eq. (12) and its
polynomial, Eq. (9), and trigonometric, Eq. (10), components for a
phenomenological EoS drawn from a Gaussian process prior Landry and Kumar
(2018); Essick  _et al._ (2020a). The polynomial fit alone is accurate to
about $\mathcal{O}(1\%)$, while the total fit is good to about one part in
$10^{5}$. For reference, we also plot $c_{s}^{2}=(1/h)dh/dz$, as a measure of
the complexity of the EoS. Even though $c_{s}^{2}$ is not globally nearly
constant, it is slowly varying and nearly monotonic.

Given the generic form of the enthalpy parametrization, there is no guarantee
that a particular fit will satisfy stability $c_{s}^{2}\geq 0$ and causality
$c_{s}^{2}\leq 1$. If $h_{t}(z)=0$, the fit is guaranteed to be stable, and a
sufficient but not necessary condition for causality is
$\gamma_{i}\leq\gamma_{i-1}/i$, which becomes necessary and sufficient in the
case of a constant sound speed. If $h_{t}(z)$ is nonzero, then $h_{*}(z)$ can
oscillate, changing on scales of order the most quickly varying Fourier mode.
Therefore, both conditions must be checked on a grid of spacing

${\delta z\lesssim\frac{1}{j_{\max}k}\,,}$ (13)

where, as above, $j_{\max}$ is the index of the fastest varying “Fourier”
mode.

While an unstable fit to the EoS cannot be tolerated in a numerical
simulation, an acausal fit may be used if it is very nearly causal (i.e. if
$\sqrt{c_{s}^{2}-1}$ is small compared to the velocity resolution of the
simulation). In practice, however, fits typically are neither acausal nor
unstable; if they are it is often a sign that the fit to the EoS is poor and
more parameters should be used.

###  III.3 Computing thermodynamic quantities

Given the expansion of Eq. (12), we can analytically compute the thermodynamic
quantities needed for GRMHD evolution as formulated in SpECTRE Deppe  _et al._
(2022a, c), or similar codes Mösta  _et al._ (2014). For example, the energy
density is

$\displaystyle\frac{de}{dz}=\rho\frac{de}{d\rho}=\rho_{0}\exp\left(z\right)h(z)\Rightarrow$
$\displaystyle
e(z)=\rho_{0}\int_{z_{0}}^{z}\exp\left(z^{\prime}\right)h(z^{\prime})dz^{\prime}+e(z_{0})\,.$
(14)

Since $h(z)$ is expressed in terms of sines, cosines, and polynomials, the
integral of Eq. (III.3) can be computed using the following identities

$\displaystyle\int\exp\left(z\right)\sin(nkz)dz$
$\displaystyle=\exp\left(z\right)\frac{\sin(nkz)-nk\cos(nkz)}{1+n^{2}k^{2}}+C\,,$
(15) $\displaystyle\int\exp\left(z\right)\frac{z^{n}}{n!}dz$
$\displaystyle=\exp\left(z\right)\frac{z^{n}}{n!}-\int\frac{z^{n-1}}{(n-1)!}\exp\left(z\right)dz=\dots\,,$
(16)

where the ellipses indicate that integration by parts can be repeated until
the integral becomes trivial. Equation (III.3) is also a gamma function, but
it is typically incomplete. Nonetheless, all integrals can be evaluated
analytically and $e(z)$ has an expansion of the form

$\displaystyle e(z)$ $\displaystyle=\exp(z)\times$
$\displaystyle\left(\sum_{i}\gamma_{i}^{\prime}z^{i}+\sum_{j}a_{j}^{\prime}\sin(kjz)+b_{j}^{\prime}\cos(kjz)\right)$
$\displaystyle+e_{*}\,,$ (17)

where the constant $e_{*}$ is determined by setting $e(z_{\min})=e_{\min}$ and
the coefficients satisfy

$\displaystyle\gamma^{\prime}_{i}$
$\displaystyle=\frac{1}{i!}\sum_{i_{\max}\geq\ell\geq
i}(-1)^{i_{\max}-\ell}\ell!\gamma_{\ell}\,,$ (18) $\displaystyle
a^{\prime}_{j}$
$\displaystyle=\frac{a_{j}}{1+j^{2}k^{2}}+\frac{b_{j}jk}{1+j^{2}k^{2}}\,,$
(19) $\displaystyle b^{\prime}_{j}$
$\displaystyle=\frac{b_{j}}{1+j^{2}k^{2}}-\frac{a_{j}jk}{1+j^{2}k^{2}}\,.$
(20)

The pressure $p(z)$ can also be evaluated analytically with a similar
expansion given that

$p(z)=\rho_{0}h(z)\exp\left(z\right)-e(z)=h\rho-e\,.$ (21)

This equation showcases the benefits of setting $\gamma_{i}\geq 0$ in Eq. (9)
to avoid cancellations in the enthalpy expansion. The pressure is computed as
the difference of two relatively large quantities, each typically 1–3 orders
of magnitude larger than the pressure itself in the relevant density interval.
If the expansion of $h(z)$ additionally had large coefficients (i.e. much
larger than the enthalpy) terms of $e(z)$ will be computed by sums of
alternating large numbers, which is numerically undesirable. However, because
$\gamma_{\ell}\sim 1/\ell!$ for EoSs with slowly varying speed of sound, the
terms in Eq. (18) are of comparable size, and about the same size as
corresponding terms of $h(z)$. Thus the terms of $p(z)$ are computed to
comparable precision as the terms of $e(z)$ and $h(z)$. We find this holds
more broadly, even when the speed of sound is not slowly varying, as
$\gamma_{\ell}$ is typically decreasing even if it is not decreasing
exponentially as in the constant-$c_{s}^{2}$ case.

Lastly, we can also analytically compute

$\frac{dp}{d\rho}=\frac{dh}{dz}\,,$ (22)

through

$\displaystyle\frac{dh}{dz}$ $\displaystyle=\sum_{i}i\gamma_{i}z^{i-1}$
$\displaystyle+\sum_{j}jk\left[a_{j}\cos(jkz)-b_{j}\sin(jkz)\right]\,.$ (23)

As can be seen from Eqs. (III.3) and (23), parameters that enter linearly in
the original expansion of $h(z)$ also appear linearly in all relevant
thermodynamic quantities.

###  III.4 Low-Density Stitching

The enthalpy parametrization is best suited for high-density regions where
pressure and energy density are comparable. Low-density regions with $p\ll
h\rho\sim e$ might be better fit by direct parametrizations of the pressure.
We therefore combine the enthalpy parametrization with a simpler low-density
parametrization below $\rho_{\min}$. Incidentally, this density region
coincides with the region of validity of nuclear theory calculations Tews
(2020); Drischler  _et al._ (2016); Essick  _et al._ (2020b) and terrestrial
experiments Adhikari  _et al._ (2021, 2022); Roca-Maza  _et al._ (2015);
Essick  _et al._ (2021a, b). The low-density EoS is therefore better
constrained and thus there is reduced need for flexibility in the EoS
parametrization. Moreover, the low-density EoS has a reduced impact on NS
observables, especially if the simulation resolution is low, such that
$(dp/dr)\Delta r>\delta p$, where $\Delta r$ is the grid spacing and $\delta
p$ is the difference induced by EoS mismodeling.

A number of options exist for the low-density EoS, including direct
parametrizations of nuclear models Tews (2020) or chiral effective field
theory($\chi$-EFT) results Tews  _et al._ (2018b); Pang  _et al._ (2021). Here
we select the existing spectral parametrization implementation Foucart  _et
al._ (2019), as it is more flexible than single-polytropes, but smoother than
piecewise-polytropes and tabulated EoSs. In certain cases, we explore
extending the spectral parametrization up to relatively high densities $\sim
2\rho_{\mathrm{nuc}}$ if it can fit the target EoS well-enough in this density
regime. Due to the low number of parameters in the spectral parametrization,
all degrees of freedom are determined by requiring differentiability of the
pressure and continuity of the internal energy at the stitching points. We
verify this stitching maintains $\mathcal{C}^{1}$ smoothness, see App. A.

###  III.5 Free parameters and fitting

The number of free parameters needed to achieve good fits of arbitrary EoSs
will impact the simulation cost. Indeed, the cost of evaluating any EoS-
dependent quantity is proportional to the number of coefficients used for the
enthalpy parametrization. Therefore it is prudent to use only as many terms as
necessary to achieve an accurate fit; accuracy in the context of numerical
simulations is measured relative to other simulation errors. There is no
definitive metric for EoS mismodeling error, as the relevant error will depend
on the application. For example, in applications to BNS inspirals, the
relevant errors are in GW phase, matter hydrodynamic variables, and magnetic
field variables. When considering the fitness of an EoS parametrization for
use in simulations, all these factors should be taken into consideration.

Nonetheless, it is pragmatically necessary to define surrogate goodness-of-fit
statistics in order to both fit the enthalpy parametrization to data and
determine approximately if such a fit is good. We describe the fitting
procedure of the parametrization to a tabulated model that we employ in App.
A. Briefly, we fit the specific enthalpy $h(z)$ on a linear grid in $z$ but
with variable precision, requiring higher precision at lower densities to
achieve equal cost across density scales. However, fitting is not the only way
to extract coefficients for use in the enthalpy parametrization; for example,
coefficients to approximate a polytropic EoS are derived in App. B using a
Taylor expansion of the specific enthalpy. Nonetheless, for realistic nuclear
models, fitting the specific enthalpy is usually necessary.

One convenient benchmark is to examine the error in radius of a typical
neutron star induced by using a enthalpy parametrization fit as compared to a
tabulated model. In order to demonstrate the general requirements for fitting,
we fit a collection of realistic nuclear theoretic EoSs, compute the error in
the radius of a $1.4M_{\odot}$ NS ($\Delta R_{\mathrm{typ}}$) and display the
results in Table 1. These fits are all carried out with $i_{\max}=12$, and
$j_{\max}\leq 5$, and have typical NS radius error of less than $70$ m. EoS
modelling error would therefore likely not be limiting in simulations with
$\sim 70$ meter resolution; this is a conservative choice of error measure as
realistic simulation errors will likely dominate static errors. Additionally,
we fit phenomenolgical EoSs drawn from a Gaussian process-mixture model priors
Landry and Essick (2019); Essick  _et al._ (2020a). We examine two cases,
first are draws from a model-agnostic prior, which are only loosely informed
by nuclear theory calculations. The second class are Gaussian process draws
conditioned on $\chi$-EFT up to $1.5\rho_{\mathrm{nuc}}$ Essick (2021); Essick
_et al._ (2021b). Both nuclear-theoretic and phenomenological EoSs show
comparable fit quality, indicating that the enthalpy parametrization is able
to reproduce a wide range of EoS models.

We list $j_{\max}$ in Table 1 as we expect that the number of trigonometric
terms is the leading-order driver of cost to evaluate the parametrization. We
quantify this further in Sec. IV.1.2. Contrarily we expect little dependence
of evaluation cost of $i_{\max}$ because evaluation of polynomials using
Horner’s method is extremely efficient. For realistic EoSs, fine-tuning of
low-density stitching and nonlinear parameters can reduce the number of
trigonometric correction terms that are required to achieve a good fit. Even
when no fine-tuning is required, typically good fits are achieved with
$j_{\max}\sim 4$. We quantify this in Fig. 2 by showing the error in the
radius of a typical star for 6 different $\chi$-EFT informed Gaussian process
draws, when no fine tuning of nonlinear or low-density parameters is
performed. The fits are better when more trigonometric correction terms are
included, all falling below $100$ m error by $j_{\max}=4$. These errors are
often due to the EoS at low-densities, and so typically fine-tuning of certain
parameters, such as the low density polytropic index, or the energy density of
EoS at the stitching density, must be carried out to achieve $\sim 10$-meter-
error fits. In practice, though, this may not be necessary as quantities such
as the tidal deformability are determined by the bulk of the matter, interior
to the crust, therefore crust modeling errors may be less significant then
predicted by using the radius as a metric. These considerations will be
especially important for BNS simulations where GW emission is predominately
determined by tidal deformability, and other sources of error may overshadow
EoS modeling error.

Table 1: A list of EoS fits with the enthalpy parametrization to nuclear theoretic and phenomenological EoS. Theoretic EoSs are listed according to the conventions of Read  _et al._ (2009a). Phenemenological EoSs are drawn from Gaussian process priors. The EoSs gp1 and gp2 are drawn from a model agnostic Gaussian process prior Landry and Essick (2019); Essick  _et al._ (2020a). EoSs gp$\chi$eft1, gp$\chi$eft3, and gp$\chi$eft5 are drawn from Gaussian process priors conditioned on $\chi$-EFT predictions at low densities. These three EoSs represent draws from hadronic, hyperonic, and quarkyonic conditioned GPs, respectively. EoS | $R_{\mathrm{typ}}\ [\mathrm{km}]$ | $\Delta R_{\mathrm{typ}}\ [\mathrm{km}]$ | $\Delta M_{\max}\ [M_{\odot}]$ | $j_{\max}$ | Ref.
---|---|---|---|---|---
alf2 | 12.968 | -0.028 | -0.003 | 3 | Alford  _et al._ (2005)
bsk19 | 10.763 | -0.006 | -0.001 | 5 | Potekhin  _et al._ (2013)
ap4 | 10.595 | -0.02 | 0.001 | 3 | Akmal  _et al._ (1998)
H4 | 12.931 | 0.01 | 0.001 | 3 | Lackey  _et al._ (2006)
bbb2 | 11.442 | -0.05 | -0.008 | 5 | Baldo  _et al._ (1997)
eng | 12.306 | -0.071 | -0.009 | 2 | Engvik  _et al._ (1996)
mpa1 | 11.696 | -0.062 | -0.003 | 5 | Müther  _et al._ (1987)
ms1 | 14.223 | -0.015 | -0.007 | 5 | Müller and Serot (1996)
qmc700 | 11.942 | -0.008 | -0.002 | 5 | Rikovska-Stone  _et al._ (2007)
sly | 11.873 | -0.053 | -0.003 | 5 | Douchin and Haensel (2001)
wff2 | 10.373 | -0.049 | -0.002 | 5 | Wiringa  _et al._ (1988)
gp1 | 12.302 | -0.02 | -0.007 | 4 | Landry and Essick (2019)
gp2 | 12.345 | -0.024 | 0.001 | 4 | Landry and Essick (2019)
gp$\chi$eft1 | 10.496 | -0.052 | 0.001 | 5 | Essick (2021)
gp$\chi$eft3 | 10.509 | -0.049 | 0.001 | 5 | Essick (2021)
gp$\chi$eft5 | 10.789 | -0.057 | -0.002 | 5 | Essick (2021)
Figure 2: Radius error in fitting Gaussian Process-generate EoSs conditioned
on $\chi$-EFT Essick (2021); Essick  _et al._ (2021b) with the enthalpy
parametriation. We plot two hadronic-conditioned draws, two quark-conditioned
draws, and two hyperonic draws. This indicates the draws are from processes
conditioned on EoS models of the given type, so that e.g. the hadronic process
is consistent with known hadronic EoSs. Nonetheless the processes use
“agnsotic” kernels which lead to very compatible distributions on EoSs for
each of the three cases Landry and Essick (2019); Essick  _et al._ (2020a). A
problem with stitching stability affected multiple of the fits at
$j_{\max}=3$, so we exclude these.

###  III.6 Use cases

The primary function of the enthalpy parametrization is to represent EoS
models for use in numerical simulations containing dense matter. Given the
wide range of models of nuclear matter, the enthalpy parametrization is
intentionally very flexible. Existing parametrizations of the nuclear EoS
typically have a handful of parameters, and extending them might be
nontrivial. In contrast, well-interpolated tables have many “parameters”, or
tabulation points, some of which we would prefer not to resolve in simulations
(such as artificially rapid changes in some pressure derivative). The enthalpy
balances these requirements in such a way that the maximal level of
flexibility can be found without introducing extraneous parameters. This
allows us to resolve EoSs from nuclear theory, Sec. IV.2, as well as EoSs
which extend or modify nuclear models, Sec. IV.3.2. Such flexibility is
crucial for determining the observational implications of new degrees of
freedom at arbitrary density scales.

Furthermore, the requirements laid out in Sec. II are tailored for a specific
application of EoS parametrizations, namely numerical simulations involving
NSs. These requirements are domain specific and need not necessarily lead to
efficient parametrizations for different applications, for example EoS
inference using astrophysical data. Besides the general faithfulness and
computational efficiency considerations, EoS parametrizations employed in
inference need to satisfy an additional requirement: they must provide a
reliable path from the observed data to the EoS constraints. Specifically, the
data must be the primary driver of inference while the impact of the EoS
parametrization itself must be either minimal or driven by first principles
and nuclear theory. Parametrizations that impose a functional form for the EoS
in terms of a finite number of parameters may fail this requirement Greif  _et
al._ (2019); Carney  _et al._ (2018). Specifically, the spectral, piecewise-
polytropic, and speed-of-sound parametrizations impose additional
phenomenological correlations between different densities that are not guided
by nuclear theory but instead by the arbitrary functional form of the
parametrization itself Legred  _et al._ (2022). Though we have not repeated
the analysis of Legred  _et al._ (2022), we expect that the enthalpy
parametrization has the same pitfall as it possesses many nearly-irrelevant
degrees of freedom that are not constrained by current observations and will
generically impart correlations between density scales. We therefore caution
against using it for inference purposes.

##  IV Parametrization Verification and simulations

In this section, we look in depth at fitting nuclear and phenomenological
models with the enthalpy parametrization and perfom numerical simulations.
First, we use SLy1.35 Foucart  _et al._ (2019), a spectral fit to the SLy EoS
Douchin and Haensel (2001); Read  _et al._ (2009a) with a low-density
polytropic exponent of 1.35962. This represents a nuclear EoS which has been
effectively simplified by being fit with a spectral EoS. Therefore, this test
allows us to analyze the performance of the enthalpy parametrization on a
problem where lower dimensional parametrizations are applicable, in terms of
both accuracy and computational performance. We next consider a tabulated DBHF
Gross-Boelting  _et al._ (1999) EoS, derived from relativistic, _ab initio_
calculations of protons and neutrons _dressed_ via interactions with one-boson
exchange potentials.444The EoS we use has employed the Bonn A potential
defined in Ref. Gross-Boelting  _et al._ (1999). It is relatively stiff, with
a typical NS radius of $\sim 13.5\mathrm{km}$. This allows us to assess the
accuracy with which we can fit realistic nuclear models. We then modify the
DBHF EoS using a constant-speed-of-sound parametrization Alford  _et al._
(2015) to construct a model with a strong phase transition, DBHF_2507. With
this we assess the ability of the enthalpy parametrization to augment
realistic low-density models with phenomenological extensions inspired by
nuclear theory.

Using the three models presented above, we study the evolution of isolated NSs
by numerical simulation. As in Ref. Deppe  _et al._ (2022c), we work with
SpECTRE within the Cowling approximation and examine NS modes that are sourced
by density perturbations due to numerical noise. We neglect spacetime dynamics
and magnetic fields, which will likely be most relevant in crust physics where
magnetic and matter pressure are comparable. We run each simulation for 40,000
CFL-limited time steps Courant  _et al._ (1967). For the DG-FD hybrid solver
of SpECTRE we use a sixth-order ($\text{P}_{5}$) discontinuous Galerkin scheme
where each element uses $6^{3}$ Gauss-Lobatto points on the mesh. If an
element switches its mesh from discontinuous Galerkin to finite difference, we
use $11^{3}$ uniformly spaced grid points for finite difference cells. The
finite difference solver needs to compute the solution (in our case $\rho$,
$p$, and $Wv^{i}$, where $W$ is the Lorentz factor and $v^{i}$ the spatial
velocity) at cell interfaces (halfway between grid points). We compute these
using two different reconstruction schemes: the widely employed monotonized
central Van Leer (1977) and a positivity-preserving adaptive order scheme
which was recently implemented in SpECTRE Deppe  _et al._ (tion). In the
$n$-th order adaptive scheme, we first try reconstructing the finite-
difference interface values with a degree $n-1$ polynomial without any
limiting procedure. If the reconstructed values are (i) not positive or (ii)
trigger a certain oscillation-detecting criterion, we repeat the
reconstruction with progressively lower-order methods. In this work we use the
fifth-order adaptive scheme which first tries reconstruction with a quartic
polynomial and switches to monotonized central if the reconstructed values
fail to satisfy the conditions described above. Finally, if the monotonized
central reconstruction did not produce positive values at the interface,
first-order reconstruction is used.

Table 2:  Analysis settings for the SpECTRE simulations with fits to the SLy1.35, DBHF, DBHF_2507 EoSs, and a polytropic EoS. _Labels_ are chosen to serve as unique identifiers for the runs. The integer suffix represents the approximate grid spacing of the run in meters. _Elts_ describes the number of computational elements used in the (three-dimensional) domain. _FD $\Delta x$_ represents the finite-difference grid spacing, in meters, of the finite-difference cells in each element when using finite-difference instead of discontinuous-Galerkin methods; this is the primary measure of resolution of the run and allows for easy comparison to other codes. $\rho_{c}$ represents the initial central (rest-mass) density of the NS being simulated. _Cost_ , in cpu-minutes per CFL-limited timestep, is the approximate cost of computing one time step in this simulation. While runtime depends on an array of factors and may not always be indicative of EoS evaluation speed, differences of $\gtrsim 20\%$ represent on otherwise identical runs likely represent EoS-induced slowdown. _Figs_ represents which figures contain plots pertaining to this run. _Info_ represents the section in which more information about the EoS can be found. _Radius_ represents the TOV radius of the NS being simulated. _Elts/D_ represents the approximate number of computational elements across the diameter of the star.  Label | Elts | FD $\Delta x$ (m)555We express resolution in finite difference grid spacing for easy comparison to finite difference codes. | $\rho_{c}(1/M_{\odot}^{2})$ | Cost(cpum/st) | Figs | Info | Radius (km) | Elts/D
---|---|---|---|---|---|---|---|---
spectral-sly-mc-220 | $\left(24\right)^{3}$ | 224 | 0.00138 | 3.5 | 5 | IV.1.2 | 11.5 | 9
enthalpy-sly-mc-220 | $\left(24\right)^{3}$ | 224 | 0.00138 | 4.1 | 5 | IV.1.2 | 11.5 | 9
spectral-dbhf-mc-130 | $\left(24\right)^{3}$ | 134 | 0.001 | 4.1 | 8 | IV.2.2 | 13.4 | 18
enthalpy-dbhf-mc-130 | $\left(24\right)^{3}$ | 134 | 0.001 | 3.7 | 8 | IV.2.2 | 13.5 | 18
enthalpy-pt-mc-130 | $\left(24\right)^{3}$ | 134 | 0.0021 | 5.1 | N/A | IV.3.2 | 11.8 | 16
enthalpy-pt-ppao-70 | $\left(48\right)^{3}$ | 67 | 0.0021 | 24.6 | 10 | IV.3.2 | 11.8 | 32
enthalpy-pt-mc-70 | $\left(48\right)^{3}$ | 67 | 0.0021 | 20.7 | 10 | IV.3.2 | 11.8 | 32
enthalpy-polytrope-mc-130 | $\left(24\right)^{3}$ | 134 | 0.00128 | 2.07 | 16 | B | 14.1 | 19
polytropic-polytrope-mc-130 | $\left(24\right)^{3}$ | 134 | 0.00128 | 2.1 | 16 | B | 14.1 | 19
enthalpy-smoothpt-170 | $\left(24\right)^{3}$ | 168 | 0.0021 | 3.5 | NA | IV.3.3 | 11.9 | 12

###  IV.1 SLy1.35

####  IV.1.1 SLy1.35: fit results

Figure 3: Fitting SLy1.35 (a spectral model of SLy) with the enthalpy
parametrization, expressed through the difference in pressure divided by the
density. The SLy1.35 EoS value for $p(\rho)/\rho$ is marked by a maroon dashed
line for comparison to the residuals. The $j_{\max}=5$ and $j_{\max}=2$ fit
residuals are marked in light blue and indigo. The vertical blue dashed line
marks the stitching density between the enthalpy and the spectral
parametrizations, while the vertical red dot-dashed line marks the central
density of the NS we simulate in Sec. IV.1.1. The solid red horizontal line
marks an error level of $3\times 10^{-3}$ for comparison with fits in Sec.
IV.2.1; errors below $3\times 10^{-3}$ at $(1,3)\rho_{\mathrm{nuc}}$ serve as
a heuristic for a good fit.

We fit SLy1.35 with the enthalpy parametrization and show the error in
pressure divided by density as a function of density in Fig. 3. We vary the
number of trigonometric terms in Eq. (10) and show results with $j_{\max}=2$
and $j_{\max}=5$. The $j_{\max}=5$ fit shows exceptional agreement; the error
measure, $\Delta p/\rho$, is near or below $1\times 10^{-4}$ over essentially
the entire domain. The $j_{\max}=5$ fit shows increased error, though
$\Delta(p/\rho)$ remains near or below $3\times 10^{-3}$ above
$\rho_{\mathrm{nuc}}$. We stitch to a spectral parametrization below
$\rho_{\mathrm{nuc}}$, marked in the Fig. 3 as a vertical dashed blue line.
Even though the low-density behavior of the EoS is a spectral EoS fitting a
spectral EoS, it is not guaranteed the low-density fit is good, because we
prioritize smooth stitching to the enthalpy solution above accurate low-
density EoS modeling, see App. A. In line with this, we see a significantly
better low-density fit for $j_{\max}=5$.

Figure 4: Radius error $\Delta R$ as a function of mass for the SLy1.35 EoS
and the $j_{\max}=10$ and $j_{\max}=5$ enthalpy parametrization fits. We mark
the mass of the stars with central density $\rho_{c}\sim
3.04\rho_{\mathrm{nuc}}$, (simulated in Sec. IV.1.2), with dashed-dot lines.
Consistently with the microscopic comparison of Fig. 3, the enthalpy fit can
reproduce macroscopic quantities with excellent agreement. The error decreases
with more trigonometric terms, but always remains small compared to $200$ m
grid resolution.

####  IV.1.2 SLy1.35: Relativistic simulations

We carry out simulations directly with SLy1.35 using the defining spectral
expansion Foucart  _et al._ (2019) as well as the $j_{\max}=5$ enthalpy
parametrization fit; details are given in Table 2. We evolve a NS with an
initial central density of $\sim 3.04\rho_{\mathrm{nuc}}$ which has a Tolman-
Oppenheimer-Volkoff (TOV) Oppenheimer and Volkoff (1939) mass of about
$1.4M_{\odot}$ and a radius of about $11.5\leavevmode\nobreak\ \mathrm{km}$,
see Figs. 3 and 4. The simulation resolution corresponds approximately to a
$220$ m finite difference grid spacing. We plot the central density as a
function of time and its spectrum

$\hat{\rho}_{\rm c}(\omega)=\int_{0}^{T}\rho_{\rm c}(t)e^{-i\omega t}dt\,,$
(24)

in Fig. 5 and find essentially identical evolution between the spectral and
enthalpy fits, in line with expectations from the static tests of Sec. IV.1.1.
This demonstrates that the enthalpy parametrization is able to faithfully
reproduce results from lower-dimensional parametrizations.

Figure 5: NS central density as a function of time (top panel) and its spectrum (bottom panel) for SpECTRE simulations with SLy1.35 (red dashed) and its $j_{\max}=10$ enthalpy fit (blue solid). These runs are labeled spectral-sly-mc-220 and enthalpy-sly-mc-220 in Table 2. In both plots the curves are nearly indistinguishable. We plot times in both milliseconds (ms), and dynamical times ($t_{\mathrm{dyn}}\equiv 1/\sqrt{\rho_{\rm c}}$) Parametrization Evaluation Cost (ns) | $p(\rho)$ | $\epsilon(\rho)$
---|---|---
enthalpy, $j_{\max}=5$ | 224 | 225
enthalpy, $j_{\max}=2$ | 120 | 120
spectral | 62 | 315
Table 3: Evaluation cost in nanoseconds for the spectral and two different
enthalpy fits to SLy1.35 for the pressure and internal energy evaluated at
$\rho=5\times 10^{-4}M_{\odot}^{2}$. The spectral parametrization has a
shorter (longer) pressure (internal energy) evaluation time. The enthalpy
evaluation cost further increases with the number of trigonometric terms
employed.

With regards to computational cost, the enthalpy parametrization results in an
overall 15% increase in total runtime compared to the spectral parametrization
on similar hardware. To isolate the EoS evaluation cost, we benchmark the
$p(\rho)$ and $\epsilon(\rho)$ evaluation in Table 3. The two parametrizations
have comparable evaluation times though exact numbers are sensitive to the
number of trigonometric terms in the enthalpy case. While $p(\rho)$ evaluation
is in general faster with the spectral parametrization, the opposite is true
for $\epsilon(\rho)$. This is because the spectral parametrization needs to
perform a quadrature to calculate $\epsilon(\rho)$, see Sec. II. In the
enthalpy parametrization, trigonometric terms cause slowdowns, although even
with $2\times j_{\max}=10$ terms the $p(\rho)$ cost does not exceed a factor
of $4$. Further studies with single-polytropic nuclear EoSs suggest that this
disadvantage effectively disappears if $j_{\max}=0$, see App. B.

###  IV.2 DBHF

####  IV.2.1 DBHF: Fitting a tabulated nuclear model

We fit the tabulated DBHF EoS with the enthalpy parametrization and further
explore the effect of low-density stitching to the spectral parametrization by
probing two different stitch densities: $\rho_{\mathrm{nuc}}$ and
$2.5\rho_{\mathrm{nuc}}$; these fits are referred to as “low-stitch” and
“high-stitch” in what follows. In the high-stitch case we use $j_{\max}=10$
trigonometric terms, while we find that $j_{\max}=5$ is enough for the low-
stitch one. See App. A for more details. We examine the microscopic and
macroscopic performance of both fits in Figs. 6 and 7. On the microscopic
side, the low-stitch fit achieves higher accuracy above
$1.1\rho_{\mathrm{nuc}}$, but worse accuracy below.

The macroscopic side presents a clearer picture. When we use the enthalpy
parametrization to describe the EoS down to a density of
$\rho_{\mathrm{nuc}}$, we obtain excellent agreement with the tabulated EoS,
with radius differences ${\cal{O}}(1)$ m for astrophysically relevant NS
masses. However, when we stitch to the spectral parametrization at
$2.5\rho_{\mathrm{nuc}}$ the radius error increases to $\mathcal{O}(100)$ m at
$1.4M_{\odot}$. The improved agreement between the $2.5\rho_{\mathrm{nuc}}$
and the $\rho_{\mathrm{nuc}}$ stitching fits can be attributed to the high
accuracy of the enthalpy parametrization in the range $\rho_{\mathrm{nuc}}$ to
$2.5\rho_{\mathrm{nuc}}$, as seen in Fig. 6. The difference in the two errors
is particularly pronounced near $2\rho_{\mathrm{nuc}}$, consistent with the
observed strong correlation between the pressure at twice saturation density
and radius of a $1.4M_{\odot}$ star Lattimer and Prakash (2001). This further
establishes the importance of the enthalpy parametrization as a flexible EoS
parametrization at nuclear saturation and above, in this case it appears
errors in $p/\rho$ must be at most $3\times 10^{-3}$ to achieve high-precision
reproduction of astrophysical observables. However, this is not an indication
that the spectral parametrization _cannot_ fit DBHF well, it just cannot fit
DBHF well while maintaining $\mathcal{C}^{1}$ pressure smoothness at the high-
density transition to the enthalpy parametrization and the low density
transition to the crust; see App. A.

Figure 6: Same as Fig. 3 but for the DBHF nuclear EoS model and two enthalpy
fits that are stitched to a spectral parametrization at $\rho_{\mathrm{nuc}}$
(low-stitch, light blue) and $2.5\rho_{\mathrm{nuc}}$ (high-stitch, indigo).
We also plot the tabulated DBHF model $p(\rho)/\rho$ in dashed-teal for
reference. The vertical dashed lines denote the stitching densities. We also
mark the value $\Delta p/\rho=3\times 10^{-3}$ as a solid red horizontal line
for reference. We mark the central density of the star we simulate in Sec.
IV.2.2 with a vertical red, dot-dash line. Figure 7: The NS mass-radius
relation but for the DBHF nuclear model and two enthalpy fits that are
stitched to a spectral parametrization at $\rho_{\mathrm{nuc}}$ (low-stitch,
light blue) and $2.5\rho_{\mathrm{nuc}}$ (high-stitch, indigo). We find
visibly improved fits to the $M$–$R$ relation when the enthalpy
parametrization extends down to lower densities. Red dots mark the NSs we
evolve in Sec. IV.2.2.

####  IV.2.2 DBHF: Relativistic simulations

We next turn to SpECTRE simulations using the DBHF fits from the previous
section. Since the low- and high-stitch fits predict $\mathcal{O}(100)$ m
differences for $R_{1.4}$, we target a resolution at that level in order to
resolve their effect. We select a NS central density of
$2.21\rho_{\mathrm{nuc}}$ that lies between the two stitching densities of
$\rho_{\mathrm{nuc}}$ and $2.5\rho_{\mathrm{nuc}}$, see Figs. 6 and 7. Run
details and settings are given in Table 2. As a consequence, the NS resulting
from the high-stitch EoS is fully described by the spectral part of the EoS.
In the low-stitch EoS, the NS is described with the enthalpy parametrization
out to $r/R\approx 7/8$, about two-thirds of the coordinate volume of the
star.

While the high-stitch EoS does represent a spectral fit to the DBHF EoS, the
spectral parameters are selected by the requirement that the spectral
parametrization reproduces the correct low-density behavior of DBHF and is
smoothly stitched to the enthalpy parametrization. It is important to note
that a better spectral fit to any particular astrophysical quantity, such as
$R_{1.4}$ may be possible, but the fit accuracy is typically lower than in the
enthalpy parametrization case. Even if the few degrees of freedom in the
spectral model are fit to minimize errors in astrophysical observables, fixing
the low-density behavior of the EoS often results in $R(M)$ deviations of 50 m
or more Foucart  _et al._ (2019). In what follows, we leverage the mismatch of
Fig. 7 to examine how well we can resolve EoSs with $\sim 100$ m radius
differences in simulations with similar resolution. Since the high-stitch NS
is fully described by the spectral parametrization, this test also serves as a
comparison of runtimes between the spectral and enthalpy parametrizations. We
do not utilize the tabulated version of DBHF because SpECTRE currently cannot
perform GRMHD simulations using tables.

We carry out simulations as detailed in Table 2 and plot the spectrum of the
central density of each star in Fig. 8. The two spectra disagree both in the
location of the NS modes and their strength, a consequence of the EoS
mismodeling shown in Figs. 6 and 7. In particular, the fundamental radial
modes disagree by nearly $3.5\%$, a difference of $\sim 130\,\mathrm{Hz}$ in
this case. Table 2 further shows the simulation runtime which is comparable in
the $130$ m resolution case; each run took about a day on $\sim 70$ processing
elements. Based on the benchmarking results of Sec. IV.1.1, total EoS
evaluation time should be comparable for the two runs, as the spectral
parametrization evalautes the pressure about twice as quickly as the enthalpy
parametrization, but evaluates the internal energy about 3 times slower. Since
the number of pressure and internal energy evaluations throughout the entire
simulation is not known _a priori_ , we cannot preemptively conclude which
should run faster, though it is likely the difference would be small. This is
reflected in the runtime; differences of $10\%$ are found, with the enthalpy
parametrization running slightly faster. Nonetheless, this could be due to an
array of confounding factors such as task allocation efficiency, and hardware
differences. We therefore conclude that the enthalpy parametrization, despite
having more flexibility, is not slower than lower dimensional parametrizations
at these resolutions for practical problems.

Figure 8: NS central density spectrum for SpECTRE simulations with enthalpy
fits to the DBHF nuclear EoS that are stitched to the spectral parametrization
at $\rho_{\mathrm{nuc}}$ (low-stitch, blue solid) and $2.5\rho_{\mathrm{nuc}}$
(high-stitch, red dashed). These runs are labeled enthalpy-dbhf-mc-130 and
spectral-dbhf-mc-130 respectively in Table 2. The simulated star has a central
density of $\approx 2.2\rho_{\mathrm{nuc}}$; it is marked in Figs. 6 and 7. In
the red case, the NS is completely described by the spectral EoS as its
central density is below $2.5\rho_{\mathrm{nuc}}$.

###  IV.3 DBHF_2507: Phase transitions

We now turn our attention to EoS with strong phase transitions and study both
smooth and non-smooth (i.e., piecewise) EoSs. We base our studies on DBHF_2507
which is constructed by combining DBHF with the constant-speed-of-sound
phenomenological parametrization for strong phase transitions Alford  _et al._
(2015). We select a transition density of $\rho_{t}=2.5\rho_{\mathrm{nuc}}$
and latent heat ratio $\Delta e/e=0.7$ Han and Steiner (2019); Chatziioannou
and Han (2020). The pressure remains constant during the phase transition,
while above that it has a constant speed of sound with $c_{s}^{2}=1$. The
induced phase transition causes a second stable branch to appear in the
$M$–$R$ relation above masses $\sim 1.6M_{\odot}$.

####  IV.3.1 DBHF_2507: piecewise parametrization

In its original form described above, DBHF_2507 is piecewise smooth, and it
can be represented effectively by a piecewise version of the enthalpy
parametrization. Below the phase transition we use either the low- or high-
stitch fits from Sec. IV.2, and transition to a new enthalpy segment after the
transition. In the high-stitch case, the hadronic part of the EoS is
completely described with the spectral parametrization. During the transition,
DBHF_2507 possesses a formally constant pressure as a function of density;
however, constructing TOV solutions using the method of Lindblom Lindblom
(1992) — the TOV method implemented in SpECTRE — requires $dh/dz=dp/d\rho$ to
be strictly positive. Therefore, in the transition region we modify the EoS to
exhibit $dh/dz=\delta$, where $\delta$ is some quantity large enough to
guarantee that $h(z)$ is numerically invertible, but still small enough to
have a small impact on the TOV solution relative to the target
resolution.666At a central density of $5.07\rho_{\mathrm{nuc}}$, the
difference in radius induced by using $\delta=1\times 10^{-4}$ instead of
$\delta=1\times 10^{-3}$ is less than $5$ m. After the end of the phase
transition, $h(z)$ is given by a constant speed of sound form, see Eq.
(III.1); this is similar to the procedure demonstrated in Gieg  _et al._
(2019).

The advantage of the enthalpy parametrization in this problem is that it is
able to model constant-speed-of-sound matter (see App. B for the polytropic
case) efficiently and with no fine-tuning. Compare this to polytropic (or
spectral) models, which can only model constant-speed-of-sound matter well
when

$\Gamma\equiv\frac{\rho}{p}\frac{dp}{d\rho}=\frac{\rho}{p}hc_{s}^{2}=\frac{p+e}{p}c_{s}^{2}$
(25)

is slowly varying. This is typically not true until some density greater than
the phase transition, where $p=p_{0}+c_{s}^{2}\Delta e\approx c_{s}^{2}\Delta
e$, especially if $c_{s}^{2}$ is small compared to $1$, such as models where
$c_{s}^{2}=1/3$ in the core Kurkela  _et al._ (2010); Tews  _et al._ (2018a).
In contrast, the enthalpy parametrization can model constant-speed-of-sound
matter to arbitrary precision, and benchmarking results demonstrate that in
such cases it can even outperform polytropic EoSs by up to 25%.

We plot the $M$–$R$ curve in Fig. 9, using the low-stitch and high-stitch fits
for the hadronic part of the EoS as discussed in Sec. IV.2.1. The low-stitch
EoS shows better agreement with DBHF_2507, consistent with previous results;
see Fig. 7. The transition mass and radius for the low-stitched model are
functionally identical to the tabulated values with errors of $\lesssim
0.01\,M_{\odot}$ and $\lesssim 1$ m. In the high-stitch case the errors
increase to $0.05\,M_{\odot}$ and $\sim 100$ m. Nonetheless, errors decrease
with increasing central density and the maximum mass $M_{\rm max}$ is
consistent to $\sim 0.01\,M_{\odot}$ for both fits. This indicates that the
enthalpy parametrization can produce effective EoS fits at high densities even
when extending a (relatively) poor low-density fit.777Such comparisons to
tabulated models might be difficult to interpret, as a $1\%$ interpolation
inconsistency in $\rho(p)$ can lead to differences of $\mathcal{O}(100\
\mathrm{m})$ on the second stable branch. This problem is more pronounced here
as the DBHF_2507 construction requires computing
$\rho(e)=2.5\rho_{\mathrm{nuc}}$ via table-based root finding, a procedure
that depends on the interpolation strategy and, in turn, affects the
transition mass and radius. For hadronic EoSs this issue is suppressed, as
differences in interpolation are smoothed over by the integration of the TOV
equations. The enthalpy parametrization, having an analytic expression for
$e(\rho)$, does not face this issue.

Figure 9: Same as Figs. 6 & 7 but for the DBHF_2507 EoS. The procedure by
which the low- and high-stitch EoSs of Sec. IV.2.1 are extended through the
phase transition is described in Sec. IV.3. Consistent with Fig. 7, the low-
stitch case can more accurately reproduce the parameters of the phase
transition, though caution must be exercised when comparing to tabulated
models as differences in interpolation in this case can be substantial. See
the text of Sec. IV.3.2. In the top panel, the black vertical dashed line
marks the onset of the phase transition. The red dot-dashed line in the top
panel and the red dots in the bottom panel mark the NSs we use in subsequent
simulations, analogous to Figs. 6 & 7 respectively.

####  IV.3.2 DBHF_2507: Relativistic simulations

We perform SpECTRE simulations with both the high- and low-stitch EoSs and NSs
with central density of $\rho_{c}=4.67\rho_{\mathrm{nuc}}$, above the
transition density from nuclear to quark matter; see the red dots in Fig. 9.
Preliminary results with low spatial resolutions demonstrated that for such
$\rho_{c}>\rho_{\mathrm{t}}$, the NS undergoes strong density oscillations
that are quickly damped. Given that the fundamental mode is long-lived
Kokkotas and Ruoff (2001), this short damping timescale is probably related to
numerical dissipation. Therefore we perform and compare simulations at various
grid resolutions to ensure convergence, increasing the number of computational
elements while the number of grid points inside each element is fixed. The
main results presented below correspond to a $\sim 70$ m resolution. We also
restrict to the low-stitch fit since the differences between the low- and
high-stitch fits are likely resolvable for $<130$ m resolution, see Fig. 9.

Figure 10: Normalized NS central density as a function of time (top panel) and
its spectrum (bottom panel) for SpECTRE simulations with enthalpy fits to the
DBHF_2507 nuclear EoS that are stitched to the spectral parametrization at
$\rho_{\mathrm{nuc}}$ for two different choices of finite-difference
reconstruction schemes. The adapative order reconstructor is marked in blue
and the monotonized central reconstructor is marked in red. Run details are
listed as enthalpy-pt-ppao-70 and enthalpy-pt-mc-70 respectively in Table 2.
The simulated star has a central density of $\sim 4.67\rho_{\mathrm{nuc}}$; it
is marked in Fig. 9. We find excellent agreement on mode frequencies but
slight differences in power distribution.

We plot the central density and the spectrum in the top and bottom panels of
Fig. 10 for both reconstruction schemes. We find good agreement in the
frequency and damping time of the density modes, though the monotonized
central scheme predicts more than double the power of the adaptive order
method below $\sim 80$ kHz. Interestingly, we find that the presence of a
quark core in the NS changes the spectrum qualitatively, c.f., Fig. 8. The
spectrum is now dominated by modes in the ${\cal{O}}(10)$ kHz range, an order
of magnitude higher than the hadronic NS case of Fig. 8. The spacing of the
modes, about $16$ kHz, is of order $c/2R_{\mathrm{core}}$, where
$R_{\mathrm{core}}\sim 6$ km is the radius of the quark core. We attribute
this to density perturbations that are confined to the quark core and are only
weakly coupled to the bulk behavior of the star across the transition. In
order to confirm this, we plot the density profile of the star extracted from
the run enthalpy-pt-ppao-70 as a function of radius in Fig. 11 at different
times. Most of the oscillation power sourced in the quark core is reflected
back into the core at the quark-hadronic boundary, with only a small fraction
getting transmitted into the hadronic region. The reflected pulse gets
inverted (fixed-end reflection); this is consistent with theoretical
expectations since the sound speed changes from $c_{s}^{2}=1$ in the quark
core to $c_{s}^{2}\sim 0.3$ in the hadronic region at the boundary.

Figure 11: Rest-mass density profile relative to the initial profile
$\Delta\rho\equiv\rho(r,t)-\rho(r,0)$ as a function of radius for a hybrid
star described with DBHF_2507 (left panel; details in Sec. IV.3.2) and a
simple polytrope with $\Gamma=2$ (right panel; details in App. B) for
different times (top to bottom). We denote the NS surface with a vertical
solid gray line in each panel and the quark-hadronic boundary with a vertical
red line in the left panel. We show snapshots of the density at four different
times in order to examine the dynamical behavior of the density oscillations.
For the hybrid star (left) density perturbations are partially transmitted and
reflected at the quark-hadronic boundary, while for the polytrope (right) the
wave smoothly propagates back and forth within the NS interior. Small black
arrows highlight the wave packet and its traveling direction. The hybrid star
snapshots are from the run enthalpy-pt-ppao-70 and the polytrope snapshots are
from a simulation with identical domain, finite-difference reconstruction
scheme, and central density, but a polytropic EoS (33) in place of DBHF_2507.
See polytropic-polytrope-mc-130 for details of a lower-resolution polytropic
simulation.

The initial perturbation needed to drive these modes is provided by numerical
noise near the transition, with $\mathcal{O}(50-100\mathrm{kHz})$ being the
scale of the sound crossing frequency of a computational element of our
domain. High frequency modes, in particular $p$-modes being primarily confined
to the core of the star is in line with expectations for radial modes Sen  _et
al._ (2022).

Simulation runtimes are provided in Table 2. We find that the adaptive-order
simulation (enthalpy-pt-ppao-70) has a $\sim 20\%$ longer runtime than the
monotonized central (enthalpy-pt-mc-70) simulation, which is expected. Turning
to the EoS and comparing enthalpy-pt-mc-130 and enthalpy-dbhf-mc-130 at
identical resolutions, number of CPUs, and reconstruction schemes, we find a
$\mathcal{O}(40\%)$ increase in runtime for DBHF_2507 as compared to DBHF. We
do not attribute this runtime slowdown to increased EoS evaluation time, as
the enthalpy EoS employed beyond the phase transition uses no trigonometric
correction terms, and thus is nearly computationally identical to the
polytrope profiled in Table 4. This indicates that individual EoS evaluations
(above the transition) are actually somewhat cheaper than in either of the
fits to DBHF, discussed in Sec. IV.2 (below the transition they are
identical). Instead, we attribute the slowdown to the non-smoothness of this
EoS; $p(\rho)$ is not analytically $\mathcal{C}^{1}$ across the phase
transition, and $\rho(p)$ is an incredibly sensitive function near the
transition. Since our default primitive recovery scheme requires root-finding
to determine $\rho(p)$ Kastaun  _et al._ (2021), this can result in
significant slowdowns.

Comparing the simulations enthalpy-pt-mc-130 and enthalpy-pt-mc-70, we also
find that refining the grid resolution from 134 m to 67 m results in only a
4-fold increase in runtime, compared to the expected $(134/67)^{3}=8$.888Since
we run simulations with the fixed number of time steps, refining the grid does
not lead to any major slowdown from time stepping. This would not be the case
if we ran to a fixed final time. We attribute this lower-than-expected
increase to the DG-FD hybrid scheme Deppe  _et al._ (2022b). At higher
resolutions each cell is smaller, therefore finite-difference cells are more
tightly concentrated in regions with discontinuities. This results in a lower
fraction of the NS reverting to the slower finite-difference scheme from the
faster discontinuous Galerkin scheme. We display a slice of the NS for
enthalpy-pt-ppao-70 in Fig. 12, and mark the cells which have reverted to
finite-difference in red. The majority of the NS interior is indeed using the
discontinuous Galerkin method, with finite difference being used only at the
phase-transition and surface interfaces.

Figure 12: NS rest mass density (colorbar, units of $M_{\odot}^{-2}$) on the
$y-z$ plane at $t=100M_{\odot}$ for the run enthalpy-pt-ppao-70. Red marks
subdomain elements where finite-difference is used. The finite-difference
cells are confined near the NS surface (outer circle, white solid line) and
the phase transition layer (inner circle) where discontinuities are expected.
The majority of the star is still evolved with the more computationally
efficient discontinuous-Galerkin method.

####  IV.3.3 Smooth transitions

Due to its flexibility, the enthalpy parametrization can also be used to model
smoother transitions in the EoS, such as those that may arise from a crossover
transition. To demonstrate this, we begin with DBHF_2507 and average the speed
of sound at nearby points over the entire EoS. This is distinct from the small
perturbation added to the EoS in in Sec. IV.3.2, as in this case we expect the
resulting EoS to be well described by a smooth interpolant. To demonstrate
this, we fit this new smoothed EoS with $j_{\max}=4$ trigonometric terms and
display the fit in Fig. 13. We find sub-$1\%$ agreement in the
$1-2\rho_{\mathrm{nuc}}$ density region which most directly affects
macroscopic observables. Relative errors are typically higher below the
enthalpy-spectral transition point (set to $\rho_{\mathrm{nuc}}$ here). This
is because smoothing the EoS is done by locally averaging the speed of sound,
so that after averaging the speed of sound is locally close to being constant.
For polytropic and nearly polytropic EoSs the speed of sound is not nearly
constant at low densities, so a spectral parametrization cannot effectively
fit the smooth $c_{s}^{2}$ EoS.

We expect simulating NSs with EoSs displaying smooth but rapidly varying
speeds of sound to be slower. This is for two primary reasons. First, fitting
EoSs to some fixed degree of precision for more complicated EoSs typically
requires adding more parameters, increasing evaluation time. Second, EoSs with
rapid changes in $dp/d\rho$ tend to slow down primitive recovery, as the
function $\rho(p)$ must be evaluated by root-finding. Even though the EoS in
this case is analytically smooth, root-finding algorithms require more
evaluations if the function is quickly varying.

Figure 13: Same as Fig. 6 but for the smoothed version of DBHF_2507
constructed in Sec. IV.3.3. The enthalpy parametrization achieves sub-$1\%$
errors in the most relevant region, $1-2\rho_{\mathrm{nuc}}$.

We perform a run with identical central density, $\rho\sim
4.67\rho_{\mathrm{nuc}}$ for 10,000 CFL limited time steps, in order to bound
performance decreases. We display the results in Table 2 as enthalpy-
smoothpt-170 . We find that the EoS presented in Fig. 13 requires a comparable
time per evolution step to enthalpy-dbhf-mc-130 indicating the EoS is
sufficiently smooth to not induce a large slowdown at this resolution. We
further plot the oscillations of the central density of both the smooth-
transition DBHF_2507 model (smooth) and enthalpy-pt-mc-130 (sharp) in Fig. 14.
We find that the smoothed fit does not lead to the characteristic decoupling
of core modes, meaning that such a model would be a poor representation of the
true DBHF_2507 EoS, even if it is able to reproduce other characteristics of
DBHF_2507, such as a small radius near $M_{\max}$.

Figure 14: NS central density as a function of time for the run enthalpy-
smoothpt-170 (blue, solid) and enthalpy-pt-mc-130 (red, dashed). The smooth
fit poorly reflects the mode structure of the true DBHF_2507 EoS, even though
the behavior of the microscopic EoS is qualitatively similar.

In general, we expect smooth EoSs to be most effectively represented by
globally smooth parametrizations, while EoSs with discontinuities will be
better modeled by piecewise parametrizations. In addition, nonsmoothness can
lead to a loss of accuracy in simulations Foucart  _et al._ (2019), so an
additional trade-off may exist between accuracy and performance in the choice
to use a smooth versus a piecewise representation. The enthalpy
parameterization is flexible enough to be effective in both the piecewise EoS
and the smooth EoS cases.

##  V Discussion

We introduced a new enthalpy-based parametrization for the cold nuclear EoS
that can capture a wide range nuclear models and their phenomenological
extensions using polynomials and trigonometric terms. The _enthalpy
parametrization_ emphasizes _flexibility_ , as it is able to effectively model
both smooth and non-smooth nuclear models, and _computational performance_ as
its evaluation cost scales with the number of parameters used. For example, it
displays comparable performance to single-polytrope parametrizations for the
case of polytropic EoSs, while the computational cost scales with the number
of fit parameters for more complex (such as non-smooth) models. This trade-off
between computational performance and flexibility, allows us to tune EoS fits
to the resolution requirements of the problem at hand.

Computational performance is achieved by inexpensive evaluation of the various
thermodynamic quantities. The $p(\rho)$ evaluation cost does not exceed
$\mathcal{O}(4)$ times that of a polytropic EoS for any case we investigated,
even when many trigonometric terms are used. In cases where this slowdown is
significant, the enthalpy parametrization may be sped up significantly by
using Clenshaw’s method Press  _et al._ (2007). We obtain faster evaluation of
$\epsilon(\rho)$ than the existing spectral parametrization in all cases, as
the latter evaluates $\epsilon$ numerically, while the enthalpy
parametrization computes all thermodynamic quantities analytically. Overall,
the additional computational cost of the enthalpy parametrization on top of
other existing parametrizations is always smaller than the cost of other
simulation components.

With the caveat that quantifying EoS fitting accuracy is subtle and depends on
the parameters one compares, we overall find that the enthalpy parametrization
is able to successfully fit nuclear models. In principle and in the context of
numerical simulations, EoS parametrizations need only fit the nuclear EoS as
well as the simulation resolution. Nonetheless, even subpercent errors in the
pressure near $\rho_{\mathrm{nuc}}-2\rho_{\mathrm{nuc}}$ can lead to $\sim
100$ m differences in NS radii. In contrast to lower-dimensional or less
flexible parametrizations, we show that the enthalpy parametrization is able
to fit tabulated and phenomenological nuclear models to effectively arbitrary
precision by using additional parameters. The optimal number of parameters is
then determined by balancing accuracy and computational cost for a given
numerical resolution.

The enthalpy parametrization’s flexibility allows us to efficiently and with
little fine tuning represent both smooth and non-smooth nuclear models. The
latter may correspond to models with strong phase transitions that we can fit
and numerically evolve using SpECTRE. Our simulations demonstrate that we can
stably evolve such stars in the Cowling approximation. However, studying the
evolution of hybrid hadronic-quark NSs away from an unstable EoS branch that
falls between the hadronic and the quark branches Espino and Paschalidis
(2022) hinges on full metric evolution coupled to GRMHD. SpECTRE’s hybrid DG-
FD scheme is crucial for the computational performance of these simulations.
The DG-FD scheme allows phase transitions to be modeled with lower-order
finite-difference methods while continuing to use higher-order discontinuous-
Galerkin methods throughout the individual hadronic and quark regions. This
leads to better computational scaling than might be expected upon mesh
refinement, as better resolution of boundaries (such as the quark-hadronic
matter boundary) within the star reduces the amount of the domain which uses
the slower finite-difference approach.

The enthalpy parametrization is a step toward ensuring that numerical
simulations can efficiently represent a wide range of nuclear phenomenology.
Accurate simulations of NSs will continue being crucial for the interpretation
of new astrophysical and experimental data. Even with current EoS constraints,
the space of potential BNS phenomenology is large, and many questions remain
regarding the impact of magnetic fields, instabilities, temperature effects
Carbone and Schwenk (2019); Raithel  _et al._ (2019) and transport physics.
Future steps include extending the applicability of the enthalpy
parametrization beyond cold, beta-equilibrated nuclear matter, and
incorporating more physical effects in SpECTRE simulations.

The simulations presented here were performed with SpECTRE commit hash
2df19579a84385b3d5ab4663e3da7e33012e0355. The earliest release of SpECTRE with
this commit is version 2023.01.13 Deppe  _et al._ (2023). Input files for the
runs performed, including enthalpy fit parameters for each EoS studied, are
available on Github SXS (2023).

###### Acknowledgements.

I.L. thanks Tianqi Zhao for helpful conversations in preparing this
manuscript. The authors thank Reed Essick, Ingo Tews, Phil Landry, and Achim
Schwink for access to $\chi$-EFT conditioned Gaussian process draws.
Charm++/Converse Kale  _et al._ (2020) was developed by the Parallel
Programming Laboratory in the Department of Computer Science at the University
of Illinois at Urbana-Champaign. This project made use of python libraries
including scipy and numpy Virtanen  _et al._ (2020); Oliphant (06). Figures
were produced using matplotlib Hunter (2007) and ParaView Ayachit (2015).
Computations were performed with the Wheeler cluster at Caltech, which is
supported by a grant from the Sherman Fairchild Foundation and Caltech. This
work was supported in part by the Sherman Fairchild Foundation at Caltech and
Cornell, as well as by NSF Grants No. PHY-2011961, No. PHY-2011968, and No.
OAC-1931266 at Caltech and by NSF Grants No. PHY-1912081 and No. OAC-1931280
at Cornell. IL and KC acknowledge support from the Department of Energy under
award number DE-SC0023101. FF gratefully acknowledges support from the
Department of Energy, Office of Science, Office of Nuclear Physics, under
contract number DE-AC02-05CH11231, from NASA through grant 80NSSC22K0719, and
from the NSF through grant AST-2107932. The authors are grateful for
computational resources provided by the LIGO Laboratory and supported by
National Science Foundation Grants PHY-0757058 and PHY-0823459.

##  Appendix A Fitting the enthalpy parametrization

Figure 15: Cost, Eq 26, in arbitrary units, of the fit to the phenomenological
EoS of Fig. 1 as a function of $k$. The minimum occurs at $k$ slightly larger
than $k_{F}$, in this case near $k=1.4k_{F}$.

In this Appendix we provide details about the procedure with which we fit some
tabulated EoS data with the enthalpy parametrization which includes the
following parameters:

  * •

The upper $\rho_{\max}$ and lower $\rho_{\min}$ density limits are chosen
based on the densities of interest. For NS simulations, reasonable values are
$\rho_{\min}=\rho_{\mathrm{nuc}}$, $\rho_{\max}=7\rho_{\mathrm{nuc}}$, but the
upper limit depends on the maximum density expected in the simulation.

  * •

The scaling parameters $\rho_{0}$ and the wavenumber of the trigonometric
correction terms $k$ are not fit, but rather fixed. When we extend a model via
a constant speed of sound, Sec. IV.3.2, the choice of $\rho_{0}$ is determined
by the modeling problem. When the EoS is fit, $\rho_{0}$ is chosen
$\rho_{0}\in(0,\rho_{\min}]$, so that
$z_{\max}=\log(\rho_{\max}/\rho_{0})\lesssim i_{\max}$, see Sec. C. We find
$\rho_{0}=\rho_{\min}/2$ is generally a robust choice. Analogous to how
$\rho_{0}$ controls the scale of polynomial terms, $k$ controls the scale of
trigonometric oscillations. As described in Sec. III, typically $k\approx
k_{F}$ is a good choice, but small perturbations $k\in[k_{F}/2,2k_{F}]$ may
improve the fit quality for certain problems, depending on the details of the
EoS. Figure 15 shows that the effect of varying $k$ is small for the
particular test problem displayed in Fig. 1.

  * •

The parameters $i_{\max}$ and $j_{\max}$ determine the number of polynomials
and trigonometric terms respectively; see Eq. (9) and Eq. (10). The quality of
the fit is a strong function of $i_{\max}$ and $j_{\max}$, but increasing
$j_{\max}$ above $\sim 10$ comes at a considerable computational cost even at
low resolutions. On the other hand the cost of increasing $i_{\max}$ is small,
typically of order $2\%$ or less of the total cost of the $p(\rho)$ evaluation
per additional polynomial term.

  * •

The coefficients of the polynomial $\gamma_{i}$, Eq. (9), and the
trigonometric $a_{j},b_{j}$, Eq. (10), expansion are fit through a linear
least-squares approach.

  * •

The energy density of the EoS at the stitching point, $e_{\min}=e(z_{\min})$.
This is the integration constant associated with solving $de/dz=\rho h$. This
parameter is constrained by $\rho_{\min}h_{\min}-e_{\min}=p_{\min}\geq 0$. In
principle $e(z_{\min})$ can be computed from EoS tables, but in practice EoS
tables may be too coarsely tabulated, or may contain violations of the first
law of thermodynamics at levels which significantly affect the computed value
of $\epsilon$. For example, a fractional error of $1\times 10^{-3}$ in
$e_{\min}$ will often translate to a fractional error of $\sim 1\times
10^{-1}$ in $\epsilon$, which therefore shifts the value of $p(z)$ by $10\%$,
as $h\rho=p+\epsilon$ is fixed by the parametrization. Therefore in certain
cases it is more effective to treat $e_{\min}$ as a free parameter, and
further use it to optimize the values of $p(\rho)$. In practice $e_{\min}$ is
set by the low-density EoS parametrization to guarantee thermodynamic
consistency.

. Given a target EoS with enthalpy $h(z_{i})$ at discrete densities $z_{i}$,
the linear fit is based on minimizing the cost function

$C(a_{j},b_{j},c_{i})=\sum_{k}\left(\frac{h_{*}(z_{k};a_{j},b_{j},c_{i})-h(z_{k})}{\sigma(z)}\right)^{2}\,,$
(26)

where $h_{*}$ is given in Eq. (12). The factor $\sigma(z)$ is the fit
tolerance which can be chosen such that the fit is optimal at different
density regions. We choose to target similar relative uncertainty on the non-
rest-mass component of the enthalpy density $(h-1)\rho=p+{\epsilon\rho}$
across density scales: $\sigma(z)\propto\rho(z)\propto\exp(z)$. Overall, the
tolerance scales as $1/\rho$, so the fit is relatively better (with respect to
$h$) at low densities. The energy density at the stitching point is then
selected; if the tabulated EoS is sufficiently high-resolution, it can be
computed by, e.g., the trapezoidal rule. Otherwise, there is no canonical
choice for this value, we choose it to maximize agreement with tabulated
$p(\rho)$ at high densities.

Finally, the EoS fit is completed by stitching to some other EoS
parametrizaton at $\rho_{\mathrm{stitch}}=\rho_{\min}$. In the majority of
cases this is the spectral parametrization, though we also explore another
enthalpy segment in Sec. IV.3.2 and a polytrope in App. B. The low-density
spectral EoS itself transitions to a lower-density polytrope at some fixed
reference density $\rho_{r}$. Following Ref. Foucart  _et al._ (2019), we
define $x\equiv\log(\rho/\rho_{r})$ and write the spectral pressure as

$p_{s}(x)=\begin{cases}p_{0}\exp\left[\Gamma_{0}x\right]&x\leq 0\,,\\\
p_{0}\exp\left[\sum_{i=0}^{3}\frac{1}{i+1}\Gamma_{i}x^{i+1}\right]&x>0\,,\end{cases}$
(27)

where $p_{0}$ controls the overall pressure, and $\Gamma_{0},\Gamma_{i}$ are
the spectral coefficients. The low-density behavior fixes $\Gamma_{0}$, while
requiring a $\mathcal{C}^{1}$ transition to the enthalpy parametrization,
i.e., continuity in pressure, energy density, and pressure derivative fixes
$3$ more parameters. In practice because $e_{\min}=e_{\text{stitch}}$ is an
integration constant in the enthalpy parametrization, we can freely set it to
the value computed for $e_{\text{stitch}}$ from the low-density
parametrization, guaranteeing exact consistency.

The remaining $1$ degree of freedom is selected by either maximizing
smoothness across the lower-density transition to the polytrope or maximizing
accuracy of the low-density EoS. Smoothness is prioritized when the stitching
density is below the core density of typical NS. Then, we set $\Gamma_{1}=0.0$
Foucart  _et al._ (2019), guaranteeing that Eq. (27) is $\mathcal{C}^{2}$
across the transition $\rho_{r}$. This typically produces good fits to the
overall $M$–$R$ curve for the entire EoS. If the spectral parametrization is
stitched to the enthalpy parametrization at a higher density (near the core
density of astrophysical NSs as is the case in the high-stitch fit of
Sec.IV.2) we instead allow $\Gamma_{1}$ to vary, choosing it to maximize the
agreement of the total parametrized EoS with the target. Both strategies
typically result in machine-precision level $\mathcal{C}^{1}$-stitching to the
enthalpy parametrization, with residuals much smaller than mismodeling in the
low-density regime.

##  Appendix B Approximating a single polytrope

As an example of the strategy for fitting a target EoS with the enthalpy
parametrization, we consider a single-polytrope. In this case, the enthalpy
coefficients can be computed analytically. The general goal is to express the
EoS in the form of Eq. (III.1), i.e., compute the enthalpy as a function of
log-density.

The polytropic exponent is defined as

$\Gamma(z)\equiv\frac{d\log
p}{d\log\rho}=\frac{\rho}{p}\frac{dh}{dz}=\frac{\rho}{p}\frac{d\left(\frac{1}{\rho}\frac{de}{dz}\right)}{dz}\,.$
(28)

For a constant polytropic exponent $\Gamma(z)=\Gamma_{0}$ and using the
identity

$p(z)=h(z)\rho(z)-e(z)=\frac{de}{dz}-e(z)\,,$ (29)

Eq. (28) becomes

$\frac{d^{2}e}{dz^{2}}-(\Gamma_{0}+1)\frac{de}{dz}+\Gamma_{0}e=0\,.$ (30)

The solution to this differential equation is

$e(z)=(e_{0}-\rho_{0})\exp(\Gamma_{0}z)+\rho_{0}\exp(z)\,,$ (31)

where we have enforced $e(z=0)=e_{0}$ and $e(z\to-\infty)\to\rho(z)$999This is
equivalent to assuming the specific internal energy $\epsilon$ is $0$ in
ordinary, low-density cold matter. This can be done by defining the baryon
“rest” mass to be the average mass of a baryon in the outer crust of a NS
(despite the fact these baryons may be bound in, e.g. iron and therefore
differ from the mass of a free neutron/proton by up to $1\%$). and the
enthalpy is

$h(z)=\frac{1}{\rho}\frac{de}{dz}=\frac{e_{0}-\rho_{0}}{\rho_{0}}\Gamma_{0}\exp{\left[(\Gamma_{0}-1)z\right]}+1\,.$
(32)

Comparing with Eq. (9) the polynomial coefficients of the enthalpy expansion
are $\gamma_{i}=h_{0}\left(\Gamma_{0}-1\right)^{i}/i!$, with
$h_{0}=\Gamma_{0}\left(e_{0}-\rho_{0}\right)/\rho_{0}$, if $i\neq 0$, and
$\gamma_{0}=h_{0}\left(\Gamma_{0}-1\right)+1$. In practice, evaluating the
polynomial expansion of Eq. (9) requires many floating point operations.
Nonetheless, this computation is not necessarily slower than evaluating a
simple polytrope if $\Gamma_{0}$ is not an integer, because floating-point
exponentiation typically at least an order of magnitude slower than
multiplication and addition.

We use this enthalpy parametrization of the polytrope model to compare against
the direct single-polytrope SpECTRE implementation and verify the predicted
Cowling-approximation NS modes Font  _et al._ (2002); Deppe  _et al._ (2022c).
The low-density EoS in the enthalpy parametrization case is stitched to the
exact polytropic expression

$P(\rho)=\frac{100}{M_{\odot}^{-2.0}}\rho^{2.0}\,.$ (33)

We evolve a NS with central density $1.28\times 10^{-3}M_{\odot}^{-2}\approx
2.84\rho_{\mathrm{nuc}}$, which is the same as the stars evolved in Refs. Font
_et al._ (2002); Deppe  _et al._ (2022c). The number of terms necessary in the
polynomial expansion depends on the desired accuracy. For a resolution of
$130$ m we find that $i_{\max}=8$ is more than sufficient. This is consistent
with theoretical expectations, as the first neglected term, is of order
$1/9!\approx 2\times 10^{-6}$, indicating errors should be of this scale or
smaller. Results are shown in Fig. 16 where the enthalpy fit to the polytrope
and the direct single-polytropic parametrization return essentially identical
results. We display the run details in Table 2, as enthalpy-polytrope-mc-130
and polytropic-polytrope-mc-130.

Figure 16: NS central density as a function of time (top panel) and its
spectrum (bottom panel) for SpECTRE simulations with an enthalpy fit to a
single polytrope with $\Gamma_{0}=2.0$ (blue) and a direct single-polytropic
parametrization (red, dashed,). Known Cowling frequencies Font  _et al._
(2002) are marked as dashed vertical lines. The spectra are identical by eye.

We find effectively no difference in the runtime for the simulations using
each of the polytropic and enthalpy parametrizations. Examining the cost of
individual EoS calls in Table 4, we find the enthalpy parametrization is
somewhat faster in both $p(\rho)$ and $\epsilon(\rho)$ evaluation, indicating
that in this case, EoS evaluation time is not a significant contribution to
runtime. This speedup is also expected to extend to constant-speed-of-sound
matter, which has an identical functional form to a polytropic EoS when
expanded in the enthalpy / parametrization, the only difference being the
values of the coefficients.

The reason the polytropic EoS evaluation is not faster despite having a very
simple analytic expression is the inefficiency of floating-point
exponentiation. In the case of our test problem the floating point exponent
2.0 is only known during runtime, and so the compiler cannot optimize EoS
calls. If the exponent is known to be an integer at compile time, the calls
can be evaluated using repeated multiplication. We implement this improvement
for this particular test problem, and find the cost of polytropic $p(\rho)$
evaluations to be $10$ ns on identical hardware, indicating a 5-fold
improvement. Nonetheless, this speedup is not reflected in the total evolution
runtime for resolutions of at least $\sim 120$ m as the EoS evaluation cost is
subdominant to other simulation components.

Parametrization Evaluation Cost (ns) | $p(\rho)$ | $\epsilon(\rho)$
---|---|---
enthalpy | 43 | 44
polytrope | 57 | 58
Table 4: Performance (in nanoseconds) of the enthalpy and single-polytropic
fits to a single polytrope in evaluating the pressure and internal energy at
$\rho=5.0\times 10^{-4}M_{\odot}^{2}$. The enthalpy parametrization
outperforms then polytrope in both cases.

##  Appendix C Numerical considerations for enthalpy coefficients

Here we expand upon numerical considerations for choices of polynomial,
coefficients $\gamma_{i}$. The choice, $\gamma_{i}\geq 0$ in Eq. (9)
effectively bounds the number of terms in the polynomial expansion which can
be practically used. To see this, consider
$z_{\max}\equiv\log\left(\rho_{\max}/\rho_{0}\right)$. Coefficients
$\gamma_{i}$ must satisfy
$\gamma_{i}z_{\max}^{i}\lesssim\mathcal{O}(\gamma_{0})\sim 1$, otherwise they
would be larger than the total enthalpy in this region, which is typically
also of this scale. With this in mind, we consider $i$ as “too large” if $i\gg
z_{\max}$, as any term which satisfies
$\gamma_{i}z_{\max}^{i}\lesssim\gamma_{0}$ has

$\frac{\gamma_{i}z^{i}}{\gamma_{i}z_{\max}^{i}}=\left(\frac{z}{z_{\max}}\right)^{i}\,,$
(34)

small except when $z$ is nearly $z_{\max}$. That is, the degree of freedom is
only relevant at the highest densities, and this density region shrinks as $i$
gets larger. For typical scales, such as $\rho_{0}=0.5\rho_{\mathrm{nuc}}$,
$\rho_{\min}=\rho_{\mathrm{nuc}}$, and $\rho_{\max}=7\rho_{\mathrm{nuc}}$,
then $z_{\max}\sim\log(14)\approx 2.6$.

One can decrease $\rho_{0}$ to increase the relevant value of $z_{\max}$, but
this requires adding more parameters, which may not be desirable, since many
of them may be irrelevant, or degenerate. One way to view this, is that in the
Taylor expansion of the exponential function (equivalently the expansion of
$h(z)$ in a constant speed of sound case), the term $z^{i}/i!$ is the largest
term on $z\in(i-1,i)$, and is generally decreasing in relevance away from this
region relative to other terms. Therefore the $i$th term of this expansion is
most relevant near $z\lesssim i$, and is unimportant far from this region.
This implies that flexibility is essentially equidistributed in $\log(\rho)$
for this approximation (the same argument applies to the spectral
parametrization), and that higher polynomial terms cannot resolve features at
low densities. Instead, we choose to switch to a new function basis at this
point, optimized to capture the largest scale features at lowest order of
approximation.

## References

  * Abbott  _et al._ (2017a) B. P. Abbott _et al._ (LIGO Scientific Collaboration, Virgo Collaboration), Phys. Rev. Lett. 119, 161101 (2017a), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ (2017b) B. P. Abbott _et al._ (LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, Insight-Hxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GW-EM, DES, DLT40, GRAWITA, Fermi-LAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, J-GEM, GROWTH, JAGWAR, CaltechNRAO, TTU-NRAO, NuSTAR, Pan-STARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKI-GW Follow-up, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR, SKA South Africa/MeerKAT), Astrophys. J. Lett. 848, L12 (2017b), arXiv:1710.05833 [astro-ph.HE] .
  * Dietrich  _et al._ (2021) T. Dietrich, T. Hinderer, and A. Samajdar, Gen. Rel. Grav. 53, 27 (2021), arXiv:2004.02527 [gr-qc] .
  * Chatziioannou (2020) K. Chatziioannou, Gen. Rel. Grav. 52, 109 (2020), arXiv:2006.03168 [gr-qc] .
  * Piekarewicz (2022) J. Piekarewicz, (2022), arXiv:2209.14877 [nucl-th] .
  * Abbott  _et al._ (2020) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 892, L3 (2020), arXiv:2001.01761 [astro-ph.HE] .
  * Cromartie  _et al._ (2019) H. T. Cromartie _et al._ , Nature Astron. 4, 72 (2019), arXiv:1904.06759 .
  * Fonseca  _et al._ (2021) E. Fonseca _et al._ , Astrophys. J. Lett. 915, L12 (2021), arXiv:2104.00880 [astro-ph.HE] .
  * Miller  _et al._ (2019) M. C. Miller, C. Chirenti, and F. K. Lamb, (2019), arXiv:1904.08907 [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._ (2019) T. E. Riley _et al._ , Astrophys. J. Lett. 887, L21 (2019), arXiv:1912.05702 [astro-ph.HE] .
  * Riley  _et al._ (2021) T. E. Riley _et al._ , Astrophys. J. Lett. 918, L27 (2021), arXiv:2105.06980 [astro-ph.HE] .
  * Antoniadis  _et al._ (2013) J. Antoniadis, P. C. Freire, N. Wex, T. M. Tauris, R. S. Lynch, _et al._ , Science 340, 1233232 (2013), arXiv:1304.6875 [astro-ph.HE] .
  * Adhikari  _et al._ (2021) D. Adhikari _et al._ (PREX), Phys. Rev. Lett. 126, 172502 (2021), arXiv:2102.10767 [nucl-ex] .
  * Adhikari  _et al._ (2022) D. Adhikari _et al._ (CREX), (2022), arXiv:2205.11593 [nucl-ex] .
  * Abbott  _et al._ (2018) B. P. Abbott _et al._ (LIGO Scientific Collaboration, Virgo Collaboration), Phys. Rev. Lett. 121, 161101 (2018), arXiv:1805.11581 [gr-qc] .
  * Essick (2021) R. Essick, (2021), arXiv:2111.04244 [astro-ph.HE] .
  * Pang  _et al._ (2021) P. T. H. Pang, I. Tews, M. W. Coughlin, M. Bulla, C. Van Den Broeck, and T. Dietrich, Astrophys. J. 922, 14 (2021), arXiv:2105.08688 [astro-ph.HE] .
  * Raaijmakers  _et al._ (2020) G. Raaijmakers _et al._ , Astrophys. J. Lett. 893, L21 (2020), arXiv:1912.11031 [astro-ph.HE] .
  * Raaijmakers  _et al._ (2021) G. Raaijmakers, S. K. Greif, K. Hebeler, T. Hinderer, S. Nissanke, A. Schwenk, T. E. Riley, A. L. Watts, J. M. Lattimer, and W. C. G. Ho, (2021), arXiv:2105.06981 [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] .
  * Legred  _et al._ (2021) I. Legred, K. Chatziioannou, R. Essick, S. Han, and P. Landry, Phys. Rev. D 104, 063003 (2021).
  * Drischler  _et al._ (2016) C. Drischler, A. Carbone, K. Hebeler, and A. Schwenk, Phys. Rev. C 94, 054307 (2016), arXiv:1608.05615 [nucl-th] .
  * Antón  _et al._ (2006) L. Antón, O. Zanotti, J. A. Miralles, J. M. Martí, J. M. Ibáñez, J. A. Font, and J. A. Pons, The Astrophysical Journal 637, 296 (2006).
  * Font (2008) J. A. Font, Living Rev. Rel. 11, 7 (2008).
  * Baumgarte and Shapiro (2010) T. W. Baumgarte and S. L. Shapiro,  _Numerical Relativity: Solving Einstein’s Equations on the Computer_ (Cambridge University Press, 2010).
  * Baiotti and Rezzolla (2017) L. Baiotti and L. Rezzolla, Rept. Prog. Phys. 80, 096901 (2017), arXiv:1607.03540 [gr-qc] .
  * Radice  _et al._ (2020) D. Radice, S. Bernuzzi, and A. Perego, Ann. Rev. Nucl. Part. Sci. 70, 95 (2020), arXiv:2002.03863 [astro-ph.HE] .
  * Foucart (2020) F. Foucart, Front. Astron. Space Sci. 7, 46 (2020), arXiv:2006.10570 [astro-ph.HE] .
  * Kyutoku  _et al._ (2021) K. Kyutoku, M. Shibata, and K. Taniguchi, Living Rev. Rel. 24, 5 (2021), arXiv:2110.06218 [astro-ph.HE] .
  * Rezzolla and Zanotti (2013) L. Rezzolla and O. Zanotti,  _Relativistic Hydrodynamics_ (Oxford University Press, 2013).
  * Margalit and Metzger (2017) B. Margalit and B. D. Metzger, Astrophys. J. 850, L19 (2017), arXiv:1710.05938 [astro-ph.HE] .
  * Radice and Dai (2019) D. Radice and L. Dai, Eur. Phys. J. A 55, 50 (2019), arXiv:1810.12917 [astro-ph.HE] .
  * Shibata  _et al._ (2019) M. Shibata, E. Zhou, K. Kiuchi, and S. Fujibayashi, Phys. Rev. D 100, 023015 (2019), arXiv:1905.03656 [astro-ph.HE] .
  * Köppel  _et al._ (2019) S. Köppel, L. Bovard, and L. Rezzolla, Astrophys. J. Lett. 872, L16 (2019), arXiv:1901.09977 [gr-qc] .
  * Annala  _et al._ (2022) E. Annala, T. Gorda, E. Katerini, A. Kurkela, J. Nättilä, V. Paschalidis, and A. Vuorinen, Phys. Rev. X 12, 011058 (2022), arXiv:2105.05132 [astro-ph.HE] .
  * Camilletti  _et al._ (2022) A. Camilletti, L. Chiesa, G. Ricigliano, A. Perego, L. C. Lippold, S. Padamata, S. Bernuzzi, D. Radice, D. Logoteta, and F. M. Guercilena,  (2022), 10.1093/mnras/stac2333, arXiv:2204.05336 [astro-ph.HE] .
  * Typel  _et al._ (2022) S. Typel _et al._ , (2022), arXiv:2203.03209 [astro-ph.HE] .
  * Siegel  _et al._ (2018) D. M. Siegel, P. Mösta, D. Desai, and S. Wu, Astrophys. J. 859, 71 (2018), arXiv:1712.07538 [astro-ph.HE] .
  * Read  _et al._ (2009a) J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, Phys. Rev. D 79, 124032 (2009a), arXiv:0812.2163 [astro-ph] .
  * Foucart  _et al._ (2019) F. Foucart, M. D. Duez, A. Gudinas, F. Hebert, L. E. Kidder, H. P. Pfeiffer, and M. A. Scheel, Phys. Rev. D 100, 104048 (2019), arXiv:1908.05277 [gr-qc] .
  * Raithel and Paschalidis (2022) C. A. Raithel and V. Paschalidis, Phys. Rev. D 106, 023015 (2022), arXiv:2204.00698 [gr-qc] .
  * Lindblom (2010) L. Lindblom, Phys. Rev. D 82, 103011 (2010), arXiv:1009.0738 [astro-ph.HE] .
  * Greif  _et al._ (2019) S. K. Greif, G. Raaijmakers, K. Hebeler, A. Schwenk, and A. L. Watts, Mon. Not. Roy. Astron. Soc. 485, 5363 (2019), arXiv:1812.08188 [astro-ph.HE] .
  * Lindblom (1992) L. Lindblom, Astrophys. J. 398, 569 (1992).
  * Deppe  _et al._ (2022a) N. Deppe, W. Throwe, L. E. Kidder, N. L. Vu, F. Hébert, J. Moxon, C. Armaza, G. S. Bonilla, Y. Kim, P. Kumar, G. Lovelace, A. Macedo, K. C. Nelli, E. O’Shea, H. P. Pfeiffer, M. A. Scheel, S. A. Teukolsky, N. A. Wittek, _et al._ , “SpECTRE v2022.08.01,” 10.5281/zenodo.6949324 (2022a).
  * Kidder  _et al._ (2017) L. E. Kidder _et al._ , J. Comput. Phys. 335, 84 (2017), arXiv:1609.00098 [astro-ph.HE] .
  * Kale  _et al._ (2020) L. Kale, B. Acun, S. Bak, A. Becker, M. Bhandarkar, N. Bhat, A. Bhatele, E. Bohm, C. Bordage, R. Brunner, R. Buch, S. Chakravorty, K. Chandrasekar, J. Choi, M. Denardo, J. DeSouza, M. Diener, H. Dokania, I. Dooley, W. Fenton, J. Galvez, F. Gioachin, A. Gupta, G. Gupta, M. Gupta, A. Gursoy, V. Harsh, F. Hu, C. Huang, N. Jagathesan, N. Jain, P. Jetley, P. Jindal, R. Kanakagiri, G. Koenig, S. Krishnan, S. Kumar, D. Kunzman, M. Lang, A. Langer, O. Lawlor, C. W. Lee, J. Lifflander, K. Mahesh, C. Mendes, H. Menon, C. Mei, E. Meneses, E. Mikida, P. Miller, R. Mokos, V. Narayanan, X. Ni, K. Nomura, S. Paranjpye, P. Ramachandran, B. Ramkumar, E. Ramos, M. Robson, N. Saboo, V. Saletore, O. Sarood, K. Senthil, N. Shah, W. Shu, A. B. Sinha, Y. Sun, Z. Sura, E. Totoni, K. Varadarajan, R. Venkataraman, J. Wang, L. Wesolowski, S. White, T. Wilmarth, J. Wright, J. Yelon, and G. Zheng, “Uiuc-ppl/charm: Charm++ version 6.10.2,”  (2020).
  * Cowling (1941) T. G. Cowling, mnras 101, 367 (1941).
  * Deppe  _et al._ (2022b) N. Deppe, F. Hébert, L. E. Kidder, and S. A. Teukolsky, Classical and Quantum Gravity 39, 195001 (2022b).
  * Deppe  _et al._ (2022c) N. Deppe _et al._ , Phys. Rev. D 105, 123031 (2022c), arXiv:2109.12033 [gr-qc] .
  * Lindblom and Indik (2014) L. Lindblom and N. M. Indik, Phys. Rev. D89, 064003 (2014), [Erratum: Phys. Rev.D93,no.12,129903(2016)], arXiv:1310.0803 [astro-ph.HE] .
  * Han and Steiner (2019) S. Han and A. W. Steiner, Phys. Rev. D 99, 083014 (2019), arXiv:1810.10967 [nucl-th] .
  * Pang  _et al._ (2020) P. T. H. Pang, T. Dietrich, I. Tews, and C. Van Den Broeck, Phys. Rev. Res. 2, 033514 (2020), arXiv:2006.14936 [astro-ph.HE] .
  * Wysocki  _et al._ (2020) D. Wysocki, R. O’Shaughnessy, L. Wade, and J. Lange, (2020), arXiv:2001.01747 [gr-qc] .
  * Özel and Psaltis (2009) F. Özel and D. Psaltis, Phys. Rev. D 80, 103003 (2009), arXiv:0905.1959 [astro-ph.HE] .
  * Blackman  _et al._ (2015) J. Blackman, S. E. Field, C. R. Galley, B. Szilágyi, M. A. Scheel, M. Tiglio, and D. A. Hemberger, Phys. Rev. Lett. 115, 121102 (2015), arXiv:1502.07758 [gr-qc] .
  * Varma  _et al._ (2019) V. Varma, S. E. Field, M. A. Scheel, J. Blackman, L. E. Kidder, and H. P. Pfeiffer, Phys. Rev. D 99, 064045 (2019), arXiv:1812.07865 [gr-qc] .
  * Wijngaarden  _et al._ (2022) M. Wijngaarden, K. Chatziioannou, A. Bauswein, J. A. Clark, and N. J. Cornish, Phys. Rev. D 105, 104019 (2022), arXiv:2202.09382 [gr-qc] .
  * Breschi  _et al._ (2022) M. Breschi, S. Bernuzzi, K. Chakravarti, A. Camilletti, A. Prakash, and A. Perego, (2022), arXiv:2205.09112 [gr-qc] .
  * Shibata and Uryu (2000) M. Shibata and K. Uryu, Phys. Rev. D 61, 064001 (2000), arXiv:gr-qc/9911058 .
  * Etienne  _et al._ (2008) Z. B. Etienne, J. A. Faber, Y. T. Liu, S. L. Shapiro, K. Taniguchi, and T. W. Baumgarte, Phys. Rev. D 77, 084002 (2008), arXiv:0712.2460 [astro-ph] .
  * Baiotti  _et al._ (2008) L. Baiotti, B. Giacomazzo, and L. Rezzolla, Phys. Rev. D 78, 084033 (2008), arXiv:0804.0594 [gr-qc] .
  * Duez  _et al._ (2008) M. D. Duez, F. Foucart, L. E. Kidder, H. P. Pfeiffer, M. A. Scheel, and S. A. Teukolsky, Phys. Rev. D 78, 104015 (2008), arXiv:0809.0002 [gr-qc] .
  * Radice  _et al._ (2014) D. Radice, L. Rezzolla, and F. Galeazzi, Mon. Not. Roy. Astron. Soc. 437, L46 (2014), arXiv:1306.6052 [gr-qc] .
  * DeBuhr  _et al._ (2018) J. DeBuhr, B. Zhang, M. Anderson, D. Neilsen, E. W. Hirschmann, T. Grenga, and S. Paolucci, Astrophys. J. 867, 112 (2018), arXiv:1512.00386 [astro-ph.IM] .
  * Yagi and Yunes (2013) K. Yagi and N. Yunes, Phys. Rev. D 88, 023009 (2013), arXiv:1303.1528 [gr-qc] .
  * Read  _et al._ (2009b) J. S. Read, C. Markakis, M. Shibata, K. Uryu, J. D. Creighton, _et al._ , Phys. Rev. D 79, 124033 (2009b), arXiv:0901.3258 [gr-qc] .
  * Ujevic  _et al._ (2022) M. Ujevic, H. Gieg, F. Schianchi, S. V. Chaurasia, I. Tews, and T. Dietrich, (2022), arXiv:2211.04662 [gr-qc] .
  * Hotokezaka  _et al._ (2011) K. Hotokezaka, K. Kyutoku, H. Okawa, M. Shibata, and K. Kiuchi, Phys. Rev. D 83, 124008 (2011), arXiv:1105.4370 [astro-ph.HE] .
  * Lackey  _et al._ (2014) B. D. Lackey, K. Kyutoku, M. Shibata, P. R. Brady, and J. L. Friedman, Phys. Rev. D 89, 043009 (2014), arXiv:1303.6298 [gr-qc] .
  * Dietrich  _et al._ (2018) T. Dietrich, S. Bernuzzi, B. Bruegmann, and W. Tichy, in  _26th Euromicro International Conference on Parallel, Distributed and Network-based Processing_ (2018) pp. 682–689, arXiv:1803.07965 [gr-qc] .
  * Dietrich  _et al._ (2017) T. Dietrich, S. Bernuzzi, and W. Tichy, Phys. Rev. D96, 121501 (2017), arXiv:1706.02969 [gr-qc] .
  * O’Boyle  _et al._ (2020) M. F. O’Boyle, C. Markakis, N. Stergioulas, and J. S. Read, Phys. Rev. D 102, 083027 (2020), arXiv:2008.03342 [astro-ph.HE] .
  * Press  _et al._ (2007) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,  _Numerical Recipes 3rd Edition: The Art of Scientific Computing_ , 3rd ed. (Cambridge University Press, USA, 2007).
  * McLerran and Reddy (2019) L. McLerran and S. Reddy, Phys. Rev. Lett. 122, 122701 (2019), arXiv:1811.12503 [nucl-th] .
  * Tews  _et al._ (2018a) I. Tews, J. Carlson, S. Gandolfi, and S. Reddy, Astrophys. J. 860, 149 (2018a), arXiv:1801.01923 [nucl-th] .
  * Kapusta and Welle (2021) J. I. Kapusta and T. Welle, Phys. Rev. C 104, L012801 (2021), arXiv:2103.16633 [nucl-th] .
  * Landry and Kumar (2018) P. Landry and B. Kumar, Astrophys. J. 868, L22 (2018), arXiv:1807.04727 [gr-qc] .
  * Essick  _et al._ (2020a) R. Essick, P. Landry, and D. E. Holz, Phys. Rev. D 101, 063007 (2020a), arXiv:1910.09740 [astro-ph.HE] .
  * Mösta  _et al._ (2014) P. Mösta, B. C. Mundim, J. A. Faber, R. Haas, S. C. Noble, T. Bode, F. Löffler, C. D. Ott, C. Reisswig, and E. Schnetter, Class. Quant. Grav. 31, 015005 (2014), arXiv:1304.5544 [gr-qc] .
  * Tews (2020) I. Tews, Frontiers in Physics 8 (2020), 10.3389/fphy.2020.00153.
  * Drischler  _et al._ (2016) C. Drischler, A. Carbone, K. Hebeler, and A. Schwenk, Phys. Rev. C 94, 054307 (2016), arXiv:1608.05615 [nucl-th] .
  * Essick  _et al._ (2020b) R. Essick, I. Tews, P. Landry, S. Reddy, and D. E. Holz, Phys. Rev. C 102, 055803 (2020b), arXiv:2004.07744 [astro-ph.HE] .
  * Roca-Maza  _et al._ (2015) X. Roca-Maza, X. Viñas, M. Centelles, B. K. Agrawal, G. Colo’, N. Paar, J. Piekarewicz, and D. Vretenar, Phys. Rev. C 92, 064304 (2015), arXiv:1510.01874 [nucl-th] .
  * Essick  _et al._ (2021a) R. Essick, P. Landry, A. Schwenk, and I. Tews, Phys. Rev. C 104, 065804 (2021a).
  * Essick  _et al._ (2021b) R. Essick, I. Tews, P. Landry, and A. Schwenk, Phys. Rev. Lett. 127, 192701 (2021b), arXiv:2102.10074 [nucl-th] .
  * Tews  _et al._ (2018b) I. Tews, J. Margueron, and S. Reddy, Phys. Rev. C 98, 045804 (2018b), arXiv:1804.02783 [nucl-th] .
  * Landry and Essick (2019) P. Landry and R. Essick, Phys. Rev. D 99, 084049 (2019), arXiv:1811.12529 [gr-qc] .
  * Alford  _et al._ (2005) M. Alford, M. Braby, M. W. Paris, and S. Reddy, Astrophys. J. 629, 969 (2005), arXiv:nucl-th/0411016 .
  * Potekhin  _et al._ (2013) A. Y. Potekhin, A. F. Fantina, N. Chamel, J. M. Pearson, and S. Goriely, Astron. Astrophys. 560, A48 (2013), arXiv:1310.0049 [astro-ph.SR] .
  * Akmal  _et al._ (1998) A. Akmal, V. R. Pandharipande, and D. G. Ravenhall, Phys. Rev. C 58, 1804 (1998), arXiv:nucl-th/9804027 .
  * Lackey  _et al._ (2006) B. D. Lackey, M. Nayyar, and B. J. Owen, Phys. Rev. D 73, 024021 (2006), arXiv:astro-ph/0507312 [astro-ph] .
  * Baldo  _et al._ (1997) M. Baldo, I. Bombaci, and G. F. Burgio, Astron. Astrophys. 328, 274 (1997), arXiv:astro-ph/9707277 .
  * Engvik  _et al._ (1996) L. Engvik, G. Bao, M. Hjorth-Jensen, E. Osnes, and E. Ostgaard, Astrophys. J. 469, 794 (1996), arXiv:nucl-th/9509016 .
  * Müther  _et al._ (1987) H. Müther, M. Prakash, and T. L. Ainsworth, Phys. Lett. B 199, 469 (1987).
  * Müller and Serot (1996) H. Müller and B. D. Serot, Nucl. Phys. A606, 508 (1996), arXiv:nucl-th/9603037 [nucl-th] .
  * Rikovska-Stone  _et al._ (2007) J. Rikovska-Stone, P. A. M. Guichon, H. H. Matevosyan, and A. W. Thomas, Nucl. Phys. A 792, 341 (2007), arXiv:nucl-th/0611030 .
  * Douchin and Haensel (2001) F. Douchin and P. Haensel, Astron. Astrophys. 380, 151 (2001), arXiv:astro-ph/0111092 .
  * Wiringa  _et al._ (1988) R. B. Wiringa, V. Fiks, and A. Fabrocini, Phys. Rev. C 38, 1010 (1988).
  * Carney  _et al._ (2018) M. F. Carney, L. E. Wade, and B. S. Irwin, Phys. Rev. D98, 063004 (2018), arXiv:1805.11217 [gr-qc] .
  * Legred  _et al._ (2022) I. Legred, K. Chatziioannou, R. Essick, and P. Landry, Phys. Rev. D 105, 043016 (2022), arXiv:2201.06791 [astro-ph.HE] .
  * Gross-Boelting  _et al._ (1999) T. Gross-Boelting, C. Fuchs, and A. Faessler, Nucl. Phys. A 648, 105 (1999), arXiv:nucl-th/9810071 .
  * Alford  _et al._ (2015) M. G. Alford, G. F. Burgio, S. Han, G. Taranto, and D. Zappalà, Phys. Rev. D 92, 083002 (2015), arXiv:1501.07902 [nucl-th] .
  * Courant  _et al._ (1967) R. Courant, K. Friedrichs, and H. Lewy, IBM Journal of Research and Development 11, 215 (1967).
  * Van Leer (1977) B. Van Leer, Journal of Computational Physics 23, 276 (1977).
  * Deppe  _et al._ (tion) N. Deppe _et al._ , (in preparation).
  * Oppenheimer and Volkoff (1939) J. Oppenheimer and G. Volkoff, Phys. Rev. 55, 374 (1939).
  * Lattimer and Prakash (2001) J. M. Lattimer and M. Prakash, Astrophys. J. 550, 426 (2001), arXiv:astro-ph/0002232 .
  * Chatziioannou and Han (2020) K. Chatziioannou and S. Han, Phys. Rev. D 101, 044019 (2020), arXiv:1911.07091 [gr-qc] .
  * Gieg  _et al._ (2019) H. Gieg, T. Dietrich, and M. Ujevic, Particles 2, 365 (2019), arXiv:1908.03135 [gr-qc] .
  * Kurkela  _et al._ (2010) A. Kurkela, P. Romatschke, and A. Vuorinen, Phys. Rev. D 81, 105021 (2010), arXiv:0912.1856 [hep-ph] .
  * Kokkotas and Ruoff (2001) K. D. Kokkotas and J. Ruoff, Astron. Astrophys. 366, 565 (2001), arXiv:gr-qc/0011093 .
  * Sen  _et al._ (2022) S. Sen, S. Kumar, A. Kunjipurayi, P. Routaray, T. Zhao, and B. Kumar, (2022), arXiv:2205.02076 [nucl-th] .
  * Kastaun  _et al._ (2021) W. Kastaun, J. V. Kalinani, and R. Ciolfi, Phys. Rev. D 103, 023018 (2021), arXiv:2005.01821 [gr-qc] .
  * Espino and Paschalidis (2022) P. L. Espino and V. Paschalidis, Phys. Rev. D 105, 043014 (2022), arXiv:2105.05269 [astro-ph.HE] .
  * Carbone and Schwenk (2019) A. Carbone and A. Schwenk, Phys. Rev. C 100, 025805 (2019), arXiv:1904.00924 [nucl-th] .
  * Raithel  _et al._ (2019) C. A. Raithel, F. Ozel, and D. Psaltis, Astrophys. J. 875, 12 (2019), arXiv:1902.10735 [astro-ph.HE] .
  * Deppe  _et al._ (2023) N. Deppe, W. Throwe, L. E. Kidder, N. L. Vu, F. Hébert, J. Moxon, C. Armaza, G. S. Bonilla, Y. Kim, P. Kumar, G. Lovelace, A. Macedo, K. C. Nelli, E. O’Shea, H. P. Pfeiffer, M. A. Scheel, S. A. Teukolsky, N. A. Wittek, I. Anantpurkar, M. Boyle, A. Carpenter, A. Ceja, H. Chaudhary, F. Foucart, N. Ghadiri, M. Giesler, J. S. Guo, D. A. B. Iozzo, I. Legred, D. Li, S. Ma, D. Melchor, M. Morales, E. R. Most, M. A. Pajkos, T. Ramirez, N. Ring, H. R. Rüter, J. Sanchez, L. C. Stein, S. Thomas, D. Vieira, T. Wlodarczyk, and D. Wu, “Spectre,”  (2023).
  * SXS (2023) SXS, “paper-2023-spectre-enthalpy-eos,” https://github.com/sxs-collaboration/paper-2023-spectre-enthalpy-eos (2023).
  * Virtanen  _et al._ (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, Nature Methods 17, 261 (2020).
  * Oliphant (06 ) T. Oliphant, “NumPy: A guide to NumPy,” USA: Trelgol Publishing (2006–), [Online; accessed ¡today¿].
  * Hunter (2007) J. D. Hunter, Computing in Science & Engineering 9, 90 (2007).
  * Ayachit (2015) U. Ayachit,  _The ParaView Guide: A Parallel Visualization Application_ (Kitware, Inc., Clifton Park, NY, USA, 2015).
  * Font  _et al._ (2002) J. A. Font, T. Goodale, S. Iyer, M. A. Miller, L. Rezzolla, E. Seidel, N. Stergioulas, W.-M. Suen, and M. Tobias, Phys. Rev. D 65, 084024 (2002), arXiv:gr-qc/0110047 .