# High-resolution numerical relativity simulations of spinning binary neutron star mergers

###### Abstract

The recent detection of gravitational waves and electromagnetic counterparts emitted during and after the collision of two neutron stars marks a breakthrough in the field of multi-messenger astronomy. Numerical relativity simulations are the only tool to describe the binary’s merger dynamics in the regime when speeds are largest and gravity is strongest.

In this work we report state-of-the-art binary neutron star simulations for irrotational (non-spinning) and spinning configurations. The main use of these simulations is to model the gravitational-wave signal. Key numerical requirements are the understanding of the convergence properties of the numerical data and a detailed error budget. The simulations have been performed on different HPC clusters, they use multiple grid resolutions, and are based on eccentricity reduced quasi-circular initial data. We obtain convergent waveforms with phase errors of accumulated over orbits to merger. The waveforms have been used for the construction of a phenomenological waveform model which has been applied for the analysis of the recent binary neutron star detection. Additionally, we show that the data can also be used to test other state-of-the-art semi-analytical waveform models.

^{†}

^{†}publicationid: pubid: 978-1-5386-4975-6 ©2018 IEEE

July 13, 2019

## I Introduction

Neutron stars (NSs) are among the most compact objects in the universe
with central densities multiple times higher than nuclear density.
Similar conditions are unreachable on earth which makes
NSs an exceptional laboratory to test nuclear physics predictions.
In particular the merger of two NSs allows the study of the high
density region of the equation of state (EOS) governing
NS matter.
In addition, NS mergers also allow us to
reveal the central engine for luminous short Gamma ray bursts (sGRBs),
to understand the origin of heavy elements in the universe,
which after their creation produce the optical and near-infrared
EM counterparts, called kilonovae, and
to test astrophysical predictions about binary populations.

The first detection of gravitational waves (GW) combined with an observation of a sGRB and a kilonova marks a breakthrough in the field of multi-messenger astronomy [1, 2]. It is expected that with the increasing sensitivity of advanced GW interferometers multiple GW detections of merging BNSs will be made in the next years [3]. In order to extract information from the data, the measured signal is cross-correlated with a GW template family to obtain a “best match”. This requires accurate GW templates to relate the source properties to the observed GW signal and consequently a detailed analysis of the compact binary coalescence close to the moment of merger.

While an analytical approach to the two-body dynamics in general
relativity is possible for the stage in which the bodies
are well-separated, a numerical solution of the
field equations, dealing with all their nonlinearities,
is needed for a faithful description of the last few orbits. However,
general relativistic simulations are computationally challenging
and expensive. The main reasons are:
(i) the nonlinearity of the equations,
(ii) the intrinsic multi-scale character of the problem
(covering the neutron star interior and the radiation zone),
(iii) no symmetries can be
exploited for generic binary simulations
(3D in space plus time),
(iv) the appearance of shocks and discontinuities in the matter fields.
Over the last years, significant progress has been made simulating BNSs,
with detailed descriptions of physical processes
as finite temperature EOSs, magnetic fields, neutrino transport, e.g. [4, 5, 6, 7, 8],
with new numerical techniques such as discontinuous
Galerkin methods [9, 10, 11]
and high-order convergent schemes [12, 13],
and with the possibility
to study a larger region of the BNS parameter space
with spinning [14], eccentric [15],
and high-mass ratio [16] configurations.

In this article we present recent state-of-the-art numerical relativity
simulations of irrotational and spinning BNSs employing high-order convergent methods.
Notice that although spin is one of the elementary
parameters of a binary system, most studies of BNSs in numerical
relativity have not considered spinning NSs,
see e.g. [17, 18, 19] for
a few exceptions of spinning BNSs described in a consistent approach.
However, it is important to include spin in simulations since
NSs in binaries are spinning objects [20].
The most famous example is the double pulsar PSR J0737-3039 for which
the orbital period and both spin periods are known [21, 22].
In addition to the novelty of describing spinning NSs,
we make use of the fact that recently it became possible to produce NS binaries with
arbitrary eccentricity, in particular, with low eccentricities.
While standard techniques lead to configurations with eccentricities of ,
astrophysical BNS systems will have significantly smaller eccentricities.
Not reducing the eccentricity leads to spurious oscillations in the phase evolution
which reduces the quality of the waveform and hampers waveform model development.

The article is structured as follows: We discuss the configurations and numerical methods in Sec. II, the code performance in Sec. III, and we assess the accuracy of the gravitational waveforms in Sec. IV. We focus in particular on the GW phase which is the most important quantity for waveform modeling and find that to date the presented simulations are the most accurate simulations of spinning BNSs, see e.g. [23, 12, 24, 25] for high-resolution simulations of irrotational BNSs. We present applications of the waveforms in Sec. V and conclude in Sec. VI.

## Ii Numerical Methods & Fiducial binaries

### Ii-a Binary configurations

Name | EOS | ||||||

MS1b | MS1b | 1.3504 | -0.099 | 288.0 | 1.8 | (64,96,128,192) | (0.291,0.194,0.145,0.097) |

MS1b | MS1b | 1.3500 | +0.000 | 288.0 | 1.7 | (64,96,128,192) | (0.291,0.194,0.145,0.097) |

MS1b | MS1b | 1.3504 | +0.099 | 288.0 | 1.9 | (64,96,128,192) | (0.291,0.194,0.145,0.097) |

MS1b | MS1b | 1.3509 | +0.149 | 288.0 | 1.8 | (64,96,128,192) | (0.291,0.194,0.145,0.097) |

H4 | H4 | 1.3717 | +0.000 | 190.0 | 0.9 | (64,96,128,192) | (0.249,0.166,0.125,0.083) |

H4 | H4 | 1.3726 | +0.141 | 190.0 | 0.4 | (64,96,128,192) | (0.249,0.166,0.125,0.083) |

SLy | SLy | 1.3500 | +0.000 | 73.5 | 0.4 | (64,96,128,192,256) | (0.234,0.156,0.117,0.078,0.059) |

SLy | SLy | 1.3502 | +0.052 | 73.5 | 0.4 | (64,96,128,192) | (0.234,0.156,0.117,0.078) |

SLy | SLy | 1.3506 | +0.106 | 73.5 | 0.7 | (64,96,128,192) | (0.234,0.156,0.117,0.078) |

In total nine different physical configurations have been simulated employing three different EOSs (MS1b, H4, SLy), see [26] for more details about the used EOSs. The intrinsic rotation of the NSs is characterized by the dimensionless spin of the stars with being the angular momentum of the stars and their masses in isolation. All spins are either aligned or anti-aligned to the orbital angular momentum. We summarize the physical parameters of the simulated BNSs in Table I.

### Ii-B Initial data construction

For construction of the initial data, the pseudospectral SGRID code [27, 28, 29] is used. SGRID allows the computation of generic neutron star binary configurations by combining the conformal thin sandwich equations [30, 31, 32] together with the constant rotational velocity (CRV) approach [33, 14, 34].

The code uses pseudospectral methods to accurately compute spatial derivatives. The computational domain is split into six patches including spatial infinity such that exact boundary conditions can be imposed. Due to the pseudospectral nature of the code, only relatively few grid points are needed to reach the precision required for our simulations. Thus SGRID does not need much memory, and we run it on a single compute node or a workstation. The most computational expensive routines of SGRID are OpenMP parallelized.

To obtain initial data with reduced eccentricity we employ
an iterative procedure varying the binary’s initial
radial velocity and the eccentricity parameter until
the eccentricity reaches a certain threshold,
see [16].
Figure 1 presents one example of the eccentricity reduction
procedure for the SLy configuration.
In most cases, three iteration steps are sufficient to
obtain final eccentricities .

In addition to its capability to compute eccentricity reduced initial configurations, SGRID’s advantage is that constraint solved initial data in hydrodynamical equilibrium can be computed for spinning NSs. Solving the coupled system of elliptic equations is a challenging task, but required to achieve the necessary accuracy to allow gravitational waveform modeling. As an indicator for the accuracy of the numerical simulations we present central density () oscillations in Fig. 2. In particular, we compute density oscillations

(1) |

for SLy (top panel) and SLy (bottom panel). The time evolution of is characterized by two effects (i) an almost linear decrease of the central density caused by the numerical discretization and (ii) oscillations caused by assumptions of the initial data. While effect (i) clearly decreases with increasing resolution, effect (ii) is almost independent of the resolution. The oscillations are of the order of . For systems not in hydrodynamical equilibrium those density oscillations can be about , see e.g. the study in [15, 16]. Important for our further consideration is that independent of the intrinsic rotation of the individual stars the magnitude of the oscillations remains unchanged, cf. bottom panel of Fig. 2 in which we show the central density oscillations for SLy.

### Ii-C Dynamical Evolution

We employ the BAM code [35, 36, 37, 13] for our simulations. The code contains state-of-the-art methods to deal with general relativistic hydrodynamics (GRHD) as well as black hole spacetimes. Finite difference stencils are used for the spatial discretization of the metric variables, and high resolution shock-capturing methods for the hydrodynamic variables. The evolution algorithm is based on the method of lines.

The code uses an adaptive mesh refinement (AMR) technique in which the domain consists of a hierarchy of nested Cartesian grids (refinement levels). The grid spacing of each refinement level is half the grid spacing of its surrounding coarser refinement level. The innermost levels move dynamically during the time evolution following the motion of the neutron stars such that the strong field region is always covered with the highest resolution. We show a sketch of the refinement strategy in Fig. 3. For the far field region, we have the option to use a “cubed-sphere” multi-patch AMR, which is particularly suited to accurately simulate the distant wave zone. However, to save computational costs we do not employ the cubed-sphere multi-patch AMR in this work. The grid configuration of the presented simulations consists of a total of 7 refinement levels labeled with increasing resolution. The specific grid setup is summarized in Table I.

The BAM code evolves spacetimes using either the BSSN [38, 39, 40] or Z4c [41, 42] formulations. For the simulations proposed in this article, the Z4c evolution system is employed, which leads to an improved accuracy compared to the BSSN formulation, see e.g. [43, 44].

The numerical fluxes for the GRHD system are built using a
flux-splitting approach based on the local Lax-Friedrich flux
and performing the reconstruction on the characteristic
fields [45, 46, 47].
For the reconstruction we use the fifth-order
WENOZ algorithm [48].
We refer to [13] for further information about the
GRHD implementation.

The code employs a hybrid OpenMP/MPI parallelization strategy. Each refinement level is divided into equally sized sub-boxes with ghost zones, whose sizes depend on the applied finite-differencing stencil. Each of the sub-boxes is then owned and evolved by a single MPI process. The ghost zones are synchronized after each evolution step. Thus, each MPI process owns exactly one sub-box of every mesh refinement box and works on the same number of grid points. This optimizes load balancing. In addition, each MPI process launches an equal number of OpenMP threads using shared memory. This helps to decrease the required memory, since fewer MPI processes with fewer ghost zones are needed.

The evolution of the Einstein equations and relativistic hydrodynamics uses approximately 150-200 grid variables (including storage of previous timesteps) of double precision data type, i.e. 8 bytes per grid point and variable. Together with additional variables not directly used in the main evolution of the Einstein equations, and taking into account (i) additional requirements for MPI buffer zones, as well as (ii) the fact that memory usage varies during evolution due the AMR regridding, we estimate that the most demanding configurations need a maximum of about TB.

## Iii Parallel Performance & Computational Resources

Let us discuss the parallel performance of the BAM code in production runs for BNS spacetime.

In the top panel of Fig. 4 we present strong scaling results of the BAM code obtained on the HPC system SuperMUC at the Leibniz Supercomputing Centre of the Bavarian Academy of Sciences and Humanities and on the Italian Tier-0 supercomputer Marconi at CINECA. We used the Intel 16 compiler on SuperMUC and on the Broadwell partition of Marconi and the Intel 17 compiler on Marconi’s Knights Landing partition. The strong scaling test is based on production-like simulations consisting of two NSs covered by a total of seven refinement levels with up to level . We consider different grid setups and employ points in each inner level, the outermost levels use as usual twice the number of points in each direction for up to points.

The plot shows speed (top panel) and efficiency (bottom panel) defined as

(2) |

The efficiency in the bottom panel of Fig. 4 refers to the smallest number of cores on which the grid configuration has been tested on.

We find that for small jobs the code speed is best on the Broadwell
partition of Marconi, although efficiency is better on SuperMUC’s Haswell
and Sandy Bridge nodes.
The largest runs have comparable performances on Broadwell and Sandy Bridge.
Efficiency drops below 75% between 2048 and 4096 cores depending on
the machine.
Interpretation of the scaling results on Knights Landing architectures
requires some care. Such processors have approximately half clock
speed than the others and peak performance can be only reached
exploiting large vector instructions and optimizing memory loads. Although they anticipate some of
the features required for Exascale Computing (e.g. power efficiency
with a large FLOP-s/watt ratios), significant code refactoring is
needed for an optimal usage.
The BAM code, in particular, is optimized for more traditional
architectures and main routines are not vectorized. In our experience the typical
speed is lower by a factor about 5-7 in most of the applications; the
performances are significantly worse for small jobs.

The simulations presented here are the largest BNS simulations performed with the BAM code so far. We used allocations at the HPC clusters SuperMUC (Germany), JURECA (Germany), Stampede (US), and Marconi (Italy), and utilized a total million CPU hours. Our largest production runs use but an optimal setup for simulations with is currently under investigation with preliminary results ongoing. The most demanding setups run on about 1000 compute cores with a total runtime of approximately 2-3 months.

## Iv Waveform Error budget

The main uncertainties in GWs obtained from numerical simulations are (I) errors caused by the discretization of the underlying continuum problem and (II) errors due to the extraction of GWs at finite radii. An additional source of error is the artificial atmosphere (necessary for the simulation of low density regions) and the related possible violation of rest mass conservation. Fortunately, this effect is subdominant and converges to zero with increasing resolution and it is therefore included in (I).

Considering error (I) any finite-differencing approximation computed at resolution can be written as

(3) |

where is the continuum solution for , and the convergence order. Although it is impossible to compute one can extrapolate from a dataset with different resolutions to obtained a more accurate result. A way to proceed is to consider Richardson extrapolation [49]. This requires a dataset at different finite resolutions and an accurate measure of the convergence order .

While BAM’s flux scheme based on reconstruction of characteristic fields [45, 46, 47] and a WENOZ scheme [48] can achieve fifth order convergence for smooth matter fields, we do not find such high convergence orders for our simulations. Ref. [13] pointed out that the surface of the neutron star, which is only continuous, leads to non-ideal weights in the WENOZ reconstruction and consequently to second order convergence.

Second order convergence is, however, achieved independently of the particular details of the
numerical simulation, i.e. configuration parameters and grid setup.
As exemplary cases we present phase differences for SLy (top panel)
and H4 (bottom panel) in Fig. 5.
No alignment of the waveforms has been performed at the beginning of the
simulations. For both setups second order convergence is achieved almost
up to the moment of
merger. We mark the moment of merger (the time where the GW amplitude peaks)
for different resolutions by vertical colored lines.
This allows us to use our datasets for Richardson extrapolation and construct a
more accurate waveform.

To mitigate errors due to extraction at a finite radius (II), we evaluate the waveform at different radii with and extrapolate the phase and amplitude to infinite radius using a polynomial of order ,

(4) |

Figure 6 present as an exemplary case the phase differences between finite radii extracted waveforms and polynomial extrapolated waveforms with (top panel) and (bottom panel) for MS1b. We obtain similar results independent of the chosen extrapolation order, cf. bottom panel of Fig. 6 for the phase differences between extrapolation order , and . However, for higher extrapolation orders we find spurious oscillations and that the extrapolation is less robust, i.e. it depends more on the chosen extraction radii. Therefore, we have decided to use extrapolation. The error is estimated by computing the phase difference with respect to a waveform extracted at finite radius of . Notice that this error is a conservative estimate since also larger extraction radii are available to validate the results and the phase difference between different extrapolation orders is basically zero. While during the inspiral the finite radii extrapolation dominates the overall error budget, close to the merger errors due to the numerical discretization dominate.

## V Applications

### V-a GW templates construction

GW data analysis pipelines generate a large number, , of waveform templates that are then cross-correlated with the signal. As a consequence, waveform models that can be produced in a fast and efficient way must be in place. In Ref. [50] we presented the first closed form tidal approximant combining Post-Newtonian (PN), effective-one-body (EOB), and numerical relativity results. [We call the approximant in the following NRTidal which is the name used in the LSC Algorithm Library Suite.]

The model can be added to any binary black hole baseline and describes the GW phase evolution due to tidal effects during the inspiral up to merger. The two main assumptions for the model are: (i) tidal effects are proportional to the effective tidal coupling constant [51, 52]

(5) |

where is the quadrupolar Love number describing the static
quadrupolar deformation of one body in the gravitoelectric field of
the companion, , and is the compactness of star ,
and (ii) tidal and spin effects can be decoupled.

We now check the validity of our assumptions on the simulation data. In the top panel of Fig. 7 we present the dimensionless GW phase as a function of the GW frequency. Additionally, we extract the tidal contribution to the phase

(6) |

with being the dimensionless GW frequency from

(7) |

The middle panel of Fig. 7
shows for all combinations of irrotational
NS configurations. This together with Eq. (6)
allows us to present an estimate for , which
is the most important quantity in the NRTidal model
(dashed black line in the middle and bottom panel of Fig. 7).
According to Fig. 7 all combinations
lead to a similar result for
which validates assumption (i).

The bottom panel of Fig. 7 shows combinations of simulations including spinning configurations. Combining spin aligned, anti-aligned, and non-spinning data of the same EOS allows us to estimate spin-spin interactions (red line). Clearly visible is that the spin-spin interactions are almost zero and not well resolved in our simulations. This is important since information about the EOS are encoded in spin-spin interaction describing the deformation of the NS due to its intrinsic rotation. Additionally, we show a combination computed with approximately the same spin magnitudes but different EOSs (orange curve). If appreciable spin-tidal coupling was present in our simulations, the curve would deviate from zero. Since both red and orange curves are so close to zero, both spin-tide coupling and the imprint of the EOS on the spin-spin interaction are small. Future simulations with even larger resolutions or possibly higher spin magnitudes might reveal these effects, but current state-of-the art simulations suggest that a decoupling of spin and tidal effects [assumption (ii)] is valid. This significantly simplifies waveform development.

Knowing the time domain behavior of the frequency-domain phase correction in the NRTidal is constructed as follows, first, fitting with a Padé approximant (black dashed line in Fig. 7), second, computing the frequency domain tidal phase numerically under the stationary phase approximation [52], and, finally, fitting the numerical result again with a Padé approximant. The resulting NRTidal model thus gives both a closed-form expression for the time domain and frequency domain tidal corrections.

### V-B GW template verification

A second usage of the simulation data is the verification of different waveform models. In the following we compare our data to the previously introduced NRTidal approximant as well as two distinct EOB models with spin aligned to the angular momentum and tidal effects SEOBNRv4T [53, 54] and TEOBResumS [55, 56, 57, 58].

We show results for the particular H4 configuration in Fig. 8. We find that all the approximants deviate significantly from NR in the representation of last few orbits (). However, the deviation is small: rad for NRTidal and SEOBNRv4T approximant and rad for TEOBResumS. Note the NR error is about up to the moment of merger. Our comparison shows that with state-of-the-art techniques, numerical relativity reaches an accuracy at which calibration of tidal EOB models to simulations becomes possible.

## Vi Conclusion

In this article we presented new high resolution simulations of binary neutron stars. We are able to simulate spinning neutron star binaries and produce waveform that convergence in rigorous self-convergence tests. Using extrapolation based on convergence tests we compute waveforms with phase uncertainties of rad, accumulated over orbits. The simulations have been performed on a variety of HPC systems and required about million CPU hours. Higher resolution runs are ongoing and will allow us to reduce numerical uncertainties even further.

With regard to numerical methods, we are exploring the feasibility of a next generation code, bamps [9, 59, 60], that implements discontinuous Galerkin (DG) methods for numerical relativity and general relativistic hydrodynamics. In principle, such methods allow high-order schemes combined with parallel AMR, with scaling to a much larger number of processors than reported here, although a full-featured implementation of DG methods for binary neutron star simulations does not exist yet.

With the numerical methods employed in BAM, second order convergence is achieved for the gravitational wave phase for all configurations. The absence of higher order convergence is caused by the presence of discontinuities at the neutron star surface. However, the clean second order convergence allows a proper error estimate and the use of Richardson extrapolation to obtain improved representations of the continuum solution. Overall, accurate numerical relativity waveforms are urgently needed to further develop, improve, and test waveform models and to be prepared for future gravitational wave detections in the new era of gravitational wave astronomy.

## Acknowledgments

T.D. acknowledges support by the European Union’s Horizon 2020 research and innovation program under grant agreement No 749145, BNSmergers. S.B. acknowledges support by the European Union’s H2020 under ERC Starting Grant, grant agreement no. BinGraSp-714626. W.T. was supported by the National Science Foundation under grants PHY-1305387 and PHY-1707227. Computations were performed on SuperMUC at the LRZ (Munich) under the project number pr48pu, JURECA (Jülich) under the project number HPO21, Stampede (Texas, XSEDE allocation - TG-PHY140019), Marconi (CINECA) under ISCRA-B the project number HP10BMAB71, and under PRACE allocation from the 14th Tier-0 call.

## References

- [1] B. P. Abbott et al., “GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral,” Phys. Rev. Lett., vol. 119, no. 16, p. 161101, 2017.
- [2] “Multi-messenger Observations of a Binary Neutron Star Merger,” Astrophys. J., vol. 848, p. L12, 2017.
- [3] B. P. Abbott et al., “Upper Limits on the Rates of Binary Neutron Star and Neutron Star–black Hole Mergers From Advanced Ligo’s First Observing run,” Astrophys. J., vol. 832, no. 2, p. L21, 2016.
- [4] K. Dionysopoulou, D. Alic, C. Palenzuela, L. Rezzolla, and B. Giacomazzo, “General-Relativistic Resistive Magnetohydrodynamics in three dimensions: formulation and tests,” Phys. Rev., vol. D88, p. 044020, 2013.
- [5] Y. Sekiguchi, K. Kiuchi, K. Kyutoku, and M. Shibata, “Dynamical mass ejection from binary neutron star mergers: Radiation-hydrodynamics study in general relativity,” Phys.Rev., vol. D91, no. 6, p. 064059, 2015.
- [6] C. Palenzuela, S. L. Liebling, D. Neilsen, L. Lehner, O. L. Caballero, E. O’Connor, and M. Anderson, “Effects of the microphysical Equation of State in the mergers of magnetized Neutron Stars With Neutrino Cooling,” Phys. Rev., vol. D92, no. 4, p. 044045, 2015.
- [7] K. Kiuchi, K. Kyutoku, Y. Sekiguchi, and M. Shibata, “Global simulations of strongly magnetized remnant massive neutron stars formed in binary neutron star mergers,” 2017.
- [8] F. Foucart, “Monte-Carlo closure for moment-based transport schemes in general relativistic radiation hydrodynamics simulations,” 2017.
- [9] M. Bugner, T. Dietrich, S. Bernuzzi, A. Weyhausen, and B. Brügmann, “Solving 3D relativistic hydrodynamical problems with weighted essentially nonoscillatory discontinuous Galerkin methods,” Phys. Rev., vol. D94, no. 8, p. 084004, 2016.
- [10] L. E. Kidder et al., “SpECTRE: A Task-based Discontinuous Galerkin Code for Relativistic Astrophysics,” J. Comput. Phys., vol. 335, pp. 84–114, 2017.
- [11] M. Dumbser, F. Guercilena, S. Koeppel, L. Rezzolla, and O. Zanotti, “A strongly hyperbolic first-order CCZ4 formulation of the Einstein equations and its solution with discontinuous Galerkin schemes,” 2017.
- [12] D. Radice, L. Rezzolla, and F. Galeazzi, “Beyond second-order convergence in simulations of binary neutron stars in full general-relativity,” Mon.Not.Roy.Astron.Soc., vol. 437, pp. L46–L50, 2014.
- [13] S. Bernuzzi and T. Dietrich, “Gravitational waveforms from binary neutron star mergers with high-order weighted-essentially-nonoscillatory schemes in numerical relativity,” Phys. Rev., vol. D94, no. 6, p. 064062, 2016.
- [14] W. Tichy, “Constructing quasi-equilibrium initial data for binary neutron stars with arbitrary spins,” Phys. Rev. D, vol. 86, p. 064024, 2012.
- [15] N. Moldenhauer, C. M. Markakis, N. K. Johnson-McDaniel, W. Tichy, and B. Brügmann, “Initial data for binary neutron stars with adjustable eccentricity,” Phys. Rev., vol. D90, no. 8, p. 084043, 2014.
- [16] T. Dietrich, N. Moldenhauer, N. K. Johnson-McDaniel, S. Bernuzzi, C. M. Markakis, B. Brügmann, and W. Tichy, “Binary Neutron Stars with Generic Spin, Eccentricity, Mass ratio, and Compactness - Quasi-equilibrium Sequences and First Evolutions,” Phys. Rev., vol. D92, no. 12, p. 124007, 2015.
- [17] S. Bernuzzi, T. Dietrich, W. Tichy, and B. Brügmann, “Mergers of binary neutron stars with realistic spin,” Phys.Rev., vol. D89, p. 104021, 2014.
- [18] N. Tacik et al., “Binary Neutron Stars with Arbitrary Spins in Numerical Relativity,” Phys. Rev., vol. D92, no. 12, p. 124012, 2015, [Erratum: Phys. Rev.D94,no.4,049903(2016)].
- [19] T. Dietrich, S. Bernuzzi, M. Ujevic, and W. Tichy, “Gravitational waves and mass ejecta from binary neutron star mergers: Effect of the stars’ rotation,” Phys. Rev., vol. D95, no. 4, p. 044045, 2017.
- [20] D. R. Lorimer, “Binary and Millisecond Pulsars,” Living Rev. Rel., vol. 11, p. 8, 2008.
- [21] M. Burgay, N. D’Amico, A. Possenti, R. Manchester, A. Lyne et al., “An Increased estimate of the merger rate of double neutron stars from observations of a highly relativistic system,” Nature, vol. 426, pp. 531–533, 2003.
- [22] A. Lyne, M. Burgay, M. Kramer, A. Possenti, R. Manchester et al., “A Double - pulsar system - A Rare laboratory for relativistic gravity and plasma physics,” Science, vol. 303, pp. 1153–1157, 2004.
- [23] S. Bernuzzi, M. Thierfelder, and B. Brügmann, “Accuracy of numerical relativity waveforms from binary neutron star mergers and their comparison with post-Newtonian waveforms,” Phys.Rev., vol. D85, p. 104030, 2012.
- [24] K. Hotokezaka, K. Kyutoku, H. Okawa, and M. Shibata, “Exploring tidal effects of coalescing binary neutron stars in numerical relativity. II. Long-term simulations,” Phys. Rev., vol. D91, no. 6, p. 064060, 2015.
- [25] K. Kiuchi, K. Kawaguchi, K. Kyutoku, Y. Sekiguchi, M. Shibata, and K. Taniguchi, “Sub-radian-accuracy gravitational waveforms of coalescing binary neutron stars in numerical relativity,” Phys. Rev., vol. D96, no. 8, p. 084060, 2017.
- [26] J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, “Constraints on a phenomenologically parameterized neutron- star equation of state,” Phys. Rev., vol. D79, p. 124032, 2009.
- [27] W. Tichy, “Black hole evolution with the BSSN system by pseudo-spectral methods,” Phys.Rev., vol. D74, p. 084005, 2006.
- [28] ——, “A New numerical method to construct binary neutron star initial data,” Class.Quant.Grav., vol. 26, p. 175018, 2009.
- [29] ——, “Long term black hole evolution with the BSSN system by pseudo-spectral methods,” Phys.Rev., vol. D80, p. 104034, 2009.
- [30] J. Wilson and G. Mathews, “Instabilities in Close Neutron Star Binaries,” Phys.Rev.Lett., vol. 75, pp. 4161–4164, 1995.
- [31] J. Wilson, G. Mathews, and P. Marronetti, “Relativistic numerical model for close neutron star binaries,” Phys.Rev., vol. D54, pp. 1317–1331, 1996.
- [32] J. York, James W., “Conformal ’thin sandwich’ data for the initial-value problem,” Phys.Rev.Lett., vol. 82, pp. 1350–1353, 1999.
- [33] W. Tichy, “Initial data for binary neutron stars with arbitrary spins,” Phys.Rev., vol. D84, p. 024041, 2011.
- [34] ——, “The initial value problem as it relates to numerical relativity,” Rept. Prog. Phys., vol. 80, p. 026901, 2017.
- [35] B. Brügmann, J. A. Gonzalez, M. Hannam, S. Husa, U. Sperhake et al., “Calibration of Moving Puncture Simulations,” Phys.Rev., vol. D77, p. 024027, 2008.
- [36] M. Thierfelder, S. Bernuzzi, and B. Brügmann, “Numerical relativity simulations of binary neutron stars,” Phys.Rev., vol. D84, p. 044012, 2011.
- [37] T. Dietrich, S. Bernuzzi, M. Ujevic, and B. Brügmann, “Numerical relativity simulations of neutron star merger remnants using conservative mesh refinement,” Phys. Rev., vol. D91, no. 12, p. 124041, 2015.
- [38] T. Nakamura, K. Oohara, and Y. Kojima, “General Relativistic Collapse to Black Holes and Gravitational Waves from Black Holes,” Prog. Theor. Phys. Suppl., vol. 90, pp. 1–218, 1987.
- [39] M. Shibata and T. Nakamura, “Evolution of three-dimensional gravitational waves: Harmonic slicing case,” Phys. Rev., vol. D52, pp. 5428–5444, 1995.
- [40] T. W. Baumgarte and S. L. Shapiro, “On the numerical integration of Einstein’s field equations,” Phys. Rev., vol. D59, p. 024007, 1999.
- [41] S. Bernuzzi and D. Hilditch, “Constraint violation in free evolution schemes: comparing BSSNOK with a conformal decomposition of Z4,” Phys. Rev., vol. D81, p. 084003, 2010.
- [42] D. Hilditch, S. Bernuzzi, M. Thierfelder, Z. Cao, W. Tichy et al., “Compact binary evolutions with the Z4c formulation,” Phys. Rev., vol. D88, p. 084057, 2013.
- [43] M. Ruiz, D. Hilditch, and S. Bernuzzi, “Constraint preserving boundary conditions for the Z4c formulation of general relativity,” Phys. Rev., vol. D83, p. 024025, 2011.
- [44] A. Weyhausen, S. Bernuzzi, and D. Hilditch, “Constraint damping for the Z4c formulation of general relativity,” Phys. Rev., vol. D85, p. 024038, 2012.
- [45] G. Jiang, “Efficient Implementation of Weighted ENO Schemes,” J. Comp. Phys., vol. 126, pp. 202–228, Jun. 1996.
- [46] A. Suresh, “Accurate Monotonicity-Preserving Schemes with Runge Kutta Time Stepping,” J. Comp. Phys., vol. 136, pp. 83–99, Sep. 1997.
- [47] A. Mignone, P. Tzeferacos, and G. Bodo, “High-order conservative finite difference GLM-MHD schemes for cell-centered MHD,” J.Comput.Phys., vol. 229, pp. 5896–5920, 2010.
- [48] R. Borges, M. Carmona, B. Costa, and W. S. Don, “An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws,” Journal of Computational Physics, vol. 227, no. 6, pp. 3191–3211, 2008. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0021999107005232
- [49] L. F. Richardson, “The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam,” Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, vol. 210, no. 459-470, pp. 307–357, 1911. [Online]. Available: http://rsta.royalsocietypublishing.org/content/210/459-470/307
- [50] T. Dietrich, S. Bernuzzi, and W. Tichy, “Closed-form tidal approximants for binary neutron star gravitational waveforms constructed from high-resolution numerical relativity simulations,” Phys. Rev., vol. D96, no. 12, p. 121501, 2017.
- [51] T. Damour and A. Nagar, “Effective One Body description of tidal effects in inspiralling compact binaries,” Phys. Rev., vol. D81, p. 084016, 2010.
- [52] T. Damour, A. Nagar, and L. Villain, “Measurability of the tidal polarizability of neutron stars in late-inspiral gravitational-wave signals,” Phys.Rev., vol. D85, p. 123007, 2012.
- [53] T. Hinderer et al., “Effects of neutron-star dynamic tides on gravitational waveforms within the effective-one-body approach,” 2016.
- [54] J. Steinhoff, T. Hinderer, A. Buonanno, and A. Taracchini, “Dynamical Tides in General Relativity: Effective Action and Effective-One-Body Hamiltonian,” Phys. Rev., vol. D94, no. 10, p. 104028, 2016.
- [55] S. Bernuzzi, A. Nagar, T. Dietrich, and T. Damour, “Modeling the Dynamics of Tidally Interacting Binary Neutron Stars up to the Merger,” Phys.Rev.Lett., vol. 114, no. 16, p. 161103, 2015.
- [56] T. Damour and A. Nagar, “New effective-one-body description of coalescing nonprecessing spinning black-hole binaries,” Phys.Rev., vol. D90, no. 4, p. 044018, 2014.
- [57] A. Nagar, T. Damour, C. Reisswig, and D. Pollney, “Energetics and phasing of nonprecessing spinning coalescing black hole binaries,” 2015.
- [58] A. Nagar et al., in prep.
- [59] D. Hilditch, A. Weyhausen, and B. Brügmann, “Pseudospectral method for gravitational wave collapse,” Phys. Rev., vol. D93, no. 6, p. 063006, 2016.
- [60] H. R. Rüter, D. Hilditch, M. Bugner, and B. Brügmann, “Hyperbolic Relaxation Method for Elliptic Equations,” 2017.