# Realistic binary neutron star initial data with Elliptica

Alireza Rashti1,2, Andrew Noe1,2 1 Institute for Gravitation & the Cosmos, The
Pennsylvania State University, University Park PA 16802, USA
2 Department of Physics, The Pennsylvania State University, University Park,
PA 16802, USA

###### Abstract

This work introduces the Elliptica pseudo-spectral code for generating initial
data of binary neutron star systems. Building upon the recent Elliptica code
update, we can now construct initial data using not only piecewise polytropic
equations of state, but also tabulated equations of state for these binary
systems. Furthermore, the code allows us to endow neutron stars within the
binary system with spins. These spins can have a magnitude close to the mass
shedding limit and can point in any direction.

Keywords: initial data, binary neutron star, elliptic solver, equation of
state

##  1 Introduction

Binary Neutron Star (BNS) systems are very common in our Universe. Estimates
suggest their merger rate falls between $250$ to $2810$ Gpc-3yr-1 [1]. The
coalescence of these BNS systems is a source of myriad phenomena such as,
among others, ejecta [2, 3], accretion disk [4], jets [5], r-process
nucleosynthesis [6], and kilonova [7]. These phenomena are treasure trove of
information revealing aspects of physics at large scales, for instance,
gravity in strong regimes and physical constants of the Universe, as well as
aspects of physics regarding small scales, like, Equation of State (EOS) in
supranuclear dense matter and formation of heavy elements in the periodic
table.

As such, current detectors such as LIGO [8], VIRGO [9], and KAGRA [10] and the
next generation ones like Cosmic Explorer [11], the DECi-hertz Interferometer
Gravitational-wave Observatory (DECIGO) [12], Einstein Telescope [13], LIGO
Voyager [14], the Laser Interferometer Space Antenna (LISA) [15], NEMO [16],
and TianQin [17] are designed to look eagerly into sky and observe the
physical signals emitted from the coalescence of compact binaries.

To unlock the wealth of information encoded in gravitational waves and their
electromagnetic counterparts, accurate theoretical models are crucial. These
models are essential for understanding events like GW170817 [18], the short
gamma-ray burst GRB170817A [19], and the kilonova transient AT2017gfo [20].
However, finding analytical solution of the governing Partial Differential
Equations, when BNS systems are coalescing, is not feasible as the PDE are in
a highly non-linear regime where no approximation is applicable [21]. In light
of this, Numerical Relativity (NR) community have put significant efforts to
solve these equations numerically and hence make sense of the observations.

Simulation of compact binary system in NR often involves two steps. The first
step is to find the solution of Einstein-Euler PDE on a hypersurface of the
spacetime manifold, namely, constructing constraint satisfying and self
consistent Initial Data (ID) that present the binary system of interest at
some time. ID codes such as, COCAL [22, 23, 24], Elliptica [25], FUKA [26],
LORENE [27, 28, 29, 30, 31], NRPyElliptic [32], SGRID [33, 34, 35], SpECTRE’s
elliptic solver [36, 37], Spells [38, 39, 40, 41], TwoPunctures [42, 43], are
developed and utilized for this step.

The next step is to use the ID and simulate the system’s evolution over time;
hence we can find the solution over the spacetime of interest. Dynamical
evolution codes, for instance, BAM [44, 45, 46], BAMPS [47, 48, 49], Dendro-GR
[50], Einstein Toolkit [51, 52], ExaHyPE [53] GR-Athena++ [54, 55, 56] GRaM-X
[57, 58], GRChombo [59, 60], Nmesh [61], NRPy+ [62], SACRA-MPI [63], Simflowny
[64], SpEC [65, 66], SpECTRE [67, 68], and SPHINCS_BSSN [69], among others,
are made and employed for this step.

Previously, Elliptica [25] was limited to the construction of ID for black
hole neutron star binary systems. Additionally, the code was only supporting
polytropic or piecewise polytropic EOS. In this work, we extend Elliptica’s
infrastructure to construction ID for BNS systems as well as supporting
tabulated EOS, for instance, the CompOSE tables [70, 71, 72].

To this aim, the Einstein-Euler equations are cast into coupled non-linear
elliptic PDE [73, 74, 34], and solved iteratively using Newton-Raphson method
by the Schur domain decomposition [25]. During the solve, the cubed spherical
patches adapt to the surface of the Neutron Stars thus separating matter and
vacuum and preventing Gibbs phenomena in the spectral method. To achieve the
target values of interest such as momenta of the system, mass and center of
NS, they are checked and adjusted during the solve. The EOS are approximated
by either (piecewise) polytropic or tabulate ones.

The remainder of paper is organized as follows. In section 2 we explain the
mathematical background and formalism we use for ID of BNS systems, in
particular, Einstein-Euler equations and EOS in Elliptica. Section 3 details
the underling algorithms and numerical techniques for construction of BNS’s
ID. We present the implementation of tabulated EOS by spline interpolant means
and diagnostics in Elliptica. Additionally, we explain the iterative procedure
for finding the physical and constraint satisfying ID for BNS systems. In
section 4, we present various convergence tests and comparison against post
Newtonian answers to showcase the proof of concept for the new version the
code. In section 5, we discuss the possible improvements and future work.

Throughout the article, we use geometric units where $G=1$ is the constant of
gravity, $c=1$ is the speed of light, and solar mass $M_{\odot}=1$.

##  2 Formalism

###  2.1 Einstein-Euler equations

To derive the equations that govern the gravity and matter on a spatial-like
hypersurface $\Sigma_{t}$ of a 4-dimensional manifold $\mathcal{M}$ with a
metric $g_{\mu\nu}$, we first write the line element of $\mathcal{M}$ as

$ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}=-\alpha^{2}dt^{2}+\gamma_{ij}(dx^{i}+\beta^{i}dt)(dx^{j}+\beta^{j}dt).$
(1)

This specific form of the line element is particularly well-suited for the
$3+1$ formalism. It provides a clear view of the key variables (fields):
$\alpha$, $\beta^{i}$, and $\gamma_{ij}$. Here, $\alpha$ is the lapse gauge
and indicates the way that a sequence of spatial slices, $\Sigma_{t}$s, are
combined to form the complete spacetime manifold $M$. $\gamma_{ij}$ is the
induced 3-metric on each spatial hypersurface and can be written as
$\gamma_{\mu\nu}=g_{\mu\nu}+n_{\mu}n_{\nu}$, where $n^{\mu}$ is the normal
vector on $\Sigma_{t}$. $\beta^{i}$ is the shift vector. The importance of
this gauge, i.e., shift vector, is that it represents the coordinate frame
being used in each hypersurface. In this work, we take the shift vector as
follows

$\beta^{i}=B^{i}+\epsilon_{ijk}{\Omega}^{j}_{\mathrm{BNS}}(r^{k}-r^{k}_{\mathrm{CM}})+\frac{v_{r}}{r_{\mathrm{BNS}}}(r^{i}-r^{i}_{\mathrm{CM}}).$
(2)

Here, $\Omega^{i}_{\mathrm{BNS}}$ is the orbital angular velocity of the BNS
system, $r^{i}_{\mathrm{CM}}$ is the position of the system’s center of mass,
$r_{\mathrm{BNS}}$ is the coordinate distance between the NS centers, $v_{r}$
is the radial velocity of the inspiraling coordinate system, and
$\epsilon_{ijk}$ is the Levi-Civita symbol. Eq. (2) proves numerically
convenient when applying Boundary Condition (BC) for $\beta^{i}$ at the edge
of the computational grid where the position vector $\vec{r}$ has large
values.

By definition the extrinsic curvature on $\Sigma_{t}$ is
$K_{\mu\nu}=-\frac{1}{2}\pounds_{n}\gamma_{\mu\nu},$ in which $\pounds_{n}$ is
the Lie derivative along the normal vector. We note that by construction
$K_{\mu\nu}n^{\mu}=0$, therefore, we can use spatial indices to describe the
extrinsic curvature. Next, by utilizing conformal decomposition, we write

$\displaystyle\gamma_{ij}=$ $\displaystyle\psi^{4}\bar{\gamma}_{ij},$ (3)
$\displaystyle K^{ij}=$ $\displaystyle A^{ij}+\frac{1}{3}K\gamma^{ij}.$ (4)

Here $\psi$ is the conformal factor, $\bar{\gamma}_{ij}$ the conformal
3-metric, $A^{ij}$ the traceless part of $K^{ij}$. Moreover, since we are
using the Extended Conformal Thin Sandwich (XCTS) formalism [74, 73], $A^{ij}$
is decomposed as

$\displaystyle A^{ij}$ $\displaystyle=\psi^{-10}\bar{A}^{ij},$ (5)
$\displaystyle\bar{A}^{ij}$
$\displaystyle=\frac{1}{2\bar{\alpha}}\left((\bar{L}\beta)^{ij}-\bar{\gamma}^{ik}\bar{\gamma}^{jl}\bar{u}_{kl}\right),$
(6)

where

$\displaystyle\bar{u}_{ij}$
$\displaystyle=\frac{\partial\bar{\gamma}_{ij}}{\partial t},$ (7)
$\displaystyle(\bar{L}\beta)^{ij}$
$\displaystyle=\bar{D}^{i}\beta^{j}+\bar{D}^{j}\beta^{i}-\frac{2}{3}\bar{\gamma}^{ij}\bar{D}_{k}\beta^{k},$
(8) $\displaystyle\alpha$ $\displaystyle=\psi^{6}\bar{\alpha},$ (9)

and $\bar{D}$ is the covariant derivative compatible with $\bar{\gamma}_{ij}$.

To complete the formulation we need to incorporate the source terms. We assume
the fluid in NS are governed by ideal fluid, hence the stress energy tensor
can be written as

$\displaystyle T_{\mu\nu}$ $\displaystyle=$
$\displaystyle(\rho_{0}+\rho_{0}\epsilon+P)u_{\mu}u_{\nu}+Pg_{\mu\nu},$
$\displaystyle=$ $\displaystyle\rho_{0}hu_{\mu}u_{\nu}+Pg_{\mu\nu},$

where, $\rho_{0}$ is the rest mass density, $\epsilon$ the specific internal
energy, $P$ the pressure, $h$ the specific enthalpy, and $u^{\mu}$ the
4-velocity of the fluid. Additionally, for $3+1$ decomposition purposes, we
project the stress energy tensor with respect to $\Sigma_{t}$ as follows

$\displaystyle E$ $\displaystyle=n_{\mu}n_{\nu}T^{\mu\nu},$ (11)
$\displaystyle S$
$\displaystyle=\gamma^{ij}\gamma_{i\mu}\gamma_{j\nu}T^{\mu\nu},$ (12)
$\displaystyle j^{i}$ $\displaystyle=-\gamma^{i}_{\mu}n_{\nu}T^{\mu\nu},$ (13)

where, $E$ is the measured energy by the Eulerian observer whose 4-velocity is
$n^{\mu}$. $S$ is the trace of matter stress tensor, and $j^{i}$ is the
momentum flux.

Unique answer to a linear elliptic equation with a source is guaranteed by the
maximum/minimum principle. This principle becomes important during solve of
constraint equations for high mass NS (for further discussion see [75]). To
maintain the maximum/minimum principle, we rescale the stress energy
projections as

$\displaystyle E$ $\displaystyle=\psi^{-6}\bar{E},$ (14) $\displaystyle S$
$\displaystyle=\psi^{-6}\bar{S},$ (15) $\displaystyle j^{i}$
$\displaystyle=\psi^{-6}\bar{j}^{i}.$ (16)

Finally, following XCTS formalism, we write Einstein’s equations in quasi
equilibrium condition

$\displaystyle\bar{D}^{2}\psi-\frac{1}{8}\psi\bar{R}-\frac{1}{12}\psi^{5}K^{2}+\frac{1}{8}\psi^{-7}\bar{A}_{ij}\bar{A}^{ij}+2\pi\psi^{-1}\bar{E}$
$\displaystyle=$ $\displaystyle 0,$ (17)
$\displaystyle\bar{D}^{2}(\bar{\alpha}\psi^{7})-(\bar{\alpha}\psi^{7})\left[\frac{1}{8}\bar{R}+\frac{5}{12}\psi^{4}K^{2}+\frac{7}{8}\psi^{-8}\bar{A}_{ij}\bar{A}^{ij}\right]$
$\displaystyle+\psi^{5}(\partial_{t}K-\beta^{k}\partial_{k}K)-2\pi\bar{\alpha}\psi^{5}(\bar{E}+2\bar{S})$
$\displaystyle=$ $\displaystyle 0,$ (18) $\displaystyle
2\bar{\alpha}\left[\bar{D}_{j}(\frac{1}{2\bar{\alpha}}(\bar{L}\beta)^{ij})-\bar{D}_{j}(\frac{1}{2\bar{\alpha}}\bar{u}^{ij})-\frac{2}{3}\psi^{6}\bar{D}^{i}K\right]-16\pi\bar{\alpha}\psi^{4}\bar{j}^{i}$
$\displaystyle=$ $\displaystyle 0,$ (19)

with the BC

$\lim_{r\to\infty}\psi=1,\ \ \ \lim_{r\to\infty}B^{i}_{0}=0,\ \ \
\lim_{r\to\infty}\alpha\psi=1.$ (20)

Additionally, we pick the free data as

$\displaystyle\bar{\gamma}_{ij}$ $\displaystyle=\delta_{ij},$ (21)
$\displaystyle K$ $\displaystyle=0,$ (22) $\displaystyle\bar{u}_{ij}$
$\displaystyle=0.$ (23)

For hydrodynamic equations, following [34], we decompose the fluid into two
parts: the rotational part of the fluid and is represented by a cross product,
and the irrotational part of the fluid that is represented by a velocity
potential (see, e.g.,[76]). In particular, the rotational part, which
represents the NS spin, reads

$w^{i}=\epsilon_{ijk}\Omega^{j}_{\mathrm{NS}}(x^{k}-x^{k}_{c}),$ (24)

here, $\Omega^{j}_{\mathrm{NS}}$ is a free parameter to set the spin level,
see sec. 4.3, and $x^{k}_{c}$ denotes the NS coordinate. The irrotational part
is shown by the potential $\phi(x,y,z)$ and obeys the following equations.
[35]

$\displaystyle\frac{c\left(\rho_{0}\right)\alpha}{h}\psi^{-4}\bar{\gamma}^{ij}\partial_{i}\partial_{j}\phi-\frac{\rho_{0}\alpha}{h}\psi^{-4}\bar{\gamma}^{ij}\bar{\Gamma}^{k}_{ij}\partial_{k}\phi+2\frac{\rho_{0}\alpha}{h}\psi^{-5}\bar{\gamma}^{ij}(\partial_{i}\psi)(\partial_{j}\phi)$
$\displaystyle+\left(D_{i}\frac{\rho_{0}\alpha}{h}\right)\left(D^{i}\phi\right)D_{i}\left[\frac{\rho_{0}\alpha}{h}w^{i}-\rho_{0}\alpha
u^{0}(\beta^{i}+\xi^{i})\right]=0,$ (25)

where, $\partial_{i}$ denotes the spatial partial derivative with respect to
Cartesian coordinate $x^{i}$ and $\xi^{i}$ is the Killing vector. Here,
following the smoothing method in [35]

$c(\rho_{0})=\rho_{0}+\bar{\epsilon}\rho_{0c}\left(\frac{\rho_{0c}-\rho_{0}}{\rho_{0c}}\right)^{4},$
(26)

in which, $\rho_{0c}$ is the values of $\rho_{0}$ at the NS center, and
$\bar{\epsilon}$ is a constant number – for which we generally use $0.1$.
Since the fluid is enclosed inside the patches that cover the NS, we further
need to impose a BC on the NS surface as

$D^{i}\phi
D_{i}\rho_{0}+w^{i}D_{i}\rho_{0}-hu^{0}(\beta^{i}+\xi^{i})D_{i}\rho_{0}=0.$
(27)

We note that since eqs. (25,27) only include derivatives of $\phi$, we can not
uniquely determine $\phi$. Hence, we further demand that the value of $\phi$
at the NS center is a constant number – like 0.

###  2.2 Equation of state

In order to close the system of equations we need an EOS. Elliptica deploys
specific enthalpy, $h$, to create a link between the macroscopic properties of
the matter and the gravity. As such, if we have a piecewise EOS we write

$\displaystyle\rho_{0}(h)$
$\displaystyle=K_{i}^{-n_{i}}\left(\frac{h-1-a_{i}}{n_{i}+1}\right)^{n_{i}},$
$\displaystyle P(h)$
$\displaystyle=K_{i}^{-n_{i}}\left(\frac{h-1-a_{i}}{n_{i}+1}\right)^{n_{i}+1},$
$\displaystyle\epsilon(h)$ $\displaystyle=\frac{a_{i}+n_{i}(h-1)}{n_{i}+1},$
(28)

where $n_{i}=\frac{1}{\Gamma_{i}-1}$ is the polytropic index and $K_{i}$’s are
specific to the given EOS; $a_{i}$’s ensure the continuity of EOS [77] and are
set as

$\displaystyle a_{0}$ $\displaystyle=0,$ $\displaystyle a_{i}$
$\displaystyle=a_{i-1}+\frac{K_{i-1}}{\Gamma_{i-1}-1}\rho_{i}^{\Gamma_{i-1}-1}-\frac{K_{i}}{\Gamma_{i}-1}\rho_{i}^{\Gamma_{i}-1}.$
(29)

Similarly, for tabulated EOS we represent $\rho_{0}\left(h\right)$,
$P\left(h\right)$, and $\epsilon\left(h\right)$ as spline interpolants given
discrete data points $\left(\rho_{0},P,\epsilon,h\right)$. Although typically
only the quantities $\left(\rho_{0},P,\epsilon\right)$ or
$\left(\rho_{0},P,e\right)$ are provided in tables, where $e$ is the total
energy density, we may convert between variables via the relations:

$\displaystyle e$ $\displaystyle=\rho_{0}\left(1+\epsilon\right),$
$\displaystyle h$ $\displaystyle=1+\epsilon+\frac{P}{\rho_{0}}.$ (30)

Having written thermodynamic variables in terms of specific enthalpy, we now
write specific enthalpy in terms of the metric and fluid variables [34]

$\displaystyle h$ $\displaystyle=$
$\displaystyle\sqrt{L^{2}-(D_{i}\phi+w_{i})(D^{i}\phi+w^{i})},$ $\displaystyle
L^{2}$ $\displaystyle=$
$\displaystyle\frac{b+\sqrt{b^{2}-4\alpha^{4}[(D_{i}\phi+w_{i})w^{i}]^{2}}}{2\alpha^{2}},$
$\displaystyle b$ $\displaystyle=$
$\displaystyle[(\xi^{i}+\beta^{i})D_{i}\phi-C]^{2}+2\alpha^{2}(D_{i}\phi+w_{i})w^{i}.$
(31)

##  3 Numerical method

###  3.1 Coordinate system

The computation grid is tiled by cubed spherical coordinate systems [78],
except at the NS centers, where we use simple Cartesian boxes to avoid
coordinate singularities at $r=0$ [79, 25]. The relation between Cartesian
coordinate, denoted by $x^{i}=(x,y,z)$, and cubed spherical coordinate,
$X^{i}=(X,Y,Z)$, is

$\displaystyle
X(x,y,z)=\frac{x}{z},\;Y(x,y,z)=\frac{y}{z},\;Z(x,y,z)=\frac{z-r_{\mathrm{in}}}{r_{\mathrm{out}}-r_{\mathrm{in}}},$
(32)

here, $X$ and $Y$ take value $\in[-1,1]$, and $Z\in[0,1]$. Additionally,
$r_{\mathrm{in}}$ and $r_{\mathrm{out}}$ are defined

$\displaystyle r_{\mathrm{in}}$
$\displaystyle=\frac{\sigma_{\mathrm{in}}(X,Y)}{\sqrt{1+X^{2}+Y^{2}}},\;r_{\mathrm{out}}$
$\displaystyle=\frac{\sigma_{\mathrm{out}}(X,Y)}{\sqrt{1+X^{2}+Y^{2}}}.$ (33)

where, the shape of the inner boundary of a patch along the radial direction
is determined by $\sigma_{\mathrm{in}}(X,Y)$ and its outer boundary by
$\sigma_{\mathrm{out}}(X,Y)$. $\sigma(X,Y)$ is related to Cartesian
coordinates by the equation $\sigma(X,Y)=\sqrt{x^{2}+y^{2}+z^{2}}$. Finally,
we note that while eq. (32) is written for patches along the z-axis, one can
generalize this along any other axes.

As the distance from NS increases, we expect the fields fall as powers of
$r^{-1}$. To account for this behavior we use a new transformation for $Z$
coordinate, denoted by $\widetilde{Z}$, for the paches covering large radii of
the grid. The transformation reads

$\widetilde{Z}=\frac{\sigma_{\mathrm{out}}}{\sigma_{\mathrm{out}}-\sigma_{\mathrm{in}}}\left(1-\frac{\sigma_{\mathrm{in}}}{r}\right),$
(34)

Here $r=\sqrt{x^{2}+y^{2}+z^{2}}$ and still $\widetilde{Z}\in[0,1]$.

Finally, Elliptica uses Chebyshev polynomials of the first kind for the basis
of the spectral expansion and deploys the extrema of Chebyshev polynomials for
its collocation points [25].

Figure 1: Grid patches for a BNS system. Left: shown is the $x-y$ plane of the
computational grid that is covered by different patches. Right: a zoomed-in
view of the $x-y$ plane focusing on the NS regions. A Cartesian box is used
around each NS center to avoid coordinate singularities. By using various
$\sigma(X,Y)$ values, cubed spherical patches can adapt to different shapes,
effectively capturing the shape of the NS surface. This allows for the
treatment of matter and vacuum in separate patches.

Fig. (1) illustrates the patches covering the $x-y$ plane of the computational
grid used in this work.

###  3.2 Elliptic solver

Given an elliptic PDE, we linearize the equation to use the Newton-Raphson
method (see [25] for a complete description). Thus, the original elliptic PDE
becomes a matrix equation, like $Ax=B$, to be solved. Now the challenge is to
solve this matrix equation efficiently. To this aim, we use the Schur domain
decomposition method [80, 25]. At the core of this method there is a divide
and conquer strategy; in this strategy, first the coupled equations are solved
and then the system decomposes into independent (uncoupled) subsystems,
amenable for parallel solving. The Schur domain decomposition allows to solve
the full system using shared memory multiprocessing in which each CPU core is
assigned to solve a matrix equations at each patch independently.

The Schur domain decomposition method arranges the Jacobian matrix in Newton-
Raphson method into two parts. The first part is associated with all uncoupled
equations and the second part comprises of the coupled equations. These
couplings between equations occur due to the inter-BC at the interfaces of
adjacent patches [25]. As such, the system of matrix equations can be seen as
two equations with two unknowns as follows

$\left(\begin{array}[]{c c}B&E\\\
F&C\end{array}\right)\left(\begin{array}[]{c}\tilde{v}\\\ \tilde{w}\\\
\end{array}\right)=\left(\begin{array}[]{c}f\\\ g\\\ \end{array}\right),$ (35)

where, the vector $\tilde{v}$ denotes all the unknowns that are uncoupled,
i.e., $\tilde{v}$ comprises all unknowns that stem from inner points of a
patch and hence uncoupled from the other patches. The vector $\tilde{w}$ is
the unknowns that are due to coupled equations, i.e., stem from inter-BC
between adjacent patches. By such an arrangement, if the unknown $\tilde{w}$
is solved, then the equation of $\tilde{v}$ can be solved consequently.

Therefore, we first solve for $\tilde{w}$ and then we use it to solve for
$\tilde{v}$. Accordingly, in Elliptica, we solve eq. 35 like:

$\displaystyle(C-FE^{\prime})\tilde{w}$ $\displaystyle=g-Ff^{\prime},$ (36)
$\displaystyle\tilde{v}$ $\displaystyle=f^{\prime}-E^{\prime}\tilde{w},$ (37)

where

$\displaystyle E^{\prime}$ $\displaystyle=B^{-1}E,$ $\displaystyle f^{\prime}$
$\displaystyle=B^{-1}f.$ (38)

As we can see eq (36) only involves the unknown $\tilde{w}$ and it can be
summarized as

$S\tilde{w}=g^{\prime},$ (39)

where, $S$ is called Schur complement matrix.

After solving eq. (39) for $\tilde{w}$, we solve for $\tilde{v}$ using eq.
(37) – hence we find the solution of the whole system.

Algorithm 1 Schur complement domain decomposition method.

1:Solve $BE^{\prime}=E$ for $E^{\prime}$;

2:Solve $Bf^{\prime}=f$ for $f^{\prime}$;

3:Compute $g^{\prime}=g-Ff^{\prime}$;

4:Compute $S=(C-FE^{\prime})$;

5:Solve $S\tilde{w}=g^{\prime}$ for $\tilde{w}$;

6:Compute $\tilde{v}=f^{\prime}-E^{\prime}\tilde{w}$;

A summary of the Schur complement domain method to solve elliptic equations is
shown in algorithm 1. For a more in-depth discussion about implementation of
Schur domain decomposition and parallelization of matrix solver the reader may
consult [25]. Finally, we use the publicly available and open-source UMFPACK
direct solver [81].

###  3.3 Diagnostics

In Elliptica we compute the baryon mass of each NS using [82]

$M_{B}=\int_{\mathrm{NS}}\rho_{0}\alpha\psi^{6}\sqrt{\bar{\gamma}}d^{3}x,$
(40)

in which $\bar{\gamma}$ is the determinant of $\bar{\gamma}_{ij}$ and the
integration is taken over the volume of the NS.

To measure the NS spins, we can use the flat space coordinate rotational
Killing vector, following [83], on the surface of NS:

$\displaystyle\vec{\phi}_{x}$
$\displaystyle=-(z-z_{c})\vec{\partial}_{y}+(y-y_{c})\vec{\partial}_{z},$
$\displaystyle\vec{\phi}_{y}$
$\displaystyle=+(z-z_{c})\vec{\partial}_{x}+(x-x_{c})\vec{\partial}_{z},$
$\displaystyle\vec{\phi}_{z}$
$\displaystyle=-(y-y_{c})\vec{\partial}_{x}+(x-x_{c})\vec{\partial}_{y},$ (41)

in which $(x_{c},y_{c},z_{c})$ is the coordinate center of the NS and
$(\vec{\partial}_{i})_{i\in\\{x,y,z\\}}$ are the basis vectors associated with
the Cartesian coordinates. Accordingly, NS spin $S_{i}$, for the $i$
direction, is measured on the surface of NS by

$S_{i}=\frac{1}{8\pi}\oint_{\mathrm{NS}}(\vec{\phi}_{i})^{j}s^{k}K_{jk}dA,$
(42)

where, $s^{k}$ is the normal vector on the NS surface. The dimensionless spin
is defined

$\chi_{i}:=\frac{S_{i}}{M^{2}_{\mathrm{NS}}}.$ (43)

Additionally, following [79], we can first compute the angular momentum
$J_{i}$, the center $R^{i}_{c}$, and linear momentum $P_{i}$ of the NS, and
then calculate $S_{i}$ as follows

$S_{i}=J_{i}-\epsilon_{ijk}R^{j}_{c}P_{k}.$ (44)

Since the chosen free data $\bar{\gamma}_{ij}$ in eq. (21) satisfies the
quasi-isotropic gauge condition [82] and $K$ in eq. (22) meets the asymptotic
maximal gauge condition [82], the ADM linear momentums and angular momentums
of the system are defined [82]

$\displaystyle P^{\infty}_{i}$
$\displaystyle=\frac{1}{8\pi}\lim_{S_{t}\to\infty}\oint_{S_{t}}(K_{jk}-K\gamma_{jk})(\vec{\partial}_{i})^{j}s^{k}dA,$
(45) $\displaystyle J^{\infty}_{i}$
$\displaystyle=\frac{1}{8\pi}\lim_{S_{t}\to\infty}\oint_{S_{t}}(K_{jk}-K\gamma_{jk})(\vec{\phi}_{i})^{j}s^{k}dA,$
(46)

here,

$\displaystyle\vec{\phi}_{x}$
$\displaystyle=-(z-z_{\mathrm{CM}})\vec{\partial}_{y}+(y-y_{\mathrm{CM}})\vec{\partial}_{z},$
$\displaystyle\vec{\phi}_{y}$
$\displaystyle=+(z-z_{\mathrm{CM}})\vec{\partial}_{x}+(x-x_{\mathrm{CM}})\vec{\partial}_{z},$
$\displaystyle\vec{\phi}_{z}$
$\displaystyle=-(y-y_{\mathrm{CM}})\vec{\partial}_{x}+(x-x_{\mathrm{CM}})\vec{\partial}_{y}.$
(47)

Lastly, to calculate the total ADM mass of the system, we use [82]

$M_{\mathrm{ADM}}=\int_{\Sigma_{t}}\left[\psi^{5}E+\frac{1}{16\pi}\left(\bar{A}_{ij}\bar{A}^{ij}\psi^{-7}-{\bar{R}}\psi-\frac{2}{3}K^{2}\psi^{5}\right)\right]\sqrt{\bar{\gamma}}d^{3}x.$
(48)

###  3.4 Tabulated equations of state

We can import tabulated EOS from the CompOSE repository [70, 71, 72]. A
general CompOSE table is first restricted to $T=0$. The baryon density $n_{b}$
is converted to the rest mass density $\rho_{0}$ via the neutron mass $m_{n}$:
$\rho_{0}=m_{n}n_{b}$.

We numerically differentiate and interpolate the EOS table to generate a
Hermite spline representation of the functions $\rho_{0}\left(h\right)$,
$P\left(h\right)$, and $\epsilon\left(h\right)$ (see sec. (A)). In practice,
the logarithms of these quantities are actually used to generate the
interpolants.

A number of slight modifications may be made to the EOS in pre-processing in
order to make it more amenable to interpolation. These modifications are
necessary for both the physical consistency of the EOS and to facilitate the
convergence of the solution.

The first modification involves adjusting specific enthalpy at the NS surface
for the EOS tables in which the specific enthalpy falls below $1$. To this
aim, we scale all $\rho_{0}$ points by the multiplicative constant
$\eta=h_{\mathrm{min}}$ corresponding to the lowest value of $h$ in the table.
Accordingly, $h$ points are scaled by the value $\eta^{-1}$, i.e:

$\displaystyle\rho_{0}\to\eta\rho_{0}$ $\displaystyle
h=\frac{e+P}{\rho_{0}}\to\frac{e+P}{\eta\rho_{0}}=\eta^{-1}h.$ (49)

This alteration can be carried through to dynamical evolution codes by e.g.
using the same EOS table or by scaling the rest-mass density (or baryon number
density) by the same factor.

In addition, several other features of tabular EOS may pose problems for both
interpolation and the use of specific enthalpy as the independent
thermodynamic variable. Among these problems is the presence of a region near
the surface of the NS where both $\frac{de}{dh}$ and $\frac{d\rho_{0}}{dh}$
diverge. This region is common to many EOS with a ‘crust’, including both
tabular EOS such as SFHo [84] and piecewise polytropics [85]. While this
region is unavoidable if we wish to accurately represent the equation of state
via specific enthalpy, cf. [86, 87], it also limits the accuracy of the
solution within the NS (see fig. 3).

Another consideration is the spacing of the data points in the CompOSE table,
which may be highly irregular especially for tables generated as a piecewise
combination of different models. To re-grid the EOS, we numerically
differentiate the table points using Fornberg’s method, explained in sec. (B),
and generate a low-degree Hermite interpolant, sec. (A), which is sampled to
produce new data points on an evenly-spaced grid. This method decreases
oscillations when the new data points are themselves interpolated.
Additionally, the interpolation of the EOS truncates the jumps in derivatives
of the thermodynamic variables, i.e., $\frac{de}{dh}$ and
$\frac{d\rho_{0}}{dh}$, that as we mentioned before, in some tables are not
well-defined. This low-degree Hermite interpolant corresponds to approximating
the derivatives as finite at these points (since infinite values would be
unphysical anyway).

Figure 2: Comparison of changes with respect to the specific enthalpy ($h$) in
total energy density ($e$) and the rest mass density ($\rho_{0}$) near the NS
surface for three EOS: K96 (single polytrope), SLy (piecewise polytrope), and
SFHo (table). A lower specific enthalpy corresponds to a point closer to the
NS surface, with $h=1$ defining the surface itself.

Fig. 2 shows the resultant tabulated SFHo EOS after re-griding using a low-
degree Hermite interpolant with comparison to EOS of K96 and SLy. The K96 EOS
is a single polytrope, corresponding to $K=96.7$ and $\Gamma=2$ [88]. SLy is
an approximate piecewise polytrope of tabulated SLy4, and therefore quantities
such as $\frac{de}{dh}$ are not continuous in certain regions corresponding to
transitions between piecewise segments [89]. SFHo, here restricted to the
zero-temperature regime, shows similar discontinuities [84].

Examining the functions relevant to the ID generation, e.g. $\rho_{0}$ and
$e$, we see that the single-polytrope model is considerably simpler and more
amenable to numerical solution than the models with discontinuities. The
difference is typically relevant to the surface of the NS (where $h=1$), as
shown in the nearly vertical profile of $\frac{d\rho_{0}}{dh}$ and
$\frac{de}{dh}$ in fig. (2). The sharp changes in $\frac{d\rho_{0}}{dh}$ and
$\frac{de}{dh}$ with the piecewise polytropic and tabular EOS are not
perfectly captured by the grid, and thus they limit the convergence of the
constraint violations (see sec. 4.1).

###  3.5 Initial data construction

Construction of ID often necessitates an iterative approach. This process
involves progressively refining the ID until a suitable solution is achieved.
However, rapid or uncontrolled updates of the fields from one step to the
next, or lack of adjustments of NS masses or their centers can lead to a
divergent answer and code crashes. Another challenge is to find the NS surface
after each update, as the true surface is unknown a-priori, and then creating
a new computational grid with patch’s surfaces adapted to the new NS surface.

We deploy the following iterative procedure to construct ID of a BNS system.
In particular, we start from a coarse resolution grid and progressively refine
the answer, while controlling diagnostics, until the Hamiltonian and momentum
constraints, eqs. (55) and (56), reach a plateau. Then, we increase the
resolution and repeat this iteration until we achieve the level of desired
accuracy.

Step 0: We superimpose two Tolman-Oppenheimer-Volko (TOV) star solutions as
the initial guess of the fields $\\{\psi,\alpha\psi,B^{i}\\}$. To initialize
the $\phi$ fields for each NS we use the approximation
$\phi=-\Omega^{z}_{\mathrm{BNS}}(y_{\mathrm{NS}}-y_{\mathrm{CM}})x$.

Step 1: We solve the elliptic equations iteratively in this specific order:
first eq. (25) for the matter field, and then eqs. (17), (18), and (19) for
the metric fields. This approach has been found to enhance the solution’s
stability, particularly when solving the matter field first. The order of
solving the metric fields themselves, however, appears to be less critical.
Moreover, during each iteration step, we focus on solving a single elliptic
equation while keeping the other fields fixed. In essence, the fixed fields
act as source terms influencing the equation being solved. Finally, within the
Newton-Raphson iterative method, we perform only one update step per field and
then incorporating the newly solved value into the source terms for the next
equation. In the following $\Xi$ denotes a field from the set
$\\{\phi,\psi,\alpha\psi,B^{i}\\}$.

Step 2: As mentioned earlier, iterative solvers are often sensitive to prompt
changes; code crashes can happen if a sudden update take place in the system.
As such, we update the field solutions that are obtained from Step 1 in a
relaxed fashion. To this aim, we use
$\Xi=\lambda\Xi_{\mathrm{new}}+(1-\lambda)\Xi_{\mathrm{old}}$, in which
$\Xi_{\mathrm{old}}$ denotes the solution before entering Step 1 and
$\Xi_{\mathrm{new}}$ is the solution after exiting that step. $\lambda$
denotes the weight of update. We generally use $\lambda=0.2$. This choice of
$\lambda$ value is proven to work for all experiments we have done.

Step 3: After updating all $\Xi$ fields, we see the baryon mass of each NS
deviates, often by a few percents, from the target value. Additionally, for
spinning NS, since at Step 0 we used TOV solution and then we added the spin
vector to the NS, it is not surprising to observe the baryon mass is different
from the prescribed target value. Moreover, the starting resolution is often
coarse and again the baryon mass may change as we go to a higher resolution.
If we do not account for these changes, the baryon mass deviates even further
at later iterations and may lead to a code crash. To adjust the baryon mass,
we note that $\rho_{0}=\rho_{0}(h)$, and as shown in eq. (31) the specific
enthalpy depends on a constant $C$, i.e., $h=h(C)$, which implies
$\rho_{0}=\rho_{0}(C)$. Therefore by using a root finder, we find the value of
$C$ in eq. (40) such that the baryon mass remains the same as the target
value.

Step 4: Since at Step 0 we began with a rough approximation of the solution,
the linear ADM momentums are initially not zero. Furthermore, for asymmetric
masses or spinning systems, we do not know in advance where is the exact
position of the system center of mass – as one needs to consider the full
general relativity effect to find it. Therefore, we iteratively find the
system center of mass $\vec{r}_{\mathrm{CM}}$, by demanding the ADM momentum
in each direction to be zero. The ADM momentum in $z$-direction proves
consistently to be small, generally,
$\frac{|P^{\infty}_{z}|}{M_{\mathrm{ADM}}}<10^{-11}$, hence we keep
$z_{\mathrm{CM}}=0$, as the original value. However, often the initial value
of ADM momenta in $x$ or $y$ direction is $\approx 10^{-5}M_{\mathrm{ADM}}$.
Therefore, we adjust $x_{\mathrm{CM}}$ and $y_{\mathrm{CM}}$ as follows

$\displaystyle x_{\mathrm{CM,new}}$
$\displaystyle=x_{\mathrm{CM,old}}+\lambda\frac{P^{\infty}_{y}}{\Omega^{z}_{\mathrm{BNS}}M_{\mathrm{ADM}}},$
$\displaystyle y_{\mathrm{CM,new}}$
$\displaystyle=y_{\mathrm{CM,old}}-\lambda\frac{P^{\infty}_{x}}{\Omega^{z}_{\mathrm{BNS}}M_{\mathrm{ADM}}}.$
(50)

Here, $\lambda$ is generally chosen $0.2$. It is worth noting that $K_{ij}$ is
a function of $\beta^{i}$, while $\beta^{i}$ through eq. (2) is a function of
$r^{i}_{\mathrm{CM}}$. Therefore, by adjusting the $r^{i}_{\mathrm{CM}}$, we
influence $K_{ij}$ – and hence all other fields as they are coupled to another
– and we can derive $P^{\infty}_{i}$ to zero iteratively.

Step 5: We update the stress energy tensor, in particular, we update the
specific enthalpy for each NS in a relaxed fashion as $h=\lambda
h_{\mathrm{new}}+(1-\lambda)h_{\mathrm{old}}$; $\lambda$ is usually set to
$0.5$. Then, we use the new value of specific enthalpy to update
$\rho_{0}(h)$, $\epsilon(h)$, and $P(h)$.

Step 6: Finding the orbital angular velocity, $\Omega^{z}_{\mathrm{BNS}}$, for
a quasi-circular orbit demands full solution of general relativity. Since we
stat with a Newtonian approximation for this value we need to refine it.
Following [90], so called the force balance method, we first compute
$\partial_{i}\ln\Gamma$ at each NS center, where $\Gamma$ computed as

$\Gamma=\frac{\alpha
u^{0}[1-(\beta^{i}+\xi^{i}+\frac{w^{i}}{hu^{0}})\frac{D_{i}\phi}{\alpha^{2}hu^{0}}-\frac{w_{i}w^{i}}{(\alpha^{2}hu^{0})^{2}}]}{\sqrt{1-(\beta^{i}+\xi^{i}+\frac{w^{i}}{hu^{0}})(\beta_{i}+\xi_{i}+\frac{w_{i}}{hu^{0}})\frac{1}{\alpha^{2}}}}.$
(51)

Then, we use a root finder to find $\Omega^{z}_{\mathrm{BNS}}$ such that

$\partial_{i}\ln[\alpha^{2}-(\beta^{i}+\xi^{i}+\frac{w^{i}}{hu^{0}})(\beta_{i}+\xi_{i}+\frac{w_{i}}{hu^{0}})]+2\partial_{i}\ln\Gamma=0,$
(52)

where, $\partial_{i}=\frac{\partial}{\partial x^{i}}$. Since in our setup the
NS centers sit on y-axis, we compute eq. (52) along this axis. We find
$\Omega^{z}_{\mathrm{BNS}}$ for each NS centers and update it accordingly.

Step 7: At this point the specific enthalpy profile and hence NS surface are
changed. We need to find the new location of NS surface for adjusting cubed
spherical patches that are covering the NS – so we can separate matter and
vacuum into different patches. Since some parts of the NS may need to extend
to the patches that are currently covering vacuum, we extrapolate $h$ into
these patches so the root finder can find where $h=1$. To this aim we use the
following formula to extrapolate specific enthalpy to the vacuum

$f\left(r\right)=\left(a+\frac{b}{r}\right)\exp\left(-c_{0}\frac{r}{r_{0}}\right),$
(53)

here, $r$ denotes the coordinate distance from the NS center, and $c_{0}$ is a
constant (typically $0.01$). We find the values of $a$, and $b$ by demanding
$C^{1}$ continuity across the NS surface. Additionally, we use the same eq. 53
to extrapolate $\phi$ field outside the NS. This step is required when the NS
surface is expanded and we want to interpolate $\phi$ from the current grid to
a new grid.

Step 8: We identify the center of each NS by locating the coordinate where the
specific enthalpy reaches its maximum value. After multiple updates to the
matter fields in previous steps, the NS centers can exhibit slight drifts from
their initial positions. These drifts can accumulate over time, potentially
leading to code crashes. To address this issue, we employ a corrective measure
that adjusts the specific enthalpy function as follows

$h_{\mathrm{new}}(\vec{r})=h_{\mathrm{old}}(\vec{r})-(\vec{r}-\vec{r}_{0})\cdot\vec{\nabla}h_{\mathrm{old}}(\vec{r_{0}}),$
(54)

where, $r_{0}$ denotes the coordinate of NS center. This adjustment ensures
that the maximum value of the specific enthalpy remains at the same location,
effectively preventing the NS centers from drifting significantly.

Step 9: We find the profile of $\sigma_{\mathrm{out}}(X,Y)$ that is necessary
for eq. (33) to have patch surfaces fitting the NS surface. To this end, for a
given angular $\theta$ and $\phi$ in spherical coordinate, we use a root
finder for $r$ to solve $h(r,\theta,\phi)-1=0$ and hence finding the new NS
surface. Here, $r$ is the coordinate distance from the NS center.

Step 10: If the NS surfaces are changed or if the resolution increases at the
next iteration, we create a new grid. To ensure a smooth transition, we use
spectral interpolation techniques to transfer data from the previous grid onto
the new one.

Step 11: We monitor the Hamiltonian and momentum constraints, eqs (55) and
(56). We restart from Step 1 unless the constraints reach their truncation
error and are level off. In this case, when the constraints reach a plateau,
we stop the iterative process if there is no higher resolution demanded;
otherwise we go to the next resolution and start from Step 1.

##  4 Results

###  4.1 Convergence test

Since Elliptica is a pseudo-spectral code, the first expectation of the code
is spectral convergence feature. As such, we calculate the Hamiltonian and
momentum constraints, respectively, using

$\displaystyle H$ $\displaystyle:=R-K_{ij}K^{ij}+K^{2}-16\pi E$ (55)
$\displaystyle M^{i}$ $\displaystyle:=D_{j}(K^{ij}-\gamma^{ij}K)-8\pi j^{i}$
(56)

For the convergence test, we generate ID for two symmetric BNS systems: one
with polytropic K96 EOS and the other with tabulated SFHo EOS. In these
systems, the NS have baryon mass $1.4$ with no spin, and their separation is
$50$ units.

(a) K96, NS interior regions

(b) SFHo, NS interior regions

(c) K96, NS exterior regions

(d) SFHo, NS exterior regions

(e) K96, all regions

(f) SFHo, all regions

Figure 3: The $L_{2}$-norm of the Hamiltonian and momentum constraint
violations for the BNS system with the K96 EOS (single polytrope) and the SFHo
EOS (table). The constraint violations at each point are summed over all
points in the specified regions. The regions not including a tabular EOS show
exponential convergence, while regions with such an EOS reach a limit at high
resolutions.

Fig. (3) shows the $L_{2}$-norm of the Hamiltonian and momentum constraints as
a function of grid resolution (focusing on the final iteration at each
resolution) for the two systems. In figs. (3(a),3(b)), the convergence test
emphasizes on the NS. We observe the constraint violations decrease
exponentially for smooth matter field, i.e., K96 EOS. For the SFHo EOS,
initially for low resolutions the constraints converge exponentially but for
higher resolutions the constraints level off. This behavior is expected; when
the grid resolution is coarse the true discontinuities of EOS variables are
not seen by the spectral expansion. For fine resolutions all features of EOS
variables emerge and while spectral convergence try resolve these features but
it is not successful. Hence we see constraints are soon level off and do not
decrease as we increase the resolution. We note that, for the tabular EOS
(SFHo), the convergence is limited due to the sharp features of
$\frac{d\rho_{0}}{dh}$ and $\frac{de}{dh}$ near the surface of the NS (see
fig. (2)). Indeed, the overall magnitude of the constraint violations is
higher and the rate of convergence is not exponential with the tabular EOS (in
contrast to the simple polytrope). Nonetheless, the profile of the convergence
tests still suggest convergence up to a limit at high resolutions.

Ensuring spectral convergence of the code for regions that fields are smooth,
we calculate $L_{2}$-norm of constraints at immediate neighboring regions of
NS, where there is no matter fields. Figs. (3(c), 3(d)) demonstrates they
convergence spectrally.

Finally, the overall converge of the constraints, calculated over all regions,
are shown in figs. (3(e), 3(f)); while both systems exhibit convergence as the
resolution increases, the rate of convergence for BNS with SFHo EOS is smaller
than K96 EOS.

###  4.2 Post-Newtonian test

We calculate the binding energy $E_{b}$ of a symmetric BNS system with NS
baryon mass $1.4$ and tabulated SLy4 EOS for varied separations of NS. Here
$E_{b}=M_{\mathrm{ADM}}-M_{\infty}$ and
$M_{\infty}=M^{1}_{\mathrm{TOV}}+M^{2}_{\mathrm{TOV}}$; we use eq. 48 to
calculate $M_{\mathrm{ADM}}$ and $M^{1/2}_{\mathrm{TOV}}$ are the
corresponding gravitational mass of NS in isolation – found by a TOV solver.

To determine $\Omega^{z}_{\mathrm{BNS}}$ for a specific separation, we employ
an iterative approach based on the force balance method described in equation
eq. 52. To validate the accuracy of the generated ID against the expected
analytical values for significant separations, particularly when the BNS
system exhibits quasi-circular motion, we calculate the binding energy $E_{b}$
corresponding to the given $\Omega^{z}_{\mathrm{BNS}}$ utilizing the post-
Newtonian formula outlined in [91].

Figure 4: Analytic post-Newtonian curve versus NR curve. The values of $E_{b}$
for given $\Omega^{z}_{\mathrm{BNS}}$ are compared for post-Newtonian method
and NR method. For binaries with a large separation post-Newtonian and NR
values are matching. For BNS of a close separation there is a slight deviation
from post-Newtonian prediction as the system is not fully in quasi-circular
status. Here, $\mu={M^{1}_{\mathrm{TOV}}M^{2}_{\mathrm{TOV}}}/{M_{\infty}}$
and the binary has a symmetric baryon mass of $1.4$ for each NS and uses
tabulate SLy4 EOS.

Figure (4) illustrates the binding energy of the BNS system across different
separations compared to the post-Newtonian data points. In cases of binaries
with significant separations, meeting the quasi-circular motion criteria, we
observe that the post-Newtonian and NR values align. However, for smaller
separations, the NR values of $E_{b}$ are marginally higher.

###  4.3 Spin

We generate ID for different $\Omega^{z}_{\mathrm{NS}}$, in eq. (24), values
pertaining to one of the NS in the BNS system, where each NS has a baryon mass
of $1.4$ and they are separated by $30$ units. In this scenario, we utilize
the K96 EOS. The maximum dimensionless spin achievable is $\approx 0.56$,
which corresponds to the mass shedding limit of a single NS as discussed in
[92]. By increasing the $\Omega^{z}_{\mathrm{NS}}$ beyond $0.02$, the spinning
NS becomes too oblate that the NS surface finding routine fails (Step 9 in
section 3.5).

Figure 5: The relation between $\Omega^{z}_{\mathrm{NS}}$ and $\chi_{z}$ is
depicted for a symmetric mass BNS system with a baryonic mass of $1.4$ and a
separation of $30$. Initially, $\chi_{z}$ exhibits a linear growth pattern
with $\Omega^{z}_{\mathrm{NS}}$, but as $\Omega^{z}_{\mathrm{NS}}$ is raised
beyond a certain point, $\chi_{z}$ starts to increase non-linearly.

Fig. (5) shows the relation between $\Omega^{z}_{\mathrm{NS}}$ in eq. 24 and
$\chi_{z}$ in eq. 43. There is a linear connection between
$\Omega^{z}_{\mathrm{NS}}$ and $\chi_{z}$ for low values of
$\Omega^{z}_{\mathrm{NS}}$. However, as $\Omega^{z}_{\mathrm{NS}}$ increases,
the relationship transitions into a nonlinear pattern, causing $\chi_{z}$ to
rise more steeply.

Moreover, fig. (5) can be utilized as an approximate reference to determine
the appropriate value of $\Omega^{z}_{\mathrm{NS}}$ corresponding to a desired
spin level.

##  5 Summary

In this work we have presented a significant upgrade to the Elliptica
infrastructure. Previously limited to black hole-neutron star system and
(piecewise) polytropic EOS for NS, Elliptica can now construct ID for spinning
BNS systems, incorporating realistic tabulated EOS for the NS matter.

To incorporate tabular EOS, we have developed a number of techniques centered
around interpolation that allow convergence in the resulting solution while
remaining as true as possible to the underlying microphysics. While the
convergence of the solution is negatively affected by the complexity of the
EOS, we maintain convergence up to a limit imposed by the surface features.

We have validated our code through convergence tests and comparisons with
established analytical results, particularly in the post-Newtonian regime.
These tests demonstrate the code’s accuracy and robustness.

For smooth EOS, such as polytropic models, we achieve spectral convergence,
indicating an optimal error reduction rate with increasing resolution.
However, for EOS tables with discontinuities, spectral convergence is
diminished. At high resolutions, we observe a convergence plateau inside NS,
though global convergence is still ensured. This behavior is expected for
discontinuous functions.

A particularly interesting area for future development lies in constructing
self-consistent ID for BNS systems that contain magnetic fields. This
capability would be especially valuable for studying systems containing
pulsars or magnetars, where these fields play a significant role in the
dynamics.

AR and AN gratefully acknowledge David Radice for valuable discussions and
insightful feedback on the manuscript. AR acknowledges support from NASA under
award No. 80NSSC21K1720. The numerical simulations were performed on Roar
Collab High Performance Computing Cluster at The Pennsylvania State
University. Elliptica is now public 111https://github.com/rashti-
alireza/Elliptica and open-source under GPL-3.0 license.

##  Appendix A Interpolation of the equation of state

To generate a suitable interpolant of a tabular EOS, we first approximate
$\frac{\partial p}{\partial h}$, $\frac{\partial\epsilon}{\partial h}$, and
$\frac{\partial\rho_{0}}{\partial h}$ by finite-difference methods on the
CompOSE data using Fornberg’s algorithm [93], explained in B, for unevenly-
spaced grids. Having the derivatives, we then generate a spline interpolant
composed of Hermite polynomials. The procedure for e.g. the pressure
$p\left(h\right)$ is:

Step 0: We are given data points $(p_{j},h_{j})$ monotonically increasing in
$h$.

Step 1: For each point $j$, we take the set of $N$ points
$X_{j}=\left\\{h_{j-\frac{N}{2}},...,h_{j},...,h_{j+\frac{N}{2}}\right\\}$
centered around $h_{j}$. We shift the indices when needed for points near the
boundaries $h_{0}$ and $h_{j_{\mathrm{max}}}$.

Step 2: At each point, approximate
$\left.\frac{dp}{dh}\right|_{h_{j}}\approx\sum_{\nu=j-\frac{N}{2}}^{j+\frac{N}{2}}\delta_{N,\nu}^{1}p_{\nu}$,
with the coefficients $\delta$ calculated by Fornberg’s algorithm applied over
the set $X_{j}$ obtained in Step 1.

Step 3: Having $p_{j}$ and its derivative $\frac{dp(h_{j})}{dh}$, we generate
the interpolating Hermite polynomial (according to [94]) of desired order.

The same procedure is applied to obtain splines for $e\left(h\right)$,
$\rho_{0}\left(h\right)$, and $\epsilon\left(h\right)$. Derivatives of these
functions are thereafter approximated by analytical derivatives of the spline
interpolant.

##  Appendix B Application of Fornberg’s finite difference method

Step 2 of sec. A entails finding derivatives such as $\frac{dp}{dh}$, which
EOS tables do not provide. we evaluate these derivatives numerically using
finite difference methods on the data points $\left(p_{j},h_{j}\right)$. Since
the data points are typically unevenly spaced, we use Fornberg’s finite
difference algorithm to generate finite difference coefficients that
approximate $\frac{dp(h_{j})}{dh}$ at each point [93]. Specifically, the
algorithm calculates the weights $\delta_{N,\nu}^{m}$ such that

$\left.\frac{d^{m}p}{{dh}^{m}}\right|_{h_{j}}\approx\sum_{\nu=j-\frac{N}{2}}^{j+\frac{N}{2}}\delta_{N,\nu}^{m}p_{\nu},$
(57)

where $N$ is the number of points used (and determines the order of the finite
difference approximation).

While Fornberg’s algorithm can provide derivatives of arbitrarily high orders
(given enough data points), we are only interested in the first derivative. In
addition, we do not need the finite difference coefficients for all orders, so
we simplify the algorithm slightly. Then given a subset $X_{j}$ centered
around $h_{j}$, from Step 1 of sec. A, we find the finite difference weights
using algorithm 2 (where $c_{1}$, $c_{2}$, and $c_{3}$ are introduced just to
simplify the notation, and $x_{k}$ is the $k$th element of $X_{j}$).

Algorithm 2 Simplification of Fornberg’s finite difference algorithm, adapted
from [93].

1:$\delta^{0}_{0,0}=c_{1}=1$

2:for $n=1,...,N$ do

3: $c_{2}=1$

4: for $\nu=0,...,n-1$ do

5: $c_{3}=x_{n}-x_{\nu}$

6: $c_{2}=c_{2}\cdot c_{3}$

7: if $n=0,1$ then

8: $\delta^{n-1}_{n-1,\nu}=0$

9: end if

10:
$\delta^{0}_{n,\nu}=\frac{\left(x_{n}-h_{j}\right)}{c_{3}}\delta^{0}_{n-1,\nu}$

11:
$\delta^{1}_{n,\nu}=\frac{1}{c_{3}}\left(\left(x_{n}-h_{j}\right)\delta^{1}_{n-1,\nu}-\delta^{0}_{n-1,\nu}\right)$

12: end for

13:
$\delta^{0}_{n,n}=\frac{c_{1}}{c_{2}}\left(h_{j}-x_{n-1}\right)\delta^{0}_{n-1,n-1}$

14:
$\delta^{1}_{n,n}=\frac{c_{1}}{c_{2}}\left(\delta^{0}_{n-1,n-1}+\left(h_{j}-x_{n-1}\right)\delta^{1}_{n-1,n-1}\right)$

15: $c_{1}=c_{2}$

16:end for

References

## References

  * [1] B. P. Abbott et al.  GW190425: Observation of a Compact Binary Coalescence with Total Mass $\sim 3.4M_{\odot}$.  Astrophys. J. Lett., 892(1):L3, 2020.
  * [2] Swami Vivekanandji Chaurasia, Tim Dietrich, Nathan K. Johnson-McDaniel, Maximiliano Ujevic, Wolfgang Tichy, and Bernd Brügmann.  Gravitational waves and mass ejecta from binary neutron star mergers: Effect of large eccentricities.  Phys. Rev. D, 98(10):104005, 2018.
  * [3] Swami Vivekanandji Chaurasia, Tim Dietrich, Maximiliano Ujevic, Kai Hendriks, Reetika Dudi, Francesco Maria Fabbri, Wolfgang Tichy, and Bernd Brügmann.  Gravitational waves and mass ejecta from binary neutron star mergers: Effect of the spin orientation.  Phys. Rev. D, 102(2):024087, 2020.
  * [4] Rodrigo Fernández, Alexander Tchekhovskoy, Eliot Quataert, Francois Foucart, and Daniel Kasen.  Long-term GRMHD simulations of neutron star merger accretion discs: implications for electromagnetic counterparts.  Mon. Not. Roy. Astron. Soc., 482(3):3373–3393, 2019.
  * [5] Lunan Sun, Milton Ruiz, Stuart L. Shapiro, and Antonios Tsokaros.  Jet launching from binary neutron star mergers: Incorporating neutrino transport and magnetic fields.  Phys. Rev. D, 105(10):104028, 2022.
  * [6] Daniel M. Siegel.  r-Process nucleosynthesis in gravitational-wave and other explosive astrophysical events.  Nature Reviews Physics, 4(5):306–318, April 2022.
  * [7] Rodrigo Fernández and Brian D. Metzger.  Electromagnetic Signatures of Neutron Star Mergers in the Advanced LIGO Era.  Ann. Rev. Nucl. Part. Sci., 66:23–45, 2016.
  * [8] J. Aasi et al.  Advanced LIGO.  Class. Quant. Grav., 32:074001, 2015.
  * [9] F. Acernese et al.  Advanced Virgo: a second-generation interferometric gravitational wave detector.  Class. Quant. Grav., 32(2):024001, 2015.
  * [10] T Akutsu, M Ando, K Arai, Y Arai, S Araki, A Araya, N Aritomi, H Asada, Y Aso, S Bae, Y Bae, L Baiotti, R Bajpai, M A Barton, K Cannon, Z Cao, E Capocasa, M Chan, C Chen, K Chen, Y Chen, C Y Chiang, H Chu, Y K Chu, S Eguchi, Y Enomoto, R Flaminio, Y Fujii, Y Fujikawa, M Fukunaga, M Fukushima, D Gao, G Ge, S Ha, A Hagiwara, S Haino, W B Han, K Hasegawa, K Hattori, H Hayakawa, K Hayama, Y Himemoto, Y Hiranuma, N Hirata, E Hirose, Z Hong, B Hsieh, G Z Huang, H Y Huang, P Huang, Y C Huang, Y Huang, D C Y Hui, S Ide, B Ikenoue, S Imam, K Inayoshi, Y Inoue, K Ioka, K Ito, Y Itoh, K Izumi, C Jeon, H B Jin, K Jung, P Jung, K Kaihotsu, T Kajita, M Kakizaki, M Kamiizumi, N Kanda, G Kang, K Kawaguchi, N Kawai, T Kawasaki, C Kim, J Kim, J C Kim, W S Kim, Y M Kim, N Kimura, N Kita, H Kitazawa, Y Kojima, K Kokeyama, K Komori, A K H Kong, K Kotake, C Kozakai, R Kozu, R Kumar, J Kume, C Kuo, H S Kuo, Y Kuromiya, S Kuroyanagi, K Kusayanagi, K Kwak, H K Lee, H W Lee, R Lee, M Leonardi, K L Li, L C C Lin, C Y Lin, F K Lin, F L Lin, H L Lin, G C Liu, L W Luo, E Majorana, M Marchio, Y Michimura, N Mio, O Miyakawa, A Miyamoto, Y Miyazaki, K Miyo, S Miyoki, Y Mori, S Morisaki, Y Moriwaki, K Nagano, S Nagano, K Nakamura, H Nakano, M Nakano, R Nakashima, Y Nakayama, T Narikawa, L Naticchioni, R Negishi, L Nguyen Quynh, W T Ni, A Nishizawa, S Nozaki, Y Obuchi, W Ogaki, J J Oh, K Oh, S H Oh, M Ohashi, N Ohishi, M Ohkawa, H Ohta, Y Okutani, K Okutomi, K Oohara, C Ooi, S Oshino, S Otabe, K Pan, H Pang, A Parisi, J Park, F E Peña Arellano, I Pinto, N Sago, S Saito, Y Saito, K Sakai, Y Sakai, Y Sakuno, S Sato, T Sato, T Sawada, T Sekiguchi, Y Sekiguchi, L Shao, S Shibagaki, R Shimizu, T Shimoda, K Shimode, H Shinkai, T Shishido, A Shoda, K Somiya, E J Son, H Sotani, R Sugimoto, J Suresh, T Suzuki, T Suzuki, H Tagoshi, H Takahashi, R Takahashi, A Takamori, S Takano, H Takeda, M Takeda, H Tanaka, K Tanaka, K Tanaka, T Tanaka, T Tanaka, S Tanioka, E N Tapia San Martin, S Telada, T Tomaru, Y Tomigami, T Tomura, F Travasso, L Trozzo, T Tsang, J S Tsao, K Tsubono, S Tsuchida, T Tsutsui, T Tsuzuki, D Tuyenbayev, N Uchikata, T Uchiyama, A Ueda, T Uehara, K Ueno, G Ueshima, F Uraguchi, T Ushiba, M H P M van Putten, H Vocca, J Wang, T Washimi, C Wu, H Wu, S Wu, W R Xu, T Yamada, K Yamamoto, K Yamamoto, T Yamamoto, K Yamashita, R Yamazaki, Y Yang, K Yokogawa, J Yokoyama, T Yokozawa, T Yoshioka, H Yuzurihara, S Zeidler, M Zhan, H Zhang, Y Zhao, and Z H Zhu.  Overview of KAGRA: Calibration, detector characterization, physical environmental monitors, and the geophysics interferometer.  Progress of Theoretical and Experimental Physics, 2021(5), 02 2021\.  05A102.
  * [11] David Reitze et al.  Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO.  Bull. Am. Astron. Soc., 51(7):035, 2019.
  * [12] Seiji Kawamura et al.  Current status of space gravitational wave antenna DECIGO and B-DECIGO.  PTEP, 2021(5):05A105, 2021.
  * [13] M. Punturo et al.  The Einstein Telescope: A third-generation gravitational wave observatory.  Class. Quant. Grav., 27:194002, 2010.
  * [14] Rana X. Adhikari et al.  Astrophysical science metrics for next-generation gravitational-wave detectors.  Class. Quant. Grav., 36(24):245010, 2019.
  * [15] Pau Amaro-Seoane et al.  Laser Interferometer Space Antenna.  2 2017.
  * [16] K. Ackley et al.  Neutron Star Extreme Matter Observatory: A kilohertz-band gravitational-wave detector in the global network.  Publ. Astron. Soc. Austral., 37:e047, 2020.
  * [17] Jun Luo et al.  TianQin: a space-borne gravitational wave detector.  Class. Quant. Grav., 33(3):035010, 2016.
  * [18] B. P. Abbott et al.  GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral.  Phys. Rev. Lett., 119(16):161101, 2017.
  * [19] B. P. Abbott et al.  Gravitational Waves and Gamma-rays from a Binary Neutron Star Merger: GW170817 and GRB 170817A.  Astrophys. J. Lett., 848(2):L13, 2017.
  * [20] D. A. Coulter et al.  Swope Supernova Survey 2017a (SSS17a), the Optical Counterpart to a Gravitational Wave Source.  Science, 358:1556, 2017.
  * [21] Francois Foucart, Pablo Laguna, Geoffrey Lovelace, David Radice, and Helvi Witek.  Snowmass2021 Cosmic Frontier White Paper: Numerical relativity for next-generation gravitational-wave probes of fundamental physics.  3 2022.
  * [22] Koji Uryu and Antonios Tsokaros.  A new code for equilibriums and quasiequilibrium initial data of compact objects.  Phys. Rev. D, 85:064014, 2012.
  * [23] Antonios Tsokaros, Kōji Uryū, and Luciano Rezzolla.  New code for quasiequilibrium initial data of binary neutron stars: Corotating, irrotational, and slowly spinning systems.  Phys. Rev. D, 91(10):104030, 2015.
  * [24] Lambros Boukas, Antonios Tsokaros, and Koji Uryu.  The Parallel Compact Object CALculator: An Efficient General Relativistic Initial Data Solver for Compact Objects.  12 2023.
  * [25] Alireza Rashti, Francesco Maria Fabbri, Bernd Brügmann, Swami Vivekanandji Chaurasia, Tim Dietrich, Maximiliano Ujevic, and Wolfgang Tichy.  New pseudospectral code for the construction of initial data.  Phys. Rev. D, 105(10):104027, 2022.
  * [26] L. Jens Papenfort, Samuel D. Tootle, Philippe Grandclément, Elias R. Most, and Luciano Rezzolla.  New public code for initial data of unequal-mass, spinning compact-object binaries.  Phys. Rev. D, 104(2):024057, 2021.
  * [27] Lorene.  Langage Objet pour la RElativité NumériquE.
  * [28] Philippe Grandclement.  Accurate and realistic initial data for black hole-neutron star binaries.  Phys. Rev. D, 74:124002, 2006.  [Erratum: Phys.Rev.D 75, 129903 (2007)].
  * [29] Keisuke Taniguchi, Thomas W. Baumgarte, Joshua A. Faber, and Stuart L. Shapiro.  Quasiequilibrium sequences of black-hole-neutron-star binaries in general relativity.  Phys. Rev. D, 74:041502, 2006.
  * [30] Keisuke Taniguchi, Thomas W. Baumgarte, Joshua A. Faber, and Stuart L. Shapiro.  Quasiequilibrium black hole-neutron star binaries in general relativity.  Phys. Rev. D, 75:084005, 2007.
  * [31] Keisuke Taniguchi, Thomas W. Baumgarte, Joshua A. Faber, and Stuart L. Shapiro.  Relativistic black hole-neutron star binaries in quasiequilibrium: Effects of the black hole excision boundary condition.  Phys. Rev. D, 77:044003, 2008.
  * [32] Thiago Assumpcao, Leonardo R. Werneck, Terrence Pierre Jacques, and Zachariah B. Etienne.  Fast hyperbolic relaxation elliptic solver for numerical relativity: Conformally flat, binary puncture initial data.  Phys. Rev. D, 105(10):104037, 2022.
  * [33] Wolfgang Tichy.  A New numerical method to construct binary neutron star initial data.  Class. Quant. Grav., 26:175018, 2009.
  * [34] Wolfgang Tichy.  Constructing quasi-equilibrium initial data for binary neutron stars with arbitrary spins.  Phys. Rev. D, 86:064024, 2012.
  * [35] Wolfgang Tichy, Alireza Rashti, Tim Dietrich, Reetika Dudi, and Bernd Brügmann.  Constructing binary neutron star initial data with high spins, high compactnesses, and high mass ratios.  Phys. Rev. D, 100(12):124046, 2019.
  * [36] Nils L. Fischer and Harald P. Pfeiffer.  Unified discontinuous Galerkin scheme for a large class of elliptic equations.  Phys. Rev. D, 105(2):024034, 2022.
  * [37] Nils L. Vu et al.  A scalable elliptic solver with task-based parallelism for the SpECTRE numerical relativity code.  Phys. Rev. D, 105(8):084027, 2022.
  * [38] Harald P. Pfeiffer, Lawrence E. Kidder, Mark A. Scheel, and Saul A. Teukolsky.  A Multidomain spectral method for solving elliptic equations.  Comput. Phys. Commun., 152:253–273, 2003.
  * [39] Francois Foucart, Lawrence E. Kidder, Harald P. Pfeiffer, and Saul A. Teukolsky.  Initial data for black hole-neutron star binaries: A Flexible, high-accuracy spectral method.  Phys. Rev. D, 77:124051, 2008.
  * [40] Nick Tacik et al.  Binary Neutron Stars with Arbitrary Spins in Numerical Relativity.  Phys. Rev. D, 92(12):124012, 2015.  [Erratum: Phys.Rev.D 94, 049903 (2016)].
  * [41] Nick Tacik, Francois Foucart, Harald P. Pfeiffer, Curran Muhlberger, Lawrence E. Kidder, Mark A. Scheel, and B. Szilágyi.  Initial data for black hole–neutron star binaries, with rotating stars.  Class. Quant. Grav., 33(22):225012, 2016.
  * [42] Marcus Ansorg, Bernd Bruegmann, and Wolfgang Tichy.  A Single-domain spectral method for black hole puncture data.  Phys. Rev. D, 70:064011, 2004.
  * [43] Marcus Ansorg.  A Double-domain spectral method for black hole excision data.  Phys. Rev. D, 72:024018, 2005.
  * [44] Bernd Bruegmann, Jose A. Gonzalez, Mark Hannam, Sascha Husa, Ulrich Sperhake, and Wolfgang Tichy.  Calibration of Moving Puncture Simulations.  Phys. Rev. D, 77:024027, 2008.
  * [45] Marcus Thierfelder, Sebastiano Bernuzzi, and Bernd Bruegmann.  Numerical relativity simulations of binary neutron stars.  Phys. Rev. D, 84:044012, 2011.
  * [46] Tim Dietrich, Sebastiano Bernuzzi, Maximiliano Ujevic, and Bernd Brügmann.  Numerical relativity simulations of neutron star merger remnants using conservative mesh refinement.  Phys. Rev. D, 91(12):124041, 2015.
  * [47] Marcus Bugner, Tim Dietrich, Sebastiano Bernuzzi, Andreas Weyhausen, and Bernd Brügmann.  Solving 3D relativistic hydrodynamical problems with weighted essentially nonoscillatory discontinuous Galerkin methods.  Phys. Rev. D, 94(8):084004, 2016.
  * [48] David Hilditch, Andreas Weyhausen, and Bernd Brügmann.  Pseudospectral method for gravitational wave collapse.  Phys. Rev. D, 93(6):063006, 2016.
  * [49] Sarah Renkhoff, Daniela Cors, David Hilditch, and Bernd Brügmann.  Adaptive hp refinement for spectral elements in numerical relativity.  Phys. Rev. D, 107(10):104043, 2023.
  * [50] Milinda Fernando, David Neilsen, Hyun Lim, Eric Hirschmann, and Hari Sundar.  Massively Parallel Simulations of Binary Black Hole Intermediate-Mass-Ratio Inspirals.  7 2018.
  * [51] Frank Löffler, Joshua Faber, Eloisa Bentivegna, Tanja Bode, Peter Diener, Roland Haas, Ian Hinder, Bruno C Mundim, Christian D Ott, Erik Schnetter, Gabrielle Allen, Manuela Campanelli, and Pablo Laguna.  The einstein toolkit: a community computational infrastructure for relativistic astrophysics.  Classical and Quantum Gravity, 29(11):115001, 2012.
  * [52] Roland Haas, Steven R. Brandt, William E. Gabella, Miguel Gracia-Linares, Beyhan Karakaş, Rahime Matur, Miguel Alcubierre, Daniela Alic, Gabrielle Allen, Marcus Ansorg, Maria Babiuc-Hamilton, Luca Baiotti, Werner Benger, Eloisa Bentivegna, Sebastiano Bernuzzi, Tanja Bode, Bernd Bruegmann, Manuela Campanelli, Federico Cipolletta, Giovanni Corvino, Samuel Cupp, Roberto De Pietri, Peter Diener, Harry Dimmelmeier, Rion Dooley, Nils Dorband, Matthew Elley, Yaakoub El Khamra, Zachariah Etienne, Joshua Faber, Toni Font, Joachim Frieben, Bruno Giacomazzo, Tom Goodale, Carsten Gundlach, Ian Hawke, Scott Hawley, Ian Hinder, Sascha Husa, Sai Iyer, Thorsten Kellermann, Andrew Knapp, Michael Koppitz, Pablo Laguna, Gerd Lanferman, Frank Löffler, Joan Masso, Lars Menger, Andre Merzky, Jonah Maxwell Miller, Mark Miller, Philipp Moesta, Pedro Montero, Bruno Mundim, Andrea Nerozzi, Scott C. Noble, Christian Ott, Ravi Paruchuri, Denis Pollney, David Radice, Thomas Radke, Christian Reisswig, Luciano Rezzolla, David Rideout, Matei Ripeanu, Lorenzo Sala, Jascha A Schewtschenko, Erik Schnetter, Bernard Schutz, Ed Seidel, Eric Seidel, John Shalf, Ken Sible, Ulrich Sperhake, Nikolaos Stergioulas, Wai-Mo Suen, Bela Szilagyi, Ryoji Takahashi, Michael Thomas, Jonathan Thornburg, Malcolm Tobias, Aaryn Tonita, Paul Walker, Mew-Bing Wan, Barry Wardell, Helvi Witek, Miguel Zilhão, Burkhard Zink, and Yosef Zlochower.  The einstein toolkit, November 2020.  To find out more, visit http://einsteintoolkit.org.
  * [53] Sven Köppel.  Towards an exascale code for GRMHD on dynamical spacetimes.  J. Phys. Conf. Ser., 1031(1):012017, 2018.
  * [54] Boris Daszuta, Francesco Zappa, William Cook, David Radice, Sebastiano Bernuzzi, and Viktoriya Morozova.  GR-Athena++: Puncture Evolutions on Vertex-centered Oct-tree Adaptive Mesh Refinement.  Astrophys. J. Supp., 257(2):25, 2021.
  * [55] William Cook, Boris Daszuta, Jacob Fields, Peter Hammond, Simone Albanesi, Francesco Zappa, Sebastiano Bernuzzi, and David Radice.  GR-Athena++: General-relativistic magnetohydrodynamics simulations of neutron star spacetimes.  11 2023.
  * [56] Alireza Rashti, Maitraya Bhattacharyya, David Radice, Boris Daszuta, William Cook, and Sebastiano Bernuzzi.  Adaptive mesh refinement in binary black holes simulations.  Class. Quant. Grav., 41(9):095001, 2024.
  * [57] Swapnil Shankar, Philipp Mösta, Steven R. Brandt, Roland Haas, Erik Schnetter, and Yannick de Graaf.  GRaM-X: A new GPU-accelerated dynamical spacetime GRMHD code for Exascale computing with the Einstein Toolkit.  10 2022.
  * [58] http://doi.org/10.5281/zenodo.6131529.  CarpetX for the Einstein Toolkit.
  * [59] Katy Clough, Pau Figueras, Hal Finkel, Markus Kunesch, Eugene A. Lim, and Saran Tunyasuvunakool.  GRChombo : Numerical Relativity with Adaptive Mesh Refinement.  Class. Quant. Grav., 32(24):245011, 2015.
  * [60] Tomas Andrade et al.  GRChombo: An adaptable numerical relativity code for fundamental physics.  J. Open Source Softw., 6(68):3703, 2021.
  * [61] Wolfgang Tichy, Liwei Ji, Ananya Adhikari, Alireza Rashti, and Michal Pirog.  The new discontinuous Galerkin methods based numerical relativity program Nmesh.  Class. Quant. Grav., 40(2):025004, 2023.
  * [62] Ian Ruchlin, Zachariah B. Etienne, and Thomas W. Baumgarte.  SENR/NRPy+: Numerical Relativity in Singular Curvilinear Coordinate Systems.  Phys. Rev. D, 97(6):064036, 2018.
  * [63] Kenta Kiuchi, Kyohei Kawaguchi, Koutarou Kyutoku, Yuichiro Sekiguchi, and Masaru Shibata.  Sub-radian-accuracy gravitational waves from coalescing binary neutron stars in numerical relativity. II. Systematic study on the equation of state, binary mass, and mass ratio.  Phys. Rev. D, 101(8):084006, 2020.
  * [64] Carlos Palenzuela, Borja Miñano, Daniele Viganò, Antoni Arbona, Carles Bona-Casas, Andreu Rigo, Miguel Bezares, Carles Bona, and Joan Massó.  A Simflowny-based finite-difference code for high-performance computing in numerical relativity.  Class. Quant. Grav., 35(18):185007, 2018.
  * [65] Maria Okounkova.  Stability of Rotating Black Holes in Einstein Dilaton Gauss-Bonnet Gravity.  Phys. Rev. D, 100(12):124054, 2019.
  * [66] Michael Boyle et al.  The SXS Collaboration catalog of binary black hole simulations.  Class. Quant. Grav., 36(19):195006, 2019.
  * [67] Lawrence E. Kidder et al.  SpECTRE: A Task-based Discontinuous Galerkin Code for Relativistic Astrophysics.  J. Comput. Phys., 335:84–114, 2017.
  * [68] Nils Deppe et al.  Simulating magnetized neutron stars with discontinuous Galerkin methods.  Phys. Rev. D, 105(12):123031, 2022.
  * [69] S. Rosswog and P. Diener.  SPHINCS_BSSN: A general relativistic Smooth Particle Hydrodynamics code for dynamical spacetimes.  Class. Quant. Grav., 38(11):115002, 2021.
  * [70] S. Typel, M. Oertel, and T. Klähn.  CompOSE CompStar online supernova equations of state harmonising the concert of nuclear physics and astrophysics compose.obspm.fr.  Phys. Part. Nucl., 46(4):633–664, 2015.
  * [71] S. Typel et al.  CompOSE Reference Manual.  Eur. Phys. J. A, 58(11):221, 2022.
  * [72] M. Oertel, M. Hempel, T. Klähn, and S. Typel.  Equations of state for supernovae and compact stars.  Rev. Mod. Phys., 89(1):015007, 2017.
  * [73] James W. York, Jr.  Conformal ’thin sandwich’ data for the initial-value problem.  Phys. Rev. Lett., 82:1350–1353, 1999.
  * [74] Harald P. Pfeiffer and James W. York, Jr.  Extrinsic curvature and the Einstein constraints.  Phys. Rev. D, 67:044022, 2003.
  * [75] Eric Gourgoulhon.  3+1 formalism and bases of numerical relativity.  3 2007.
  * [76] Luciano Rezzolla and Olindo Zanotti.  Relativistic Hydrodynamics.  Oxford University Press, 09 2013.
  * [77] Jocelyn S. Read, Benjamin D. Lackey, Benjamin J. Owen, and John L. Friedman.  Constraints on a phenomenologically parameterized neutron-star equation of state.  Phys. Rev. D, 79:124032, 2009.
  * [78] C. Ronchi, R. Iacono, and P.S. Paolucci.  The “cubed sphere”: A new method for the solution of partial differential equations in spherical geometry.  Journal of Computational Physics, 124(1):93–114, 1996.
  * [79] Wolfgang Tichy, Alireza Rashti, Tim Dietrich, Reetika Dudi, and Bernd Brügmann.  Constructing binary neutron star initial data with high spins, high compactnesses, and high mass ratios.  Phys. Rev. D, 100(12):124046, 2019.
  * [80] Yousef Saad.  Iterative Methods for Sparse Linear Systems.  Society for Industrial and Applied Mathematics (SIAM, 3600 Market Street, Floor 6, Philadelphia, PA 19104, Philadelphia, Pa, 2003.
  * [81] Timothy A. Davis.  Algorithm 832: UMFPACK V4.3 - An unsymmetric-pattern multifrontal method.  ACM Transactions on Mathematical Software, 30(2):196–199, jun 2004\.
  * [82] Eric Gourgoulhon.  3 + 1 Formalism in General Relativity.  Number July. Springer, Berlin New York, 2011.
  * [83] Manuela Campanelli, Carlos O. Lousto, Yosef Zlochower, Badri Krishnan, and David Merritt.  Spin Flips and Precession in Black-Hole-Binary Mergers.  Phys. Rev. D, 75:064030, 2007.
  * [84] Andrew W. Steiner, Matthias Hempel, and Tobias Fischer.  Core-collapse supernova equations of state based on neutron star observations.  Astrophys. J., 774:17, 2013.
  * [85] Michael F. O’Boyle, Charalampos Markakis, Nikolaos Stergioulas, and Jocelyn S. Read.  Parametrized equation of state for neutron star matter with continuous sound speed.  Phys. Rev. D, 102(8):083027, 2020.
  * [86] Lee Lindblom.  Spectral Representations of Neutron-Star Equations of State.  Phys. Rev. D, 82:103011, 2010.
  * [87] Francois Foucart, Matthew D. Duez, Alana Gudinas, Francois Hebert, Lawrence E. Kidder, Harald P. Pfeiffer, and Mark A. Scheel.  Smooth Equations of State for High-Accuracy Simulations of Neutron Star Binaries.  Phys. Rev. D, 100(10):104048, 2019.
  * [88] Thomas W. Baumgarte and Stuart L. Shapiro.  Numerical Relativity: Solving Einstein’s Equations on the Computer.  Cambridge University Press, 2010.
  * [89] P. Haensel, J. L. Zdunik, and F. Douchin.  Equation of state of dense matter and the minimum mass of cold neutron stars.  Astron. Astrophys., 385:301, 2002.
  * [90] Wolfgang Tichy.  The initial value problem as it relates to numerical relativity.  Rept. Prog. Phys., 80(2):026901, 2017.
  * [91] Luc Blanchet.  Gravitational Radiation from Post-Newtonian Sources and Inspiralling Compact Binaries.  Living Rev. Rel., 17:2, 2014.
  * [92] Marcus Ansorg, A. Kleinwachter, and R. Meinel.  Highly accurate calculation of rotating neutron stars: detailed description of the numerical methods.  Astron. Astrophys., 405:711, 2003.
  * [93] Bengt Fornberg.  Generation of finite difference formulas on arbitrarily spaced grids.  Mathematics of Computation, 51(184):699–706, 1988.
  * [94] Richard L. Burden and Douglas Faires, J.  Numerical Analysis 9th Edition.  Brooks/Cole, 2011.

  *[BNS]: Binary Neutron Star
  *[EOS]: Equation of State
  *[PDE]: Partial Differential Equation
  *[NR]: Numerical Relativity
  *[ID]: Initial Data
  *[NS]: Neutron Star
  *[BC]: Boundary Condition
  *[XCTS]: Extended Conformal Thin Sandwich
  *[TOV]: Tolman-Oppenheimer-Volko