CoRe database of binary neutron star merger waveforms and its application in waveform development

We present the Computational Relativity CoRe collaboration's public database of gravitational waveforms from binary neutron star mergers. The database currently contains 367 waveforms from numerical simulations that are consistent with general relativity and that employ constraint satisfying initial data in hydrodynamical equilibrium. It spans 164 physically distinct configuration with different binary parameters (total binary mass, mass-ratio, initial separation, eccentricity, and stars' spins) and simulated physics. Waveforms computed at multiple grid resolutions and extraction radii are provided for controlling numerical uncertainties. We also release an exemplary set of 18 hybrid waveforms constructed with a state-of-art effective-one-body model spanning the frequency band of advanced gravitational-wave detectors. We outline present and future applications of the database to gravitational-wave astronomy.

Since 2009 our team has contributed to some of the research lines mentioned above. Here, we present the largest-to-date public database of BNS waveforms composed of new simulations and those published in [34,37,43,45,52,54,56,58,82,[87][88][89][90][91][92][93][94][95]. The combined set of simulations required about 150 million CPU-hours on supercomputers in Europe and the United States. We publicly release these data with the goal of supporting researchers and further developments in the field of GW astronomy (www.computational-relativity.org). We plan to extend the database with waveforms from upcoming simulations and from other groups/codes. This article describes the simulation methods in 3 + 1 NR, summarizes the quality of the computed waveforms and the key parameters that characterize the GW, and concludes outlining some of the many applications. We use geometrized units c = G = 1 and express results in terms of solar masses (M = 1.9889 × 10 33 g) if not otherwise stated. Conversion factors to CGS are [L] = GM /c 2 1.47670×10 5 cm and [T ] = GM /c 3 4.92549 × 10 −6 s.

SIMULATION METHODS
Initial data. Initial data are constructed by solving the Einstein constraint equations in the conformal thin sandwich formalism and by imposing hydrodynamical equilibrium for the star fluid [96][97][98]. The fluid's , the number of orbits [note that for highly eccentric orbits close to head-on, the number of orbits can drop below 1], and the employed resolution ∆x of the finest level covering the entire NS for different configurations. Different colored markers refer to different EOS, see top color bar. In the last panel we also include simulations with different grid resolutions and numerical methods (fluxes, mesh refinement strategies etc.); simulations of a fixed configuration performed at the same resolutions but using different methods are marked with vertical bars in this panel.
flow is chosen to be either irrotational [99], or prescribed according to the constant rotational velocity formalism [100,101]. Binaries in quasi-circular orbits are built imposing a helical Killing vector [102], whereas for eccentric orbits an approximate "helliptical" Killing vector is used [56,103]. We use either the public Lorene [104] or the SGRID [56,101,105] code. Both codes use multidomain pseudo-spectral methods with surface fitting coordinates [102,106].
Evolutions. Dynamical simulations are performed using free-evolution schemes for the Einstein equations and general relativistic hydrodynamics (GRHD). For the spacetime, we employ either the BSSNOK [107][108][109] formalism or the Z4c formalism [110][111][112][113][114]. The latter has improved constraint propagation and damping properties with respect to BSSNOK, especially in matter simulations [110,114]. We use the moving puncture gauge [115][116][117][118][119], which can handle automatically the gravitational collapse without the need for excision [120][121][122]. GRHD is solved in flux-conservative form [123]. Some mergers are simulated with microphysical EOS and neutrino cooling is taken into account with a leakage scheme. Viscous effects in GR are also simulated in a few cases using the large eddy scheme (GRLES) developed in [82].
We use two different NR codes: BAM [124][125][126] and THC [42,[127][128][129]. Both codes use a simple mesh refinement scheme whereby the grid hierarchy is composed of nested Cartesian boxes, some of which can be moved to track the orbital motion of the stars [125,130,131]. The grid setup is controlled by the resolution ∆x in the finest levels. The finest refinement levels cover entirely the NSs during the inspiral. The other levels are constructed by progressively coarsening the resolution by factors of two and extend to the wave-extraction zone. Discretization is based on fourth (or higher) order finitedifferencing stencils and GRHD is handled with either standard finite-volume or high-order finite-differencing high-resolution shock-capturing methods [42,43]. The THC code also implements a neutrino leakage scheme and the GRLES [51,82].
Wave extraction. GWs are extracted on coordinate spheres with radius r using the spin-weighted s = −2 spherical harmonics decomposition of the Weyl scalar Ψ 4 , e.g., [125]. Some of the THC simulations employ the Cauchy characteristic extraction technique to obtain Ψ 4 at future null infinity [132]. The metric multipoles are reconstructed using the fixed frequency method [133]. We release the = m = 2 metric multipole h 22 as a function of the coordinate time t and of the retarded time u = t − r * , where r * (r) is the tortoise coordinate defined by assuming r is the isotropic radius and using the binary total mass for the Schwarzschild spacetime. We also release the GW energy and angular momentum emitted during the simulation that are computed as in [134,135]. Our waveforms are extracted at different extraction radii and can be further extrapolated to obtain null-infinity estimates, e.g. [40,136,137].

INPUT PHYSICS
BNS simulations require several assumptions on the NS fluid and input models describing the matter interactions. The yet unknown EOS is among the most important quantities determining the NS properties and the binary dynamics. It determines the tidal deformations and interactions of the stars during the inspiral [138][139][140][141][142], the lifetime and rotation frequency of the merger remnant [36,91,[143][144][145], and the amount of unbound matter ejected during the merger process, e.g. [38,73,74,88]. We release data for 16 different EOSs. Two EOSs are polytropic models with adiabatic index Γ = 2 [52]. Nine EOSs are zero-temperature nuclear physics model represented by piecewise polytropic fits [146]. They are augmented with a Γ-law pressure component during the simulation to approximate temperature effects [147]. Five EOSs are tabulated finite-temperature microphysical models developed in [148][149][150][151][152], which we also release. Finite-temperature effects are crucial during and after merger, when compressional and shock heating are present, e.g. [38,67,70,153,154].
The role of magnetic fields on the post-merger dynamics is currently a key open question [60,66,82,84,86,155,156]. Large magnetic field instabilities might cause turbulence and induce viscosity, potentially affecting the merger remnant, mass outflows and the GW emission. Also, while not directly relevant for GW emission on the dynamical timescale of our simulations, neutrino transport plays a crucial role in the merger remnant, e.g. [35,38,51,64,67,68,70,77,79]. We plan to include more data from simulations with advanced radiation transport schemes and magnetic field effects as robust NR results become available. Binary mass. In contrast to BBHs, BNS dynamics cannot be rescaled by the binary total mass (A, B label the NSs) since M enters the description of tidal interactions during the inspiral and determines the merger remnant. Formation scenarios and the constraints from GW170817 indicate that NS masses lie within ∼ 1.0−2.3M [65,85,[157][158][159][160]. Current observations range from ∼ 1.0M [161,162] to ∼ 2.0M [163,164] or possibly even ∼ 2.3M [165], with BNS masses varying in ∼ 2.5 − 2.9 M [166,167]. The wide mass range in our database fully covers the observational and a large fraction of the theoretical limits.
GWs from a M = 3.2M merger are shown in the top panel of Fig. 2. High-mass mergers likely result in a prompt BH formation; while high-mass BNS emit strong GWs, they are EM faint due to smaller ejecta and disk masses [74,144].
Mass ratio. The mass ratio has a clear imprint on the GW/EM signals: BNS with larger q are less luminous in GWs [54,94,141], but their larger mass ejecta can power bright EM transients [38,54,168]. The isolated NS mass distribution implies mass ratios up to q max 2.3, but population synthesis models predicts lower values q max 1.8−1.9 [56,169]. The largest observed mass ratio in BNSs is q ∼ 1.3 [170,171]. The CoRe database contains data with mass ratios up to q = 2.1, which is the largest simulated so far [37,56]. In this simulation the companion NS is tidally disrupted during the merger leading to postmerger GWs with small amplitude (second panel in Fig. 2).
Spins. The dimensionless spin of a NS in a binary can be defined as where the angular momentum S A is computed from the isolated NS with the same EOS, rotational velocity, and baryonic mass as the constituents of the binary [52,54,56]. The maximum NS spin is not precisely known, since it depends on the EOS, but existing EOS models predict breakup spins below |χ| ∼ 0.7, corresponding to spin periods of less than 1 ms [172]. The fastest spinning NS in a BNS system is PSR J1946+2052 [173] which will have |χ| ∼ 0.05 at merger. For spins parallel to orbital angular momentum (say z-direction), the effective spin [174] is the quantity determining the leading-order spin-orbit effects on the phase evolution of the binary. Spin-orbit interactions quantitatively change the inspiral-merger and remnant dynamics [52,175,176]. Neglecting their effect can bias the GW parameter estimation [177][178][179]. Spin-precession effects in BNS have been first simulated in Refs. [56,58]; the computed GW signal is shown in Fig. 2 (third panel).
Eccentricity. The emission of GWs causes field binaries to circularize to eccentricities e 10 −5 by the time they enter the LIGO-Virgo band [180]. Therefore, for an accurate modeling of the GWs, it is important to simulate small eccentricity binaries. The residual (numerical) eccentricity of the initial data can be reduced to e 10 −3 − 10 −4 using an iterative procedure [56,181,182].
On the contrary, dynamically assembled BNSs or those belonging to hierarchical triplets could be highly eccentric even at the time of merger [183,184]. An example of a highly eccentric merger with e = 0.6 is shown in Fig. 2 (forth panel). The bursts in the GW amplitude are caused by the close encounters of the two stars. These encounters also induce f -mode oscillations which allow an independent constraint of the EOS for upcoming 3rd generation GW detectors [47,185,186].
Tidal parameters. Tidal interactions in the post-Newtonian (PN) formalism are described by a multipolar set of parameters proportional to the relativistic Love numbers [138][139][140]187]. The dominant effect depends on the gravitoelectric quadrupolar Love numbers k A 2 and the compactness C A of the NS through the expression Λ A 2 = 2k 2 /(3C 5 A ). Tidal interactions are attractive and enter at leading PN order in the GW phasing evolution through the combination [142,[188][189][190], The tidal parameterΛ is a key quantity to characterize the non-perturbative regime of the merger dynamics as shown in [34,94,191] and discussed below. Furthermore, it provides a simple but effective parameterization of the characteristic GW post-merger frequencies [38,89,192] and of the disk mass [92]. The value ofΛ for GW170817 is constrained to be 630 on the basis of the analysis of the GW signal alone [6,193]. In addition, Refs. [92,194] suggested that the observation of an EM counterpart to GW170817 allows to place a lower bound on the tidal deformability ofΛ 400, Ref. [92], orΛ 200, Ref. [194]. Further constraints arise from the theoretical modeling of matter near nuclear density, e.g. [195], other astrophysical observations, e.g. [158,[196][197][198][199], and from the combination of all these constraints by considering a large set of possible nuclear physics EOSs, e.g. [200,201].

DATA QUALITY
Waveforms' error budgets based on convergence tests and finite radius extraction have been presented in [40,42,43,45,95,128,129]. Phase convergence is typically observed for about 10 − 15 orbits at sufficiently high resolutions, corresponding to about 96 grid points per NS diameter. The error due to finite-radius extraction dominates in the early part of the simulations, but truncation errors increase towards merger and afterwards where the uncertainty is the largest [40]. Typical accumulated phase errors up to merger are estimated as δφ ∼ 0.2 − 1.5 rad for simulations in which convergence can be proven. We stress, however, that the lowest resolutions employed in our runs are not convergent and do not give quantitatively reliable results for multiple orbits. Post-merger GWs are typically less accurate, but monotonic behavior with grid resolution can be observed at sufficiently high resolutions, e.g. [91]. Our post-merger data are sufficiently robust to infer the energy and frequency content, e.g. [89,91,94].
We assessed systematic errors due to the use of nonlinear numerical schemes used for GRHD [128,135] with extensive testing of different algorithms and/or extensive code comparisons. We have tested consistency between BAM and THC for datasets: BAM:0097 and THC:0036, BAM:0063 and THC:0029, BAM:0064 and THC:0028, using exactly the same initial data. We found that phase differences are below the estimated uncertainties. A simple polytropic EOS setup has also been compared to results obtained with the SpEC code [202][203][204] with similar results [205].
We stress that all of our waveforms are computed using constraint-satisfying initial data in hydrostatic equilibrium. Constraint violating and/or non-hydrostatic initial data exhibits large unphysical fluid oscillations that contaminate the GW signals. These oscillations are significantly reduced and converge to zero if equilibrium is imposed [95]. Systematic errors generated by the initial data were studied by comparing the evolution of a binary produced by SGRID and Lorene using the same evolution setup (BAM:0026, BAM:0027) [52]. Differences in the GW phase and collapse time to black-hole were found to be compatible with those expected from finite grid resolutions effects.

APPLICATIONS
The CoRe waveform database has wide applicability to the study of strong-field BNS dynamics and for GW astronomy.
Our simulations showed that, despite the complexity of the physics involved, the main quantities characterizing the merger dynamics, like the mass-rescaled GW frequency and the binding energy per unit mass, are determined by parameters likeΛ, emerging from perturbative (PN and effective-one-body, EOB) analysis [34,89]. About 100 simulations of the CoRe database were used to compute the total GW luminosity in terms of tidal parameters and the mass-ratio for all BNS with aligned spins |χ z | 0.14, and to set upper limits to the total emitted energy [94].
A related application is the study of the merger outcome. NR data are crucial to understand the formation of massive NS remnant [35,156,206] and prompt black hole formation at merger [144,207,208].
The data we provide can be used to verify and develop inspiral-merger waveform models for LIGO-Virgo analysis. BAM simulations have been already used in the development of the TEOBResum model [87]. Further analytical-numerical comparisons showed that state-ofart tidal EOB models might underestimate tidal effects at merger for stiff EOS and small M [93]. Our spinning BNS are currently used to test the performances of the TEOBResumS model [20,209].
High-resolution BAM simulations (BAM:0037, BAM:0064, BAM:0095) were employed to construct the tidal phase model NRtidal [45,95,179]. The latter is a closedform expression fitting the inspiral-merger GW composed of PN, TEOBResum, and NR data used to augment any BBH waveform model with tidal effects [179]. Notably, NRtidal was used in the LIGO-Virgo analysis of GW170817 [6,193,210,211], and other groups are using similar approaches for GW modeling [212].
A main open challenge is the modeling of GWs from merger remnants [144,145,192,[213][214][215][216][217][218][219]. Several features of the signal are understood, but quantitative models are missing. We anticipate that CoRe data will be used to develop new post-merger models to be employed for the analysis of current and third-generation detectors. The latter are the most promising observatories to capture high-frequency GW signals, e.g. [90,91,220].
Our data can also be injected in synthetic detector noise to test parameter estimation pipelines, similarly to what was done for BBHs [30]. For BNSs, however, complete waveforms spanning thousands of GW cycles during the inspiral and tens of GW post-merger cycles would be needed. To address the problem, we generate hybrid waveforms combining analytical models and NR data and covering the frequency range of ground-based interferometers [221] (see also [191,222,223]). We release 18 of these hybrids corresponding to equal, unequal masses and spinning BNSs.
The CoRe database will have a reach beyond the applications we have just discussed. In the future, we plan to include more quantities from our simulations. For example, mass outflows ejected during merger, e.g. [38,51,70,73,74,79,86,88,168], and disk masses and profiles [92]. These data will be crucial for the interpretation of EM counterparts.
Marconi (PRACE proposal 2016153522 and ISCRA-B project number HP10B2PL6K), on the Hydra and Draco clusters of the Max Planck Computing and Data Facility, the compute cluster Minerva of the Max-Planck Institute for Gravitational Physics, and on the ARA cluster of the University of Jena.