Second MAYA Catalog of Binary Black Hole Numerical Relativity Waveforms

Gravitational wave observations have become increasingly frequent, with the LIGO Scientific, Virgo and KAGRA Collaborations (LVK) reporting 90 detections through the end of their third observing run, including 83-85 likely binary black holes (BBHs) LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration (2021a); Abbott et al. (2021a, b, 2019), and more are expected in the fourth observing run. The BBH systems observed by the LVK span a large range of parameters, including mass ratios from equal mass to the possible 9:1 mass ratio of GW190814 LIGO Scientific Collaboration and Virgo Collaboration (2020). While most of the events are consistent with nonspinning black holes (BHs), some show evidence of strong precession LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration (2021a), and some events even show signs of eccentricity Gayathri et al. (2022); Romero-Shaw et al. (2021). As we prepare for next-generation gravitational-wave detectors such as the Laser Interferometer Space Antenna (LISA) Baker et al. (2019); Amaro-Seoane et al. (2017); Robson et al. (2019), the Einstein Telescope (ET) Punturo et al. (2010); Abernathy et al. (2011), and Cosmic Explorer (CE) Evans et al. (2021), we anticipate a significant increase in the volume of signals as well as a more diverse parameter space. Understanding the properties of the systems these signals come from and their underlying populations LIGO Scientific Collaboration and Virgo Collaboration (2021) and using the signals to test General Relativity (GR) LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration (2021b) relies upon having accurate predictions of the anticipated signals.

Due to the complex nature of Einstein’s theory of GR, analytic solutions don’t exist for the full two body problem of merging compact objects. A number of techniques have been developed to solve for the dynamics and radiation of such systems, and the most appropriate technique depends upon the parameters of the binary system in question. When the objects are far apart, post-Newtonian (PN) approximations can be used with high accuracy and confidence Isoyama et al. (2020); Blanchet (2014); Schäfer and Jaranowski (2018); Futamase and Itoh (2007); Blanchet et al. (2008); Faye et al. (2012, 2015). For systems where the masses of the black holes are highly unequal, small mass ratio approximations allow for efficient simulations Hinderer and Flanagan (2008); Miller and Pound (2021); Barack and Pound (2019). Great progress has also been made towards using small mass ratio approximations for the inspiral of near equal mass ratio systems van de Meent and Pfeiffer (2020); Wardell et al. (2021); Warburton et al. (2021); Albertini et al. (2022a, b).

In the highly nonlinear regime of the merger of comparably massed compact objects, numerical relativity (NR) is the most accurate approach. NR computationally solves Einstein’s equations by expressing them in a 3+1 formalism and separating the constraint and evolution equations in order to create an initial value problem that can be evolved on a computer Baumgarte and Shapiro (1998a, 2010, b). For the case of BBHs in vacuum, there is no missing physics within the simulations, as long as GR is correct, and the only limitations are computational in nature. A NR simulation will solve Einstein’s equations for a single binary system described by the masses and spin vectors of the two component BHs. For BBHs in vacuum, the results of the simulation can be scaled by the total mass, so we set the total ADM mass of the system to 1, and the masses can then simply be described by the mass ratio, $q=m_{1}/m_{2}\geq 1$. By creating template banks of NR simulations, NR waveforms have been used directly in gravitational wave (GW) detection and parameter estimation Lange et al. (2017); Abbott et al. (2016); Gayathri et al. (2022); Calderon Bustillo et al. (2022).

These simulations are placed discretely throughout the parameter space and are time consuming and computationally expensive. Therefore, while NR waveforms can be used directly as template banks for GW analysis, for many analyses, analytic or semi-analytic models are created, using information from NR for calibration Varma et al. (2019a); Bohé et al. (2017); Khan et al. (2016); Blackman et al. (2017); Husa et al. (2016); Taracchini et al. (2014); Hannam et al. (2014); Santamaria et al. (2010); Nakano et al. (2011); Pratten et al. (2021). However, creating reliable models and template banks relies upon NR waveforms being sufficiently accurate and densely covering the complete parameter space Ferguson et al. (2021); Pürrer and Haster (2020).

There are several NR codes and public waveform catalogs created by the NR community. The SXS collaboration has a public waveform catalog generated using their SpEC code Mroué et al. (2013); Boyle et al. (2019) and have been developing their next-generation code, SpECTRE Deppe et al. (2021). There are several waveform catalogs generated with codes derived from an open source software used by many NR groups Loffler et al. (2012), the Einstein Toolkit Jani et al. (2016); Healy et al. (2017a, 2019); Healy and Lousto (2022, 2020). The waveforms in the catalog presented in this paper are also generated with a code derived from the Einstein Toolkit, called Maya.

We previously released the Georgia Tech waveform catalog Jani et al. (2016) using our Maya code including 452 waveforms covering the spinning parameter space up to $q=8$ and including nonspinning waveforms up to $q=15$. These waveforms have been used to study exceptional LVK events Lange et al. (2017), the detectability of intermediate-mass BHs Jani et al. (2019); Chandra et al. (2020); Abbott et al. (2022), as well as post-merger chirps Calderón Bustillo et al. (2019), to highlight a few recent studies. They have also been used to study the different stages of the merger, including finding a relationship between the frequency at merger and the spin of the remnant BH Ferguson et al. (2019); Healy et al. (2014). Using the Georgia Tech catalog, it was also found that the BH recoil can be measured from GW observations Calderón Bustillo et al. (2018), a method which has been applied to real events Calderón Bustillo et al. (2022). This catalog has since been rebranded to the MAYA Catalog.

This paper introduces the second MAYA Catalog, expanding beyond the parameter space coverage of the first catalog and introducing a new format. This format can be read using the mayawaves python library Ferguson et al. (2023a, b), which is described further in Section  II.2. The new catalog includes a suite of eccentric simulations as well as improved coverage of the high mass ratio parameter space and the precessing parameter space. In Section  II, we describe the Maya code used to create this catalog as well as the mayawaves python library. We then dive into the description of the new catalog in Section  III. Section  IV compares our simulations to other NR waveforms and waveform models, including an eccentric model, and Section  V compares the models for the remnant BH parameters to the results obtained from our catalog. Finally, in Section  VI, we include a discussion of the impact of center-of-mass (COM) drift corrections.

In this section, we’ll describe the code used to perform these NR simulations (Maya) as well as the python library used to interact with and analyze the simulations (mayawaves).

The MAYA catalog of NR waveforms is accessible at https://cgp.ph.utexas.edu/waveform. At that location, you will see a table with the metadata describing each simulation in the catalog. From that table, you can download each of the simulations. They are stored in an h5 format that is designed to be read using the mayawaves python library. Any simulations that contain all the necessary data are also provided in the format detailed in  Schmidt et al. (2017) for use with PyCBC LIGO Scientific Collaboration (2018); Nitz et al. (2020).

This catalog introduces 181 new simulations that expand our coverage of the BBH parameter space, for a total of 635 simulations. The catalog website has also been updated to include the original catalog Jani et al. (2016) in the new mayawaves compatible format. Figure  2 shows the coverage of the previous catalog and this catalog, highlighting the many more eccentric simulations, precessing simulations, and high mass ratio simulations introduced by this catalog. Figure  3 shows the spin coverage of the combined catalogs. Figure  4 shows the distribution of the number of inspiral GW cycles after junk radiation for each simulation in the original and new catalogs. Simulations in the new catalog have a median of $\sim 16.5$ GW cycles and the overall catalog has a median of $\sim 11$ GW cycles.

This catalog is an amalgamation of several suites of simulations performed for various studies. This includes a suite of 80 non-precessing, eccentric simulations with mass ratios up to $q=4$. It includes 38 simulations for a study of secondary spin for mass ratios up to $q=15$. We also include 9 simulations performed for LVK followup studies. Finally we have 31 simulations placed to optimally fill in gaps in our parameter space using the network and method described in  Ferguson (2022). Each of these suites is described in more detail below. The remaining 23 simulations are an assortment of simulations performed over the last few years.

Within the NR community, there are several different methods and techniques for solving Einstein’s equations directly. While each method has its benefits and drawbacks, all NR waveforms are considered to the be the gold standard for studying the merger of comparable-mass BBH systems. To ensure accuracy and consistency in analyses that use NR waveforms of each technique, it is important to perform cross-code comparisons. Several of these have been done over the past couple decades, and the methods used by the Einstein Toolkit and the Simulating eXtreme Spacetimes (SXS) Collaboration have proven to be reliable Ajith et al. (2012); Hinder et al. (2013); Healy et al. (2018). Here, we provide a comparison between this catalog and comparable waveforms available in the SXS catalog.

Similarly, there are many analytic and semi-analytic models for computing GWs. These models have the benefit of being faster than NR allowing for continuous sampling of the parameter space. However, given that these models are trained or calibrated to NR waveforms, they are generally not considered to be as accurate as NR waveforms. In this section we will also perform a comparison between this catalog and waveforms generated using such models.

For each of these comparisons, we compute the mismatch on a flat noise curve, defined as:

where

with $S_{n}$ being a flat PSD, and $*$ denoting the complex conjugate. Mismatches are very sensitive to systematic effects such as the amount of tapering and the starting frequency, and we find that these effects can change the mismatches by up to $\sim 10^{-3}$. For consistency and to reduce Gibbs oscillations, we consider only waveforms with more than 10 inspiral cycles and taper 6 of them. We compute all mismatches using a total mass of $M_{tot}=100M_{\odot}$ and a starting frequency of 30 Hz.

For the SXS comparison, we consider all equivalent non-precessing systems that exist in our catalog and in the SXS LVK catalog. Figure  5(a) shows the mismatch as a function of inclination, using all available modes and a flat noise curve. The grey lines show each individual system and the black line shows the median. Figure  5(b) shows the distribution of mismatches (including all 20 inclinations) with a median mismatch of $10^{-3}$.

We also compare our waveforms to the current state-of-the-art waveform models. For the quasicircular case, we compare to IMRPhenomXPHM Pratten et al. (2021) for both precessing and non-precessing systems, and for the eccentric case, we compare to TEOBResumS-DALI Chiaramello and Nagar (2020); Nagar et al. (2021); Nagar and Rettegno (2021); Placidi et al. (2022). Figure  6(a) shows the mismatch between the quasicircular NR waveforms and the IMRPhenomXPHM waveforms as a function of inclination, including all available modes on a flat noise curve. Each of the precessing systems appears as a faint orange line with the vibrant orange line showing the median value for all precessing simulations. All non-precessing simulations appear as faint blue lines with the vibrant blue line showing the median value for all non-precessing simulations. The black line shows the median mismatch of all quasicircular simulations. Figure  6(b) shows the distribution of the mismatches for precessing (blue) and aligned (orange) systems including 20 uniformly distributed inclinations. For non-precessing systems, the median mismatch is $\sim 10^{-3}$, and for the precessing systems, the median mismatch is $\sim 2\times 10^{-2}$.

For the eccentric case, we use TEOBResumS-DALI as the waveform model and consider only aligned spin systems. We optimize over the initial eccentricity and frequency of the model waveforms to obtain the best mismatch, since the model’s eccentricity parameter is not physically meaningful. Figure  7 shows the mismatches between TEOBResumS-DALI waveforms and 20 of our waveforms which fall within the validity range of the model ($e\leq 0.2$, $q\leq 3$, and $a_{1},a_{2}\leq 0.7$). The median mismatch over all systems and inclinations is $\sim 5\times 10^{-3}$.

In this section, we compare the properties of the remnant BH computed from NR simulations to those predicted using NRSur7dq4Remnant Varma et al. (2019a). We use the SurfinBH python library  Varma et al. (2019b, 2018) to interact with NRSur7dq4Remnant and compute the predicted remnant properties. We provide the masses and spins of the initial BHs at a time $100M$ before merger, rotated into the same frame used by NRSur7dq4Remnant. We only consider systems which fall within the range of validity for NRSur7dq4Remnant ($q\leq 4$ and dimensionless spin parameters $a_{1},a_{2}\leq 0.8$). The masses and spins of the initial and remnant BHs are computed from the BH apparent horizons, and the magnitude of the kick velocity is computed from the linear momentum radiated by the GWs.

Figure  8 shows the remnant properties obtained from NR simulations compared to the remnant properties obtained from NRSur7dq4Remnant. For each of the remnant properties, the non-precessing systems have much lower errors than the precessing systems. However, for remnant spin, all residuals fall within $~{}0.0125$, and for remnant mass, all residuals fall within $~{}0.0012$. The recoil velocity is more challenging since it is very sensitive to spin angles. The system with the largest residual has an error of $29\%$.

During the NR simulation of the inspiral of BBH systems, the center of mass experiences a wobble and a drift. These are caused by a combination of physical effects due to the asymmetrical emission of GWs and numerical error. Since the GWs are decomposed with the spherical harmonics centered on the initial COM, this can cause mode mixing. In recent years there has been much discussion of the impact that these COM drifts can cause to the waveform modes. We include here a brief analysis of the impact of correcting for such drifts.

Using the Scri python package Boyle et al. (2020); Boyle (2013); Boyle et al. (2014); Boyle (2016), we apply corrections for a translation and boost in the COM, computed using equations 6 and 7 of  Woodford et al. (2019). We then perform a mismatch analysis comparing the COM corrected waveforms and the uncorrected waveforms for the overall strain and for individual modes. Figure  9 shows the mismatches between raw and COM corrected strains. In Figure  9(a), we see the mismatch with all modes combined as a function of inclination. The median mismatch is fairly constant for all inclinations, with a value of $~{}\sim 4\times 10^{-5}$. In Figure  9(b), we see the mismatch values for several modes as a function of mass ratio. While there does not appear to be a strong correlation between the impact of COM corrections and mass ratio, the figure does reveal that certain modes are more impacted than others. In particular, the $\ell=4$, $m=2$ mode has mismatches up to $\sim 10^{-2}$.

Using the mayawaves library, our waveforms can be read in their raw form or with COM corrections. The waveforms stored in the PYCBC compatible format LIGO Scientific Collaboration (2018); Nitz et al. (2020); Schmidt et al. (2017) have been corrected for COM drift.

This paper introduces the second MAYA catalog of NR waveforms. It contains 181 waveforms that fill in and expand our coverage of the BBH parameter space. This includes 80 from an eccentric suite ranging up to eccentricities of 0.1, as well as 7 systems with both precession and eccentricity. It also includes 38 waveforms studying secondary spin, 9 simulations as NR followup for LVK events, and 31 simulations placed to optimize parameter space coverage.

As the most accurate way of studying the merger of comparably massed objects, NR waveforms are crucial for understanding the observations of current and next-generation GW detectors. We must, therefore, ensure we have sufficient NR coverage of the anticipated parameter space. This catalog makes significant strides towards this goal by pushing and filling in considerable gaps in existing NR coverage of the BBH parameter space. With mass ratios up to $q=15$, highly precessing simulations, and an extensive eccentric suite, this catalog greatly improves the coverage of public NR catalogs.

By performing a thorough error analysis and comparison to other NR codes, we are confident in the accuracy and reliability of the waveforms presented in this catalog. Through a comparison to state-of-the-art GW models, we show consistency but also highlight the continuing need for NR waveforms particularly in the less well understood precessing space.

The waveforms are provided at https://cgp.ph.utexas.edu/waveform and can be read using the mayawaves python package. We have also updated the previous Georgia Tech catalog to include more data by using the new format, leading to a total of 635 public waveforms. All waveforms that fit the requirements are also provided in the format described in  Schmidt et al. (2017) for use with PYCBC.