The hydrodynamic evolution equations in SPHINCS_BSSN are modeled via a high-accuracy version of smooth particle hydrodynamics (SPH), see [42–46] for reviews of the method. The basics of the relativistic SPH equations have been derived very explicitly in Section 4.2 of Rosswog [43], 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 studies [38, 47, 48] 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⁻⁶ s for time and to 6.17797 × 10¹⁷ gcm⁻³ for densities. We further use the metric signature (-,+,+,+), and we measure all energies in units of m0c2, where m₀ is the average baryon mass¹. 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 wai, 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., [49]):

where α is the lapse function, βⁱ 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 [49–52]. 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

and W_ab(h_k) is a shorthand for W(|r→a-r→b|/hk). As a numerical energy variable, we use the canonical energy per baryon

which is evolved according to

with

and

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 ξⁱ which tends to a Cartesian basis vector at spatial infinity, the ADM linear momentum along the direction of ξⁱ is defined as [53, Equation (8.40)] (see also [54, Sec. II])

where ∂Σ is the boundary of a spacelike hypersurface Σ that extends up to spatial infinity, dσ_j is its surface element, K^j_i is the extrinsic curvature, and K = Kⁱ_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 ξⁱ. 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 -g=αγ [49, 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 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 [47]. Here, we only provide the details of our preferred kernels. The first of these favorites is the Wendland C6-smooth kernel [55]:

where we use exactly 300 constributing neighbor particles. The normalization constant is σ = 1, 365/(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 [56]:

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 neighbor particles. For this kernel, the normalization constant is σ₈ = 1.17851074088357 in 3D. Contrary to Wendland kernels [55], 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 [57], the Wendland C2, C4, and C6 kernels (see [58] for the explicit expressions) and the W8H 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 (“neighbors”) for our cubic lattice. For neighbor numbers just beyond 200, the density accuracy in this experiment is about three orders of magnitude better for W8H 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 [47], we choose in this study, the W8H-kernel, and we found 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 neighbor 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 [59], for more details on the exact procedure, we refer to the MAGMA2 code study [38] where the same approach was used. A large number of experiments confirms that we get very similar results when using W^WC6 and W8H.

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 [38]. 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.1). Since no changes compared to our previous study have been made, we refer the interested reader for the equations, implementation details and tests to Rosswog et al. [39].

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ⁱ, 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 [33], and the modifications needed for piecewise polytropic equations of states (EOSs) are laid out explicitly in Appendix A of Rosswog et al. [39].

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 [40], 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. [39]. 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 study, we restrict ourselves to results obtained for the APR3 EOS [60] only. This EOS allows for a maximum mass of MTOVmax=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_⊙ [61–67] for the “best educated guess” of the maximum neutron star mass and a recent Bayesian study [68] suggests a maximum TOV mass of 2.52 -0.29+0.33 M_⊙, all broadly consistent with our choice of the APR3 EOS. More sophisticated treatments of high-density nuclear matter physics will be addressed in future.

We evolve the spacetime according to the (“Φ-version” of the) BSSN equations [69, 70]. We have written wrappers around code extracted from the well tested McLachlan thorn of the Einstein Toolkit [71, 72]. 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 βⁱ (shift vector) and they read

where γ = det(γ_ij), Γ~jki 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, Rij=R~ij+Rijϕ, where

The derivatives on the right-hand side of the BSSN equations are evaluated via standard finite differencing techniques and, unless mentioned otherwise, we use the 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. [34], to which we refer the interested reader. For the gage 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.

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 third order Runge–Kutta algorithm [73], three times per numerical time step.

The mesh-to-particle mapping step is performed via fifth order Hermite interpolation, that we developed [33] extending the work of Timmes and Swesty [74]. 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 [33] 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 [34, 39], 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 [75] 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).

In Rosswog et al. [39], we used kernel functions of different orders borrowed from vortex methods [76, 77], 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.

The remnant of a binary NS merger can, depending on the EOS, the total mass and spin (and further processes which are not modeled 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 ~0.3, see e.g., Figures 14 and 16 in Rosswog and Diener [33]. 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 a. [79], 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 neighbors 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 [71] 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, Equations (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₀ 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⁰ 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ⁱ 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 neighbor particles in any negative way until they can be safely removed once their lapse value has dropped below the removal threshold.

As described in Diener et al. [34, Sec. 2.2.2], the initial neutron stars are modeled 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 [38], Rosswog and Diener [33], and Diener et al. [34, Sec. 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 ν₀ = ν_tot/N_part, where ν_tot is the total baryon number for the star and N_part is 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 NID(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 NID(r→a) while using SPH particles of the same baryon number². 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 NID(r→a) do not change significantly anymore³.

While we want to construct initial conditions with densities N_a as close as possible to NID(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 ρ = ρ₀ + δρ, where ρ₀ is the true solution, the resulting pressure is

With P₀ = P(ρ₀), the relative error ϵ_P in the pressure reads

Neutron star equations of state can be approximated by piecewise polytropes [40], 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 Equation (62) to

so 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, ~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 study were produced using (67).

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. [34, Sec. 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. [34, 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 APM iteration uses “ghost” or “boundary” particles outside the star that prevent the particles modeling the star from being pushed outside the stellar surface, see Diener et al. [34, Sec. 2.2.2] and [34, 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 with 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 [34], we placed the boundary particles on a lattice, between two ellipsoidal surfaces, and 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. [34, 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 the 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 <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 toward 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.

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 [36, 80]. In coorbiting coordinates of Cartesian type, the conformal flatness condition can be written as Gourgoulhon et al. [36, Sec. IV.A], Papenfort et al. [80, 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(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(j)i 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⁴. 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_⊙ LORENE ID with 2 million SPH particles, see Section (4.1), the two estimates are as follows:

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 as follows:

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 [34, Sec. 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 cancelation (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.

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 [84], see [85] 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¹⁴ g cm⁻³, then the remnant bounces back and, after several more oscillations, the peak density settles near a value of 9.5 × 10¹⁴ g cm⁻³. As expected, the lapse is the lowest, where the density is the 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 ψ₄ (black), and how these are calculated in detail can be found in Appendix A of Diener et al. [34]. All Ψ₄-based GW waveforms were analyzed using kuibit [86]. Again, we find a rather good agreement of the quadrupole result with the more sophisticated ψ₄ 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 Figure 16. More than 2% of the initial ADM mass and more than 20% of the initial ADM angular momentum are radiated.

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 three 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, ~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⁻³ M_⊙. Only ~6 × 10⁻⁴ M_⊙ of this material is unbound from the black hole.

In the left panel of Figure 11, 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 toward 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 is below 0.02.

In the right panel of Figure 11, we show a comparison between the maximum GW amplitude times at the distance to the observer as extracted from the quadrupole formula and from Ψ₄. 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, ~3.8 × 10⁻² M_⊙ of energy and 1.18 M_⊙² of angular momentum is radiated away by GWs. By analyzing the properties of the final horizon using the tools from the Einstein Toolkit[87], we find that the black hole that formed has a dimensionless spin parameter of 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.,[88]). 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 study.

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 [89, 90] 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 [91] 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. [92], 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 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. [92]. (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. [93]. This relation is

where R_1.6 is the circumferential radius of a 1.6M_⊙ star with the given EOS and Mchirp=(m1m2)3/5/(m1+m2)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₁ = m₂ = 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.