# The Lagrangian Numerical Relativity code SPHINCS_BSSN_v1.0

Stephan Rosswog 1,2, Francesco Torsello 2 and Peter Diener 3,4,∗

###### Abstract

We present version 1.0 of our Lagrangian Numerical Relativity code
SPHINCS_BSSN. This code evolves the full set of Einstein equations, but
contrary to other Numerical Relativity codes, it evolves the matter fluid via
Lagrangian particles in the framework of a high-accuracy version of Smooth
Particle Hydrodynamics (SPH). The major new elements introduced here are: i) a
new method to map the stress–energy tensor (known at the particles) to the
spacetime mesh, based on a local regression estimate; ii) additional measures
that ensure the robust evolution of a neutron star through its collapse to a
black hole; and iii) further refinements in how we place the SPH particles for
our initial data. The latter are implemented in our code SPHINCS_ID which now,
in addition to LORENE, can also couple to initial data produced by the initial
data library FUKA. We discuss several simulations of neutron star mergers
performed with SPHINCS_BSSN_v1.0, including irrotational cases with and
without prompt collapse and a system where only one of the stars has a large
spin ($\chi=0.5)$.

##  1 Keywords:

Numerical Relativity, gravitational waves, neutron stars, Smooth Particle
Hydrodynamics, Initial Data

##  2 Introduction

The first neutron star merger event, GW170817, opened the era of gravitational
wave-based multi-messenger astrophysics with a bang. The inspiral stages were
recorded in the detectors for about one minute [Abbott et al., 2017c], and
subsequently a firework was observed all across the electromagnetic spectrum
[Abbott et al., 2017d]. Starting with a (special) short gamma-ray burst (sGRB)
detected $1.7~{}$s after the peak gravitational wave (GW) emission [Goldstein
et al., 2017; Savchenko et al., 2017], followed by an initially blue and
subsequently rapidly reddening kilonova [Arcavi et al., 2017; Evans et al.,
2017; Cowperthwaite et al., 2017]. The event also displayed a rising X-ray
flux starting nine days after merger [Troja et al., 2017], peaking after 160
days, that then started to decline steeply afterwards [Hallinan et al., 2017;
Margutti et al., 2017; Lyman et al., 2018; Troja et al., 2018]. This was
interpreted as the imprints of a structured jet observed at an angle of $\sim
25^{\circ}$ from the jet core [Lamb and Kobayashi, 2017; Fong et al., 2017;
Lamb and Kobayashi, 2018; Hotokezaka et al., 2018]. Several years later,
broad-band synchroton afterglow was detected [Hajela et al., 2022; Troja et
al., 2020; Hajela et al., 2022] that was interpreted as a “kilonova afterglow”
due to a mildly relativistic ejection component interacting with the
interstellar medium [Nakar and Piran, 2011; Hotokezaka et al., 2015, 2018].

These observations allowed for a number of remarkable conclusions to be drawn.
For example, the pre-merger gravitational wave signal placed constraints on
the neutron star tidal deformability and therefore on the nuclear matter
equation of state [Abbott et al., 2017c]. The time delay between the GW peak
and the sGRB showed that GWs propagate at the speed of light to an enormous
precision [Abbott et al., 2017b]. The event also allowed for a determination
of the Hubble parameter [Abbott et al., 2017a] completely independent of
previous approaches. The bolometric light curve evolution of the kilonova was
remarkably consistent with a broad range of decaying r-process elements
[Metzger et al., 2010; Kasen et al., 2017; Rosswog et al., 2018], thereby
confirming the long-held suspicion that neutron star mergers are major sources
of r-process elements in the cosmos [Lattimer and Schramm, 1974; Symbalisty
and Schramm, 1982; Eichler et al., 1989; Rosswog et al., 1999; Freiburghaus et
al., 1999]. While one would naively expect a red kilonova due to the extreme
neutron-richness of the original neutron star material and the related large
opacities [Kasen et al., 2013], the early blue kilonova emission underlined
that about 0.01 M⊙ of the ejecta contained light (nucleon numbers $A<130$)
r-process material, which is also consistent with the identification of
strontium lines [Watson et al., 2019]. While underlining that a broad range of
heavy elements has been produced, these observations also stress the
importance of weak interactions that have transformed a substantial fraction
of the neutrons into protons to produce the light r-process material. The
“kilonova afterglow” in turn, hints at a broad velocity distribution within
the ejecta, extending to at least mildly relativistic velocities ($\gtrsim
0.6c$).

The phenomena described above illustrate the richness of the properties of the
ejected material and they stress the importance of understanding the detailed
properties of this $\sim 1\%$ of the total neutron star binary mass. Numerical
simulations play an integral part in understanding and interpreting multi-
messenger observations of compact binary systems. The initial generations of
simulation models had a strong focus on the strong-field spacetime dynamics,
the motion of the neutron star fluid in it and the resulting gravitational
wave emission, often with highly idealized microphysics. Today’s frontiers,
however, have shifted more towards a complex multi-physics modelling in which
General Relativity/strong-field gravity is only one ingredient out of many.
Apart from including more physical processes such as magnetic field evolution
or neutrino transport, the now observationally established connection with the
electromagnetic emission also places higher demands on the length and time
scales that need to be modelled in a neutron star merger event.

The vast majority of today’s Numerical Relativity codes makes use of Eulerian
hydrodynamics. While these codes have produced a multitude of important
results, a larger methodological variety would be desirable, for independent
checks of results but also to potentially address problems where established
methods struggle, as, for example, in the long-term evolution of merger
ejecta. The SPHINCS_BSSN code is the first Lagrangian hydrodynamics code that
solves the full set of Einstein equations. The first results, restricted to
standard relativistic hydrodynamics tests and oscillating and collapsing
neutron stars, were presented in Rosswog and Diener [2021]. In Diener et al.
[2022] the scope was extended to neutron star mergers, at that stage using
simple polytropic equations of state and LORENE-based initial conditions
[Grandclément et al., 2001; Gourgoulhon et al., 2001; LORENE, 2001] produced
with an extension of the “Artificial Pressure Method”, originally proposed in
Rosswog [2020b], to the case of neutron star binaries. Subsequently, further
technical improvements were introduced [Rosswog et al., 2022] and nuclear
matter properties were approximated in terms of piecewise polytropic equations
of state [Read et al., 2009]. We have further improved our simulation
technology and in this paper, we describe the ingredients of what we have
tagged as “version 1.0” of our code, SPHINCS_BSSN_v1.0. We emphasize the
latest improvements while still giving a broad overview of the complete
methodology. The new elements include a further refined method to map the
particle properties (specifically their stress–energy tensor) to the spacetime
mesh and additional measures to ensure that we can robustly evolve a neutron
star through its collapse to a black hole. We also describe further
refinements in the particle setup of our initial data, produced with the code
SPHINCS_ID SPHINCS_ID [2023].

Our article is structured as follows. In Section 3 we describe the
methodology, with Section 3.1 focusing on the hydrodynamics, Section 3.2 on
the equation of state, Section 3.3 on the spacetime evolution and Section 3.4
on how spacetime and matter evolution are coupled. In Section 3.5 we describe
measures to robustly evolve the collapse of a neutron star to a black hole and
in Section 4 we summarize our latest improvements in constructing SPH initial
conditions in full General Relativity. In Section 5 we discuss several
examples of neutron star mergers, and we summarize our results in Section 6

##  3 Methodology of SPHINCS_BSSN_v1.0

Here we describe the methodological elements of the code SPHINCS_BSSN_v1.0. We
only concisely summarize those parts that have been laid out already elsewhere
in the literature, but we describe in detail the elements that are presented
here for the first time. These are in particular: a) a more sophisticated
coupling between the spacetime and the fluid (via a local polynomial
regression estimate), b) specific measures (enhancement of grid resolution and
the potential transformation of fluid into “dust” particles) that enable us to
robustly simulate the formation of black holes.

###  3.1 Hydrodynamics

The hydrodynamic evolution equations in SPHINCS_BSSN are modelled via a high-
accuracy version of Smooth Particle Hydrodynamics (SPH), see Monaghan [2005];
Rosswog [2009]; Springel [2010]; Price [2012] and Rosswog [2015b] for reviews
of the method. The basics of the relativistic SPH equations have been derived
very explicitly in Section 4.2 of Rosswog [2009] and we will only present the
final equations here. Several accuracy-enhancing elements such as kernels,
gradient estimators and dissipation steering strategies (for either Newtonian
or relativistic cases) have been explored in a recent series of papers Rosswog
[2015a, 2020b, 2020a] and most of them are also implemented in
SPHINCS_BSSN_v1.0.
We use units in which $G=c=1$ and masses are measured in solar units. These
“code units” approximately correspond in physical units to $1.47663~{}$km for
lengths, to $4.92549\times 10^{-6}$ s for time and to $6.17797\times 10^{17}$
gcm-3 for densities. We further use the metric signature (-,+,+,+) and we
measure all energies in units of $m_{0}c^{2}$, where $m_{0}$ is the average
baryon mass111This quantity depends on the actual nuclear composition, but
simply using the atomic mass unit $m_{u}$ gives a precision of better than 1%.
We therefore use the approximation $m_{0}\approx m_{u}$ in the following. A
more detailed discussion can be found in Section 2.1.1 of Diener et al.
[2022]. Greek indices run from 0 to 3 and latin indices from 1 to 3.
Contravariant spatial indices of a vector quantity $\vec{w}$ at particle $a$
are denoted as $w_{a}^{i}$, while covariant ones will be written as
$(w_{i})_{a}$.
To discretize our fluid equations we choose a “computing frame” in which the
computations are performed. Quantities in this frame usually differ from those
calculated in the local fluid rest frame. The line element in a 3+1 split of
spacetime reads (e.g. Baumgarte and Shapiro [2010])

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

where $\alpha$ is the lapse function, $\beta^{i}$ the shift vector and
$\gamma_{ij}$ the spatial 3-metric. We use a generalized Lorentz factor

$\Theta\equiv\frac{1}{\sqrt{-g_{\mu\nu}v^{\mu}v^{\nu}}}\quad{\rm with}\quad
v^{\alpha}=\frac{dx^{\alpha}}{dt}.$ (2)

The coordinate velocities $v^{\alpha}$ are related to the four-velocities,
normalized as $U_{\alpha}U^{\alpha}=-1$, by

$v^{\alpha}=\frac{dx^{\alpha}}{dt}=\frac{U^{\alpha}}{\Theta}=\frac{U^{\alpha}}{U^{0}}.$
(3)

We choose the computing frame baryon number density $N$ as density variable,
which is related to the baryon number density as measured in the local fluid
rest frame, $n$, by

$N=\sqrt{-g}\,\Theta\,n.$ (4)

Here $g$ is the determinant of the spacetime metric. Note that this density
variable is very similar to what is used in Eulerian approaches [Alcubierre,
2008; Baumgarte and Shapiro, 2010; Rezzolla and Zanotti, 2013; Shibata, 2016].
We keep the baryon number of each SPH particle, $\nu_{a}$, constant so that
exact numerical baryon number conservation is hard-wired. At every
(Runge–Kutta sub-)step, the computing frame baryon number density at the
position of a particle $a$ is calculated via a weighted sum (actually very
similar to Newtonian SPH)

$N_{a}=\sum_{b}\nu_{b}\,W(|\vec{r_{a}}-\vec{r}_{b}|,h_{a}),$ (5)

where the smoothing length $h_{a}$ characterizes the support size of the SPH
smoothing kernel $W$, see below. As numerical momentum variable, we choose the
canonical momentum per baryon

$(S_{i})_{a}=(\Theta\mathcal{E}v_{i})_{a},$ (6)

where $\mathcal{E}=1+u+P/n$ is the relativistic enthalpy per baryon with $u$
being the internal energy per baryon and $P$ the gas pressure. The quantity
$S_{i}$ evolves according to

$\frac{d(S_{i})_{a}}{dt}=\left(\frac{d(S_{i})_{a}}{dt}\right)_{\rm
hyd}+\left(\frac{d(S_{i})_{a}}{dt}\right)_{\rm met},$ (7)

where the hydrodynamic part is given by

$\left(\frac{d(S_{i})_{a}}{dt}\right)_{\rm
hyd}=-\sum_{b}\nu_{b}\left\\{\frac{P_{a}}{N_{a}^{2}}D^{a}_{i}+\frac{P_{b}}{N_{b}^{2}}D^{b}_{i}\right\\}$
(8)

and the gravitational part by

$\left(\frac{d(S_{i})_{a}}{dt}\right)_{\rm
met}=\left(\frac{\sqrt{-g}}{2N}T^{\mu\nu}\frac{\partial g_{\mu\nu}}{\partial
x^{i}}\right)_{a}.$ (9)

Here we have used the abbreviations

$D^{a}_{i}\equiv\sqrt{-g_{a}}\;\frac{\partial W_{ab}(h_{a})}{\partial
x_{a}^{i}}\quad{\rm and}\quad D^{b}_{i}\equiv\sqrt{-g_{b}}\;\frac{\partial
W_{ab}(h_{b})}{\partial x_{a}^{i}}$ (10)

and $W_{ab}(h_{k})$ is a shorthand for $W(|\vec{r}_{a}-\vec{r}_{b}|/h_{k})$.
As numerical energy variable we use the canonical energy per baryon

$e_{a}=\left(S_{i}v^{i}+\frac{1+u}{\Theta}\right)_{a}=\left(\Theta\mathcal{E}v_{i}v^{i}+\frac{1+u}{\Theta}\right)_{a},$
(11)

which is evolved according to

$\frac{de_{a}}{dt}=\left(\frac{de_{a}}{dt}\right)_{\rm
hyd}+\left(\frac{de_{a}}{dt}\right)_{\rm met},$ (12)

with

$\left(\frac{de_{a}}{dt}\right)_{\rm
hyd}=-\sum_{b}\nu_{b}\left\\{\frac{P_{a}}{N_{a}^{2}}\;v_{b}^{i}\;D^{a}_{i}+\frac{P_{b}}{N_{b}^{2}}\;v_{a}^{i}\;D^{b}_{i}\right\\}$
(13)

and

$\left(\frac{de_{a}}{dt}\right)_{\rm
met}=-\left(\frac{\sqrt{-g}}{2N}T^{\mu\nu}\frac{\partial g_{\mu\nu}}{\partial
t}\right)_{a}.$ (14)

It is instructive to make the connection between our numerical momentum
variable $S_{i}$, Equation (6), and the Arnowitt–Deser–Misner (ADM) momentum
of the fluid. Given a spatial vector field $\xi^{i}$ which tends to a
Cartesian basis vector at spatial infinity, the ADM linear momentum along the
direction of $\xi^{i}$ is defined as [Gourgoulhon, 2012, Equation (8.40)] (see
also [Krishnan et al., 2007, Section II])

$\displaystyle P_{\xi}^{\rm
ADM}\coloneqq\dfrac{1}{8\pi}\int_{\partial\Sigma}\mathrm{d}\sigma_{j}\left(K^{j}{}_{i}-\delta^{j}{}_{i}\,K\right)\xi^{i},$
(15)

where $\partial\Sigma$ is the boundary of a spacelike hypersurface $\Sigma$
that extends up to spatial infinity, $\mathrm{d}\sigma_{j}$ its surface
element, $K^{j}{}_{i}$ is the extrinsic curvature and $K=K^{i}{}_{i}$ its
trace. Using Gauss’ theorem, the momentum constraint, and some geometry, the
integral in Equation (15) can be written as (see Appendix B)

$\displaystyle P_{\xi}^{\rm ADM}$
$\displaystyle=\int_{\Sigma}\mathrm{d}^{3}x\sqrt{\gamma}\left[\xi^{i}s_{i}+\dfrac{1}{16\pi}\left(K^{ij}-\gamma^{ij}\,K\right)\mathscr{L}_{\xi}\gamma_{ij}\right],$
(16)

where $\gamma_{ij}$ is the spatial metric, $s_{i}$ is the spatial part of the
momentum density measured by the Eulerian observer
$s_{\rho}\coloneqq-n^{\mu}\,T_{\mu\nu}\,\gamma^{\nu}{}_{\rho}$, and
$\mathscr{L}_{\xi}$ is the Lie derivative along $\xi^{i}$. The two terms in
the integrand of Equation ((16)) are the parts of the ADM linear momentum
determined by the fluid and the spacetime, respectively. The expression in
Equation ((16)) allows us to write the ADM momentum _of the fluid_ in terms of
the SPH canonical momentum $S_{i}$, after we relate the latter with the
Eulerian spatial momentum density $s_{i}$. We do this explicitly in Appendix
B. Here we show the final result and its SPH approximation:

$\displaystyle P_{\xi}^{\rm ADM,fluid}$
$\displaystyle=\int_{\Sigma}\mathrm{d}^{3}x\sqrt{\gamma}\,\xi^{i}s_{i}=\int_{\Sigma}\mathrm{d}^{3}x\,N\,\xi^{i}S_{i}\simeq\sum_{b}\nu_{b}\left(\xi^{i}S_{i}\right)_{b}$
(17)

with the index $b$ running over all the particles, $\nu_{b}$ being the baryon
number of particle $b$, $N$ defined in Equation (4),
$s_{i}=\Theta\,n\,\alpha\,S_{i}$, $\alpha$ is the lapse function, and
$\sqrt{-g}=\alpha\sqrt{\gamma}$ [Baumgarte and Shapiro, 2010, Equation
(2.124)]. The rightmost formula in Equation (17) can be used to compute an
estimate of the ADM momentum of the fluid using only SPH fields. Such an
estimate can then be compared with another one computed using Equation (15);
this comparison tells us about the error on the ADM momentum of the fluid
introduced when we model the fluid with SPH particles. We make this comparison
at the level of the initial data (ID) in Section 4.4.
For the kernel function that is needed for the SPH approximation, we have
implemented a large variety of choices. We have performed many numerical
experiments similar to the one shown in Figure 1, some of which are documented
in detail in Rosswog [2015a]. Here we only provide the details of our
preferred kernels. The first of these favorites is the Wendland C6-smooth
kernel [Wendland, 1995]

$W(q)=W^{WC6}(q)=\frac{\sigma}{h^{3}}(1-q)^{8}_{+}(32q^{3}+25q^{2}+8q+1),$
(18)

where we use exactly 300 constributing neighbour particles. The normalization
constant is $\sigma=1365/(64\pi)$ in 3D and the symbol $(.)_{+}$ denotes the
cutoff function max$(.,0)$. Other favorites include (some) members of the
family of the harmonic-like kernels (Cabezon et al. [2008])

$\displaystyle W(q)=W^{\rm
H}_{n}(q)=\frac{\sigma_{n}}{h^{3}}\left\\{\begin{array}[]{llr}1,&\;q=0,\\\\[11.38109pt]
\left(\frac{\sin\left[\frac{\pi}{2}q\right]}{\frac{\pi}{2}q}\right)^{n}&\;0<q<2,\\\\[11.38109pt]
0,&\;{\rm else},\end{array}\right.$ (22)

namely those with $n=7$ and 8. Out of this family, we chose, after ample
experiments, the kernel with $n=8$ for which we use exactly 220 contributing
neighbour particles. For this kernel the normalization constant is
$\sigma_{8}=1.17851074088357$ in 3D. Contrary to Wendland kernels [Wendland,
1995], this kernel is not immune against the (benign) pairing instability, but
it provides an excellent density estimate. We show in Figure 1 the result of a
density measurement experiment. Particles are placed on a cubic lattice and
masses are assigned so that the mass density is exactly unity. We then use
several kernel functions, the commonly used cubic spline kernel [Monaghan,
1992], the Wendland C2, C4, and C6 kernels (see Dehnen and Aly [2012] for the
explicit expressions) and the $W^{\rm H}_{8}$ kernel, Equation (22).

Figure 1: Error on the density estimate for different kernels as a function of
support size for particles arranged on a cubic lattice. The parameter $\eta$
determines the smoothing length $h$ as a multiple of the typical particle
spacing via $h=\eta(\nu/N)^{1/3}$, the support radius of all kernels is $2h$.
We have also indicated the corresponding number of neighbour particles.

The smoothing lengths $h$ in this experiment are set as multiples of the
typical particle separation $(\nu/N)^{1/3}$ , $h=\eta(\nu/N)^{1/3}$, and each
of the kernels has a support radius of $2h$. We have also indicated some
approximate numbers of contributing particles (“neighbours”) for our cubic
lattice. For neighbour numbers just beyond 200, the density accuracy in this
experiment is about three orders of magnitude better for $W^{\rm H}_{8}$
compared to the Wendland C6 kernel. While it is a priori not entirely clear
how to weigh the excellent performance in this idealized experiment against
the desirable property of the Wendland kernels to maintain a very regular
particle distribution during dynamical evolution [Rosswog, 2015a], we choose
in this work the $W^{\rm H}_{8}$-kernel and we find very satisfactory results.
On inspecting the simulations presented here, we do not see any significant
pairing among the particles and the density distributions appear noise-free
while exhibiting sharp features.
To keep the numerical noise at a minimum, we choose at each time step the
smoothing length of each particle $a$, $h_{a}$, so that there are exactly 220
contributing neighbour particles within the support radius $2h_{a}$. In other
words, at each time step the smoothing length of particle $a$, $h_{a}$, is set
so that $2h_{a}$ equals the distance to the 221-th closest SPH particle. This
particle sits by definition at the radius where the kernel becomes zero, so
that exactly 220 particles have a non-zero contribution. This approach ensures
a very smooth and subtle evolution of the smoothing length and avoids the
introduction of noise through the update of the smoothing length. In practice,
this is achieved by using our fast RCB-tree [Gafton and Rosswog, 2011], for
more details on the exact procedure, we refer to the MAGMA2 code paper
[Rosswog, 2020b] where the same approach was used. A large number of
experiments confirms that we get very similar results when using $W^{WC6}$ and
$W^{H}_{8}$.
To robustly handle relativistic shocks, our momentum and energy equations are
augmented with dissipative terms. These terms consist of an artificial
dissipation and an artificial conductivity part. In both of these terms we
apply a slope-limited reconstruction to the mid-point of each particle pair
and this reconstruction approach has been shown to massively reduce unwanted
dissipative effects [Rosswog, 2020b]. In order to further reduce dissipation
where it is not needed, we make our dissipation parameters time-dependent;
they are increased when a shock or numerical noise is detected, but otherwise
they decay exponentially to a small value (here chosen as $\alpha_{0}=0.1$).
Since no changes compared to our previous work have been made, we refer the
interested reader for the equations, implementation details and tests to
Rosswog et al. [2022].
The quantities that we evolve numerically, $S_{i}$ and $e$, together with the
density $N$, see Equation(5), are numerically very convenient, but they are
not the physical quantities that we are interested in. We therefore have to
recover the physical quantities $n$, $u$, $v^{i}$, $P$ from $N$, $S_{i}$, $e$
at every integration (sub-)step. This “recovery” step is done in a very
similar way as in Eulerian approaches. For polytropic equations of state, our
strategy is described in detail in Section 2.2.4 of Rosswog and Diener [2021],
and the modifications needed for piecewise polytropic EOSs are laid out
explicitly in Appendix A of Rosswog et al. [2022].

###  3.2 Equations of state (EOS)

To close the system of hydrodynamic equations we need an equation of state.
Currently, we are using piecewise polytropic approximations to cold nuclear
matter equations of state [Read et al., 2009], that are enhanced with an ideal
gas-type thermal contribution to both pressure and specific internal energy, a
common approach in Numerical Relativity simulations. For explicit expressions
please see Appendix A of Rosswog et al. [2022]. To date, we have implemented
14 piecewise polytropic equations of state, but since the effects coming from
different EOSs are not the topic of this code paper, we restrict ourselves to
results obtained for the APR3 EOS [Akmal et al., 1998] only. This EOS allows
for a maximum mass of $M_{\rm TOV}^{\rm max}=2.39$ M⊙ and a 1.4 M⊙ neutron
star has a dimensionless tidal deformability of $\Lambda_{1.4}=390$. Indirect
arguments and the statistics of the radio pulsar/X-ray neutron star
distribution point to values in the range of $\sim 2.2-2.4$ M⊙ [Fryer et al.,
2015; Antoniadis et al., 2016; Margalit and Metzger, 2017; Bauswein et al.,
2017; Shibata et al., 2017; Rezzolla et al., 2018; Alsing et al., 2018] for
the “best educated guess” of the maximum neutron star mass and a recent
Bayesian study [Biswas and Datta, 2021] suggests a maximum TOV mass of
2.52${}^{+0.33}_{-0.29}$ M⊙, all broadly consistent with our choice of the
APR3 EOS. More sophisticated treatments of high-density nuclear matter physics
will be addressed in the future.

###  3.3 Spacetime evolution

We evolve the spacetime according to the (“$\Phi$-version” of the) BSSN
equations [Shibata and Nakamura, 1995; Baumgarte and Shapiro, 1999]. We have
written wrappers around code extracted from the well tested McLachlan thorn of
the Einstein Toolkit Löffler et al. [2012]; Einstein Toolkit web page [2020].
The dynamical variables used in this method are related to the
Arnowitt–Deser–Misner (ADM) variables $\gamma_{ij}$ (3-metric), $K_{ij}$
(extrinsic curvature), $\alpha$ (lapse function) and $\beta^{i}$ (shift
vector) and they read

$\displaystyle\phi$ $\displaystyle=\frac{1}{12}\log(\gamma),$ (23)
$\displaystyle\tilde{\gamma}_{ij}$ $\displaystyle=e^{-4\phi}\gamma_{ij},$ (24)
$\displaystyle K$ $\displaystyle=\gamma^{ij}K_{ij},$ (25)
$\displaystyle\tilde{\Gamma}^{i}$
$\displaystyle=\tilde{\gamma}^{jk}\tilde{\Gamma}^{i}_{jk},$ (26)
$\displaystyle\tilde{A}_{ij}$
$\displaystyle=e^{-4\phi}\left(K_{ij}-\frac{1}{3}\gamma_{ij}K\right),$ (27)

where $\gamma=\operatorname{det}(\gamma_{ij})$, $\tilde{\Gamma}^{i}_{jk}$ are
the Christoffel symbols related to the conformal metric $\tilde{\gamma}_{ij}$
and $\tilde{A}_{ij}$ is the conformally rescaled, traceless part of the
extrinsic curvature. The corresponding evolution equations read

$\displaystyle\partial_{t}\phi$ $\displaystyle=-\frac{1}{6}\left(\alpha
K-\partial_{i}\beta^{i}\right)+\beta^{i}\bar{\partial}_{i}\phi,$ (28)
$\displaystyle\partial_{t}\tilde{\gamma}_{ij}$
$\displaystyle=-2\alpha\tilde{A}_{ij}+\tilde{\gamma}_{ik}\partial_{j}\beta^{k}+\tilde{\gamma}_{jk}\partial_{i}\beta^{k}-\frac{2}{3}\tilde{\gamma}_{ij}\partial_{k}\beta^{k}+\beta^{k}\bar{\partial}_{k}\tilde{\gamma}_{ij},$
(29) $\displaystyle\partial_{t}K$
$\displaystyle=-e^{-4\phi}\left(\tilde{\gamma}^{ij}\left[\partial_{i}\partial_{j}\alpha+2\partial_{i}\phi\partial_{j}\alpha\right]-\tilde{\Gamma}_{\mathrm{(n)}}^{i}\partial_{i}\alpha\right)+\alpha\left(\tilde{A}^{i}_{j}\tilde{A}^{j}_{i}+\frac{1}{3}K^{2}\right)+\beta^{i}\bar{\partial}_{i}K+4\pi\alpha(\rho+s),$
(30) $\displaystyle\partial_{t}\tilde{\Gamma}^{i}$
$\displaystyle=-2\tilde{A}^{ij}\partial_{j}\alpha+2\alpha\left(\tilde{\Gamma}^{i}_{jk}\tilde{A}^{jk}-\frac{2}{3}\tilde{\gamma}^{ij}\partial_{j}K+6\tilde{A}^{ij}\partial_{j}\phi\right)$
$\displaystyle+\tilde{\gamma}^{jk}\partial_{j}\partial_{k}\beta^{i}+\frac{1}{3}\tilde{\gamma}^{ij}\partial_{j}\partial_{k}\beta^{k}-\tilde{\Gamma}_{\mathrm{(n)}}^{j}\partial_{j}\beta^{i}+\frac{2}{3}\tilde{\Gamma}_{\mathrm{(n)}}^{i}\partial_{j}\beta^{j}+\beta^{j}\bar{\partial}_{j}\tilde{\Gamma}^{i}-16\pi\alpha\tilde{\gamma}^{ij}s_{j},$
(31) $\displaystyle\partial_{t}\tilde{A}_{ij}$
$\displaystyle=e^{-4\phi}\left[-\partial_{i}\partial_{j}\alpha+\tilde{\Gamma}^{k}_{ij}\partial_{k}\alpha+2\left(\partial_{i}\alpha\partial_{j}\phi+\partial_{j}\alpha\partial_{i}\phi\right)+\alpha
R_{ij}\right]^{\mathrm{TF}}+\alpha(K\tilde{A}_{ij}-2\tilde{A}_{ik}\tilde{A}^{k}_{j})$
$\displaystyle+\tilde{A}_{ik}\partial_{j}\beta^{k}+\tilde{A}_{jk}\partial_{i}\beta^{k}-\frac{2}{3}\tilde{A}_{ij}\partial_{k}\beta^{k}+\beta^{k}\bar{\partial}_{k}\tilde{A}_{ij}-e^{-4\phi}\alpha
8\pi\left(T_{ij}-\frac{1}{3}\gamma_{ij}s\right),$ (32)

where

$\displaystyle\rho$
$\displaystyle=\frac{1}{\alpha^{2}}(T_{00}-2\beta^{i}T_{0i}+\beta^{i}\beta^{j}T_{ij}),$
(33) $\displaystyle s$ $\displaystyle=\gamma^{ij}T_{ij},$ (34) $\displaystyle
s_{i}$ $\displaystyle=-\frac{1}{\alpha}(T_{0i}-\beta^{j}T_{ij}),$ (35)

and $\beta^{i}\bar{\partial}_{i}$ denote partial derivatives that are upwinded
based on the shift vector. The superscript “TF” in the evolution equation of
$\tilde{A}_{ij}$ denotes the trace-free part of the bracketed term. Finally
$R_{ij}=\tilde{R}_{ij}+R^{\phi}_{ij}$, where

$\displaystyle\tilde{\Gamma}_{ijk}$
$\displaystyle=\frac{1}{2}\left(\partial_{k}\tilde{\gamma}_{ij}+\partial_{j}\tilde{\gamma}_{ik}-\partial_{i}\tilde{\gamma}_{jk}\right),$
(36) $\displaystyle\tilde{\Gamma}_{ij}^{\;\;\;k}$
$\displaystyle=\tilde{\gamma}^{kl}\tilde{\Gamma}_{ijl},$ (37)
$\displaystyle\tilde{\Gamma}^{i}_{jk}$
$\displaystyle=\tilde{\gamma}^{il}\tilde{\Gamma}_{ljk},$ (38)
$\displaystyle\tilde{\Gamma}_{\mathrm{(n)}}^{i}$
$\displaystyle=\tilde{\gamma}^{jk}\tilde{\Gamma}^{i}_{jk},$ (39)
$\displaystyle\tilde{R}_{ij}$
$\displaystyle=-\frac{1}{2}\tilde{\gamma}^{kl}\partial_{k}\partial_{l}\tilde{\gamma}_{ij}+\tilde{\gamma}_{k(i}\partial_{j)}\tilde{\Gamma}^{k}+\tilde{\Gamma}_{\mathrm{(n)}}^{k}\tilde{\Gamma}_{(ij)k}+\tilde{\Gamma}^{k}_{il}\tilde{\Gamma}_{jk}^{\;\;\;l}+\tilde{\Gamma}^{k}_{jl}\tilde{\Gamma}_{ik}^{\;\;\;l}+\tilde{\Gamma}^{k}_{il}\tilde{\Gamma}_{kj}^{\;\;\;l},$
(40) $\displaystyle R^{\phi}_{ij}$
$\displaystyle=-2\left(\partial_{i}\partial_{j}\phi-\tilde{\Gamma}^{k}_{ij}\partial_{k}\phi\right)-2\tilde{\gamma}_{ij}\tilde{\gamma}^{kl}\left(\partial_{k}\partial_{l}\phi-\tilde{\Gamma}^{m}_{kl}\partial_{m}\phi\right)+4\partial_{i}\phi\partial_{j}\phi-4\tilde{\gamma}_{ij}\tilde{\gamma}^{kl}\partial_{k}\phi\partial_{l}\phi.$
(41)

The derivatives on the right-hand side of the BSSN equations are evaluated via
standard Finite Differencing techniques and, unless mentioned otherwise, we
use sixth order differencing as a default. We have recently implemented a
fixed mesh refinement for the spacetime evolution, which is described in
detail in Diener et al. [2022], to which we refer the interested reader. For
the gauge choices, we use a variant of “1+log-slicing,” where the lapse is
evolved according to

$\partial_{t}\alpha=-2\alpha K$ (42)

and a variant of the “$\Gamma$-driver” shift condition with

$\partial_{t}\beta^{i}=\frac{3}{4}(\tilde{\Gamma}^{i}-\eta\beta^{i}),$ (43)

where $\eta$ is the “$\beta$-driver” parameter.

###  3.4 Coupling between spacetime and matter

The SPHINCS_BSSN approach of evolving the spacetime on a mesh and the matter
fluid via particles requires a continuous information exchange: the gravity
part of the particle evolution is driven by derivatives of the metric, see
Eqs. (9) and (14), which are known on the mesh, while the stress–energy
tensor, needed in Equation (28)-(35), is known on the particles. This
bidirectional information exchange is needed at every Runge–Kutta substep; in
our case, with an optimal 3rd order Runge–Kutta algorithm [Gottlieb and Shu,
1998], three times per numerical time step.
The mesh-to-particle mapping step is performed via 5th order Hermite
interpolation, that we developed [Rosswog and Diener, 2021] extending the work
of Timmes and Swesty [2000]. Contrary to a standard Lagrange polynomial
interpolation, the Hermite interpolation guarantees that the metric remains
twice differentiable as particles pass from one grid cell to another and
therefore avoids the introduction of additional noise. Our approach is
explained in detail in Section 2.4 of Rosswog and Diener [2021] to which we
refer the interested reader.
Here, we provide a further improvement to the more challenging of the two
steps, the particle-to-mesh mapping. As in our previous work [Diener et al.,
2022; Rosswog et al., 2022], we follow a “Multi-dimensional Optimal Order
Detection” (MOOD) strategy. The mapping is performed simultaneously with
different orders, and the most accurate solution is selected out of the
possible results (according to some error measure, see below), provided that
the solution meets additional physical admissibility conditions.
This is somewhat similar in spirit to ENO reconstruction schemes in mesh-based
hydrodynamics, see Zhang and Shu [2016] for a concise summary of ENO and WENO
schemes, in the sense that one does not pretend to know a priori which
reconstruction order (or here: interpolation order) is ”good enough”. Instead
one tries a variety of options and selects the best one. Our scheme is similar
to ENO methods where one tries different stencils and selects the smoothest
one according to a suitable criterion. (WENO actually goes one step further by
using a weighted sum of the different options to combine them into potentially
higher order.)
Whereas in Rosswog et al. [2022] we used kernel functions of different orders
borrowed from vortex methods [Cottet and Koumoutsakos, 2000; Cottet et al.,
2014], we determine here the best fit to the particle properties at a grid
point by using the local polynomial regression estimate that we describe in
the next section.

####  3.4.1 Local Regression Estimate (LRE)

The task is to find the function that best fits the values of the particles
surrounding the grid point $\vec{r}^{\rm G}$, see Figure 2.

Figure 2: Geometry of the particle-to-mesh mapping: function values that are
known at particle positions (blue circles) are to be mapped to a grid point
(black square).

If we assume that the values at the particle positions $\\{f^{p}\\}$ are
described by a function $f(\vec{r})$, we can Taylor-expand this function
around the desired grid point $\vec{r}^{\rm G}$:

$\displaystyle f(\vec{r})=$ $\displaystyle f(\vec{r}^{\rm
G})+(\partial_{i}f)_{\vec{r}^{\rm G}}(\vec{r}-\vec{r}^{\rm
G})^{i}+\frac{1}{2!}(\partial_{ij}f)_{\vec{r}^{\rm G}}\;(\vec{r}-\vec{r}^{\rm
G})^{i}(\vec{r}-\vec{r}^{\rm G})^{j}$ (44) $\displaystyle+$
$\displaystyle\frac{1}{3!}(\partial_{ijk}f)_{\vec{r}^{\rm
G}}\;(\vec{r}-\vec{r}^{\rm G})^{i}(\vec{r}-\vec{r}^{\rm
G})^{j}(\vec{r}-\vec{r}^{\rm G})^{k}$ $\displaystyle+$
$\displaystyle\frac{1}{4!}(\partial_{ijkl}f)_{\vec{r}^{\rm
G}}\;(\vec{r}-\vec{r}^{\rm G})^{i}(\vec{r}-\vec{r}^{\rm
G})^{j}(\vec{r}-\vec{r}^{\rm G})^{k}(\vec{r}-\vec{r}^{\rm G})^{l}$
$\displaystyle+$ higher order terms.

This Taylor expansion can be interpreted as a polynomial approximation of a
given order where the basis functions have been shifted to the point of
interest $\vec{r}^{\rm G}$. The local approximation $\tilde{f}(\vec{r})$ of
$f(\vec{r})$ around the point $\vec{r}^{\rm G}$then reads

$\tilde{f}^{G}(\vec{r})=\vec{\beta}^{G}\cdot\vec{P}^{G}(\vec{r}),$ (45)

where the coefficient vector reads

$\displaystyle\vec{\beta}^{G}=\Big{[}$ $\displaystyle f(\vec{r}^{\rm
G}),\;\;(\partial_{x}f)_{\vec{r}^{\rm G}},\,(\partial_{y}f)_{\vec{r}^{\rm
G}},\,(\partial_{z}f)_{\vec{r}^{\rm G}},\;\;$
$\displaystyle\frac{1}{2}(\partial_{xx}f)_{\vec{r}^{\rm
G}},\,(\partial_{xy}f)_{\vec{r}^{\rm G}},\,(\partial_{xz}f)_{\vec{r}^{\rm
G}},\,\frac{1}{2}(\partial_{yy}f)_{\vec{r}^{\rm
G}},\,(\partial_{yz}f)_{\vec{r}^{\rm
G}},\,\frac{1}{2}(\partial_{zz}f)_{\vec{r}^{\rm G}},\;\;\cdots\Big{]},$ (46)

and the “shifted basis functions” read

$\displaystyle\vec{P}^{G}(\vec{r})=\Big{[}$ $\displaystyle
1,\;\;\Delta^{G}x,\,\Delta^{G}y,\,\Delta^{G}z,\;\;$ (47)
$\displaystyle\Delta^{G}x\Delta^{G}x,\,\Delta^{G}x\Delta^{G}y,\,\Delta^{G}x\Delta^{G}z,\,\Delta^{G}y\Delta^{G}y,\,\Delta^{G}y\Delta^{G}z,\,\Delta^{G}z\Delta^{G}z,\;\;\cdots\Big{]}^{T},$
(48)

where $\Delta^{G}\vec{r}=\vec{r}-\vec{r}^{\rm
G}=[\Delta^{G}x,\Delta^{G}y,\Delta^{G}z]^{T}$, and $\Delta^{G}x=x-x^{G}$
(similarly for the other components). The degrees of freedom ($\mathrm{DoF}$)
for a given number of dimensions $d$ and a maximum polynomial order of the
basis $m$ are given by

$\mathrm{DoF}=\frac{(d+m)!}{d!\,m!},$ (49)

that is, in 3D, we have 1, 4, 10, 20 and 35 DoF for constant, linear,
quadratic, cubic and quartic polynomials.
The optimal coefficients at $\vec{r}^{G}$, $\vec{\beta}^{G}$, are found by
minimizing the error functional

$\epsilon^{G}\equiv\sum_{p}[f_{p}-\tilde{f}^{G}(\vec{r}_{p})]^{2}\;W_{pG}=\sum_{p}\left[f_{p}-{\sum_{i=1}^{\mathrm{DoF}}\beta^{G}_{i}P^{G}_{i}(\vec{r}_{p})}\right]^{2}\;W_{pG},$
(50)

where $W_{pG}=W(\vec{r}^{p}-\vec{r}^{G},l_{p})$ is a smooth weighting function
that ensures that particles further away from the grid point are weighted less
in the error functional, and the $p$-sum runs over all contributing particles.
The kernel width is set by $l_{p}=(\nu_{p}/N_{p})^{1/3}$. One could choose,
for example, compactly supported SPH kernels for this weighting function; we
choose simple tensor products of 1D $M_{4}$-kernels. The exact form of the
weight function does not have a strong influence on the resulting error
measure.
The optimal coefficients that minimize the error at $\vec{r}^{G}$ are
determined via

$\left(\frac{\partial\epsilon^{G}}{\partial\beta_{i}^{G}}\right)_{\vec{r}^{\rm
G}}\stackrel{{\scriptstyle!}}{{=}}0$ (51)

and, with a few steps of algebra, one finds

$\beta^{G}_{i}=\left(M_{ik}\right)^{-1}\;B_{k},$ (52)

where

$M_{ik}=\sum_{p}P^{G}_{i}(\vec{r}_{p})\;P^{G}_{k}(\vec{r}_{p})\;W_{pG}$ (53)

is the “moment matrix” and

$B_{k}=\sum_{p}f_{p}\;P^{G}_{k}(\vec{r}_{p})\;W_{pG}$ (54)

is a vector depending on the function values at the particles. Since the
moment matrix does not depend on the function values themselves (only on the
relative positions), it can be used for several function vectors (here one for
each $T_{\mu\nu}$ component). With increasing polynomial order the condition
number of the moment matrix $\mathcal{C}\equiv||MM^{-1}||$, a measure for how
close a matrix is to being singular, can become very large, therefore we use
the singular value decomposition [Press et al., 1992] to solve the system in
Equation (52). The condition numbers, however, can be massively reduced by
using re-scaled basis functions which we illustrate in Appendix D and Figure
A2. In practice, we use re-scaled basis functions.
The function value estimate at the grid point (see Equation (46)),
$\tilde{f}(\vec{r}^{G})$ is the first component of $\beta^{G}$, the
derivatives are
$\partial_{x}\tilde{f}(\vec{r}^{G}),\partial_{y}\tilde{f}(\vec{r}^{G})$,
$\partial_{z}\tilde{f}(\vec{r}^{G})$ are the components two to four, and so
on. For the case where we only allow the lowest polynomial order (i.e., a
constant polynomial), the moment matrix has only one element

$M=\sum_{p}W_{pG}\Longrightarrow M^{-1}=\frac{1}{\sum_{p}W_{pG}},$ (55)

and the function vector becomes

$B=\sum_{p}f_{p}W_{pG}.$ (56)

In this case the function value at the grid point is estimated as

$\tilde{f}^{G}=\frac{\sum_{p}f_{p}W_{pG}}{\sum_{p}W_{pG}}.$ (57)

In other words, this is just the straightforward kernel-weighted average of
the values at the contributing particles (with an exact partition of unity
enforced). We show an example of the LRE function approximation in Appendix D.

####  3.4.2 An LRE-based MOOD approach

While in a well-sampled region, as for example the stellar interior, a high-
order LRE approximation likely provides the most accurate function estimate,
this is not necessarily true near the stellar surface. There, due to the Gibbs
phenomenon, spurious oscillations may occur. Therefore, we calculate
$T_{\mu\nu}$ estimates for different polynomial orders $m$,
$\tilde{T}_{\mu\nu}^{{\rm G},m}$, and then select the “optimal order” that
best represents the particle values and is physically admissible. We use the
following error measure:

$\displaystyle E^{{\rm
G},m}\equiv\sum_{p}W_{pG}\left[\sum_{\mu,\nu}\left\\{{\tilde{T}}_{\mu\nu}^{{\rm
G},m}(\vec{r}_{p})-T_{\mu\nu,p}\right\\}^{2}\right]=\sum_{p}W_{pG}\left[\sum_{\mu,\nu}\left\\{\left(\vec{\beta}^{G,m}_{\mu,\nu}\cdot
P^{G}(\vec{r}_{p})\right)-T_{\mu\nu,p}\right\\}^{2}\right].$ (58)

In other words, based on the optimal coefficients at the grid point, we
estimate the function values at each particle position that contributes and
calculate the weighted quadratic deviation as error measure. This approach
differs from our earlier one [Rosswog et al., 2022] by not using pre-defined
kernel functions, but instead applying the LRE-approximation, and by
considering all $T_{\mu\nu}$ components in the error measure, rather than just
$T_{00}$ as before.
The most straightforward MOOD estimate would be to calculate estimates for
different polynomial orders $m$ and to select the solution with the smallest
value of $E^{{\rm G},m}$. Unfortunately, this only works in nearly all cases.
In very few cases, where a grid point lies outside the neutron star surface,
but the finite-size SPH particles still contribute to this point, we found
that the approximation with the smallest error may deliver values much larger
than those at the contributing particles, or even unphysical values such as a
negative $T_{00}$ (energy density). For these reasons, additional measures
need to be taken near the stellar surface. However, the question to be
answered is: how does an SPH particle know that it is near the surface?

####  3.4.3 Detecting not well-embedded grid points

Here we use a simple, yet, as it turns out, very robust method to detect
whether a grid point is well engulfed by particles or not. In a first step, at
each particle position we numerically calculate an estimate for an expression
that has an analytically known result and where the deviation from the exact
result can be used to identify surface particles. For this purpose we chose
$\nabla\cdot\vec{r}=3$ and the numerical estimate given by

$\left(\nabla\cdot\vec{r}\right)_{a}=\sum_{b}\frac{\nu_{b}}{N_{b}}(\vec{r}_{b}-\vec{r}_{a})\cdot\nabla_{a}W_{ab}(h_{a}).$
(59)

This expression is one of the standard SPH discretizations (similar to the
commonly used expression for $\nabla\cdot\vec{v}$, see Equation (31) in
Rosswog [2009]). To avoid another neighbour-loop over all particles,
expression (59) can be conveniently calculated alongside the SPH derivatives.
This means that the update of $\nabla\cdot\vec{r}$ is lagging behind by one
third of a time step. The property of being at the surface changes on a much
longer time scale and only averages of the deviations are used, see below, so
that using a value of $\nabla\cdot\vec{r}$ calculated a third of a time step
earlier is well suited for our purposes. In deriving SPH expressions, surface
terms are usually neglected and therefore the expression Equation (59) only
yields an accurate approximation to the exact value of 3 if it is embedded
from all sides with particles. If instead the expression is evaluated near a
surface, contributing particles are missing on one side and therefore the
numerical estimate is substantially smaller than the theoretical value. From
the relative error

$\delta_{a}\equiv\frac{|(\nabla\cdot\vec{r})_{a}-3|}{3}$ (60)

we calculate the average deviation
$\langle\delta\rangle_{G}=(\sum_{b=1}^{n}\delta_{b})/n$ over the particles
that contribute to the error measure Equation (58). If
$\langle\delta\rangle_{G}$ is above a given threshold, the corresponding grid
point is identified as being outside the particle surface. We show an example
from the inspiral of two neutron stars in Figure 3. In the stellar interior,
the values of $\langle\delta\rangle_{G}$ are $\approx 0.005$, while those at
the surface reach $\sim 0.1$. After some experimentation, we have chosen a
threshold of 0.05 for $\langle\delta\rangle_{G}$; for grid points at which
$\langle\delta\rangle_{G}$ exceed this threshold, we only use the lowest order
mapping, $m=0$, see below. This approach robustly avoids all ”outlier points”
in the mapping, and the mapped values of $T_{\mu\nu}$ accurately reflect the
matter distribution.

Figure 3: The average value of $\langle\delta\rangle_{G}$ for all the
particles that contribute to the mapping of the stress–energy tensor at grid
points in the $xy$-plane. For grid points inside the stars, the value is
typically $\delta\approx 0.005$.

Figure 4: At the top, we show a cut of the density at $t=2.95\mathrm{ms}$ for
the case of a equal mass binary system with 2 1.5 M⊙neutron stars with the
APR3 equation of state. At the bottom, we show the corresponding mapping of
the $tt$-component of the stress–energy tensor onto the grid.

####  3.4.4 Summary of the particle-to-mesh mapping algorithm

The first step in the algorithm consists in identifying the particles that
contribute to a given grid point $G$. This list of particles is identified via
a hash-grid as described in detail in Section 2.1.3 of Diener et al. [2022].
Subsequently, we perform the following steps:

  * •

calculate the LRE estimates $\tilde{T}_{\mu\nu}^{{\rm G},m}$ for the
polynomial orders $m=0,1,2,3$ using Equation (45)

  * •

if $\langle\delta\rangle_{G}>0.05$ choose $\tilde{T}_{\mu\nu}^{{\rm G},m=0}$
since this is a grid point outside the particle surface

  * •

if more than 40 particles (= twice the number of degrees of freedom for cubic
polynomials) contribute and the error $E^{G,3}$ is smallest and
$\tilde{T}_{00}^{{\rm G},m=3}>0$, choose $\tilde{T}_{\mu\nu}^{{\rm G},m=3}$

  * •

if more than 20 particles (= twice the number of degrees of freedom for
quadratic polynomials) contribute and the error $E^{G,2}$ is smallest and
$\tilde{T}_{00}^{{\rm G},m=2}>0$, choose $\tilde{T}_{\mu\nu}^{{\rm G},m=2}$

  * •

if more than 8 particles (= twice the number of degrees of freedom for linear
polynomials) contribute and the error $E^{G,1}$ is smallest and
$\tilde{T}_{00}^{{\rm G},m=1}>0$, choose $\tilde{T}_{\mu\nu}^{{\rm G},m=1}$

  * •

in all remaining cases choose the robust “parachute method”, that is,
polynomial order 0 and $\tilde{T}_{\mu\nu}^{{\rm G},m=0}$.

The additional conditions on the number of contributing particles have been
introduced to avoid the inversion of poorly conditioned matrices. In Figure 4
we illustrate how well this procedure works. The top plot shows a cut of the
density of particles at t = 2.95 ms for the case of an equal mass binary
system with two 1.5 M⊙ neutron stars with the APR3 equation of state. The
bottom plot shows the resulting mapped $tt$-component of the stress energy
tensor in the $xy$-plane of the grid. As can be seen, all the features, that
are visible in the particle density profile, are clearly reproduced in the
mapped stress-energy tensor.

###  3.5 Stably simulating the formation of black holes

The remnant of a binary NS merger can, depending on the EOS, the total mass
and spin (and further processes which are not modelled here), undergo a
collapse to a black hole (BH). This can happen either “promptly” on the
dynamical timescale of the remnant (typically $\sim 1~{}$ms), or it can be
“delayed” for several dynamical timescales or, for binaries at the low-mass
end, it may not occur at all. If a BH does form, extra care is required in
order to avoid numerical problems.
The first potential problem we noticed when we initially attempted to simulate
a collapse was simply that at some point in time, the particles become packed
so close together, that the grid resolution is insufficient to evolve the
metric accurately enough to maintain a physically valid solution. This can
result in failures in the recovery of the primitive variables.
To cure this problem, we allow for the addition of more refinement levels.
After some experiments, we decided to use the ratio of the hydrodynamic and
the BSSN time step as an indicator of when it is time to add another
refinement level. As the particles move closer together, their Courant time
step, $\Delta t_{\rm SPH}=\xi_{\rm SPH}(h/c_{s})$, decreases (here $c_{s}$ is
the speed of sound), whereas the Courant time step for the mesh, $\Delta
t_{\rm BSSN}=\xi_{\rm G}(\Delta x/c)$, stays constant. Note that we have
omitted particle and grid labels for readability, and for clarity we have
written the expression using the speed of light $c$ explicitly. $\Delta x$ is
the resolution on the finest grid. For the dimensionless pre-factors, our
default values are $\xi_{\rm SPH}=0.2$ and $\xi_{\rm G}=0.35$. Whenever the
ratio of the two time steps grows beyond a threshold,

$\frac{\Delta t_{\mathrm{BSSN}}}{\Delta t_{\mathrm{SPH}}}>C_{\mathrm{refine}}$
(61)

we add a refinement level. In numerical experiments, we found that
$C_{\mathrm{refine}}=2.5$ works well.
Our current mesh refinement hierarchy is very simple: our finer grids are
always half the size of the next coarser grid, and they are centered at the
origin. We follow this strategy also when we create a new refinement level:
the new level has the same number of points and half the size of the
previously finest grid, and it is again centered at the origin. The data on
the new grid are calculated from the data on the previously finest grid using
the same interpolation operators that we use in the prolongation operators to
fill the ghost cells of the refinement levels, i.e. via cubic Lagrange
interpolation.
After adding a new refinement level, $\Delta t_{\mathrm{BSSN}}$ is half its
previous value and the time step ratio will again be less than
$C_{\mathrm{refine}}$. As the collapse proceeds, $\Delta t_{\mathrm{SPH}}$
will continue to decrease and may eventually trigger the addition of another
refinement level. This can in principle continue indefinitely, but would
eventually slow the simulation to a halt due to very small time steps.
Therefore, at some point in time, we start to remove particles in the
innermost, collapsing core.
The best criterion would of course be to remove particles when they are deep
enough inside the forming black hole. Since we have not yet implemented an
apparent horizon finder, we do not know precisely when and where the black
hole forms at run time. Instead, we rely on the value of the lapse, $\alpha$,
at the location of the particles as a proxy. With the slicing and shift
conditions used in the code, it is observed that the value of $\alpha$ at the
horizon of a single static black hole, evolved to numerical stationarity, is
around 0.3, see e.g. Figs. 14 and 16 in Rosswog and Diener [2021]. The value
of $\alpha$ at the actual horizon during the dynamical phase of the collapse
will of course vary slightly, but will not differ too much from the value of
0.3. Thus, we remove particles when they enter regions that have a
substantially lower lapse than this threshold value.
Based on the turduckening idea of Brown et al. [2009], we should be able to
safely change the interior of a black hole as long as it is done sufficiently
deep inside. In particular removing the source (the particles) of the stress
energy tensor, should not affect how the continuing collapse is seen from the
outside. To be safe we want to wait as long as possible before starting to
remove particles, but as the particles pile up inside the black hole the
Courant time step decreases dramatically, potentially making the collapse
process very computationally expensive. However, we are saved by the
observation that particles are essentially in free fall when they are that
deep inside a black hole. Hence, the energy momentum tensor is completely
dominated by the rest mass and velocity with the pressure and internal energy
only providing negligible small corrections. Therefore we can convert
particles into “dust” by setting their pressure and internal energy to zero
once the lapse at their position is less than $\alpha_{\mathrm{dust}}=0.05$.
The dust particles will no longer affect the evolution of their neighbours and
will no longer contribute to $\Delta t_{\mathrm{SPH}}$. They will simply
evolve along geodesics until the lapse at their positions falls below
$\alpha_{\mathrm{cut}}=0.02$ at which point we simply remove them completely
from the simulation. With this two stage process, converting particles first
to dust and then later removing them, we manage to have the particles
contribute to the stress energy tensor for significantly longer without
adversely affecting the time step. We found that removing particles as soon as
the lapse dropped below $\alpha=0.05$ could lead to a delay in the collapse of
the lapse in the center of the black hole.
Post-processing our data using the apparent horizon finder from the Einstein
Toolkit Löffler et al. [2012], showed that this delay in the lapse did not
significantly affect the horizon properties, but we still prefer to avoid it.
It turns out, that even with extra refinement levels, once we start converting
particles to dust and later removing particles, it is still possible to
eventually get failures in the recovery of primitive variables. However, this
only happens when the particles are essentially in free fall and the solution
is to convert them to dust before they reach $\alpha_{\mathrm{dust}}$. This is
only necessary in very few cases, if at all.
Once particles have been converted to dust, we still have to recover the
primitive variables from the evolved variables, but this is very simple. The
relations between the evolved and primitive variables, Eqs.(4), (6) and (11),
for a dust particle with vanishing $u$ and $P$ reduce to (omitting for
simplicity the particle label)

$\displaystyle N$ $\displaystyle=\sqrt{-g}\,\Theta\,n$ (62) $\displaystyle
S_{i}$ $\displaystyle=\Theta\,v_{i}$ (63) $\displaystyle e$
$\displaystyle=\Theta\,v_{i}v^{i}+\frac{1}{\Theta}.$ (64)

That is, $S_{i}$ reduces to the spatial part of the covariant 4-velocity,
$U_{\mu}$. Therefore, we can find $U_{0}$ so that the 4-velocity is properly
normalized, $U_{\mu}U^{\mu}=-1$. Raising the index on the covariant
4-velocity, we can simply read off $\Theta=U^{0}$ and then find $n$. In case a
fluid particle is transformed to dust, we also adjust $e$ so that it is
consistent with the values of $v^{i}$ and $\Theta$. In summary, the recovery
of primitive from evolved variables is always possible for dust particles in a
straightforward way, and the particles can keep evolving, contributing to the
stress–energy tensor, but do not affect any of their neighbour particles in
any negative way, until they can be safely removed once their lapse value has
dropped below the removal threshold.

##  4 Improvements to the initial data setup

In this section, we describe recent improvements in setting up the SPH
particles in the initial data (ID) code SPHINCS_ID SPHINCS_ID [2023]. It can
now also be linked to our fork of FUKA Papenfort et al. [2021]; Frankfurt
University/Kadath Initial Data solver [2023]—extended to comply with our
needs—to produce BSSN and SPH ID for neutron star binaries. The FUKA codes are
built on an extended version of the KADATH library Grandclément [2010]. In
this section, we refer generically to LORENE Gourgoulhon et al. [2001];
Grandclément et al. [2001]; LORENE [2001] and FUKA with the term “ID solver,”
when the discussion applies to both solvers.

###  4.1 Modeling neutron stars with the Artificial Pressure Method

As described in [Diener et al., 2022, Section 2.2.2], the initial neutron
stars are modelled by placing the SPH particles according to the “Artificial
Pressure Method” (APM) which uses the solutions found by the ID solver. We
briefly summarize the original method below, and refer the reader to Rosswog
[2020b]; Rosswog and Diener [2021] and [Diener et al., 2022, Section 2.2.2]
for more details, before we describe an additional improvement.
First, particles are placed according to a freely specified geometry (lattice,
spherical surfaces, etc.), and then each particle $a$ is assigned the same
baryon number $\nu_{0}=\nu_{\rm tot}/N_{\rm part}$, where $\nu_{\rm tot}$ is
the total baryon number for the star and $N_{\rm part}$ the number of SPH
particles used to model it. Subsequently an ”artificial pressure” is defined
as

$\displaystyle\pi_{a}={\rm
max}\left(1+\dfrac{N_{a}-N_{\mathrm{ID}}(\vec{r}_{a})}{N_{\mathrm{ID}}(\vec{r}_{a})},0.1\right).$
(65)

Here $N_{a}$ is the SPH estimate of the density variable defined in Equation
(4) on particle $a$ and $N_{\mathrm{ID}}(\vec{r}_{a})$ is the result from the
ID solver. The lower bound of 0.1 is imposed to avoid negative values. The
major goal of the original APM is to minimize the difference between $N_{a}$
and $N_{\mathrm{ID}}(\vec{r}_{a})$ while using SPH particles of the same
baryon number.222Another goal of the APM is to produce a locally isotropic
particle distribution Rosswog [2020b], that is, a distribution such that for
every particle, there exist a small enough neighborhood around it such that
the particle distribution inside such neighborhood is isotropic. Therefore, at
each iteration of the APM, the particle positions are updated in order to
achieve vanishing (artificial) pressure gradients. The corresponding position
update reads:

$\displaystyle{\delta\vec{r}_{a}}^{\rm
APM}=-\frac{1}{2}h_{a}^{2}\,\nu_{a}\sum_{b}\frac{\pi_{a}+\pi_{b}}{N_{b}}\,\nabla_{a}W_{ab}(h_{a}),$
(66)

see Section 3.1 for the meaning of the involved quantities. The iteration
stops when the differences between $N_{a}$ and $N_{\mathrm{ID}}(\vec{r}_{a})$
do not change significantly anymore.333Currently, we exit the APM iteration if
the baryon number ratio $\nu_{\mathrm{max}}/\nu_{\mathrm{min}}$ does not
change more than $0.25\%$ for 300 iterations. The maximum and minimum are
taken over all the particles.

While we want to construct initial conditions with densities $N_{a}$ as close
as possible to $N_{\mathrm{ID}}(\vec{r}_{a})$, it is in the end the (physical)
pressure gradients that, apart from gravity, drive the physical fluid motion.
Therefore, it may be advantageous to construct the artificial pressure,
$\pi_{a}$, from the physical pressures rather than the densities as in
Equation (65) For a short motivation as to why to use the pressure, we will
briefly switch to a Newtonian description (the GR case with our conventions
follows in a straight-forward way) and we define, for a general EOS, the
quantity

$\bar{\gamma}\equiv\left(\frac{\partial
P}{\partial\rho}\right)_{s}\frac{\rho}{P},$ (67)

which, for the special case of a polytropic EOS, simply reduces to the
polytropic exponent $\gamma$. If we assume that we have found a numerical
solution for the density that is $\rho=\rho_{0}+\delta\rho$, where $\rho_{0}$
is the true solution, the resulting pressure is

$P(\rho_{0}+\delta\rho)\approx P(\rho_{0})+\left(\frac{\partial
P}{\partial\rho}\right)_{s}\delta\rho.$ (68)

With $P_{0}=P(\rho_{0})$, the relative error $\epsilon_{P}$ in the pressure
reads

$\epsilon_{P}=\frac{P(\rho_{0}+\delta\rho)-P_{0}}{P_{0}}=\frac{(\partial
P/\partial\rho)_{s}\delta\rho}{P_{0}}=\bar{\gamma}\frac{\delta\rho}{\rho}=\bar{\gamma}\;\epsilon_{\rho}.$
(69)

Neutron star equations of state can be approximated by piecewise polytropes
[Read et al., 2009], which even at the lowest density piece (as is the case
for all other equations of state that we are aware of) have a polytropic
exponent value of $\gamma>1.3$. The higher density pieces have substantially
larger values, often larger than $3$. Thus, we expect that for all physically
relevant EOSs (not just polytropes) the relative error in $P$ will be larger
than in $\rho$. This motivates us to change the definition of the artificial
pressure from (65) to

$\displaystyle\pi_{a}={\rm
max}\left(1+\dfrac{p(\rho_{a})-p(\rho_{\mathrm{ID}}(\vec{r}_{a}))}{p(\rho_{\mathrm{ID}}(\vec{r}_{a}))},0.1\right),$
(70)

so that the APM iteration minimizes the error on the pressure directly. In
Figure 5 we see a comparison between the errors on the pressure when using the
definitions (65) and (70). With definitions (70) the errors decrease in the
inner 93% of stellar radius, but do not improve significantly in the outermost
layers. This is the case because at finite numerical resolution, the extremely
steep (physical) gradients in the stellar surface cannot be resolved. However,
these outer layers only constitute a very small fraction of the baryonic mass
of the star, about $0.01\%$ for the star in Figure 5, and are therefore not a
matter of concern here.

Figure 5: Relative errors on the pressure
$p_{\mathrm{err}}\coloneqq|p(\rho_{a})-p(\rho_{\mathrm{ID}}(\vec{r}_{a}))|/p(\rho_{\mathrm{ID}}(\vec{r}_{a}))$
at the end of an APM iteration that uses definition (65) (black) and another
APM iteration that uses definition (70) (red). The errors are evaluated on
$2.5\times 10^{6}$ particles modeling the same neutron star, which belongs to
a $1.9M_{\odot}$ equal-mass, $47.5$ km separation, MS1b, irrotational BNS
produced with LORENE. (Top) $p_{\mathrm{err}}(r)$, with $r$ radial coordinate
of each particle measured from the barycenter of the star. The left panel
shows the data for $r\in[0,9.5]$ km; the right panel for $r\in[9.5,11]$ km.
(Bottom) $p_{\mathrm{err}}$ binned in log scale; each bin is one order of
magnitude wide. $m$ is the sum of the baryonic masses of the particles in each
bin, in units of the baryonic mass of the star
$M_{\mathrm{b}}=2.016M_{\odot}$. Except for the outermost particles with
$r\gtrsim 10$ km $\simeq 93\%R$ ($R$ being the larger radius of the star),
whose errors do not decrease, the errors for the particles in the bulk of the
star decrease using definition (70). For example, $\sim 10\%$ of the baryonic
mass of the star improves the error by roughly one order of magnitude, moving
from bin $[10^{-2},10^{-1}]$ to bin $[10^{-3},10^{-2}]$.

In summary, the errors do not reduce very significantly, but the change in the
definition of the artificial pressure is nevertheless a welcome enhancement
that complements the other improvements to SPHINCS_ID and SPHINCS_BSSN_v1.0.
The SPHINCS_ID code lets the user choose which definition to use, (65) or
(70). The ID used for the runs described in this paper were produced using
(70).

###  4.2 The initial particle distribution on surface-conforming ovals

Until recently we have prepared the initial condition for the APM on each star
by placing particles on spherical surfaces with radii in the interval $(0,R)$,
with $R$ being the larger radius of the star, see [Diener et al., 2022,
Section 2.2.2]. This algorithm places close-to-equal-mass particles on each
spherical surface, taking into account the mass of the spherical shell bounded
by a spherical surface and the next. Therefore, some information on the
density profile of the star is considered already at the level of the initial
condition for the APM iteration. We have improved our algorithm by placing the
initial particles on ovals that conform to the (scaled) surface of the star,
otherwise we follow the same steps as described in [Diener et al., 2022,
Appendix B.1], more details can be found in Appendix A. A comparison between
particle distributions, produced using ovals and ellipsoids, is shown in
Figure 6 for a star whose geometry deviates significantly from spherical. The
placement of particles on ovals improves the accuracy of the particle model of
the outer layers of the star. This is because the APM starts with initial
conditions having a smoother particle distribution on the outer layers,
compared to the distributions obtained using spheres or ellipsoids. The latter
include cuts on the outer layers when the spheres or ellipsoids cut through
the surface of the star. This improvement is even more important for
rotationally flattened stars.

Figure 6: Cuts along the $xy$ plane, with $|z|<\,\sim\\!\\!0.89$ km, of
particle distributions with $10^{6}$ particles modeling the same star. The
left panel shows particles placed on surface-conforming ovals; the right panel
shows particles placed on ellipsoids. In both cases, a small random
displacement is applied, since it leads to a better initial condition for the
APM iteration, as already noted in Diener et al. [2022]. The star belongs to a
1.9 $M_{\odot}$ equal-mass, $47.5$ km separation, MS1b, irrotational BNS
produced with LORENE. The blue points have $|z|<\,\sim\\!\\!0.16$ km and trace
the surface of the star (for more details, see Appendix A). The particles
placed on surface-conforming ovals model very accurately the shape of the star
and produce a smooth surface, thus providing a more accurate initial condition
for the APM.

###  4.3 The boundary particles used during the APM

The APM iteration uses “ghost” or “boundary” particles outside the star that
prevent the particles modelling the star from being pushed outside the stellar
surface, see [Diener et al., 2022, Section 2.2.2] and [Diener et al., 2022,
Appendix B.2]. At each step of the iteration, the boundary particles are
assigned an artificial pressure that increases linearly with their distance
from the center of the star. The boundary particles closest to the star are
assigned an artificial pressure equal to $3\pi_{\mathrm{max}}$, those farthest
away a value of $30\pi_{\mathrm{max}}$ and for those in between, the
artificial pressure varies linearly between these bounds, where
$\pi_{\mathrm{max}}$ be the maximum artificial pressure of the real particles.
For these bounds we empirically found the best APM results. This artificial
pressure gradient makes the real particles feel a stronger repulsive force,
the closer they approach the stellar surface. Previously [Diener et al.,
2022], we placed the boundary particles on a lattice, between two ellipsoidal
surfaces, now we place them on a lattice between two surface-conforming ovals
instead. This, analogously to the initial placement of real particles on
surface-conforming ovals described in Section 4.2, makes it easier for the APM
iteration to model the outer layers and the overall geometry of the star.
In [Diener et al., 2022, Appendix B.2], the parameter $\delta$ was introduced,
which is the distance between the surface of the star and the boundary
particle closest to it, along the direction of the star’s largest radius. The
modeling of the outer layers turns out to somewhat depend on the value of
$\delta$, and it would be desirable to remove this dependence. If $\delta$ is
too small, the real particles cannot approach the surface of the star; if
$\delta$ is too large, the particle distribution on the outer layers can
become non-smooth, leaving a few isolated particles outside the otherwise
smooth stellar surface. In order to reduce the dependence on $\delta$, we set
a small value of $\delta$ initially, and then let the boundary particles move
outwards (effectively increasing $\delta$) very slowly until the condition
$|r_{\mathrm{av}}-R|<h_{\mathrm{av}}/3$ is met, where $r_{\mathrm{av}}$ and
$h_{\mathrm{av}}$ are the average values of radius and smoothing lengths of
the particles in the outer layers and $R$ is largest radius of the star. The
particles modeling the outer layers are defined as those having a radial
coordinate (measured from the center of the star) larger than $99.5\%$ of $R$.
In case these should be less than 10 particles, this fraction is reduced in
steps of $0.5\%$ until this number is exceeded. This algorithm starts out
producing a smooth particle distribution on the outer layers since $\delta$ is
initially small, and preserves the smoothness when it gently allows the real
particles to move towards the surface. In addition, it reduces the dependence
on the initial value of $\delta$ since the latter is increased during the
iteration.
A comparison between a particle distribution obtained using the latest methods
described in this section and in Section 4.2, and a particle distribution
obtained with older methods, is shown in Figure 7. Note how the model of the
geometry of the star and its outer layers have improved with the new methods.

Figure 7: Projections of particle distributions with $2.5\times 10^{6}$
particles modeling the same neutron star. Real particles are in black,
boundary particles with $|z|<\,\sim\\!\\!0.74$ km are in red, and the points
lying on the surface of the star with $|z|<\,\sim\\!0.16$ km are in blue. The
star belongs to a $1.9M_{\odot}$ equal-mass, $47.5$ km separation, MS1b,
irrotational BNS produced with LORENE. (Top-left) The real and boundary
particles are placed on surface-conforming ovals, see Section 4.2 and Section
4.3. (Top-right) The real and boundary particles are placed on ellipsoids.
(Bottom-left) Particles placed with our latest APM algorithm, with the
particles in the top-left panel as initial condition. (Bottom-right) Particles
placed with our older APM algorithm, with the particles in the top-right panel
as initial condition. See main text for details. The particles in the bottom-
left panel better model the outer layers and the geometry of the star,
compared to those in the bottom-right panel.

###  4.4 The ADM linear momentum for the initial data

When considering ID produced with LORENE and FUKA, we can simplify the SPH
estimate of the ADM momentum of the fluid, Equation (17). The ID solvers
assume asymptotic flatness, conformal flatness and maximal slicing on the
initial spacelike hypersurface Gourgoulhon et al. [2001]; Papenfort et al.
[2021]. In coorbiting coordinates of Cartesian type, the conformal flatness
condition can be written as [Gourgoulhon et al., 2001, Section
IV.A],[Papenfort et al., 2021, Equation (6)]:

$\displaystyle\gamma_{ij}=A^{2}\delta_{ij},$ (71)

with $A$ being the conformal factor and $\delta_{ij}=\mathrm{diag}(1,1,1)$ the
Euclidean metric. Hence, a global, orthogonal, non-orthonormal frame
$e^{i}_{(j)}=\delta^{i}{}_{j}$ exists on the initial spacelike hypersurface.
This frame also becomes orthonormal, hence Cartesian, at spatial infinity
where $A\rightarrow 1$ due to the asymptotic flatness condition. Therefore, we
can set $\xi^{i}=e^{i}_{(j)}$ in Equation (16) to compute the $j^{\rm th}$
Cartesian component of the ADM linear momentum. Doing so, and imposing the
maximal slicing condition, the part of the ADM linear momentum determined by
the spacetime—i.e., the second term in the squared parenthesis in Equation
(16)—is zero. We show this explicitly in Appendix C.
Hence, for the ID the total ADM linear momentum is equal to the ADM momentum
of the fluid, which can be estimated with (17) for the LORENE and FUKA ID. In
addition, we can compare this estimate with the one obtained using (15). It is
possible to compute the linear ADM momentum within LORENE as a surface
integral at infinity.444We added this feature to our fork of LORENE, since we
could not find a function that does it in the original fork. However, a
function that computes the Bowen–York angular momentum as described in
[Gourgoulhon et al., 2001, Section IIID], was included in the original fork
LORENE Reference guide. Definition of Lorene::Binaire::angu_mom [2023], and we
used it as a template. LORENE can easily handle this computation thanks to its
compactified coordinates. FUKA also provides an estimate of the ADM momentum.
Hence, for each ID, we have two independent estimates of the ADM linear
momentum: one computed by the ID solver as a surface integral using (15), and
the other computed by SPHINCS_ID as an SPH estimate of a volume integral using
(17). For the irrotational, equal-mass $1.3M_{\odot}$ LORENE ID with 2 million
SPH particles, see Section (5.1), the two estimates are:

$\displaystyle\texttt{LORENE}:$ $\displaystyle\quad P^{\mathrm{ADM}}_{i}$
$\displaystyle\simeq\big{(}<10^{-15}$ $\displaystyle,\,6.336\times 10^{-11}$
$\displaystyle,\,<10^{-15}\big{)},$ (72a)
$\displaystyle\texttt{SPHINCS\\_ID}:$ $\displaystyle\quad
P^{\mathrm{ADM}}_{i}$ $\displaystyle\simeq\big{(}-1.393\times 10^{-7}$
$\displaystyle,\,6.318\times 10^{-11}$ $\displaystyle,\,<10^{-15}\big{)}.$
(72b)

Only the $x$ component is affected by a substantial error after the particles
are placed, but it stays very small nonetheless. We obtain similar results for
the other equal-mass, non-spinning systems we consider in Section 5.2. For the
equal-mass $1.3M_{\odot}$ FUKA ID used in the simulation described in Section
5.3, where one star is spinning with $\chi\simeq 0.5$, the estimates are:

$\displaystyle\texttt{FUKA}:$ $\displaystyle\quad P^{\mathrm{ADM}}_{i}$
$\displaystyle\simeq\big{(}-2.804\times 10^{-13}$ ,
$\displaystyle\,\,1.683\times 10^{-10}$ , $\displaystyle\,<10^{-15}\big{)},$
(73a) $\displaystyle\texttt{SPHINCS\\_ID}:$ $\displaystyle\quad
P^{\mathrm{ADM}}_{i}$ $\displaystyle\simeq\big{(}\hphantom{lin}6.030\times
10^{-6}$ , $\displaystyle\,\,2.809\times 10^{-4}$ ,
$\displaystyle\,<10^{-15}\big{)}.$ (73b)

It makes sense that the SPH estimate is much better for the non-spinning
systems since the particles modeling the second star are placed mirroring
those modeling the first star, with respect to the $yz$ plane [Diener et al.,
2022, Section 2.2]. Triggered by questions from the referee, we realized that
it would be even better to reflect the particle positions of star one through
the origin to set up star two. In this case particles would be placed with
complete symmetry for all components of the coordinates (including $x$) and
the helical symmetry present in the initial data would guarantee that the
velocity of a particle in star two is exactly opposite the corresponding
particle in star one. There would then be exact cancellation (to roundoff
error) of their contribution to the momentum, allowing us to reduce the small
initial value for the $x$-component of the ADM momentum for irrotational,
equal mass systems. This procedure will be used in future particle setups. In
the system with a spinning star, however, neither mirror nor reflection
symmetry can be enforced, and the terms in the sum (17) do not compensate to
the same degree.

##  5 Numerical Results for Neutron star mergers

Here we show astrophysical examples of neutron star mergers with
SPHINCS_BSSN_v1.0. Standard hydrodynamics tests such as shock tube tests are
not impacted by any of the new elements introduced here, therefore, we refer
the interested reader to our previous papers [Rosswog and Diener, 2021;
Rosswog et al., 2022]. In Section 5.1 we show a binary neutron star merger
where a remnant survives (for at least several dynamical time scales), Section
5.2 shows an example where the merger remnant promptly collapses to form a
black hole and in Section 5.3 we show results for a binary system where only
one of the neutron stars has a large spin, whereas the other has none. All
systems are of equal mass, the simulations start from an initial separation of
45 km and are performed with slightly more than 2 million SPH particles
(except for the case where collapse to a black hole happens promptly; here 1
million SPH particles are used), the APR3 EOS, initially 7 mesh refinement
levels out to $\approx 2268$ km in each coordinate direction and a minimum
initial grid spacing of $\Delta x=369$ m. Keep in mind that new refinement
levels are added dynamically, when the criterion described in Section 3.5 is
met. In the run presented in Section 5.2 where there is a collapse to a black
hole, the number of refinement levels dynamically increases up to 11.

###  5.1 Neutron star merger with surviving remnant

We show here the merger of two 1.3 M⊙ neutron stars, with initial conditions
produced by LORENE. After a few orbits of inspiral, the stars merge and remain
initially close to perfect symmetry, see panel 1 in Figure 8. The strong shear
at the interface between the stars is Kelvin–Helmholtz unstable and as matter
from this region is sprayed out, deviations from perfect symmetries emerge
(panel 2), as also frequently seen in Eulerian neutron star merger
simulations. A few milliseconds later, the remnant settles into what seems a
stationary state with a bar-like central object shedding mass via a spiral-
wave into the surrounding torus. This spiral wave ejection channel might have
played an important role in the early blue kilonova signal after the first
observed neutron star merger GW170817 [Nedora et al., 2019], see Rosswog and
Korobkin [2022] for a review.
In the left panel of Figure 9 we show the evolution of the maximum density
(red curve, right axis) together with the minimum lapse function value (black
curve, left axis). In an initial very deep compression the density reaches a
value close to $9.4\times 10^{14}$ g cm-3, then the remnant bounces back and,
after several more oscillations, the peak density settles near a value of
$9.5\times 10^{14}$ g cm-3. As expected, the lapse is lowest where the density
is highest and vice versa. The right panel shows the value of the maximum GW
amplitude times the distance to the observer, as calculated via the quadrupole
approximation (orange) and as extracted from the spacetime via the Weyl scalar
$\psi_{4}$ (black), how these are calculated in detail can be found in
Appendix A of Diener et al. [2022]. All $\Psi_{4}$ based GW waveforms were
analyzed using kuibit Bozzola [2021]. Again, we find a rather good agreement
of the quadrupole result with the more sophisticated $\psi_{4}$ method. The
radiated energy (red curve, left axis) and angular momentum (black curve,
right axis) are plotted as a percentage of the initial ADM values in the left
panel of Fig 16. More than 2% of the initial ADM mass and more than 20% of the
initial ADM angular momentum are radiated.

Figure 8: Density distribution in the orbital plane of an irrotational binary
system with 2 $\times 1.3$ M⊙ with the APR3 EOS.

Figure 9: (Left) Maximum mass density (red curve, right axis) together with
the minimum value of the lapse function $\alpha$ (black curve, left axis).
(Right) Maximum gravitational wave amplitude extracted via the Weyl scalar
$\Psi_{4}$ (black curve) and the quadrupole formula (yellow curve). Both
panels refer to the simulation of an irrotational 2 $\times$ 1.3 M⊙ neutron
star binary shown in Figure 8. For convenience with comparison with other
plots, $t=0$ correspond to the time of the maximum amplitude in the
gravitational waveform.

###  5.2 Neutron star merger with black hole formation

We also show the merger of two 1.5 M⊙ neutron stars, with initial conditions
produced by LORENE. Only 1 million particles are used here (runs with higher
resolution become very slow due to the timestep requirements) with a
corresponding initial grid spacing of $\Delta x=499$ m. In this case, the
merged object is massive enough that it undergoes a prompt collapse. During
the collapse additional grid refinement levels are added when needed and at
the end we have a total of 11 refinement levels with a finest grid resolution
of $\Delta x=32$ m.
In Figure 10 we show 3 snapshots of the equatorial density. The first, at
$t=3.28$ ms, is from well before particles start to be removed and the density
is still increasing. At $t=4.29$ ms, about 90% of the particles have already
been removed, but the maximum density of the remaining particles is still
close to the initial central density of the stars. At $5.29$ ms, matter has
been drained down to $4\times 10^{-3}$ M⊙. Only about $6\times 10^{-4}$ M⊙ of
this material is unbound from the black hole.
In the left panel of Figure11 we show the evolution of the maximal density
(red curve, right axis) and the minimum lapse (black curve, left axis). In
this plot, particles that have been converted to dust, do not count towards
the maximum density and minimum lapse. The rapid drop in the maximum density
is completely due to the conversion of particles to dust and their eventual
removal at lapse values below 0.02.
In the right panel of Figure11 we show a comparison between the maximum GW
amplitude times the distance to the observer as extracted from the quadrupole
formula and from $\Psi_{4}$. As expected, the quadrupole waveform shuts off
too early as it is sourced by the matter motion and does not know about the
quasinormal ringdown of the spacetime itself. In this simulation about
$3.8\times 10^{-2}$ M⊙ of energy and $1.18$ M⊙2 of angular momentum is
radiated away by GWs. By analyzing the properties of the final horizon using
the tools from the Einstein Toolkit (Löffler et al. [2012]), we find that the
black hole that formed has a dimensionless spin parameter of about $a/M\approx
0.8$, consistent with earlier findings that binary neutron star mergers leads
to faster spinning black holes than binary black hole mergers (see e.g.
Baiotti et al. [2008]). These last numbers should be taken by a grain of salt,
as the finite resolution may still have an impact. A more detailed analysis is
left for future work.

Figure 10: Density distribution in the orbital plane of an irrotational binary
system with 2 $\times 1.5$ M⊙ with the APR3 EOS.

Figure 11: (Left) Maximum mass density (red curve, right axis) together with
the minimum value of the lapse function $\alpha$ (black curve, left axis).
(Right) Maximum gravitational wave amplitude extracted via the Weyl scalar
$\Psi_{4}$ (black curve) and the quadrupole formula (yellow curve). Both
panels refer to the simulation of an irrotational 2 $\times$ 1.5 M⊙ neutron
star binary. For convenience with comparison with other plots, $t=0$
correspond to the time of the maximum amplitude in the gravitational waveform.

###  5.3 Neutron star merger with a single spinning star

Most commonly, irrotational binary systems are studied, and they are
considered as most realistic since dissipative effects cannot spin up neutron
stars to substantial spin values [Bildsten and Cutler, 1992; Kochanek, 1992]
and, at merger, any residual stellar spin is likely small compared to the huge
orbital angular momentum. Nature, however, likely can produce neutron star
binary systems in several ways [Tauris and van den Heuvel, 2023] and, likely
at smaller rates, more extreme systems may be produced. As one such example,
we study here a binary system where only one of the two neutron stars is
rapidly spinning while the other is irrotational. Such systems have hardly
been explored before, we are only aware of one such study by Papenfort et al.
[2022] where authors study extreme mass ratio and spin spin configurations.
We evolve a binary system with $2\times 1.3$ M⊙ stars, where one of the stars
is spinning. The chosen value of the spin parameter, $\chi=0.5$, corresponds
to a spin period of about 1.2 ms. Since LORENE cannot construct such a case,
we use the FUKA library instead. As can be seen from Figure 12, this initial
data produces a matter distribution that is substantially different from the
case shown in Section 5.1. During merger, a massive tidal tail forms and our
evolution here is, qualitatively, similar to panel 1 in Figure 1 of Papenfort
et al. [2022]. (Note however that their system has a different spin value, a
different EOS and a different mass.) In panels two and three of Figure 12 one
sees how the rapidly spinning central remnant is punching shock waves into the
remnant. This strong shock compression in the torus drives polar outflows at
$\sim 40^{\circ}$ from the polar axis, see the volume rendering in the left
panel of Figure 13, with velocities up to $\sim 0.4$ c (right panel same
figure). The density compression at merger is much milder, see left panel of
Figure 14, and the post-merger gravitational wave amplitudes (right panel) are
substantially lower than in the non-spinning case. This effect is also
reflected in the Fourier power spectra shown in Figure 15. Here the non-
spinning case (with higher power) is shown in blue and the spinning case in
orange. It can also be seen that the frequency of the main peak after merger
shifts to lower frequency in the spinning case compared to the non-spinning
case, consistent with the less compact remnant.
As a quick test, we can compare the peak frequency with the prediction of the
empirical quasi-universal relation given by Equation (4) in Vretinaris et al.
[2020]. This relation is

$f_{\mathrm{peak}}/M_{\mathrm{chirp}}=13.822-0.576M_{\mathrm{chirp}}-1.375R_{1.6}+0.479M_{\mathrm{chirp}}^{2}-0.073R_{1.6}M_{\mathrm{chirp}}+0.044R_{1.6}^{2},$
(74)

where $R_{1.6}$ is the circumferential radius of a 1.6M⊙ star with the given
EOS and $M_{\mathrm{chirp}}=(m_{1}m_{2})^{3/5}/(m_{1}+m_{2})^{1/5}$ is the
chirp mass of the binary. Using LORENE we can calculate $R_{1.6}$ for a star
with the APR3 equation of state to be $R_{1.6}=11.75$ km. With
$m_{1}=m_{2}=1.3$ M⊙ we find a chirp mass of $M_{\mathrm{chirp}}=1.13$ M⊙.
Inserting these numbers into Equation (74) we find that the relation predicts
a value of $f_{\mathrm{peak}}=3.09$ kHz, in excellent agreement with the
location of the peak of the spectrum for the non-spinning case (blue curve).
Finally, in Figure 16 we compare the radiated energy (red curve, left axis)
and angular momentum (black curve, right axis) for the non-spinning (left
panel) and spinning (right panel) case. In both cases, the values are given as
a percentage of the initial ADM value and the solid line only includes the
contribution form $\ell=2$, whereas the dashed line includes the contribution
from all modes up to $\ell=4$. The non-spinning case radiates about twice as
much energy as the spinning case and does it predominantly in the $\ell=2$
modes. The spinning case does show a bit more contribution from the higher
$\ell$ modes, consistent with a more asymmetric merger.
Since our main aim here is to demonstrate that challenging astrophysical
problems can be robustly addressed by SPHINCS_BSSN_v1.0, we leave a further
discussion of the astrophysical implications of such problems to future
publications.

Figure 12: Density distribution in the orbital plane of a 2 $\times 1.3$ M⊙
binary with the APR3 EOS. One of the stars no spin, while the other has
$\chi\simeq 0.5$.

Figure 13: (Left) Volume rendering of the density distribution at $t=9.93$ ms
of the $2\times 1.3$ M⊙ merger with a single spinning star. The rapidly
spinning central object compresses the torus by shock waves which results in
polar bulk outflows with velocities reaching $\sim 0.4c$ (Right).

Figure 14: (Left) Maximum mass density (red curce, right axis) together with
the minimum value of the lapse function $\alpha$ (black curve, left axis.
(Right) $\ell=2,m=2$ mode of the $h_{+}$ polarization of the gravitational
wave strain extracted from the Weyl scalar $\Psi_{4}$. Both panels refer to
the simulation of a 2 $\times$ 1.3 M⊙ neutron star binary where one star has a
spin of $\chi\simeq 0.5$. For convenience with comparison with other plots
$t=0$ correspond to the time of the maximum amplitude in the gravitational
waveform.

Figure 15: Spectra of the $\ell=2$, $m=2$ gravitational waveforms for the non-
spinning case (blue curve) and the case with one spinning star (orange curve).

Figure 16: The radiated energy (red curves) and angular momentum (black
curves) as function of time with $t=0$ corresponding to the peak amplitude of
the gravitational waveform. Both are shown as percantages of the initial ADM
values. The left axis is for energy and the right axis for angular momentum.
The solid lines only includes the contribution from the $\ell=2$ modes, while
the dashed lines includes all modes up to $\ell=4$. The plot on the left is
for the non-spinning case while the plot on the right is for the case with one
spinning star.

##  6 Summary

In this work, we have presented version 1.0 of our Lagrangian Numerical
Relativity code SPHINCS_BSSN. Some of the methodological elements have been
published before [Rosswog and Diener, 2021; Diener et al., 2022; Rosswog et
al., 2022], others are introduced here for the first time.
First, a new way to map the stress–energy tensor $T_{\mu\nu}$ (known at the
particle positions) to our spacetime mesh is introduced. The new method sets
up a polynomial bases of a given order at each grid point and then computes
expansion coefficients that are optimal (for the given order) in the sense
that they minimize an error functional. We do this for polynomial orders from
0 to 3 and out of those possibilities we select the one that best represents
the surrounding particle values and meets some admissibility criteria. Our
procedure is described in detail in Section 3.4 and we show an instructive
example of the method in Appendix D.
Second, we have introduced measures that make the simulation of a collapse to
a black hole more robust. We realized that in some cases our originally chosen
spacetime evolution was not resolved well enough and we now add additional
refinement levels when the hydrodynamically allowed time step drops
substantially below the time step that is admissible for the space time
evolution. Once the lapse value at a particle position has dropped to a very
small value (here $\alpha_{\rm cut}=0.02$), we remove the particle to avoid
the time step shrinking towards zero. The lapse value $\alpha_{\rm cut}$ is
well below the value where an apparent horizon forms ($\sim 0.3$). While
evolving towards this very low lapse value, the recovery of the physical
variables from the numerical ones can fail. While this happens well inside the
horizon and thus should not affect the spacetime outside of it, we
nevertheless need to keep the particle evolution going until the threshold
lapse for removal is reached. to avoid this problem, we transform the
corresponding fluid particle into “dust” with vanishing pressure and internal
energy when the lapse at a particle drops below $a_{\mathrm{dust}}$. This
allows for a simple and robust recovery of the physical variables, and the
particle’s contribution to the stress–energy tensor is counted until it is
finally removed. For more details on the procedure, see Section 3.5.
The third improvement concerns the placement of the SPH particles modeling the
fluid at the level of the initial data, and is implemented in the code
SPHINCS_ID. This code can now use initial data produced with the FUKA library,
in addition to those produced with the LORENE library. In the latest version
of the code, the particles—both physical particles modeling the stars, and
boundary particles used in the “Artificial Pressure Method”—are placed so that
they model the geometry of the stars more accurately than before. This allows
for a better approximation of hydrodynamical equilibrium with SPH particles.
After their initial placement, the particles are iterated into optimal
positions according to a variant of the Artificial Pressure Method. In the
original version of this method, the relative error between the density
provided by the ID solver and the SPH estimate, was used to define an
“artificial pressure.” The latter’s gradient pushes the particles in positions
where they reduce the error on the density. In the latest version of this
method, we instead use the relative error of the physical pressure (rather
than the density) to compute the artificial pressure. This minimizes the error
on the physical pressure directly, and leads (with everything else being the
same) to lower errors in the physical pressure and thus to more accurate
initial data.
To illustrate the working and robustness of SPHINCS_BSSN_v1.0 we have
performed three simulations: one irrotational binary merger ($2\times 1.3$ M⊙)
that remains stable on the simulation timescale, one irrotational system
($2\times 1.5$ M⊙) that collapses “promptly” (i.e. without any bounce) and one
extreme binary system where only one of the stars has a (large) spin,
$\chi=0.5$. All these simulations use the APR3 equation of state, the first
two simulations are produced using LORENE, the latter using FUKA.
Not too surprisingly, for the stable irrotational case we find an anti-
correlation between the maximum density and the minimum lapse value, see
Figure 9, left panel. Concerning the GW emission, we have rather good
agreement between the quadrupole waveform (for more details see Rosswog et al.
[2022]) and the waveform extracted from the Weyl scalar $\psi_{4}$, see Figure
9, right panel. We show the agreement for the first two cases only, but it is
similarly good for the third case. Again expected, we find that the GW
emission is strongly dominated by the $l=2,m=2$ mode. We find that GWs carry
away about 2% of the initial ADM mass and about 20% of the initial ADM angular
momentum.
For the collapsing system, we find that very little mass $<6\times 10^{-4}$M⊙
escapes the fate of falling into the BH and that the final BH is spinning
fairly fast. The dimensionless spin parameter of $a/M\approx 0.8$ is
significantly larger than the end result of an irrotational binary BH merger
where $a/M\approx 0.68$.
Last, but not least, we performed a simulation of an extreme case with only
one rapidly spinning star that has been produced using the FUKA library. We
find that the neutron star spin has a very large impact on the merger
morphology. Similar to cases with extreme mass ratios, a single puffed up
tidal tail forms. Overall, the collision is less violent in the sense that the
high-density regions become not as much compressed as in the equal mass case
and the minimum lapse values remain larger. We also observe that less energy
and angular momentum are radiated by GWs in the post-merger phase, likely
because the central regions are less perturbed in the less violent collision
and thus deviate less from rotational symmetry.
Clearly, SPHINCS_BSSN_v1.0 would benefit from the inclusion of more
microphysics and its computational performance needs to be further improved.
These issues will be addressed in future work.

## Data availablity statement

The data underlying this article will be shared on reasonable request to the
corresponding author.

## Author Contributions

SR has developed the methods in and coded the fluid part of SPINCS_BSSN v1.0.
He has further designed various versions of mapping particle properties to the
mesh, the latest of which is described here in Section 3.4 and exemplified in
Appendix 4. All of this has happened in close coordination with PD. SR has
further written the first draft of this paper.
FT developed SPHINCS_ID, produced the initial data used in the simulations
using LORENE, FUKA and SPHINCS_ID, and performed the computation of the ADM
momentum of the fluid in SPH. FT wrote the part of Section 3.1 describing the
SPH estimate of the ADM momentum, Section 4, Appendix A, Appendix B, Appendix
C. All authors contributed to manuscript revision, read, and approved the
submitted version.
PD has developed the methods necessary for grid structures (including
refinement) and coded the interface to the needed routines from McLachlan from
the Einstein Toolkitand implemented them in SPINCS_BSSN v1.0. PD derived and
implemented the Hermite polynomial routines used for mapping metric
information from the grid to the particles and worked closely with SR on
developing and testing the methods to map the stress energy tensor from the
particles to the grid.

## Funding

SR has been supported by the Swedish Research Council (VR) under grant number
2020-05044, by the research environment grant “Gravitational Radiation and
Electromagnetic Astrophysical Transients” (GREAT) funded by the Swedish
Research Council (VR) under Dnr 2016-06012, by which also FT has been
supported, and by the Knut and Alice Wallenberg Foundation under grant Dnr.
KAW 2019.0112, by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy – EXC 2121 ”Quantum Universe”
– 390833306 and by the European Research Council (ERC) Advanced Grant
INSPIRATION under the European Union’s Horizon 2020 research and innovation
programme (Grant agreement No. 101053985). SR’s calculations were performed on
the facilities of the North-German Supercomputing Alliance (HLRN), on the
resources provided by the Swedish National Infrastructure for Computing (SNIC)
in Linköping, partially funded by the Swedish Research Council through Grant
Agreement no. 2016-07213, and at the SUNRISE HPC facility supported by the
Technical Division at the Department of Physics, Stockholm University.

## Acknowledgments

We thank Ian Hawke for insightful comments on our LRE method and we are very
grateful to Sam Tootle for his help with the FUKA library. We also thank the
referees for insightful and useful comments. Special thanks go to Holger
Motzkau and Mikica Kocic for their excellent support in upgrading and
maintaining SUNRISE.

## Conflict of interest

The authors declare that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest.

## Publisher’s note

All claims expressed in this article are solely those of the authors and do
not necessarily represent those of their affiliated organizations, or those of
the publisher, the editors and the reviewers. Any product that may be
evaluated in this article, or claim that may be made by its manufacturer, is
not guaranteed or endorsed by the publisher.

## References

  * Abbott et al. [2017a] Abbott, B. P., Abbott, R., Abbott, T. D., Acernese, F., Ackley, K., Adams, C., et al. (2017a).  A gravitational-wave standard siren measurement of the Hubble constant.  _Nature_ 551, 85–88.  10.1038/nature24471
  * Abbott et al. [2017b] Abbott, B. P., Abbott, R., Abbott, T. D., Acernese, F., Ackley, K., Adams, C., et al. (2017b).  Gravitational Waves and Gamma-Rays from a Binary Neutron Star Merger: GW170817 and GRB 170817A.  _ApJL_ 848, L13.  10.3847/2041-8213/aa920c
  * Abbott et al. [2017c] Abbott, B. P., Abbott, R., Abbott, T. D., Acernese, F., Ackley, K., Adams, C., et al. (2017c).  GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral.  _Physical Review Letters_ 119, 161101.  10.1103/PhysRevLett.119.161101
  * Abbott et al. [2017d] Abbott, B. P., Abbott, R., Abbott, T. D., Acernese, F., Ackley, K., Adams, C., et al. (2017d).  Multi-messenger Observations of a Binary Neutron Star Merger.  _ApJL_ 848, L12.  10.3847/2041-8213/aa91c9
  * Akmal et al. [1998] Akmal, A., Pandharipande, V. R., and Ravenhall, D. G. (1998).  Equation of state of nucleon matter and neutron star structure.  _Phys. Rev. C_ 58, 1804–1828.  10.1103/PhysRevC.58.1804
  * Alcubierre [2008] Alcubierre, M. (2008).  _Introduction to 3+1 Numerical Relativity_ (Oxford University Press)
  * Alsing et al. [2018] Alsing, J., Silva, H. O., and Berti, E. (2018).  Evidence for a maximum mass cut-off in the neutron star mass distribution and constraints on the equation of state.  _MNRAS_ 478, 1377–1391.  10.1093/mnras/sty1065
  * Antoniadis et al. [2016] Antoniadis, J., Tauris, T. M., Ozel, F., Barr, E., Champion, D. J., and Freire, P. C. C. (2016).  The millisecond pulsar mass distribution: Evidence for bimodality and constraints on the maximum neutron star mass.  _arXiv e-prints_ , arXiv:1605.0166510.48550/arXiv.1605.01665
  * Arcavi et al. [2017] Arcavi, I., Hosseinzadeh, G., Howell, D. A., McCully, C., Poznanski, D., Kasen, D., et al. (2017).  Optical emission from a kilonova following a gravitational-wave-detected neutron-star merger.  _Nature_ 551, 64–66.  10.1038/nature24291
  * Baiotti et al. [2008] Baiotti, L., Giacomazzo, B., and Rezzolla, L. (2008).  Accurate evolutions of inspiralling neutron-star binaries: Prompt and delayed collapse to a black hole.  _Phys. Rev. D_ 78, 084033.  10.1103/PhysRevD.78.084033
  * Baumgarte and Shapiro [1999] Baumgarte, T. W. and Shapiro, S. L. (1999).  Numerical integration of Einstein’s field equations.  _Phys. Rev. D_ 59, 024007.  10.1103/PhysRevD.59.024007
  * Baumgarte and Shapiro [2010] Baumgarte, T. W. and Shapiro, S. L. (2010).  _Numerical Relativity: Solving Einstein’s Equations on the Computer_ (Cambridge: Cambridge University Press)
  * Bauswein et al. [2017] Bauswein, A., Just, O., Janka, H.-T., and Stergioulas, N. (2017).  Neutron-star Radius Constraints from GW170817 and Future Detections.  _ApJL_ 850, L34.  10.3847/2041-8213/aa9994
  * Bildsten and Cutler [1992] Bildsten, L. and Cutler, C. (1992).  Tidal interactions of inspiraling compact binaries.  _ApJ_ 400, 175–180.  10.1086/171983
  * Biswas and Datta [2021] Biswas, B. and Datta, S. (2021).  Constraining neutron star properties with a new equation of state insensitive approach.  _arXiv e-prints_ , arXiv:2112.10824
  * Bozzola [2021] Bozzola, G. (2021).  kuibit: Analyzing Einstein Toolkit simulations with Python.  _The Journal of Open Source Software_ 6, 3099.  10.21105/joss.03099
  * Brown et al. [2009] Brown, J. D., Diener, P., Sarbach, O., Schnetter, E., and Tiglio, M. (2009).  Turduckening black holes: an analytical and computational study.  _Phys. Rev. D_ 79, 044023.  10.1103/PhysRevD.79.044023
  * Cabezon et al. [2008] Cabezon, R. M., Garcia-Senz, D., and Relano, A. (2008).  A one-parameter family of interpolating kernels for smoothed particle hydrodynamics studies.  _Journal of Computational Physics_ 227, 8523–8540.  10.1016/j.jcp.2008.06.014
  * Cottet et al. [2014] Cottet, G., Etancelin, J., Perignon, F., and Picard, C. (2014).  High order Semi-Lagrangian particle methods for tranport equations.  _ESAIM: Mathematical Modelling and Numerical Analysis, EDP Sciences_ 48, 1029 – 1060.  10.1086/185939
  * Cottet and Koumoutsakos [2000] Cottet, G.-H. and Koumoutsakos, P. D. (2000).  _Vortex Methods_
  * Cowperthwaite et al. [2017] Cowperthwaite, P. S., Berger, E., Villar, V. A., and Metzger, B. D. (2017).  The Electromagnetic Counterpart of the Binary Neutron Star Merger LIGO/Virgo GW170817. II. UV, Optical, and Near-infrared Light Curves and Comparison to Kilonova Models.  _ApJL_ 848, L17.  10.3847/2041-8213/aa8fc7
  * Dehnen and Aly [2012] Dehnen, W. and Aly, H. (2012).  Improving convergence in smoothed particle hydrodynamics simulations without pairing instability.  _MNRAS_ 425, 1068–1082.  10.1111/j.1365-2966.2012.21439.x
  * Diener et al. [2022] Diener, P., Rosswog, S., and Torsello, F. (2022).  Simulating neutron star mergers with the Lagrangian Numerical Relativity code SPHINCS_BSSN.  _European Physical Journal A_ 58, 74.  10.1140/epja/s10050-022-00725-7
  * Eichler et al. [1989] Eichler, D., Livio, M., Piran, T., and Schramm, D. N. (1989).  Nucleosynthesis, Neutrino Bursts and $\gamma$-Ray from Coalescing Neutron Stars.  _Nature_ 340, 126
  * Einstein Toolkit web page [2020] [Dataset] Einstein Toolkit web page (2020).  https://einsteintoolkit.org/.  [Online; accessed 9-December-2020]
  * Evans et al. [2017] Evans, P. A., Cenko, S. B., Kennea, J. A., Emery, S. W. K., Kuin, N. P. M., Korobkin, O., et al. (2017).  Swift and NuSTAR observations of GW170817: Detection of a blue kilonova.  _Science_ 358, 1565–1570.  10.1126/science.aap9580
  * Fong et al. [2017] Fong, W., Berger, E., Blanchard, P. K., Margutti, R., Cowperthwaite, P. S., Chornock, R., et al. (2017).  The Electromagnetic Counterpart of the Binary Neutron Star Merger LIGO/Virgo GW170817. VIII. A Comparison to Cosmological Short-duration Gamma-Ray Bursts.  _ApJ_ 848, L23.  10.3847/2041-8213/aa9018
  * Frankfurt University/Kadath Initial Data solver [2023] [Dataset] Frankfurt University/Kadath Initial Data solver (2023).  https://kadath.obspm.fr/.  Accessed 28-04-2023
  * Freiburghaus et al. [1999] Freiburghaus, C., Rosswog, S., and Thielemann, F.-K. (1999).  R-process in neutron star mergers.  _ApJ_ 525, L121
  * Fryer et al. [2015] Fryer, C. L., Belczynski, K., Ramirez-Ruiz, E., Rosswog, S., Shen, G., and Steiner, A. W. (2015).  The Fate of the Compact Remnant in Neutron Star Mergers.  _ApJ_ 812, 24.  10.1088/0004-637X/812/1/24
  * Gafton and Rosswog [2011] Gafton, E. and Rosswog, S. (2011).  A fast recursive coordinate bisection tree for neighbour search and gravity.  _MNRAS_ 418, 770–781.  10.1111/j.1365-2966.2011.19528.x
  * Goldstein et al. [2017] Goldstein, A., Veres, P., Burns, E., Briggs, M. S., Hamburg, R., Kocevski, D., et al. (2017).  An Ordinary Short Gamma-Ray Burst with Extraordinary Implications: Fermi-GBM Detection of GRB 170817A.  _ApJL_ 848, L14.  10.3847/2041-8213/aa8f41
  * Gottlieb and Shu [1998] Gottlieb, S. and Shu, C. W. (1998).  Total variation diminishing Runge-Kutta schemes.  _Mathematics of Computation_ 67, 73–85
  * Gourgoulhon [2012] Gourgoulhon, É. (2012).  _3+1 Formalism in General Relativity: Bases of Numerical Relativity_.  Lecture Notes in Physics (Springer Berlin Heidelberg)
  * Gourgoulhon et al. [2001] Gourgoulhon, E., Grandclément, P., Taniguchi, K., Marck, J.-A., and Bonazzola, S. (2001).  Quasiequilibrium sequences of synchronized and irrotational binary neutron stars in general relativity: Method and tests.  _Physical Review D_ 63\.  10.1103/physrevd.63.064029
  * Grandclément et al. [2001] Grandclément, P., Bonazzola, S., Gourgoulhon, E., and Marck, J.-A. (2001).  A multidomain spectral method for scalar and vectorial poisson equations with noncompact sources.  _Journal of Computational Physics_ 170, 231–260.  10.1006/jcph.2001.6734
  * Grandclément [2010] Grandclément, P. (2010).  Kadath: A spectral solver for theoretical physics.  _Journal of Computational Physics_ 229, 3334–3357.  https://doi.org/10.1016/j.jcp.2010.01.005
  * Hajela et al. [2022] Hajela, A., Margutti, R., Bright, J. S., Alexander, K. D., Metzger, B. D., Nedora, V., et al. (2022).  Evidence for X-Ray Emission in Excess to the Jet-afterglow Decay 3.5 yr after the Binary Neutron Star Merger GW 170817: A New Emission Component.  _ApJL_ 927, L17.  10.3847/2041-8213/ac504a
  * Hallinan et al. [2017] Hallinan, G., Corsi, A., Mooley, K. P., Hotokezaka, K., Nakar, E., Kasliwal, M. M., et al. (2017).  A radio counterpart to a neutron star merger.  _Science_ 358, 1579–1583.  10.1126/science.aap9855
  * Hotokezaka et al. [2018] Hotokezaka, K., Kiuchi, K., Shibata, M., Nakar, E., and Piran, T. (2018).  Synchrotron Radiation from the Fast Tail of Dynamical Ejecta of Neutron Star Mergers.  _ApJ_ 867, 95.  10.3847/1538-4357/aadf92
  * Hotokezaka et al. [2015] Hotokezaka, K., Piran, T., and Paul, M. (2015).  Short-lived 244Pu points to compact binary mergers as sites for heavy r-process nucleosynthesis.  _Nature Physics_ 11, 1042.  10.1038/nphys3574
  * Kasen et al. [2013] Kasen, D., Badnell, N. R., and Barnes, J. (2013).  Opacities and Spectra of the r-process Ejecta from Neutron Star Mergers.  _ApJ_ 774, 25.  10.1088/0004-637X/774/1/25
  * Kasen et al. [2017] Kasen, D., Metzger, B., Barnes, J., Quataert, E., and Ramirez-Ruiz, E. (2017).  Origin of the heavy elements in binary neutron-star mergers from a gravitational-wave event.  _Nature_ 551, 80–84.  10.1038/nature24453
  * Kochanek [1992] Kochanek, C. S. (1992).  Coalescing binary neutron stars.  _ApJ_ 398, 234–247.  10.1086/171851
  * Krishnan et al. [2007] Krishnan, B., Lousto, C. O., and Zlochower, Y. (2007).  Quasi-Local Linear Momentum in Black-Hole Binaries.  _Physical Review D_ 76, 081501.  10.1103/PhysRevD.76.081501
  * Lamb and Kobayashi [2017] Lamb, G. P. and Kobayashi, S. (2017).  Electromagnetic counterparts to structured jets from gravitational wave detected mergers.  _MNRAS_ 472, 4953–4964.  10.1093/mnras/stx2345
  * Lamb and Kobayashi [2018] Lamb, G. P. and Kobayashi, S. (2018).  GRB 170817A as a jet counterpart to gravitational wave triggerGW 170817\.  _MNRAS_ 478, 733–740.  10.1093/mnras/sty1108
  * Lattimer and Schramm [1974] Lattimer, J. M. and Schramm, D. N. (1974).  Black-Hole-Neutron-Star collisions.  _ApJ, (Letters)_ 192, L145
  * Löffler et al. [2012] Löffler, F., Faber, J., Bentivegna, E., Bode, T., Diener, P., Haas, R., et al. (2012).  The Einstein Toolkit: a community computational infrastructure for relativistic astrophysics.  _Classical and Quantum Gravity_ 29, 115001.  10.1088/0264-9381/29/11/115001
  * Löffler et al. [2012] Löffler, F., Faber, J., Bentivegna, E., Bode, T., Diener, P., Haas, R., et al. (2012).  The Einstein Toolkit: A Community Computational Infrastructure for Relativistic Astrophysics.  _Class. Quantum Grav._ 29, 115001.  doi:10.1088/0264-9381/29/11/115001
  * LORENE [2001] [Dataset] LORENE (2001).  http://lorene.obspm.fr.  Accessed 28-04-2023
  * LORENE Reference guide. Definition of Lorene::Binaire::angu_mom [2023] [Dataset] LORENE Reference guide. Definition of Lorene::Binaire::angu_mom (2023).  https://lorene.obspm.fr/Refguide/classLorene_1_1Binaire.html#a3e33ff8bbf0552580707577de0937ba4.  Accessed 06-04-2023
  * Lyman et al. [2018] Lyman, J. D., Lamb, G. P., Levan, A. J., Mandel, I., Tanvir, N. R., Kobayashi, S., et al. (2018).  The optical afterglow of the short gamma-ray burst associated with GW170817.  _Nature Astronomy_ 2, 751–754.  10.1038/s41550-018-0511-3
  * Margalit and Metzger [2017] Margalit, B. and Metzger, B. D. (2017).  Constraining the Maximum Mass of Neutron Stars from Multi-messenger Observations of GW170817.  _ApJ_ 850, L19.  10.3847/2041-8213/aa991c
  * Margutti et al. [2017] Margutti, R., Berger, E., Fong, W., Guidorzi, C., Alexander, K. D., Metzger, B. D., et al. (2017).  The Electromagnetic Counterpart of the Binary Neutron Star Merger LIGO/Virgo GW170817. V. Rising X-Ray Emission from an Off-axis Jet.  _ApJ_ 848, L20.  10.3847/2041-8213/aa9057
  * Metzger et al. [2010] Metzger, B. D., Martinez-Pinedo, G., Darbha, S., Quataert, E., Arcones, A., Kasen, D., et al. (2010).  Electromagnetic counterparts of compact object mergers powered by the radioactive decay of r-process nuclei.  _MNRAS_ 406, 2650–2662.  10.1111/j.1365-2966.2010.16864.x
  * Monaghan [1992] Monaghan, J. J. (1992).  Smoothed particle hydrodynamics.  _Ann. Rev. Astron. Astrophys._ 30, 543
  * Monaghan [2005] Monaghan, J. J. (2005).  Smoothed particle hydrodynamics.  _Reports on Progress in Physics_ 68, 1703–1759
  * Nakar and Piran [2011] Nakar, E. and Piran, T. (2011).  Detectable radio flares following gravitational waves from mergers of binary neutron stars.  _Nature_ 478, 82–84.  10.1038/nature10365
  * Nedora et al. [2019] Nedora, V., Bernuzzi, S., Radice, D., Perego, A., Endrizzi, A., and Ortiz, N. (2019).  Spiral-wave Wind for the Blue Kilonova.  _ApJ Letters_ 886, L30.  10.3847/2041-8213/ab5794
  * Papenfort et al. [2022] Papenfort, L. J., Most, E. R., Tootle, S., and Rezzolla, L. (2022).  Impact of extreme spins and mass ratios on the post-merger observables of high-mass binary neutron stars.  _MNRAS_ 513, 3646–3662.  10.1093/mnras/stac964
  * Papenfort et al. [2021] Papenfort, L. J., Tootle, S. D., Grandclément, P., Most, E. R., and Rezzolla, L. (2021).  New public code for initial data of unequal-mass, spinning compact-object binaries.  _Physical Review D_ 104\.  10.1103/physrevd.104.024057
  * Press et al. [1992] Press, W. H., Flannery, B. P., Teukolsky, S. A., and Vetterling, W. T. (1992).  _Numerical Recipes_ (New York: Cambridge University Press)
  * Price [2012] Price, D. J. (2012).  Smoothed particle hydrodynamics and magnetohydrodynamics.  _Journal of Computational Physics_ 231, 759–794.  10.1016/j.jcp.2010.12.011
  * Read et al. [2009] Read, J. S., Lackey, B. D., Owen, B. J., and Friedman, J. L. (2009).  Constraints on a phenomenologically parametrized neutron-star equation of state.  _Phys. Rev. D_ 79, 124032.  10.1103/PhysRevD.79.124032
  * Rezzolla et al. [2018] Rezzolla, L., Most, E. R., and Weih, L. R. (2018).  Using Gravitational-wave Observations and Quasi-universal Relations to Constrain the Maximum Mass of Neutron Stars.  _ApJ_ 852, L25.  10.3847/2041-8213/aaa401
  * Rezzolla and Zanotti [2013] Rezzolla, L. and Zanotti, O. (2013).  _Relativistic Hydrodynamics_
  * Rosswog [2009] Rosswog, S. (2009).  Astrophysical smooth particle hydrodynamics.  _New Astronomy Reviews_ 53, 78–104
  * Rosswog [2015a] Rosswog, S. (2015a).  Boosting the accuracy of SPH techniques: Newtonian and special-relativistic tests.  _MNRAS_ 448, 3628–3664.  10.1093/mnras/stv225
  * Rosswog [2015b] Rosswog, S. (2015b).  SPH Methods in the Modelling of Compact Objects.  _Living Reviews of Computational Astrophysics (2015)_ 1\.  10.1007/lrca-2015-1
  * Rosswog [2020a] Rosswog, S. (2020a).  A Simple, Entropy-based Dissipation Trigger for SPH.  _ApJ_ 898, 60.  10.3847/1538-4357/ab9a2e
  * Rosswog [2020b] Rosswog, S. (2020b).  The Lagrangian hydrodynamics code MAGMA2.  _MNRAS_ 498, 4230–4255.  10.1093/mnras/staa2591
  * Rosswog and Diener [2021] Rosswog, S. and Diener, P. (2021).  SPHINCS_BSSN: a general relativistic smooth particle hydrodynamics code for dynamical spacetimes.  _Classical and Quantum Gravity_ 38, 115002.  10.1088/1361-6382/abee65
  * Rosswog et al. [2022] Rosswog, S., Diener, P., and Torsello, F. (2022).  Thinking Outside the Box: Numerical Relativity with Particles.  _Symmetry_ 14, 1280.  10.3390/sym14061280
  * Rosswog and Korobkin [2022] Rosswog, S. and Korobkin, O. (2022).  Heavy elements and electromagnetic transients from neutron star mergers.  _Annalen der Physik_
  * Rosswog et al. [1999] Rosswog, S., Liebendörfer, M., Thielemann, F.-K., Davies, M., Benz, W., and Piran, T. (1999).  Mass ejection in neutron star mergers.  _A & A_ 341, 499–526
  * Rosswog et al. [2018] Rosswog, S., Sollerman, J., Feindt, U., Goobar, A., Korobkin, O., Wollaeger, R., et al. (2018).  The first direct double neutron star merger detection: Implications for cosmic nucleosynthesis.  _A &A_ 615, A132.  10.1051/0004-6361/201732117
  * Savchenko et al. [2017] Savchenko, V., Ferrigno, C., Kuulkers, E., Bazzano, A., Bozzo, E., Brandt, S., et al. (2017).  INTEGRAL Detection of the First Prompt Gamma-Ray Signal Coincident with the Gravitational-wave Event GW170817.  _ApJL_ 848, L15.  10.3847/2041-8213/aa8f94
  * Shibata [2016] Shibata, M. (2016).  _Numerical Relativity_.  10.1142/9692
  * Shibata et al. [2017] Shibata, M., Fujibayashi, S., Hotokezaka, K., Kiuchi, K., Kyutoku, K., Sekiguchi, Y., et al. (2017).  Modeling GW170817 based on numerical relativity and its implications.  _Phys. Rev. D_ 96, 123012.  10.1103/PhysRevD.96.123012
  * Shibata and Nakamura [1995] Shibata, M. and Nakamura, T. (1995).  Evolution of three-dimensional gravitational waves: Harmonic slicing case.  _Phys. Rev. D_ 52, 5428–5444.  10.1103/PhysRevD.52.5428
  * SPHINCS_ID [2023] [Dataset] SPHINCS_ID (2023).  https://sphincsid.bitbucket.io.  Accessed 28-04-2023
  * Springel [2010] Springel, V. (2010).  Smoothed Particle Hydrodynamics in Astrophysics.  _ARAA_ 48, 391–430.  10.1146/annurev-astro-081309-130914
  * Symbalisty and Schramm [1982] Symbalisty, E. and Schramm, D. N. (1982).  Neutron Star Collisions and the r-Process.  _Astrophys. Lett._ 22, 143
  * Tauris and van den Heuvel [2023] Tauris, T. M. and van den Heuvel, E. P. J. (2023).  _Physics of Binary Star Evolution. From Stars to X-ray Binaries and Gravitational Wave Sources_.  10.48550/arXiv.2305.09388
  * Timmes and Swesty [2000] Timmes, F. X. and Swesty, F. D. (2000).  The Accuracy, Consistency, and Speed of an Electron-Positron Equation of State Based on Table Interpolation of the Helmholtz Free Energy.  _ApJS_ 126, 501–516
  * Troja et al. [2018] Troja, E., Piro, L., Ryan, G., van Eerten, H., Ricci, R., Wieringa, M. H., et al. (2018).  The outflow structure of GW170817 from late-time broad-band observations.  _MNRAS_ 478, L18–L23.  10.1093/mnrasl/sly061
  * Troja et al. [2017] Troja, E., Piro, L., van Eerten, H., Wollaeger, R. T., Im, M., Fox, O. D., et al. (2017).  The X-ray counterpart to the gravitational-wave event GW170817.  _Nature_ 551, 71–74.  10.1038/nature24290
  * Troja et al. [2020] Troja, E., van Eerten, H., Zhang, B., Ryan, G., Piro, L., Ricci, R., et al. (2020).  A thousand days after the merger: Continued X-ray emission from GW170817.  _MNRAS_ 498, 5643–5651.  10.1093/mnras/staa2626
  * Vretinaris et al. [2020] Vretinaris, S., Stergioulas, N., and Bauswein, A. (2020).  Empirical relations for gravitational-wave asteroseismology of binary neutron star mergers.  _Physical Review D_ 101\.  10.1103/physrevd.101.084039
  * Watson et al. [2019] Watson, D., Hansen, C. J., Selsing, J., Koch, A., Malesani, D. B., Andersen, A. C., et al. (2019).  Identification of strontium in the merger of two neutron stars.  _Nature_ 574, 497–500.  10.1038/s41586-019-1676-3
  * Wendland [1995] Wendland, H. (1995).  Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree.  _Advances in Computational Mathematics_ 4, 389–296
  * Zhang and Shu [2016] Zhang, Y. and Shu, W. (2016).  ENO and WENO schemes.  _Handbook of Numerical Ananlysis_ 10.1016/bs.hna.2016.09.009

## Appendices

##  Appendix A Surface-conforming ovals

In order to use the surface-conforming ovals as discussed in Sec. 4.2, we
first need to find the surface of the star. We do this in SPHINCS_ID in the
following way.

  1. 1.

Choose a threshold density that marks the end of the star. Currently, this
value is $\rho_{\mathrm{thres}}=10^{-12}\,\mbox{code units}\simeq 6.2\times
10^{5}\,\mathrm{g/cm^{3}}$ (see the beginning of Sec. 3.1 for the definition
of our code units). Also, choose a tolerance value $\mathsf{tol}=10^{-6}$
(code units), the choices of these values are motivated in step 3.

  2. 2.

Select a direction $(\theta,\phi)$ centered at the center of the star, and
choose two points lying on it, one inside and one outside the star. Read the
density provided by the ID solver at these points.

  3. 3.

Apply the bisection algorithm around $\rho_{\mathrm{thres}}$ until the
distance between the two points is equal to or lower than $\mathsf{tol}$.
Then, the innermost of the two points is declared to be the radius $r$ along
the considered $(\theta,\phi)$ direction. Arguably, the values of
$\rho_{\mathrm{thres}}$ and $\mathsf{tol}$ are somewhat arbitrary. We
determined them by comparing the resulting values for the 4 radii along the
positive and negative $x$ direction, and along the $y$ and $z$ direction, with
the values provided by LORENE.

  4. 4.

Repeat steps 2 and 3 for the desired number of directions, to find a set of
points $\\{\theta,\phi,r\\}$ that lie on the surface of the star. Currently,
we sample the surface every degree in both $\theta$ and $\phi$; we thus have
$180\times 360$ points.

Once the sample $\\{\theta,\phi,r\\}$ lying on the surface is found, the value
of $r(\theta,\phi)$ for any $(\theta,\phi)\in[0,\pi)\times[0,2\pi)$ is found
with bilinear interpolation.

##  Appendix B The SPH estimate of the ADM linear momentum of the fluid

Here we explicitly compute the ADM linear momentum of the fluid in terms of
the SPH canonical momentum, as mentioned in Sec. 4.4. The integral in (15) can
be written as an integral over the spacelike hypersurface $\Sigma$ using
Gauss’ theorem,

$\displaystyle P_{\xi}^{\rm ADM}$
$\displaystyle\coloneqq\dfrac{1}{8\pi}\int_{\partial\sigma}\mathrm{d}S_{j}\left(K^{j}{}_{i}-\delta^{j}{}_{i}\,K\right)\xi^{i}$
$\displaystyle=\dfrac{1}{8\pi}\int_{\Sigma}\mathrm{d}^{3}x\sqrt{\gamma}\left[\xi^{i}D_{j}\left(K^{j}{}_{i}-\delta^{j}{}_{i}\,K\right)+\left(K^{j}{}_{i}-\delta^{j}{}_{i}\,K\right)D_{j}\xi^{i}\right],$
(75)

with $\gamma$ determinant of the spatial metric $\gamma_{ij}$, and $D_{j}$
covariant derivative compatible with $\gamma_{ij}$. The second term in the
square parenthesis can be rewritten as

$\displaystyle\left(K^{j}{}_{i}-\delta^{j}{}_{i}\,K\right)D_{j}\xi^{i}$
$\displaystyle=\left(K^{ij}-\gamma^{ij}\,K\right)D_{j}\xi_{i}$
$\displaystyle=\left(K^{ij}-\gamma^{ij}\,K\right)\dfrac{D_{j}\xi_{i}+D_{i}\xi_{j}}{2}=\left(K^{ij}-\gamma^{ij}\,K\right)\dfrac{\mathscr{L}_{\xi}\gamma_{ij}}{2},$
(76)

where the first equality uses the compatibility between $\gamma_{ij}$ and
$D_{j}$, the second equality uses the symmetry of $K^{ij}$ and $\gamma^{ij}$
in their indices, and the last equality uses the definition of the Lie
derivative $\mathscr{L}_{\xi}\gamma_{ij}$. The first term in the square
parenthesis in (B) can be rewritten using the momentum constraint [Baumgarte
and Shapiro, 2010, eq.(2.126)],

$\displaystyle D_{j}\left(K^{j}{}_{i}-\delta^{j}{}_{i}\,K\right)=8\pi s_{i},$
(77)

where $s_{i}$ is the spatial part of the momentum density measured by the
Eulerian observer
$s_{\rho}\coloneqq-n^{\mu}\,T_{\mu\nu}\,\gamma^{\nu}{}_{\rho}$, with
$T_{\mu\nu}$ being the stress–energy tensor, $n^{\mu}=(1,-\beta^{i})/\alpha$
vector normal to $\Sigma$ (and 4-velocity of the Eulerian observer), $\alpha$
lapse function, $\beta^{i}$ shift vector, and
$\gamma^{\mu}{}_{\nu}=\delta^{\mu}{}_{\nu}+n^{\mu}n_{\nu}$ projector onto
$\Sigma$. Hence, the ADM linear momentum in (B) can be rewritten as

$\displaystyle P_{\xi}^{\rm ADM}$
$\displaystyle=\int_{\Sigma}\mathrm{d}^{3}x\sqrt{\gamma}\left[\xi^{i}s_{i}+\dfrac{1}{16\pi}\left(K^{ij}-\gamma^{ij}\,K\right)\mathscr{L}_{\xi}\gamma_{ij}\right].$
(78)

Now we use the definition of $s_{\mu}$ to write

$\displaystyle s_{\mu}\coloneqq-n^{\rho}\,T_{\rho\nu}\,\gamma^{\nu}{}_{\mu}$
$\displaystyle=-n^{\rho}T_{\rho\mu}-\left(n^{\rho}T_{\rho\nu}n^{\nu}\right)n_{\mu}.$
(79)

The spatial part of $s_{\mu}$ is (see Eq. (35))

$\displaystyle s_{i}$ $\displaystyle=-n^{\nu}T_{\nu
i}-\left(n^{\nu}T_{\nu\rho}n^{\rho}\right)n_{i}=-n^{\nu}T_{\nu
i}=\dfrac{\beta^{j}T_{ji}-T_{0i}}{\alpha}$ (80)

where we used $n_{\mu}=(-\alpha,0,0,0)$ and the expression for the components
of $n^{\mu}$, given above.
The SPH canonical momentum per baryon, for a perfect fluid, is defined as
(see, e.g., [Rosswog, 2009, eq. 198])

$\displaystyle
S_{i}\coloneqq\Theta\left(1+u+\dfrac{P}{n}\right)g_{i\mu}v^{\mu},$ (81)

where $v^{\mu}=U^{\mu}/U^{0}$ is the fluid velocity in the SPH computing
frame, $U^{\mu}=\mathrm{d}x^{\mu}/\mathrm{d}\tau$ is the 4-velocity of the
fluid, $\tau$ is the proper time of the fluid, $\Theta$ is the generalized
Lorentz factor
$\Theta\coloneqq\left(-g_{\mu\nu}v^{\mu}v^{\nu}\right)^{-1/2}=U^{0}$, with
$g_{\mu\nu}$ spacetime metric constructed from $\gamma_{ij}$, $\alpha$ and
$\beta^{i}$; $u$ is the specific internal energy, $P$ the pressure, and $n$
the baryon density of the fluid. The last three hydrodynamical quantities are
measured in the local rest frame of the fluid. It holds:

$\displaystyle U_{i}$
$\displaystyle=g_{i\mu}U^{\mu}=g_{i0}U^{0}+g_{ij}U^{j}=U^{0}\beta_{i}+\gamma_{ij}U^{j}\Longrightarrow\gamma_{ij}U^{j}=U_{i}-U^{0}\beta_{i},$
$\displaystyle g_{i\mu}v^{\mu}$
$\displaystyle=g_{i\mu}\dfrac{U^{\mu}}{U^{0}}=\dfrac{g_{i0}U^{0}+g_{ij}U^{j}}{U^{0}}=g_{i0}+\dfrac{\gamma_{ij}U^{j}}{U^{0}}=\beta_{i}+\dfrac{U_{i}-U^{0}\beta_{i}}{U^{0}}=\dfrac{U_{i}}{U^{0}}.$
(82)

Using (B) in (81),

$\displaystyle S_{i}=\left(1+u+\dfrac{P}{n}\right)U_{i}.$ (83)

Since the SPH canonical momentum is a momentum per baryon, we can try to build
it from the stress–energy tensor in the following way:

$\displaystyle q_{i}\coloneqq\dfrac{T_{0i}}{N},$ (84)

where $N\coloneqq\sqrt{-g}\,\Theta\,n$ is the baryon number density in the
computing frame, with $g$ determinant of $g_{\mu\nu}$. Dimensionally, $q_{i}$
is a momentum per baryon. The expression (81) holds for a perfect fluid, hence
we now specialize to the stress–energy tensor of a perfect fluid,

$\displaystyle q_{i}=\dfrac{(e+P)U_{0}U_{i}+P\,g_{0i}}{N},$ (85)

with $e=n(1+u)$ energy density of the fluid measured in its local rest frame.
We use

$\displaystyle
U_{0}=g_{0\mu}U^{\mu}=g_{00}U^{0}+g_{0i}U^{i}=g_{00}\Theta+g_{0i}\Theta
v^{i}=\Theta\left(-\alpha^{2}+\beta^{j}\beta_{j}+\beta_{j}v^{j}\right),$ (86)

in expanding (85),

$\displaystyle q_{i}$
$\displaystyle=\dfrac{e+P}{\sqrt{-g}\Theta\,n}\Theta\left(-\alpha^{2}+\beta^{j}\beta_{j}+\beta_{j}v^{j}\right)U_{i}+\dfrac{P}{\sqrt{-g}\Theta\,n}\beta_{i}.$
(87)

We substitute $e=n(1+u)$ into (87),

$\displaystyle q_{i}$
$\displaystyle=\dfrac{1}{\sqrt{-g}}\left[\left(1+u+\dfrac{P}{n}\right)U_{i}\right]\left(-\alpha^{2}+\beta^{j}\beta_{j}+\beta_{j}v^{j}\right)+\dfrac{P}{\sqrt{-g}\Theta\,n}\beta_{i}$
$\displaystyle=\dfrac{1}{\sqrt{-g}}\left[S_{i}\left(-\alpha^{2}+\beta^{j}\beta_{j}+\beta_{j}v^{j}\right)+\dfrac{P}{\Theta\,n}\beta_{i}\right].$
(88)

Using (84), the expression in (B) becomes

$\displaystyle
T_{0i}=\dfrac{N}{\sqrt{-g}}\left[S_{i}\left(-\alpha^{2}+\beta^{j}\beta_{j}+\beta_{j}v^{j}\right)+\dfrac{P}{\Theta\,n}\beta_{i}\right]=\Theta\,n\left[S_{i}\left(-\alpha^{2}+\beta^{j}\beta_{j}+\beta_{j}v^{j}\right)+\dfrac{P}{\Theta\,n}\beta_{i}\right],$
(89)

that is, we related the SPH momentum per baryon $S_{i}$ with the $T_{0i}$
components of the perfect fluid stress–energy tensor. Note that, if the
expression for the SPH momentum per baryon was known in general, the same
computation could in principle be generalized to any type of fluid.
We can now insert (89) into (80), and obtain the relation between the Eulerian
momentum density $s_{i}$ and the SPH momentum per baryon $S_{i}$,

$\displaystyle s_{i}$ $\displaystyle=\dfrac{\beta^{j}T_{ji}-T_{0i}}{\alpha}$
$\displaystyle=\dfrac{\beta^{j}T_{ji}}{\alpha}-\dfrac{\Theta\,n}{\alpha}\left[S_{i}\left(-\alpha^{2}+\beta^{j}\beta_{j}+\beta_{j}v^{j}\right)+\dfrac{P}{\Theta\,n}\beta_{i}\right]$
$\displaystyle=\dfrac{\beta^{j}}{\alpha}\left[n\left(1+u+\dfrac{P}{n}\right)U_{i}U_{j}+P\,\gamma_{ij}\right]-\dfrac{\Theta\,n}{\alpha}\left[S_{i}\left(-\alpha^{2}+\beta^{j}\beta_{j}+\beta_{j}v^{j}\right)+\dfrac{P}{\Theta\,n}\beta_{i}\right]$
$\displaystyle=\dfrac{\beta^{j}U_{j}\,n}{\alpha}S_{i}-\dfrac{\Theta\,n}{\alpha}S_{i}\,\beta_{j}v^{j}+\dfrac{P\,\beta_{i}}{\alpha}-\dfrac{P\,\beta_{i}}{\alpha}-\dfrac{\Theta\,n}{\alpha}\left(-\alpha^{2}+\beta^{j}\beta_{j}\right)S_{i}$
$\displaystyle=\dfrac{S_{i}\,n}{\alpha}\left(\beta^{j}U_{j}+\Theta\alpha^{2}-\Theta\beta^{j}\beta_{j}-\Theta\beta_{j}v^{j}\right).$
(90)

Using (B), we rewrite

$\displaystyle\beta^{j}U_{j}=\beta^{j}\Theta
g_{j\mu}v^{\mu}=\Theta\left(\beta^{j}g_{j0}v^{0}+\beta^{j}g_{jk}v^{k}\right)=\Theta\left(\beta^{j}\beta_{j}+\beta_{k}v^{k}\right),$
(91)

where we used $v^{0}=U^{0}/U^{0}=1$. Substituting (91) into (B),

$\displaystyle
s_{i}=\dfrac{S_{i}\,\Theta\,n}{\alpha}\left(\beta^{j}\beta_{j}+\beta_{k}v^{k}+\alpha^{2}-\beta^{j}\beta_{j}-\beta_{j}v^{j}\right)=\Theta\,n\,\alpha\,S_{i}.$
(92)

Using (92), the expre/ssion for the ADM momentum determined by the fluid
becomes [see (16)]

$\displaystyle P_{\xi}^{\rm ADM,fluid}$
$\displaystyle=\int_{\Sigma}\mathrm{d}^{3}x\sqrt{\gamma}\,\xi^{i}S_{i}=\int_{\Sigma}\mathrm{d}^{3}x\left(\sqrt{\gamma}\,\Theta\,n\alpha\right)\xi^{i}S_{i}=\int_{\Sigma}\mathrm{d}^{3}x\,N\,\xi^{i}S_{i}.$
(93)

where we used $\sqrt{-g}=\alpha\sqrt{\gamma}$ and
$N\coloneqq\sqrt{-g}\,\Theta\,n$.
The integral in (93) is over the entire spacelike hypersurface $\Sigma$, but
it has support only where $S_{i}$ (or $N$) is nonzero, that is, where the
fluid is. Therefore, we can approximate the integral with an SPH summation
over the particles, to obtain an estimate of the ADM linear momentum of the
fluid,

$\displaystyle P_{\xi}^{\rm ADM,fluid}$
$\displaystyle=\int_{\Sigma}\mathrm{d}^{3}x\,N\,\xi^{i}S_{i}\simeq\sum_{a}\nu_{a}\left(\xi^{i}S_{i}\right)_{a},$
(94)

with the index $a$ running over all the particles, $\nu_{a}$ being the baryon
number of particle $a$.

##  Appendix C Computation of the ADM linear momentum for the ID

Both LORENE and FUKA assume a spacetime free of gravitational waves which
implies that the whole momentum should be carried by the fluid alone. We
briefly show this here explicitly.
As discussed in (4.4), on the initial spacelike hypersurface and with the
assumptions of asymptotic flatness, conformal flatness and maximal slicing,
used by LORENE and FUKA to solve for the ID, $e^{i}_{(j)}=\delta^{i}{}_{j}$ is
an orthogonal frame which becomes orthonormal, hence Cartesian, at spatial
infinity. We can then set $\xi^{i}=e^{i}_{(j)}$ in (16), and obtain a formula
for the Cartesian components of the ADM linear momentum at $t=0$:

$\displaystyle P_{j}^{\rm ADM}$
$\displaystyle=\int_{\Sigma}\mathrm{d}^{3}x\sqrt{\gamma}\left[e^{i}_{(j)}s_{i}+\dfrac{1}{16\pi}\left(K^{mi}-\gamma^{mi}\,K\right)\mathscr{L}_{e_{(j)}}\gamma_{mi}\right].$
(95)

We now use (71) to note that

$\displaystyle\mathscr{L}_{e_{(j)}}\gamma_{mi}$
$\displaystyle=\mathscr{L}_{e_{(j)}}(A^{2}\delta_{mi})=\delta_{mi}\,\mathscr{L}_{e_{(j)}}A^{2}+A^{2}\,\mathscr{L}_{e_{(j)}}\delta_{mi}$
$\displaystyle=\delta_{mi}e^{k}_{(j)}\partial_{k}A^{2}+A^{2}\big{(}e^{k}_{(j)}\partial_{k}\delta_{mi}+\delta_{ik}\partial_{m}e^{k}_{(j)}+\delta_{mk}\partial_{i}e^{k}_{(j)}\big{)}=\delta_{mi}e^{k}_{(j)}\partial_{k}A^{2},$
(96)

were the last equality follows from $\delta_{mi}$ and $e^{k}_{(j)}$ having
constant components. Inserting (C) into (95),

$\displaystyle P_{j}^{\rm ADM}$
$\displaystyle=\int_{\Sigma}\mathrm{d}^{3}x\sqrt{\gamma}\left[e^{i}_{(j)}s_{i}+\dfrac{1}{16\pi}\left(K^{mi}-\gamma^{mi}\,K\right)\delta_{mi}e^{k}_{(j)}\partial_{k}A^{2}\right].$
(97)

Finally, using the maximal slicing condition
$K=K^{i}{}_{i}=K^{mi}\gamma_{mi}=A^{2}K^{mi}\delta_{mi}=0$, which implies
$K^{mi}\delta_{mi}=0$ because $A\neq 0$,

$\displaystyle P_{j}^{\rm ADM}$
$\displaystyle=\int_{\Sigma}\mathrm{d}^{3}x\sqrt{\gamma}\,e^{i}_{(j)}s_{i}=\int_{\mathrm{stars}}\mathrm{d}^{3}x\sqrt{\gamma}\,s_{j}.$
(98)

Hence, for the LORENE and FUKA binary neutron star ID, it suffices to
integrate (98) over a compact domain (the stars), to obtain the total linear
ADM momentum.

##  Appendix D An example of function approximation via a Local Regression
Estimate (LRE)