Second release of the CoRe database of binary neutron star merger waveforms

The first observation of gravitational waves (GWs) from a binary neutron star (BNS) coalescence accompanied by electromagnetic (EM) signals marked a milestone in GW astronomy. Numerical relativity (NR) simulations are the main tool to explore the merger dynamics in the strong-field regime and aid the development of BNS gravitational waveforms that are necessary for GW detection and parameter estimation. The largest NR waveforms public catalogs contain data from thousands of binary black holes (BBHs) simulations, covering a significant portion of the mass-ratio and spin parameter space for quasi-circular mergers [1–7] and explore mergers from eccentric and generic orbits [4–6]. Public waveforms from simulations of binaries with NSs are more limited, and include the Computational Relativity (CoRe) database [8] (164 binaries at different resolutions for a total of 367), the SACRA-MPI [9, 10] (46, total 276), the SXS (2, total 6) [3, 11], among others [12, 13]. These waveforms are crucial for developing accurate inspiral-merger GW templates with tidal effects [14–26] and postmerger emission [27–38] with direct applications to equation of state (EOS) constraints [39–42]. The NR simulations performed for these waveforms are also key to determine the properties of the remnants from the binary parameters and the input physics (EOS, mass, spins, etc), e.g. [27, 43–51] (see also [52, 53] for recent reviews). Consequently, new and extended data releases are necessary to support research in the field of GW astronomy.

Here, we present a new release of the CoRe database that comprises 90 new physically distinct BNS configurations at multiple resolutions, for a total of 254 binaries and 590 simulations. The new release includes GW strains and Weyl multipoles information up to the mode. The new data were computed in simulations presented in references [54–66] and include BNS waveforms consistent with the GW events GW170817 [56, 58, 59, 67] and GW190425 [66, 68].

The paper is organized as follows. Section 2 summarizes the employed simulation techniques. Section 3 describes the physics content of the database and the impact of the binary parameters on the waveforms. Section 4 presents a full merger waveform error analysis for a case study, and gives an overview of the average accuracy of the data. Section 5 presents, as a first application of the database, ready-to-use EOS-insensitive fitting formulas for the GW frequency, amplitude and the peak luminosity that characterize the merger as well as analogous relations for the post-merger GW spectra.

The CoRe database is hosted on the public gitlab server

         https://core-gitlfs.tpi.uni-jena.de/core_database

Associated code repositories and resources can be accessed from

           www.computational-relativity.org/

In particular, we provide the python package watpy to ease the checkout of the data and perform standard waveform analyses [69].

1.1. Notation

NR data are computed in geometrized units and solar masses ; we use these units also in this paper unless explicitly indicated. We recall that s and  km. The binary mass is , where are the gravitational masses of the two stars. The mass ratio is defined as , and the symmetric mass ratio is , where corresponds to the equal-mass case, whereas for very unequal masses. The dimensionless, mass-rescaled spin vectors are denoted with for and the spin components aligned with the orbital angular momentum L are labeled as . The effective spin parameter is defined as the mass-weighted aligned spin,

Similarly, one can define the spin parameter [70],

The quadrupolar tidal polarizability parameters are defined as for [71], where and C_i are respectively the gravito-electric Love numbers and the compactness of the ith neutron star (NS). The tidal coupling constant is [71] ¹³

that, similarly to the reduced tidal parameter [72]:

parametrizes the leading-order tidal contribution to the binary interaction potential and waveform phase (note that for q = 1).

The radiated GW (polarizations and ) is decomposed in multipoles as:

where D_L is the luminosity distance, are the spin-weighted spherical harmonics and are respectively the polar and azimuthal angles that define the orientation of the binary with respect to the observer. Each mode can be decomposed in amplitude and phase as:

with a related GW frequency:

A dimensionless frequency relates to the frequency in Hz according to the formula:

The GW strain modes are related to the Weyl Ψ₄ curvature modes by:

CoRe simulations compute only at different extraction radii R. However, the above equation can be integrated to obtain the strain, either by using the fix-frequency integration method [73] or directly in the time-domain and performing a polynomial correction, e.g. [14, 74, 75]. Comparisons between analytical and NR data often use the Regge–Wheeler–Zerilli normalized multipolar waveforms ,

The radiated energy is obtained as [46]:

whereas the angular momentum is computed as:

The data released are computed with . The binary dynamics can be characterized by the binding energy and the orbital angular momentum, we therefore work with the binding energy per reduced mass, obtained by substracting the GW energy loss from the initial ADM mass, and the dimensionless rescaled angular momentum , see [46, 76] for details. The GW luminosity peak is computed: as

The moment of merger is defined as the time of the peak of , and referred simply as 'merger' when it cannot be confused with the coalescence/merger process. Waveforms are often shown in terms of the retarded time:

where r is the coordinate extraction radius in the simulations (assumed close to the isotropic Schwarzschild radius), is the associated tortoise Schwarzschild coordinate, and is the Schwarzschild radius.

2.1. Initial data

Initial data for CoRe simulations are generated solving Einstein's constraint equations in the conformal thin sandwich (CTS) formalism [77] assuming a helical Killing vector and imposing hydrodynamical equilibrium for the NS's fluid [78, 79]. It is assumed that the fluid is either irrotational or in a quasi-equilibrium state with constant rotational velocity, which allows for consistent simulations of NS with spin [80, 81]. In the latter formalism, the rotational part w^a of the fluid's velocity is determined by an angular velocity parameter ωⁱ as:

where are the coordinates of the NS center. Possible definitions of spin for a star in a binary are discussed in references [46, 82, 83]. The spin parameters given in the CoRe database are defined to be those of single NSs in isolation with the same rest mass and the same ωⁱ as the BNS components [46, 82, 84]. To construct initial data with abritrary eccentricities, we use an extension of the helical symmetry condition that is based on approximate instantaneous first integrals of the Euler equations and a self-consistent iteration of the CTS equations [85]. This method also allows us to create low-eccentricity initial data in quasi-circular orbits using an iterative procedure that combines initial data and evolution codes [82] (see also [86–88]).

CoRe initial data are calculated using either Lorene [89–91] or SGRID [80–82, 84, 92, 93]. Both codes use multi-domain pseudospectral methods to solve the CTS equations and surface-fitting coordinates that minimize spurious stellar oscillations at the beginning of the evolutions and guarantee accurate determination of the initial binary global quantities. Lorene can construct irrotational binaries with either piecewise polytropic or tabulated EOS. In the latter case, they are often obtained as cold, β-equilibrated slices of finite-temperature, composition dependent EOS. SGRID can generate irrotational or spinning binaries with piecewise polytropic EOS and arbitrary (or reduced) eccentricity. In particular, SGRID can simulate BNS in which the individual stars rotate close to the breakup spin and have masses which are of the maximum supported NS mass allowed by the EOS [84]. Evolutions of initial data generated with SGRID and Lorene were compared in [46], where we found them to be in good agreement.

Initial data in quasi-circular orbits are characterized by the following global quantities of the 3+1 hypersurfaces: the total baryon mass M_b (a conserved quantity along the evolution); the total binary gravitational mass M, i.e. the sum of the two gravitational masses of the bodies in isolation; the initial orbital frequency and the corresponding ADM mass and angular momentum .

2.2. Evolution codes

2.3. EOS models

CoRe simulations currently employ 18 different EOS models for the neutron star matter. BAM data are computed using analytical EOS in the form:

where is a given piecewise politropic EOS model [117]. It prescribes also a value for the specific internal energy given the rest mass density ρ, augmented with a γ-law 'thermal' pressure term (usually, [45, 118]). The specific parameters we employ for the piecewise polytropic EOS mimic well-established zero-temperature EOS models [117]; tables of these parameters are available on the CoRe website ¹⁴ .

The current release significantly extends the data computed with finite-temperature EOS over the first release. The release includes data from seven finite-temperature EOS, used in the calculation of references [54–59, 63, 66]. The finite-temperature EOS include the following models: BHB [119], BLh [120, 121], BLQ [63, 121], DD2 [122, 123], LS220 [124], SFHo [125], SLy4/SRO [126, 127].

All these EOS include neutrons, protons, nuclei, electrons, positrons, and photons as relevant thermodynamics degrees of freedom. The ALF2 [128] and BLQ EOS [63, 121] are hybrid models accounting for deconfined quark matter. BHB is a hadronic model that includes Λ hyperons [119, 129].

Cold, neutrino-less β-equilibrated matter described by these microphysical EOS predicts NS maximum masses and radii within a larger range than that allowed by current astrophysical constraints, including GW170817 [40, 41, 130]. Figure 1 shows the mass-radius diagram and the quadrupolar tidal polarizability parameter-mass diagram of these EOS. The largest radius of a NS is  km (EOS 2 H) and the smallest radius  km (EOB 2B). The smallest maximum mass is (EOS H3), whereas the largest is (EOS 2 H).

Zoom In Zoom Out Reset image size

Most of these EOS can be found on the CoRe website in tabulated form. In the simulations, the EOS is called during the hydrodynamics evolution in order to compute the pressure from the rest-mass density, the temperature, and the electron fraction, i.e. in the form . Any relevant thermodynamical quantity is evaluated by multi-linear interpolating the tabulated values in , , and Y_e . As common in relativistic hydrodynamics, the EOS is called during the transformation from conservative to primitive variables. The latter takes place at each timestep and grid point and it involves a numerical root finding of the function , where the specific internal energy is implicitly given by the temperature T [134]. Hence, each root-finding step includes another root finder for the function (see [135] for a discussion on computational efficiency and a non-standard approach based on neural networks.)

2.4. Microphysics

Most of the THC simulations account for the loss of energy and lepton number due to the net emission of neutrinos using a leakage scheme [107, 134]. Accordingly, effective neutrino leakage rates are computed as a physically motivated interpolation from the emission and diffusion rates. The latter require the knowledge of the optical depth of each computational zone in such a way as to recover the correct cooling time scale. Neutrino reabsorption is included in some simulations using the M0 scheme [107]. This scheme splits neutrinos in an optically thick component, treated with the leakage scheme, and a free-streaming component. The free-streaming neutrinos and their average energies are obtained by solving the radiative transfer equations on a set of radial rays (the so-called ray-by-ray approach) fully-implicitly in time. More recently, we have implemented an energy-integrated M1 scheme in THC [108]. The new scheme can self-consistently capture the diffusion of neutrinos from the merger remnant and its reabsorption in the ejecta. M1 simulations are not included in the current release of the database, but will be made public as soon as the associated publications have been accepted. Table 1 summarizes all neutrino reactions currently included in THC together with the reference in which the form of the rates we use are derived.

Table 1. Weak reaction rates and references for their implementation. We use the following notation denotes a neutrino, ν_x denote any heavy-lepton neutrino, denotes a nucleon, and A denotes a nucleus.

Reaction	References
	[136]
	[136]
	[137]
	[137]
	[138]
	[137]
	[139]

CoRe simulations are performed for various binary masses, mass ratios, NS spins and EOS as summarized by figure 2. They cover a significant portion of the BNS parameter space and allow to quantitatively explore the connection between the GW morphology and the binary parameters in some detail. Figure 3 illustrates the variety of waveforms contained in the database. In the following, we give an overview of the database content and outline the connections between physics and waveform morphology.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The database contains waveforms from binaries with total masses ranging from to around with 45 datasets reaching mass ratios larger than and up to q = 2.1 [47, 58, 64]. EOS effect can be summarized to some extent ¹⁵ by the quadrupolar tidal polarizability parameters [71], where larger (smaller) values of are associated to stiffer (softer) EOSs. The most compact NSs (and most massive binaries) are associated with the smallest values of (and ), see the right panel of figure 1. The CoRe data encompasses well the mass and EOS variation for realistic BNSs. Waveforms from both irrotational and spinning (using the formalism outlined in section 2.1) quasi-circular mergers are included [46, 48, 60]. For aligned spins, the individual dimensionless components range in ; about seven datasets are from simulations with precession effects [48, 60]. The distribution of key parameters among the CoRe simulations is shown in figure 4.

Zoom In Zoom Out Reset image size

Most of the CoRe waveforms are produced from quasi-circular mergers. The residual eccentricities of non-iterated quasi-circular initial data is usually –10⁻¹, see the bottom right panel of figure 4. About 13 datasets have an initial eccentricity that was reduced through an iterative procedure employing the formalism described in section 2.1. A subset of waveforms refer instead to eccentric mergers with initial eccentricity values as high as [107, 140, 141]. In particular, the simulation in the bottom panels of figure 3 has an initial eccentricity of 0.55.

The effects of mass-ratio, spin, and tides on the orbital dynamics can be studied by means of gauge-invariant energy curves , that are also publicly released. We illustrate this in figure 5 for a few examples. In the inspiral, the binary's angular momentum j decreases due to GW emission and the system becomes more bound (E_b stays negative and increases). Equal mass () non-spinning BBH systems merge with , indicating that about 3% of the binary mass was radiated in GWs to the moment of merger (marker in the figure). Tidal effects in BNS make the potential governing the relative dynamics more attractive. The tidal constribution to the potential at leading order is , i.e. it is stronger for larger tidal polarizabilities and it is short-ranged thus affecting the motion mostly at high frequencies (small separations, r) towards merger. Consequently, the inspiral of an equal mass non-spinning BNS is faster than a BBH inspiral. The binding energy at the moment of merger is , which is smaller than the black hole case because the BNS system is less compact. Mass-ratio effects make the potential more repulsive, but are less effective than tides at high frequencies. The q = 2 BNS shown in figure 5 merges at smaller values than the equal mass because of tidal disruption. The remnant has also larger angular momentum [59].

Zoom In Zoom Out Reset image size

Spin-effects are dominated by the leading-order spin-orbit interaction; their character is thus repulsive or attractive depending on the projection of the spin on the orbital angular momentum [142]. This is analogous to what happens to corotating/counter-rotating circular orbits in Kerr spacetimes that move outwards (inwards) for antialigned (aligned) spin configurations with respect to the nonspinning case. In BBH simulations this effect has been named as 'hang-up' effect [143]. In figure 5, the spinning BNS with is more bound than the non-spinning BNS at the moment of merger with . Note that j in this case includes the spin contribution. Moreover, the eccentric equal-mass case, contrary to the previous ones, shows brief moments of constant E_b indicating the times when the NSs are apart (see inset of figure 5). Energy curves for BNS have been studied in detail in [46, 48, 60] to which we refer for more details. We stress that the properties of BNS systems at the moment of merger can be captured by EOS-insensitive (quasi-universal) relations [16]. The latter can be helpful in waveform modelling and used to estimate the properties of the remnant. We refer to section 5 for further discussion.

High-mass BNS produce a remnant that promptly collapses to a black hole shortly after the moment of merger and before the two cores can bounce [52, 53]. Prompt collapse implies negligible shocked dynamical ejecta, because the bulk of the mass ejection comes from the first core bounce after their collision [54]. Prompt collapse can be characterized by a threshold mass , that mainly depends on the maximum mass of cold equilibria supported by the EOS [45, 144]. The recent analysis of [145], based on 227 finite-temperature EOS and CoRe data ¹⁶ , found that the prompt collapse mass threshold for equal-mass non-spinning BNS is well described by an EOS-insensitive threshold:

where is the compactness of the maximum NS mass, and , . A prompt collapse waveform has a rapidly damped black hole ringdown after the moment of merger as shown in the top panels of figure 3. Consequently, the postmerger GW signal is practically negligible for the sensitivities of both current and next-generation detectors. The lack of shocked ejecta and of a massive disc also implies that equal-mass prompt-collapse mergers have dim EM emission. However, for very asymmetric BNS with , it is the tidal disruption of the secondary NS and its accretion onto the primary to trigger the gravitational collapse [58]. Thus, asymmetric mergers can be electromagnetically bright because they produce massive tidal dynamical ejecta and remnants with accretion discs of mass . This prompt collapse process is mainly controlled by the incompressibility parameter of nuclear matter around the TOV maximum density [146]. A robust, EOS-insensitive criterion is not known in these conditions [58, 146–149], but tidal disruption effects are subdominant to the mass effect; they produce maximal variations from the equal-mass criterion of [145, 146].

Without prompt collapse, the evolution of a NS remnant is driven by the GWs emission of lasting  ms (GW-driven phase) [150, 151]. During this phase, a remnant that collapses to a black hole is called short-lived, while a remnant that settles to an approximately axisymmetric rotating NS is called long-lived. Examples of postmerger signals from these remnants are shown in the last two panels on the right of figure 3, for the equal-mass case. The GW-driven phase is associated to a luminous GW transient that peaks at frequencies –4 kHz [27, 28, 152–155]. The spectrum of this transient is rather complex but has robust and well-studied features at a few characteristic frequencies. Most of the GW power is emitted in the mode at a nearly constant frequency ; the more compact and close to collapse the remnant is, the higher and more varying the emission frequency is. The postmerger dynamics is primarily controlled by the masses of the two stars and the bulk properties of the zero-temperature EOS, in particular maximum TOV mass and compactness [52, 156]. Finite temperature and neutrinos do not produce qualitative differences, other than possibly on the time of gravitational collapse of the remnant [157]. Quantitative differences in the GW signal introduced by finite-temperature and neutrino effects are typically subdominant compared to finite-resolution uncertainties [108, 158]. On the other hand, microphysics plays a crucial role in the EM counterparts and nucleosynthesis from mergers, e.g. [107, 159–162].

The remnant's signal from asymmetric binaries with mass ratio carries the imprint of the tidal disruption during merger [47, 58, 163]. An example of such a waveform is shown in the second panels (top to bottom) of figure 3. Comparing to the equal-mass long-lived case, the postmerger amplitude is significantly smaller because the asymmetric remnant does not experience the violent bounces of the symmetric remnant. For the same reason, the early-times modulations in frequency and amplitude present in the equal-mass case are significantly suppressed in the asymmetric case.

The evolution of a NS remnant beyond the GW-driven phase is uncertain at present. Explorations of the viscous phase using NR simulations have started [59, 164, 165], but they are still incomplete in many ways. While GW emission is expected to be significantly weaker than during merger, remnant's instabilities might enhance GW emission. Current NR results suggest that BNS remnants have an excess of both gravitational mass and angular momentum after the GW-driven phase and when compared to equilibrium configuration with the corresponding baryon mass [166, 167]. Possible mechanisms to shed (part of) this energy are CFS [168, 169] and one-arm instabilities [170–172] that would lead to potentially detectable, long GW transients at  kHz. Example of such waveforms are THC:0028, THC:0029, and THC:0036 [171].

Finally, CoRe data are available for multiple grid resolutions as discussed in section 2 and shown by figure 2. Most of the newly released data contain high resolution simulations with a minimum grid spacing as low as , e.g. the NS are resolved with a uniform mesh of spacing  m. Notably, simulations of more than 20 orbits or up to hundreds milliseconds postmerger and with microphysics were performed at these resolutions. Simulations at multiple resolutions are a crucial aspect for data quality that is discussed next.

Waveform accuracy depends on several aspects of the simulations. Within the CoRe data the largest sources of uncertainty are (i) the truncation error of the numerical scheme, that is regulated by the mesh resolution employed in the simulations, and (ii) the finite extraction radius for the GW data, e.g. [19, 65, 104, 173]. Other aspects are relevant for waveform modeling, as for example, the length of the simulation (number of orbits/GW cycles), the residual eccentricity in quasi-circular initial data, and the simulation of realistic physics (star rotation, EOS, etc).

Waveform accuracy should be studied by the user case-by-case considering amplitude and phase plots with datasets of simulations at different resolutions and extraction radii. This analysis typically requires a minimum of three simulations of the same BNS at different grid resolutions (a 'convergent series') and has been performed by the authors in references [9, 10, 19, 65, 104–106, 173]. We give below in section 4.1 a complete example of error analysis of a orbit inspiral-merger waveform.

In GW astronomy, the quality of a waveform template is commonly assessed using the Wiener product (overlap) between two waveforms for a given detector [174],

where is the power spectral density (PSD) of the detector noise and the Fourier transform of . The inner product allows to define 'accuracy standards' for either detectability or measurements (parameter estimation), e.g. [175–181]. In the former case, one is interested in quantifying the fractional loss of signal-to-noise ratio (SNR) due the use of a sub-optimal, discrete match filter. Since the number of GW events is proportional to the observable volume, and the distance is inversely proportional to the observed SNR, the fractional loss of potential events scales like the cube of the minimum overlap in the discrete template bank [175, 176, 181]. In the latter case, one is interested in quantifying the bias (or the maximum knowledge) on the GW parameters given the noise in the detector (statistical errors) [177, 178, 181]. In practice, one proceeds by defining the faithfulness between two waveforms:

where t₀ and φ₀ are a reference initial time and phase, and its complementary, the unfaithfulness, . By demanding that, at worst, the systematics biases become of the same order as the statistical ones when the noise level is doubled, it is possible to establish the condition [181]:

where ρ is the SNR and . This condition is necessary for unbiased parameter estimation (faithful waveforms); its violation does not imply that an analysis has biases [173, 178, 182, 183]. The above criterion can be used to quantify the accuracy of NR data, for example by calculating the faithfulness between data at different resolutions [65, 173, 183]. We will use the faithfulness measure in section 4.2 to discuss the average accuracy of the data of the CoRe database.

4.1. Example of NR waveform analysis

In this section we present a waveform error analysis for BAM:0066 [20]. This example effectively represents data that exhibit second order convergence. Figure 6 shows the strain Rh₂₂ at the lowest extraction radius available for this simulation, , its amplitude and frequency . Note that in this section we use R instead of D_L .

Zoom In Zoom Out Reset image size

In order to test self-convergence, we compare amplitude and phase differences of between the different resolutions. For this case, we consider the simulation at resolutions grid points on the highest refined AMR level; hereafter low (L), medium (M), and high (H). The convergence rate p is found experimentally by rescaling these differences using the scaling factor [173],

where is the grid spacing at resolution x. We show the self-convergence test in figure 7. The differences decrease with increasing resolution, as one would expect from convergent data. They also increase with increasing simulation time because truncation errors accumulate during the simulation. The optimal scaling is found for p = 2 with , thus indicating second order convergence. In presence of convergence, a measure of the error to be assigned to the (highest resolution) data is given simply by the difference between the two highest resolutions. This is a conservative estimate because (for convergent data) the truncation error is certainly smaller. Alternatively, the experimental convergence factor can be in principle used in a Richardson extrapolation of the data to provide an improved dataset and error estimate [19, 173]. Note that in this procedure the waveforms are not shifted by a relative time and phase shift because the simulations of the convergent series are run using the same initial data with a fixed initial phase.

Zoom In Zoom Out Reset image size

To assess the uncertainties originated from the waveform obtained at finite-extraction radii, R_i , we compare the phase differences between consecutive radii [173]:

and similarly for the relative amplitudes, . In figure 8 we show the differences at the extraction radii . The phase differences decrease at progressively large radii, thus indicating the numerical waveforms are converging towards their true morphology at null infinity. The phase differences are larger at early times and decrease towards merger; note this behaviour has the opposite sign of that of resolution effects [19]. The relative differences in amplitude are for all radii, indicating robust results are obtained already with relatively close extraction sphere. The waveforms can be extrapolated to null infinity using either a polynomial in of order K [173] or the method outline in [184]. The two methods give comparable results; the former is more general and can be applied to the curvature multipoles , the latter is a simpler method for the strain modes. An error due to finite extraction can be then assigned to the data at finite extraction as the difference with the extrapolated data (or viceversa). Another method is to post-process simulations using Cauchy characteristic extraction [185] and to simulate the waveform at future null infinity. This technique was used for some of the CoRe data.

Zoom In Zoom Out Reset image size

The total error budget can be computed as the sum in quadrature of the truncation and finite extraction errors, and it is shown in figure 9 for both the curvature and strain (2,2) modes. As mentioned above, the truncation phase error is typically a factor larger than the finite extraction error (for ) at merger and in simulations with tens of orbits.

Zoom In Zoom Out Reset image size

Finally, we obtain the unfaithfulness of the waveforms between the different resolutions (M–H and L–M). The Wiener integral is evaluated in the frequency range and employing the Advanced LIGO PSD P1200087 [186] from bajes [187]. Here corresponds to the initial GW frequency, and to the frequency at the moment of merger. For the faithfulness threshold in equation (20), we consider as the strict requirement, and , corresponding to the number of intrinsic parameters of a BNS. Similarly to [65], the SNR values are chosen to be . Figure 10 shows the computed values. The smallest unfaithfulness (M–H, ) passes five out of the six accuracy tests, whereas the other one (L–M, ) passes only two, namely and . However, the unfaithfulness value lies closely (or on top) of the threshold .

Zoom In Zoom Out Reset image size

Analyses similar to the one above are necessary to determine the quality of the NR data for GW modeling. Convergence of the data is a necessary requisite for robust error estimates. Other diagnostic quantities used to verify convergence in simulations are constraint violation, baryon mass conservation and the stars oscillations during the first orbits, e.g. [60, 82, 99, 102, 173]. Achieving waveform convergence in long-term evolutions of BNS is a nontrivial result and, in our experience, requires at least fifth order finite-differencing schemes or finite volume schemes with fifth order reconstructions (at the current resolutions) [19, 65, 105]. Second [19], approximately third [105] and clear fourth order convergence [65] has been demonstrated up to merger in some data using these finite-differencing conservative schemes. Extreme mass ratios and NS rotation close to the breakup limit remain challenging to simulate as well as to obtain clean convergence in GW higher (subdominant) modes like and . Work in these directions is ongoing [58, 62, 64, 65]. For example, clear fourth order convergence in the subdominant (3, 2) and (4, 4) modes for q = 1 has been shown in [65]. Postmerger waveforms typically show slower convergence due to shock formation at merger and the complex fluid dynamics in the remnant. Nonetheless, GW spectra have remarkably robust features that can be accurately quantified with NR data, as we shall discuss in section 5. We refer the reader to [33, 38] for recent work on the accuracy of CoRe postmerger waveform.

4.2. Faithfulness analysis

In an attempt to give an overview of the accuracy of the waveform database, we compute the unfaithfulness of the mode waveforms h₂₂ between the highest and second highest resolutions, for the whole database. We use again the zero-detuned, high-power Advanced LIGO PSD [186]. The minimum frequency employed in the integral of equation (18) corresponds to the initial frequency of each individual simulation.

The result of this analysis is summarized in figure 11, where is shown as a function of the number of orbits and different colors mark the microphysics scheme employed in each simulation. The unfaithfulness values are scattered on a wide range, but about 65% of the waveforms lay below the 1% level which is conventionally considered the accuracy threshold for detection purposes. Importantly, the dependence on the number of orbits (simulation length) is very weak and most of the simulations with ten or more orbits have . Several waveforms from multiple-orbits have ; according to the analysis in the previous section, these data can be considered faithful (suitable for parameter estimation) up to signal SNR of 30–80. We note that data with very few orbits (e.g. THC:0019, BAM:0029, and BAM:0082) show a remarkably low unfaithfulness. These simulations have a short inspiral and rather focus on the postmerger signal, which is not considered in this analysis. Hence, small is not necessarily an indication that these simulations are suitable for waveform modeling.

Zoom In Zoom Out Reset image size

A faithfulness analysis for postmerger signals was recently presented by some of us in [33, 38]. There, we found average mismatches of –0.4. The main source of uncertainty in the postmerger waveforms is the numerical resolution (see the above section) and the impact of the resolution on the remnant's collapse.

As a first application of the database, we present in this section new EOS-insensitive relations for the merger and postmerger waveforms. Previous work found that several key quantities characterizing the merger dynamics depend on the unknown EOS mainly throughout the tidal parameters and have a very weak dependence on other details of the matter model, e.g. [29, 58, 151, 154, 188–190]. Similarly, the GW postmerger spectrum has robust features that can be captured within a few percent accuracy by tidal parameters and/or other properties of NS equilibria in EOS-insensitive way [27, 28, 33, 154–156]. These relations have some practical use in GW astronomy because they deliver accurate estimates for the peak luminosity [53, 151] and for the remnant properties [191–193] (see also [53] for a detailed review) and because they are the building blocks to develop NR-informed waveform models.

First, we consider the mass-rescaled GW amplitude and frequency at the moment of merger, and , and update the fits developed in [33, 38, 189]. Following closely the fitting procedure of [38], we represent any quantity by the factorized function:

where each factor Q^M , Q^S , Q^T accounts for the mass ratio in terms of , spin corrections in terms of , and tidal effects in terms of . The first two factors are given by the linear polynomial expressions , and , with . The last factor is instead a rational polynomial:

with . The best fit parameters are shown in table 2. The amplitude and frequency have 1σ errors of and respectively. We also obtain a χ² of for the former and for the latter.

Table 2. Updated fit coefficients for relevant merger and PM quantities.

	a0	k						χ2	Error	R2
	0.55	1	5.27	0.31	−39.21			0.113	2.6%	0.949
		2					−2.00
		3				0.12	11.09
		4					9.72
	0.22	1	0.80	0.25	−1.99		1.80	0.329	4.5%	0.925
		2					599.99
		3				0.1	7.80
		4					84.76
Mf2		1	31.02		29.99		1.13	0.067	3.6%	0.958
		2					−0.99
		3					39.99
		4					27.77

Next, we use the public CoRe data on the emitted GW energy and extract the peak luminosity using equation (13). For BBHs, this quantity does not depend on the mass scale and it is accurately described by the fits of [194]. For BNS, it has been studied in [151]. We propose the ansatz:

where are the mass and spin dependent fits from [194] and:

Note the scaling factor for . By construction, the fit reduces to the BBH case for . The luminosity peak is calculated in geometric units; the conversion factor to CGS units is given by the Planck luminosity . Figure 12 shows the best fit for and the CoRe data; the best fitting coefficients are reported in table 3. The average 1σ deviation is about over the entire dataset with less than a dozen of outliers. The peak luminosities for BNS are the least accurately modelled (4 BNS configurations). The figure shows that the largest peak luminosities are reached by BNS with that correspond to high-mass binaries and prompt collapse mergers. These events can reach peak luminosities of , about an order of magnitude less than BBHs (of any mass). BNS mass ratios can lower of about an order of magnitude, while spins of magnitude do not significantly affect . We stress that BNS with the largest peak luminosity do not correspond in general to the BNS that radiate the largest amount of energy because postmerger emission can radiate further energy [150, 151] (see also figure 5). We can set an upper limit to the total radiated GW energy from our dataset, obtaining .

Zoom In Zoom Out Reset image size

Table 3. Best fit coefficients for the luminosity peak. The last columns show the χ², the fit's relative standard deviation and the coefficient of determination R².

	k							χ2	Error	R2
	1	2.28		−17.74	−0.57	−17.47	4.58	2.23	12%	0.961
	2					13.91	10.10
	3		14.64	−5.35	−50.54	11.61	−29.96

Finally, we illustrate the use of CoRe data to model postmerger GWs by discussing a fit of the postmerger's spectrum peak frequency f₂, e.g. [28, 29, 33, 38]. This peak frequency is a robust feature found in all NR simulations. Direct GW inference on f₂ can be used to constrain NS properties [32, 34, 155, 156, 191, 192, 195]. The peak frequency also enters as one of the central parameters in postmerger waveform models that will be employed in the future for more sophisticated matched-filter analyses [196]. Following [38] we employ again equation (23) to fit the mass-rescaled Mf₂. The best fitting coefficients are presented in table 2 and have a . Figure 13 shows Mf₂ as a function of for selected values of mass ratio and spin. The 1σ error is below ; this precision is in principle sufficient for informative measurements of the NS mass-radius sequence. For example, using the EOS-insensitive relation between f₂ and the maximum density of an equilibrium non-rotating NS put forward in [156], it would be possible to determine the maximum density of an equilibrium non-rotating NS to and the maximum mass to with a single signal at the detectability threshold.

Zoom In Zoom Out Reset image size

As a further illustration, we calibrate the EOS-insensitive relations (mass-rescaled) between f₂ and the NS radius [38, 197, 198]:

where R_X is the equilibrium radius corresponding to a NS with mass . Figure 14 shows Mf₂ as function of R_1.4, R_1.8 and . Best fit parameters are given in table 4. Other features of the postmerger spectrum can be quantified in a similar way. We release reduced postmerger data and analysis scripts on the CoRe website.

Zoom In Zoom Out Reset image size

Table 4. Updated fit coefficients for Mf₂ as a function of the NS radii R_1.4 and R_1.8.

	a0	a1	a2	a3	χ2	Error	R2
	0.24	−0.10		—	0.55	5.9	0.901
	0.23	−0.10		—	0.31	4.5	0.949
	0.15	−0.11			0.31	4.5	0.949
	0.20	−0.10			0.30	4.4	0.952

We presented a new set of BNS simulations for the second release of the CoRe database, expanding it to 254 different binary configurations covering a wide parameter space. The new data includes BNS consistent with the GW events GW170817 [59] and GW190425 [66, 68]. Simulations were performed with a large number of EOSs, including several microphysical models [59, 63, 66]. Some simulations include the effects of neutrinos, either through the leakage scheme [107, 134, 199], or using the M0 approach [54, 107]. Turbulent viscosity is included in some models using the GRLES formalism [61, 109]. Finally, we include simulations produced using a new hybrid numerical flux scheme, EFL, that was introduced in [65] showing fourth order convergence and smaller phase errors than previous simulations using WENO schemes in BAM.

We described in detail the methodology we used to assess the overall accuracy of the waveforms and presented results for all the configurations in the database. The CoRe database waveform have typical unfaithfulness of less than 10⁻², some have unfaithfulness of less than 10⁻⁴, so they are suitable for precision waveform modeling applications. However, to ensure the convergence and usability of the simulations, more extensive analysis is needed. As an example, we showed a full analysis of one of our simulations, BAM:0066, which showed a clear second order convergence and passed several accuracy tests.

Finally, as a first application of the CoRe database, we fitted phenomenological formulas for the merger amplitude, frequency, and GW luminosity. These fits are able to model the CoRe data with high accuracy ( for the merger amplitude and frequency and for the peak luminosity). We also recalibrated various quasi-universal relations between the post-merger peak frequency and the binary parameters, again finding deviations from the universal relations of only a few percent. These were used in [38] to construct the first complete inspiral, merger, and post-merger waveform model for BNS.

We release the CoRe database to the community with the hope that it will enable future discoveries in GW astronomy. Potential applications include the development of new waveform models, the validation of data analysis pipelines and new NR codes, and the planning of future GW experiments. In the future, we plan to release new simulation data on a rolling basis, with data releases taking place at the publication time of the corresponding paper.

A G, M B acknowledge partial support from the Deutsche Forschungsgemeinschaft (DFG) under Grant No. 406116891 within the Research Training Group RTG 2522/1. F Z, M B, S B, G D acknowledge support from the EU H2020 under ERC Starting Grant, No. BinGraSp-714626. D R acknowledges funding from the U.S. Department of Energy, Office of Science, Division of Nuclear Physics under Award Number(s) DE-SC0021177 and from the National Science Foundation under Grant Nos. PHY-2011725, PHY-2020275, PHY-2116686, and AST-2108467. S B acknowledges support from the Deutsche Forschungsgemeinschaft, DFG, Project MEMI Number BE 6301/2-1. B B acknowledges support from the Deutsche Forschungsgemeinschaft, DFG, Grant BR 2176/5-1. W T acknowledges funding from the National Science Foundation under Grants PHY-2011729 and PHY-2136036.

Numerical relativity simulations were performed at various supercomputing centers. ARA, a resource of Friedrich-Schiller-Universtät Jena supported in part by DFG grants INST 275/334-1 FUGG, INST 275/363-1 FUGG and EU H2020 BinGraSp-714626. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (www.lrz.de, Gauss Projects pn29ba, pn56zo, pn68wi). The authors acknowledge the national High Performance Computing Center Stuttgart (HLRS) for providing access to the supercomputer HPE Apollo Hawk under the Grant Numbers GWanalysis/44189 and INTRHYGUE/44215. The authors gratefully acknowledge the computing time granted by the Resource Allocation Board and provided on the supercomputer Lise and Emmy as part of the NHR infrastructure, where resources were granted through the Project bbp00049. Joliot-Curie at GENCI@CEA (PRACE-ra5202). Marconi at CINECA (ISCRA-B Project HP10BMHFQQ, INF20_teongrav and INF21_teongrav allocation). Bridges, Bridges2, Comet, Expanse, Stampede2 (NSF XSEDE allocation TG-PHY160025), NSF/NCSA Blue Waters (NSF AWD-1811236), supercomputers. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

The CoRe data are hosted on the gitlab server at TPI Jena. Data postprocessing was performed on the Virgo 'Tullio' server at Torino supported by INFN.

The data that support the findings of this study are openly available at the following URL/DOI: https://doi.org/10.5281/zenodo.7253784.

The simulation data discussed in this work is publicly available at

www.computational-relativity.org/gwdb/

The database metadata are summarized in the repo core-database-index, which contains a json file with the main properties of the available simulations and the different runs. A repository is associated to one distinct physical binary and contains folders for the different runs performed. For each run, we release a complete metadata file and a HDF5 file with the multipolar waveform for both and at different extraction radii and the energetics.

Access to our private codes and data is possible upon reasonable request. Any use of the simulation data must be done in accordance with the terms of use contained here: www.computational-relativity.org/terms/.

The repository

    https://git.tpi.uni-jena.de/core/watpy provides classes to work with the CoRe waveforms and tutorials. It is also available via PyPI

    https://pypi.org/project/core-watpy/.

The code includes two main modules. The coredb module contains tools to download and upload NR simulation data, menage the metadata of the simulations, visualize statistics of the database, and work with the HDF5 files provided in the CoRe website. The wave module provides methods for the visualization and the analysis of (multipolar) NR outputs, i.e. Weyl curvature and GW strain. watpy is compatible with NR files from BAM, Cactus/Einstein Toolkit (WhiskyTHC/FreeTHC) and the CoRe database.

The data and scripts employed for the development of the fits presented in this work, can be found in [69].

Any reuse of the merger and postmerger fit data must be done in accordance with the Creative Commons Attribution 4.0 International license which applies to the data files.

THC is open source and publicly available at

    https://bitbucket.org/FreeTHC/workspace/projects/THC

A tutorial, microphysics tables, and example parfiles are available at

    http://personal.psu.edu/~dur566/whiskythc.html