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

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

Tim Dietrich1 2, Sebastiano Bernuzzi3 4, Bernd Brügmann5, Wolfgang Tichy6 1Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am Mühlenberg 1, Potsdam-Golm, 14476, Germany
2Nikhef, Science Park 105, 1098 XG Amsterdam, The Netherlands
3Dipartimento di Scienze Matematiche Fisiche ed Informatiche, Universitá di Parma, I-43124 Parma, Italia
4Istituto Nazionale di Fisica Nucleare, Sezione Milano Bicocca, gruppo collegato di Parma, I-43124 Parma, Italia
5Theoretical Physics Institute, University of Jena, D-07743 Jena, Germany
6Department of Physics, Florida Atlantic University, Boca Raton, FL 33431 US

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)
TABLE I: BNS configurations. The first column defines the name of the configuration with the notation: EOS. The subsequent columns describe: the EOS [26], the NS’ individual masses , the stars’ dimensionless spins , the effective dimensionless coupling constant of the binary, the residual eccentricity [16], the resolutions employed for the evolution grid, and the grid resolutions in the finest level, where we set .

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.

Fig. 1: Eccentricity reduction procedure for SLy. We present the proper distance of the system as a function of time. The artificial eccentricity is clearly visible as oscillations in the proper distance.

Ii-B Initial data construction

Fig. 2: Central density oscillations for the non-spinning configuration SLy (top panel) and the spinning configuration SLy (bottom panel). The density oscillations are independent of the intrinsic spin of the NSs which shows the validity of the constructed initial data employing the CRV approach.

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


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

Fig. 3: Example of the adaptive mesh refinement in BAM. We also show isocontour surfaces of the NS density and the metric’s lapse function, that roughly indicates the gravitational potential.

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.

Fig. 4: Strong scaling (top) and efficiency (bottom) of the BAM code for different resolutions on Marconi and SuperMUC. The test has been performed on the phase1 of SuperMUC with Sandy Bridge-EP Xeon E5-2680 8C processors (solid lines), phase 2 of SuperMUC using Haswell Xeon Processor E5-2697v3 (dotted lines), on the Intel Knights Landing partition of Marconi (dashed lines), and on the Broadwell partition of Marconi (dot-dashed lines) with Intel Xeon E5-2697 v4. For reference ideal scaling is shown for the simulations on phase 1 of SuperMUC.

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


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

Fig. 5: Phase differences between different resolutions versus retarded time for configurations SLy (top panel) and H4 (bottom panel). The rescaling factor is chosen such that each rescaled (i.e. dashed) line would coincide with the line from the next lower resolution phase difference for exact second order convergence. The vertical lines mark the moment of merger (peak of the GW amplitude) for different resolutions.

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


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 ,


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.

Fig. 6: Phase difference between waveforms extracted at various finite radii and extrapolated waveforms assuming (top panel) and (bottom panel). In the bottom panel also the phase difference between different extrapolation orders is shown.

V Applications

Fig. 7: Top panel: Accumulated phase/number of orbits of the numerical relativity simulations with respect to a dimensionless frequency of , which corresponds to about 2 orbits after the start of the simulations. Middle panel: Estimates for , Eq. (6), for combinations of different simulations. The dashed black line is the estimate of used in the NRTidal model. Bottom panel: Combinations of different simulations to extract the spin-spin contribution to the phase in the late inspiral (red) as well as to estimate possible tidal-spin couplings (orange). For comparison, we include again the estimate of of the NRTidal model as a black dashed line.

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]


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


with being the dimensionless GW frequency from


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

Fig. 8: Comparison of the H4 configuration with the NRTidal (blue), SEOBNRv4T (red), TEOBResumS (green) waveform models. Numerical uncertainty is shown as a blue shaded region. As a gray curve we include the real part of the numerical waveform. The waveform approximants are aligned to the numerical relativity waveform in the gray shaded region.

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.


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.


  • [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.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description