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

We present the second data release of gravitational waveforms from binary neutron star merger simulations performed by the Computational Relativity (CoRe) collaboration. The current database consists of 254 different binary neutron star configurations and a total of 590 individual numerical-relativity simulations using various grid resolutions. The released waveform data contain the strain and the Weyl curvature multipoles up to $\ell=m=4$. They span a significant portion of the mass, mass-ratio,spin and eccentricity parameter space and include targeted configurations to the events GW170817 and GW190425. CoRe simulations are performed with 18 different equations of state, seven of which are finite temperature models, and three of which account for non-hadronic degrees of freedom. About half of the released data are computed with high-order hydrodynamics schemes for tens of orbits to merger; the other half is computed with advanced microphysics. We showcase a standard waveform error analysis and discuss the accuracy of the database in terms of faithfulness. We present ready-to-use fitting formulas for equation of state-insensitive relations at merger (e.g. merger frequency), luminosity peak, and post-merger spectrum.

The paper is organized as follows. Sec. 2 summarizes the employed simulation techniques. Sec. 3 describes the physics content of the database and the impact of the binary parameters on the waveforms. Sec. 4 presents a full merger waveform error analysis for a case study, and gives an overview of the average accuracy of the data. Sec. 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 http://www.computational-relativity.org/ In particular, we provide the python package watpy to ease the checkout of the data and perform standard waveform analyses.

Notation
NR data are computed in geometrized units c = G = 1 and solar masses M = 1; we use these units also in this paper unless explicitly indicated. We recall that GM /c 3 4.925490947 µs and GM /c 2 1.476625038 km. The binary mass is M = m 1 + m 2 , where m 1,2 are the gravitational masses of the two stars. The mass ratio is defined as q = m 1 /m 2 ≥ 1, and the symmetric mass ratio is ν = m 1 m 2 /M 2 ∈ [0, 1/4], where ν = 1/4 corresponds to the equal-mass case, whereas ν → 0 for very unequal masses. The dimensionless, mass-rescaled spin vectors are denoted with χ i for i = 1, 2 and the spin components aligned with the orbital angular momentum L are labeled as χ i = χ i ·L/|L|. The effective spin parameter χ eff is defined as the mass-weighted aligned spin, Similarly, one can define the spin parameter [69], The quadrupolar tidal polarizability parameters are defined as Λ i = (2/3) k 2,i C −5 i for i = 1, 2 [70], where k 2,i and C i are respectively the = 2 gravito-electric Love numbers and the compactness of the i-th neutron star (NS). The tidal coupling constant is [70] that, similarly to the reduced tidal parameter [71] Λ = 16 13 (m 1 + 12m 2 )m 4 parametrizes the leading-order tidal contribution to the binary interaction potential and waveform phase (note that κ T 2 = 3 16Λ for q = 1.) The radiated GW (polarizations h + and h × ) is decomposed in ( , m) multipoles as where D L is the luminosity distance, −2 Y m are the s = −2 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 h m (t) can be decomposed in amplitude A m (t) and phase φ m (t) as with a related GW frequency A dimensionless frequencyω = GM ω relates to the frequency in Hz according to the formula The GW strain modes are related to the Weyl Ψ 4 curvature modes ψ m bÿ CoRe simulations compute only ψ m at different extraction radii R. However, the above equation can be integrated to obtain the strain, either by using the fix-frequency integration method [72] or directly in the time-domain and performing a polynomial correction, e.g. [14,73,74]. Comparisons between analytical and NR data often use the Regge-Wheeler-Zerilli normalized multipolar waveforms Ψ m , The radiated energy is obtained as [46] whereas the angular momentum is computed as The data released are computed with max = 4. 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, [46,75] for details. The GW luminosity peak is computed as The moment of merger is defined as the time of the peak of A 22 (t), 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), r * is the associated tortoise Schwarzschild coordinate, and r S = 2M is the Schwarzschild radius.

Initial data
Initial data for CoRe simulations are generated solving Einstein's constraint equations in the conformal thin sandwich (CTS) formalism [76] assuming a helical Killing vector and imposing hydrodynamical equilibrium for the NS's fluid [77,78]. 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 [79,80]. In the latter formalism, the rotational part w a of the fluid's velocity is determined by an angular velocity parameter ω i as where x k * are the coordinates of the NS center. Possible definitions of spin for a star in a binary are discussed in Refs. [46,81,82]. 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 ω i as the BNS components [46,81,83]. 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 selfconsistent iteration of the CTS equations [84]. 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 [81] (see also [85][86][87]).
CoRe initial data are calculated using either Lorene [88][89][90] or SGRID [79-81, 83, 91, 92]. 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 ∼ 98% of the maximum supported NS mass allowed by the EOS [83]. Evolutions of initial data generated with SGRID and Lorene were compared in Ref. [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 Ω ω 22 /2 and the corresponding ADM mass M ADM and angular momentum J ADM .

Evolution codes
CoRe simulations evolve initial data using a free-evolution approach to 3+1 Einstein field equations based on the hyperbolic conformal formulations BSSNOK [93][94][95] or Z4c [96][97][98]. The latter is used for all of the newly released data (Ref. [65] also uses BSSNOK). The (1+log)-lapse and gamma-driver shift conditions are used for the gauge sector. The general relativistic hydrodynamics is solved in first-order conservative form [99]. Wave extraction is typically performed on coordinate spheres at finite radius placed in the wave zone of the computational domain (typically R ∼ 500 − 1000 M ) and calculating the Weyl pseudoscalar Ψ 4 , see e.g. [100] for details.
Simulations are performed with two independent mesh-based codes: BAM [100, 101] and THC [102][103][104], both developed and maintained within our collaboration. These codes employ adaptive mesh refinement (AMR) techniques in which the domain consists of a hierarchy of nested Cartesian grids (refinement levels). The grid spacing of each refinement level in each direction is half the grid spacing of its surrounding coarser refinement level. Finite difference stencils are used for the spatial discretization of the metric variables (usually at fourth order accuracy), and high resolution shock-capturing methods for the hydrodynamics. The Berger-Oliger or Berger-Colella algorithm is employed during the explicit mesh evolution. The latter is performed with the method of lines and Runge-Kutta schemes of third or fourth order accuracy in time. The innermost levels move dynamically during the time evolution following the motion of the NS such that the strong field region around a NS is always covered with the highest resolution. Both codes employ a hybrid OpenMP/MPI parallelization strategy and show good parallel scaling up to thousands of cores.
BAM implements high-order finite-differencing WENO schemes [19] and, more recently, an entropy-flux-limited (EFL) scheme [65], that is better adapted to the treatment of the NS surface, to accurately simulate multiple orbits and GWs from inspiral-mergers. The typical grid configurations for these simulations consist of seven refinement levels, where the innermost level split into two boxes covering each of the NSs. Standard grid parameters for resolution studies are chosen with n m ∈ [96,256] points per direction in each inner (moving) level and n ∈ [144,512] for the outer levels. The minimal grid spacing in each direction is ∆ ∼ [0.059, 0.321] M and the maximal resolution reached in the released simulation is ∆ ∼ 0.059 M . Symmetries can be imposed to reduce the computational cost of certain problems. For example, aligned-spin BNS are often simulated in bitant symmetry (z > 0 volume). The simulation parameters can vary for each simulation; the relevant ones are reported in the CoRe metadata.
THC implements both high-order finite-differencing schemes [105] and Kurganov-Tadmor-type central schemes. The latter are preferentially used with simulations with microphysics. THC can make use of microphysical EOS, and implements various neutrino transport schemes [54,106,107] (see below) and subgrid-scale treatment of turbulent mixing and dissipation (GRLES) [61,108] to accurately simulate remnants and postmerger dynamics. Most of the GRLES data in the current release employ an effective model for turbulence based on the high-resolution magnetohydrodynamics simulation of Ref. [109], where the viscosity parameter is set to ν T = mix c 2 s and c s is the sound speed of the fluid. mix is typically defined to be a function of the rest-mass density calibrated with the general-relativistic magneto-hydrodynamics simulations of [109] (see [61]). THC builds on the Cactus framework [110] and the Einstein Toolkit [111,112]. THC simulations use the Carpet adaptive mesh refinement driver for Cactus [113], which implements both vertex centered and cell-centered adaptive mesh refinement with flux correction [114,115]. The grid structure used in the THC simulations is similar to that used in BAM. The grid structure is specified by the resolution at the coarsest refinement level and at the location of the center of the neutron stars. The refinement levels on the grid hierarchy do not have to be connected and Carpet can merge different regions to create grids with complex topology. The standard resolution setup of the THC simulations uses a resolution of ∆ = 0.125 M in every direction on the finest refinement level. The maximal resolution reached in the released simulation is ∆ 0.08 M . The typical CFL is 0.125. However, an even lower CFL of 0.0625 is used on the coarsest grid to handle the gamma-driver source term in the shift evolution equation. All THC simulations included in the current release of the CoRe database use bitant symmetry.

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 P pwp (ρ) is a given piecewise politropic EOS model [116]. It prescribes also a value pwp for the specific internal energy given the rest mass density ρ, augmented with a γ-law "thermal" pressure term (usually, γ th = 1.75 [45,117]). The specific parameters we employ for the piecewise polytropic EOS mimic well-established zero-temperature EOS models [116]; 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 Refs. [54-59, 63, 66]. The finite-temperature EOS include the following models: BHBΛφ [118], BLh [119,120], BLQ [63,120], DD2 [121,122], LS220 [123], SFHo [124], SLy4/SRO [125,126].
All these EOS include neutrons, protons, nuclei, electrons, positrons, and photons as relevant thermodynamics degrees of freedom. The ALF2 [127] and BLQ EOS [63,120] are hybrid models accounting for deconfined quark matter. BHBΛφ is a hadronic model that includes Λ hyperons [118,128].
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 p = P (ρ, T, Y e ). Any relevant thermodynamical quantity is evaluated by multi-linear interpolating the tabulated values in log ρ, log T , 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 f (p) := p − P (ρ, , Y e ), where the specific internal energy is implicitly given by the temperature T [133]. Hence, each root-finding step includes another root finder for the function g(T ) = − E(T ) (see [134] for a discussion on computational efficiency and a non-standard approach based on neural networks.)

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 [106,133]. 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 [106]. This scheme splits neutrinos in an optically thick component, treated with the leakage Table 1. Weak reaction rates and references for their implementation. We use the following notation ν ∈ {ν e ,ν e , ν x } denotes a neutrino, ν x denote any heavy-lepton neutrino, N ∈ {n, p} denotes a nucleon, and A denotes a nucleus.

Reaction
Reference 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 [107]. 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.

Overview
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 gravitationalwave 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. The database contains waveforms from binaries with total masses ranging from 2.4 M to around ∼3.4 M with 45 datasets reaching mass ratios larger than q 1.4 and up to q = 2.1 [47,58,64]. EOS effect can be summarized to some extent by the quadrupolar tidal polarizability parameters Λ 1,2 [70], where larger (smaller) values of Λ i are associated to stiffer (softer) EOSs. The most compact NSs (and most massive binaries) are associated with the smallest values of Λ 1,2 (andΛ), see the right panel A NS spacetime is characterized by an infinite number of multipolar Love numbers of gravitoeletric and gravimagnetic type; Λ 2 parametrizes only the (gravitoelectric) leading order term in the Lagrangian. of Fig. 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 Sec. 2.1) quasi-circular mergers are included [46,48,60]. For aligned spins, the individual dimensionless components range in χ z ∈ [−0.25, 0.5); about 7 datasets are from simulations with precession effects [48,60]. The distribution of key parameters among the CoRe simulations is shown in Fig. 4.
Most of the CoRe waveforms are produced from quasi-circular mergers. The residual eccentricities of non-iterated quasi-circular initial data is usually e ∼ 10 −2 − 10 −1 , see the bottom right panel of Fig. 4. About 13 datasets have an initial eccentricity e 10 −3 that was reduced through an iterative procedure employing the formalism described in Sec. 2.1. A subset of waveforms refer instead to eccentric mergers with initial eccentricity values as high as ∼0.7 [106,139,140]. In particular, the simulation in the bottom panels of Fig. 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 E b (j), that are also publicly released. We illustrate this in Fig. 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 |E b | increases). Equal mass (ν = 1/4) non-spinning BBH systems merge with E b −0.12, 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 ∼ −κ T 2 /r 6 , i.e. it is stronger for larger tidal polarizabilities Λ 1,2 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 binary black hole inspiral. The binding energy at the moment of merger is |E b | ∼ 0.064, 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 Fig. 5 merges at smaller values |E b | ∼ 0.055 than the equal mass because of tidal disruption. The remnant has also larger angular momentum j ∼ 3.6 [59].
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 [141]. 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 binary black hole simulations this effect has been named as "hang-up" effect [142]. In Fig. 5, the spinning BNS withŜ = 0.1 is more bound than the non-spinning BNS at the moment of merger with E b ∼ −0.068. 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 Fig. 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 Sec. 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 m thr = k thr M TOV max , that mainly depends on the maximum mass of cold equilibria M TOV max supported by the EOS [45,143]. The recent analysis of Ref. [144], based on 227 finite-temperature EOS and CoRe data P, found that the prompt collapse mass threshold for equal-mass non-spinning BNS is well described by an EOS-insensitive threshold where C TOV max is the compactness of the maximum NS mass, and a = −3.36 ± 0.20, b = 2.35 ± 0.06. 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 q 1.4, 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 P These data are not released in the database since the waveforms are rather short and extracted at close radii. ejecta and remnants with accretion discs of mass ∼0.1 M . This prompt collapse process is mainly controlled by the incompressibility parameter of nuclear matter around the TOV maximum density [145]. A robust, EOS-insensitive criterion is not known in these conditions [58,[145][146][147][148], but tidal disruption effects are subdominant to the mass effect; they produce maximal variations from the equal-mass criterion of ∼8% [144,145].
Without prompt collapse, the evolution of a NS remnant is driven by the GWs emission of ∼10 53 erg lasting 20 milliseconds (GW-driven phase) [149,150]. 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 ∼2 − 4 kHz [27,28,[151][152][153][154]. 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 (2, 2) mode at a nearly constant frequency ω 22 (t) ≈ 2πf 2 ; the more compact and close to collapse the remnant is, the higher and more varying the ω 22 (t) 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,155]. Finite temperature and neutrinos do not produce qualitative differences, other than possibly on the time of gravitational collapse of the remnant [156]. Quantitative differences in the GW signal introduced by finite-temperature and neutrino effects are typically subdominant compared to finite-resolution uncertainties [107,157]. On the other hand, microphysics plays a crucial role in the EM counterparts and nucleosynthesis from mergers, e.g., [106,[158][159][160][161].
The remnant's signal from asymmetric binaries with mass ratio q 1.4 carries the imprint of the tidal disruption during merger [47,58,162]. 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,163,164], 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 [165,166]. Possible mechanisms to shed (part of) this energy are CFS [167,168] and one-arm instabilities [169][170][171] that would lead to potentially detectable, long GW transients at 1 kHz. Example of such waveforms are THC:0028, THC:0029, and THC:0036 [170]. Finally, CoRe data are available for multiple grid resolutions as discussed in Sec. 2 and shown by Fig. 2. Most of the newly released data contain high resolution simulations with a minimum grid spacing as low as ∆ ∼ 0.06 M , e.g., the NS are resolved with a uniform mesh of spacing ∼88.4 meters. 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
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,103,172]. Other aspects are relevant for waveform modelling, 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 Refs. [9,10,19,65,[103][104][105]172]. We give below in Sec. 4.1 a complete example of error analysis of a ∼10 orbit inspiral-merger waveform.
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 [174,175,180]. 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) [176,177,180]. In practice, one proceeds by defining the faithfulness between two waveforms where t 0 and φ 0 are a reference initial time and phase, and its complementary, the unfaithfulness,F := 1 − F. 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 [180] where ρ is the SNR and 2 1. This condition is necessary for unbiased parameter estimation (faithful waveforms); its violation does not imply that an analysis has biases [172,177,181,182]. 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,172,182]. We will use the faithfulness measure in Sec. 4.2 to discuss the average accuracy of the data of the CoRe database.

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 22 at the lowest extraction radius available for this simulation, R = 700 M , its amplitude |Rh 22 | and frequency M ω 22 . Note that in this section we use R instead of D L .
In order to test self-convergence, we compare amplitude and phase differences of Rψ 22 between the different resolutions. For this case, we consider the simulation at resolutions n m = 120, 160, 240 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 SF [172], where ∆ x 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 SF(2) = 1.4, 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,172]. 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.
To assess the uncertainties originated from the waveform obtained at finite-extraction radii, R i , we compare the phase differences between consecutive radii [172] and similarly for the relative amplitudes, ∆ * A 22 /A 22 . In Fig. 8 we show the differences at the extraction radii R = 700, 750, 800, 850, 900 M . 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 ∼10 −4 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 1/R of order K [172] or the method outline in [183]. The two methods give comparable results; the former is more general and can be applied to the curvature multipoles ψ m , 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 (CCE) [184] and to simulate the waveform at future null infinity. This technique was used for some of the CoRe data. The total error budget can be computed as the sum in quadrature of the truncation and finite extraction errors, and it is shown in Fig. 9 for both the curvature and strain (2,2) modes. As mentioned above, the truncation phase error is typically a factor ∼2 larger than the finite extraction error (for R 500M ) at merger and in simulations with tens of orbits.
Finally, we obtain the unfaithfulnessF of the waveforms between the different resolutions (M-H and L-M). The Wiener integral is evaluated in the frequency range f ∈ [f min , f mrg ] and employing the Advanced LIGO PSD P1200087 [185] from bajes [186]. Here f min corresponds to the initial GW frequency, and f mrg to the frequency at the moment of merger. For the faithfulness threshold F thr in Eq. (20), we consider 2 = 1 as the strict requirement, and 2 = 6, corresponding to the number of intrinsic parameters of a BNS. Similarly to [65], the SNR values are chosen to be ρ = 14, 30, 80. Figure 10 shows the computed values. The smallest unfaithfulness (M-H, n m = 240, 160) passes five out of the six accuracy tests, whereas the other one (L-M, n m = 160, 120) passes only two, namely F 14,6 thr and F 30,6 thr . However, the unfaithfulness value lies closely (or on top) of the threshold F 14,1 thr . Analyses similar to the one above are necessary to determine the quality of the NR data for GW modelling. 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,81,98,101,172]. Achieving waveform convergence in long-term evolutions of BNS is a nontrivial result and, in our experience, requires at least fifth order finitedifferencing schemes or finite volume schemes with fifth order reconstructions (at the current resolutions) [19,65,104]. Second [19], approximately third [104] and clear fourth order convergence [65] has been demonstrated up to merger in some data using these finite-differencing conservative schemes. Extreme mass ratios q ∼ 2 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 ( , m) = (2, 1), (3,3) and (4,4). 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 Sec. 5. We refer the reader to Ref. [33,38] for recent work on the 160 180 200 220 240 n accuracy of CoRe postmerger waveform.

Faithfulness analysis
In an attempt to give an overview of the accuracy of the waveform database, we compute the unfaithfulness of the (2,2) mode waveforms h 22 between the highest and second highest resolutions, for the whole database. We use again the zero-detuned, high-power Advanced LIGO PSD [185]. The minimum frequency f min employed in the integral of Eq. (18) corresponds to the initial frequency of each individual simulation. The result of this analysis is summarized in Fig. 11, whereF 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 haveF < 0.01. Several waveforms from multiple-orbits haveF 10 −4 ; 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, smallF is not necessarily an indication that these simulations are suitable for waveform modelling.
A faithfulness analysis for postmerger signals was recently presented by some of us in [33,38]. There, we found average mismatches of ∼0.01 − 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.

Quasi-universal relations
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,150,153,[187][188][189]. 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,[153][154][155]. These relations have some practical use in GW astronomy because they deliver accurate estimates for the peak luminosity [53,150] and for the remnant properties [190][191][192] (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, A mrg 22 /M and M f mrg 22 /ν, and update the fits developed in Ref. [33,38,188]. Following closely the fitting procedure of Ref. [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 X = 1 − 4ν, spin corrections in terms ofŜ, and tidal effects in terms of κ T 2 . The first two factors are given by the linear polynomial expressions Q M = 1 + a M 1 X, and Q S = 1 + p S 1Ŝ , with p S 1 = a S 1 (1 + b S 1 X). The last factor is instead a rational polynomial . The best fit parameters are shown in Tab. 2. The amplitude and frequency have 1σ errors of 2.6% and 4.6% respectively. We also obtain a χ 2 of ∼ 0.126 for the former and ∼ 0.329 for the latter.
Next, we use the public CoRe data on the emitted GW energy and extract the peak luminosity L peak using Eq. (13). For binary black holes, this quantity does not depend on the mass scale and it is accurately described by the fits of Ref. [193]. For BNS, it has been studied in Ref. [150]. We propose the ansatz where L BBH peak are the mass and spin dependent fits from Ref. [193] and p k (ν,Ŝ) = p k1 (Ŝ)ν + p k2 (Ŝ)ν 2 + p k3 (Ŝ)ν 3 p kj (Ŝ) = p kj0Ŝ + p kj1 .
Note the scaling factor 1/ν for L peak . By construction, the fit reduces to the BBH case for κ T 2 → 0. The luminosity peak is calculated in geometric units; the conversion factor to CGS units is given by the Planck luminosity L P = c 5 /G ≈ 3.63 × 10 59 erg s −1 . Figure 12 shows the best fit for L peak /ν and the CoRe data; the best fitting coefficients are reported in Tab. 4. The average 1σ deviation is about 12% over the entire dataset with less than a dozen of outliers. The peak luminosities for q ∼ 2 BNS are the least accurately modelled (4 BNS configurations). The figure shows that the largest peak luminosities are reached by BNS with κ T 2 80 that correspond to high-mass binaries and prompt collapse mergers. These events can reach peak luminosities of ∼10 55 erg s −1 , about an order of magnitude less than binary black holes (of any mass). BNS mass ratios q 1.5 can lower L peak of about an order of magnitude, while spins of magnitude ∼0.1 do not significantly affect L peak . 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 [149,150] (see also Fig. 5). We can set an upper limit to the total radiated GW energy from our dataset, obtaining E tot GW 0.676 M c 2 . Finally, we illustrate the use of CoRe data to model postmerger GWs by discussing a fit of the postmerger's spectrum peak frequency f 2 , e.g. [28,29,33,38]. This peak frequency is a robust feature found in all NR simulations. Direct GW inference on f 2 can be used to constrain NS properties [32,34,154,155,190,191,194]. 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 [195]. Following Ref. [38] we employ again Eq. (23) to fit the mass-rescaled M f 2 . The best fitting coefficients are presented in Table 2 and have a χ 2 ∼ 0.07. Figure 13 shows M f 2 as a function of κ T 2 for selected values of mass ratio and spin. The 1σ error is below 4%; this precision is in principle sufficient for informative measurements of the NS mass-radius sequence. For example, using the EOS-insensitive relation between f 2  and the maximum density of an equilibrium non-rotating NS put forward in [155], it would be possible to determine the maximum density of an equilibrium non-rotating NS to ∼15% and the maximum mass M TOV max to ∼12% with a single signal at the detectability threshold.
As a further illustration, we calibrate the EOS-insensitive relations (mass-rescaled) between f 2 and the NS radius [38,196,197] where R X is the equilibrium radius corresponding to a NS with mass X = 1.4, 1.8 M . Figure 14 shows M f 2 as function of R 1.4 , R 1.8 and R 1.4 /R 1.8 . Best fit parameters are given in Table 3. 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.

Conclusion
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. 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 [106,133,198], or using the M0 approach [54,106]. Turbulent viscosity is included in some models using the GRLES formalism [61,108]. 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 −2 , some have unfaithfulness of less than 10 −4 , so they are suitable for precision waveform modelling 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 (< 5% for the merger amplitude and frequency and 17% for the peak luminosity). We also recalibrated various quasiuniversal 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 numerical relativity 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.