The Lagrangian Numerical Relativity code SPHINCS BSSN v1.0

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 ( χ = 0 . 5) .

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 (≳ 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 ∼ 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 2 we describe the methodology, with Section 2.1 focusing on the hydrodynamics, Section 2.2 on the equation of state, Section 2.3 on the spacetime evolution and Section 2.4 on how spacetime and matter evolution are coupled.In Section 2.5 we describe measures to robustly evolve the collapse of a neutron star to a black hole and in Section 3 we summarize our latest improvements in constructing SPH initial conditions in full General Relativity.In Section 4 we discuss several examples of neutron star mergers, and we summarize our results in Section 5

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.

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 [2015aRosswog [ , 2020b,a] ,a] 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 × 10 −6 s for time and to 6.17797 × 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 mass 1 .Greek indices run from 0 to 3 and latin indices from 1 to 3. Contravariant spatial indices of a vector quantity ⃗ w at particle a are denoted as w i a , 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]) where α is the lapse function, β i the shift vector and γ ij the spatial 3-metric.We use a generalized Lorentz factor The coordinate velocities v α are related to the four-velocities, normalized as U α U α = −1, by 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 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, ν 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) 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 where 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 where the hydrodynamic part is given by and the gravitational part by Here we have used the abbreviations As numerical energy variable we use the canonical energy per baryon which is evolved according to 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 ξ i which tends to a Cartesian basis vector at spatial infinity, the ADM linear momentum along the direction of ξ i is defined as [Gourgoulhon, 2012, Equation (8.40)] (see also [Krishnan et al., 2007, Section II]) where ∂Σ is the boundary of a spacelike hypersurface Σ that extends up to spatial infinity, dσ 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) where γ ij is the spatial metric, s i is the spatial part of the momentum density measured by the Eulerian observer s ρ := −n µ T µν γ ν ρ , and L ξ is the Lie derivative along ξ 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: with the index b running over all the particles, ν b being the baryon number of particle b, N defined in Equation ( 4), s i = Θ n α S i , α is the lapse function, and and Shapiro, 2010, Equation (2.124)].The rightmost formula in Equation ( 17) can be used to compute an estimate of the ADM Frontiers 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 3.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] where we use exactly 300 constributing neighbour particles.The normalization constant is σ = 1365/(64π) 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]) 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 σ 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 H 8 kernel, Equation ( 19).The smoothing lengths h in this experiment are set as multiples of the typical particle separation (ν/N ) 1/3 , h = η(ν/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 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 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], 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 η determines the smoothing length h as a multiple of the typical particle spacing via h = η(ν/N ) 1/3 , the support radius of all kernels is 2h.We have also indicated the corresponding number of neighbour particles.
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 W C6 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 α 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].

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 max TOV = 2.39 M ⊙ and a 1.4 M ⊙ neutron star has a dimensionless tidal deformability of Λ 1.4 = 390.Indirect arguments and the statistics of the radio pulsar/X-ray neutron star distribution point to values in the range of ∼ 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.

Spacetime evolution
We evolve the spacetime according to the ("Φ-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 γ ij (3-metric), K ij (extrinsic curvature), α (lapse function) and β i (shift vector) and they read where γ = det(γ ij ), Γi jk are the Christoffel symbols related to the conformal metric γij and Ãij is the conformally rescaled, traceless part of the extrinsic curvature.The corresponding evolution equations read where and β i ∂i denote partial derivatives that are upwinded based on the shift vector.The superscript "TF" in the evolution equation of Ãij denotes the trace-free part of the bracketed term.Finally 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 and a variant of the "Γ-driver" shift condition with where η is the "β-driver" parameter.

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 ( 25)-( 32), is known on the particles.This bidirectional information exchange is needed at every Runge-Kutta substep; in our case, with an optimal 3 rd order Runge-Kutta algorithm [Gottlieb and Shu, 1998], three times per numerical time step.The mesh-to-particle mapping step is performed via 5 th 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 "Multidimensional 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.

Local Regression Estimate (LRE)
The task is to find the function that best fits the values of the particles surrounding the grid point ⃗ r G , see Figure 2. If we assume that the values at the particle positions {f p } are described by a function f (⃗ r), we  can Taylor-expand this function around the desired grid point ⃗ r G : This Taylor expansion can be interpreted as a polynomial approximation of a given order where the basis functions have been shifted to the point where the coefficient vector reads and the "shifted basis functions" read where The degrees of freedom (DoF) for a given number of dimensions d and a maximum polynomial order of the basis m are given by 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 ⃗ r G , ⃗ β G , are found by minimizing the error functional where ) 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 = (ν 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 ⃗ r G are determined via and, with a few steps of algebra, one finds where is the "moment matrix" and 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 µν component).With increasing polynomial order the condition number of the moment matrix C ≡ ||M M −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 ( 49).
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 ( 43) ) 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 and the function vector becomes In this case the function value at the grid point is estimated as 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.

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 µν estimates for different polynomial orders m, T G,m µν , and then select the "optimal order" that best represents the particle values and is physically admissible.We use the following error measure: 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 µν 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 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?

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 ∇ • ⃗ r = 3 and the numerical estimate given by This expression is one of the standard SPH discretizations (similar to the commonly used expression for ∇ • ⃗ v, see Equation (31) in Rosswog [2009]).To avoid another neighbour-loop over all particles, expression (56) can be conveniently calculated alongside the SPH derivatives.This means that the update of ∇ • ⃗ 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 ∇ • ⃗ 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 ( 56) 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 we calculate the average deviation ⟨δ⟩ G = ( n b=1 δ b )/n over the particles that contribute to the error measure Equation ( 55).If ⟨δ⟩ 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 ⟨δ⟩ G are ≈ 0.005, while those at the surface reach ∼ 0.1.After some experimentation, we have chosen a threshold of 0.05 for ⟨δ⟩ G ; for grid points at which ⟨δ⟩ 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 µν accurately reflect the matter distribution.

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 T G,m µν for the polynomial orders m = 0, 1, 2, 3 using Equation (42) µν since this is a grid point outside the particle surface • in all remaining cases choose the robust "parachute method", that is, polynomial order 0 and T 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.

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 ∼ 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, ∆t SPH = ξ SPH (h/c s ), decreases (here c s is the speed of sound), whereas the Courant time step for the mesh, ∆t BSSN = ξ G (∆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.∆x is the resolution on the finest grid.For the dimensionless pre-factors, our default values are ξ SPH = 0.2 and ξ G = 0.35.Whenever the ratio of the two time steps grows beyond a threshold, we add a refinement level.In numerical experiments, we found that C 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, ∆t BSSN is half its previous value and the time step ratio will again be less than C refine .As the collapse proceeds, ∆t 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, α, 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 α 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 α 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 α dust = 0.05.
The dust particles will no longer affect the evolution of their neighbours and will no longer contribute to ∆t SPH .They will simply evolve along geodesics until the lapse at their positions falls below α 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 α = 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 α 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) That is, S i reduces to the spatial part of the covariant 4-velocity, U µ .Therefore, we can find U 0 so that the 4-velocity is properly normalized, U µ U µ = −1.Raising the index on the covariant 4-velocity, we can simply read off Θ = 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 Θ.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.

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.

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 ν 0 = ν tot /N part , where ν tot is the total baryon number for the star and N part the number of SPH particles used to model it.Subsequently an "artificial pressure" is defined as Here N a is the SPH estimate of the density variable defined in Equation ( 4) on particle a and N ID (⃗ 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 ID (⃗ r a ) while using SPH particles of the same baryon number. 2 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: see Section 2.1 for the meaning of the involved quantities.The iteration stops when the differences between N a and N ID (⃗ r a ) do not change significantly anymore. 3hile we want to construct initial conditions with densities N a as close as possible to N ID (⃗ 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, π a , from the physical pressures rather than the densities as in Equation ( 62) 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 which, for the special case of a polytropic EOS, simply reduces to the polytropic exponent γ.If we assume that we have found a numerical solution for the density that is ρ = ρ 0 + δρ, where ρ 0 is the true solution, the resulting pressure is With P 0 = P (ρ 0 ), the relative error ϵ P in the pressure reads 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 γ > 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 ρ.This motivates us to change the definition of the artificial pressure from (62) to 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 ( 62) and ( 67).With definitions (67) 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.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, ( 62) or (67).The ID used for the runs described in this paper were produced using (67).

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.

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π max , those farthest away a value of 30π max and for those in between, the artificial pressure varies linearly between these bounds, where π 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 3.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 δ 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 δ, and it would be desirable to remove this dependence.If δ is too small, the real particles cannot approach the surface of the star; if δ 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 δ, we set a small value of δ initially, and then let the boundary particles move outwards (effectively increasing δ) very slowly until the condition |r av − R| < h av /3 is met, where r av and h 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 Figure 6.Cuts along the xy plane, with |z| < ∼ 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 ⊙ equal-mass, 47.5 km separation, MS1b, irrotational BNS produced with LORENE.The blue points have |z| < ∼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.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 δ 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 δ 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 3.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.

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)]: with A being the conformal factor and δ ij = diag(1, 1, 1) the Euclidean metric.Hence, a global, orthogonal, non-orthonormal frame e i (j) = δ i j exists on the initial spacelike hypersurface.This frame also becomes orthonormal, hence Cartesian, at spatial infinity where A → 1 due to the asymptotic flatness condition.Therefore, we can set ξ i = e i (j) in Equation ( 16) to compute the j 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. 4LORENE 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 ⊙ LORENE ID with 2 million SPH particles, see Section (4.1), the two estimates are: , 6.336 × 10 −11 , < 10 −15 , (69a) 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 4.2.For the equal-mass 1.3M ⊙ FUKA ID used in the simulation described in Section 4.3, where one star is spinning with χ ≃ 0.5, the estimates are: SPHINCS ID : 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.

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 4.1 we show a binary neutron star merger where a remnant survives (for at least several dynamical time scales), Section 4.2 shows an example where the merger remnant promptly collapses to form a black hole and in Section 4.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 ≈ 2268 km in each coordinate direction and a minimum initial grid spacing of ∆x = 369 m.Keep in mind that new refinement levels are added dynamically, when the criterion described in Section 2.5 is met.In the run presented in Section 4.2 where there is a collapse to a black hole, the number of refinement levels dynamically increases up to 11.

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 × 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 × 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 ψ 4 (black), how these are calculated in detail can be found in Appendix A of Diener et al. [2022].All Ψ 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 ψ 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.

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 ∆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 ∆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 × 10 −3 M ⊙ .Only about 6 × 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 Ψ 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 × 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 ≈ 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.

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 × 1.3 M ⊙ stars, where one of the stars is spinning.The chosen value of the spin parameter, χ = 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 4.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 ∼ 40 • from the polar axis, see the volume rendering in the left panel of Figure 13, with velocities up to ∼ 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 where R 1.6 is the circumferential radius of a 1.6M ⊙ star with the given EOS and M 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 chirp = 1.13 M ⊙ .Inserting these numbers into Equation ( 71) we find that the relation predicts a value of f 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 ℓ = 2, whereas the dashed line includes the contribution from all modes up to ℓ = 4.The non-spinning case radiates about twice as much energy as the spinning case and does it predominantly in the ℓ = 2 modes.The spinning case does show a bit more contribution from the higher ℓ 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.

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 µν (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 2.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 α cut = 0.02), we remove the particle to avoid the time step shrinking towards zero.The lapse value α cut is well below the value where an apparent horizon forms (∼ 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 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 2.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 × 1.3 M ⊙ ) that remains stable on the simulation timescale, one irrotational system (2 × 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, χ = 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 ψ 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 × 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 ≈ 0.8 is significantly larger than the end result of an irrotational binary BH merger where a/M ≈ 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.

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.3.4.The integral in (15) can be written as an integral over the spacelike hypersurface Σ using Gauss' theorem, with γ determinant of the spatial metric γ ij , and D j covariant derivative compatible with γ ij .The second term in the square parenthesis can be rewritten as where the first equality uses the compatibility between γ ij and D j , the second equality uses the symmetry of K ij and γ ij in their indices, and the last equality uses the definition of the Lie derivative L ξ γ ij .The first term in the square parenthesis in (72) can be rewritten using the momentum constraint [Baumgarte and Shapiro, 2010, eq.(2.126)], where s i is the spatial part of the momentum density measured by the Eulerian observer s ρ := −n µ T µν γ ν ρ , with T µν being the stress-energy tensor, n µ = (1, −β i )/α vector normal to Σ (and 4-velocity of the Eulerian observer), α lapse function, β i shift vector, and γ µ ν = δ µ ν + n µ n ν projector onto Σ.Hence, the ADM linear momentum in (72) can be rewritten as Now we use the definition of s µ to write The spatial part of s µ is (see Eq. ( 32)) where we used n µ = (−α, 0, 0, 0) and the expression for the components of n µ , given above.The SPH canonical momentum per baryon, for a perfect fluid, is defined as (see, e.g., [Rosswog, 2009, eq. 198]) that is, we related the SPH momentum per baryon S i with the T 0i components of the perfect fluid stressenergy 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 (86) into (77), and obtain the relation between the Eulerian momentum density s i and the SPH momentum per baryon S i , Using ( 79), we rewrite where we used v 0 = U 0 /U 0 = 1.Substituting ( 88) into (87), Using (89), the expre/ssion for the ADM momentum determined by the fluid becomes [see ( 16)] where we used √ −g = α √ γ and N := √ −g Θ n.The integral in (90) is over the entire spacelike hypersurface Σ, 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, with the index a running over all the particles, ν a being the baryon number of particle a.

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 (3.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) = δ i j is as error measure and (∇f ) num is the numerically calculated gradient estimate while (∇f ) theo is the exact result.Also here, the overall scale between the lowest (linear) and next-to-lowest order (quadratic) is similar, but the quadratic results are much less noisy.Again we find a substantial improvement going to the next order (cubic), but no more substantial gain when further increasing the polynomial order.We also want to briefly illustrate the effect of re-scaling the moment matrix with appropriate powers of the length l p = (ν p /N p ) 1/3 to ensure that all matrix elements have the same dimension.In practice, this is achieved by de-dimensionalizing the shifted basis functions, e.g.∆ G x → ∆ G x/l p , ∆ G x∆ G y → ∆ G x∆ G y/l 2 p etc. We show in Figure A2 the condition numbers, C ≡ ||M M −1 ||, a measure for how close a matrix is to being singular, for the above experiment, once for the straight forward and once for the re-scaled version.The condition numbers in the re-scaled case improve dramatically, e.g. by five orders of magnitude in the case of quartic basis functions (orange curves in Figure A2).We do, however, not see any noteworthy improvement of the accuracy which we interpret as a success of the well-working Singular Value Decomposition.Nevertheless, since the re-scaling is essentially free of computational cost, we always use re-scaled basis functions.

Figure 2 .
Figure2.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 more than 40 particles (= twice the number of degrees of freedom for cubic polynomials) contribute and the error E G,3 is smallest and T G,m=3 00 > 0, choose T 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 T G,m=2 00 > 0, choose T 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 T G,m=1 00 > 0, choose T G,m=1 µν

Figure 3 .
Figure3.The average value of ⟨δ⟩ 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 δ ≈ 0.005.

Figure 4 .
Figure 4.At the top, we show a cut of the density at t = 2.95ms 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.

Figure 5 .
Figure 5. Relative errors on the pressure p err := |p(ρ a ) − p(ρ ID (⃗ r a ))|/p(ρ ID (⃗ r a )) at the end of an APM iteration that uses definition (62) (black) and another APM iteration that uses definition (67) (red).The errors are evaluated on 2.5 × 10 6 particles modeling the same neutron star, which belongs to a 1.9M ⊙ equal-mass, 47.5 km separation, MS1b, irrotational BNS produced with LORENE.(Top) p err (r), with r radial coordinate of each particle measured from the barycenter of the star.The left panel shows the data for r ∈ [0, 9.5] km; the right panel for r ∈ [9.5, 11] km.(Bottom) p 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 b = 2.016M ⊙ .Except for the outermost particles with r ≳ 10 km ≃ 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 (67).For example, ∼ 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 ].

Figure 7 .
Figure 7. Projections of particle distributions with 2.5 × 10 6 particles modeling the same neutron star.Real particles are in black, boundary particles with |z| < ∼0.74 km are in red, and the points lying on the surface of the star with |z| < ∼ 0.16 km are in blue.The star belongs to a 1.9M ⊙ 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 3.2 and Section 3.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.

Figure 8 .Figure 9 .
Figure 8. Density distribution in the orbital plane of an irrotational binary system with 2 ×1.3 M ⊙ with the APR3 EOS.

Figure 10 .Figure 11
Figure 10.Density distribution in the orbital plane of an irrotational binary system with 2 ×1.5 M ⊙ with the APR3 EOS.

Figure 12 .
Figure 12.Density distribution in the orbital plane of a 2 ×1.3 M ⊙ binary with the APR3 EOS.One of the stars no spin, while the other has χ ≃ 0.5.

Figure 13
Figure 13.(Left) Volume rendering of the density distribution at t = 9.93 ms of the 2 × 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 ∼ 0.4c (Right).

Figure 14
Figure 14.(Left) Maximum mass density (red curce, right axis) together with the minimum value of the lapse function α (black curve, left axis.(Right) ℓ = 2, m = 2 mode of the h + polarization of the gravitational wave strain extracted from the Weyl scalar Ψ 4 .Both panels refer to the simulation of a 2 × 1.3 M ⊙ neutron star binary where one star has a spin of χ ≃ 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 .
Figure15.Spectra of the ℓ = 2, m = 2 gravitational waveforms for the non-spinning case (blue curve) and the case with one spinning star (orange curve).

Figure 16 . 4 .
Figure16.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 ℓ = 2 modes, while the dashed lines includes all modes up to ℓ = 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.

Figure A1 .
Figure A1.Relative errors of the function approximation (left) and the gradient of the function approximation (right) as function of radius of the grid points.