The third Gravitational-Wave Transient Catalog (GWTC-3) contains 90 binary coalescence candidates detected by the LIGO-Virgo-KAGRA Collaboration (LVK). We provide a re-analysis of binary black hole (BBH) events using a recently developed numerical relativity (NR) waveform surrogate model, NRSur7dq4, that includes all $\ell \leq 4$ spin-weighted spherical harmonic modes as well as the complete physical effects of precession. Properties of the remnant black holes' (BH's) mass, spin vector, and kick vector are found using an associated remnant surrogate model NRSur7dq4Remnant. Both NRSur7dq4 and NRSur7dq4Remnant models have errors comparable to numerical relativity simulations and allow for high-accuracy parameter estimates. We restrict our analysis to 47 BBH events that fall within the regime of validity of NRSur7dq4 (mass ratios greater than 1/6 and total masses greater than $60 M_{\odot}$). While for most of these events our results match the LVK analyses that were obtained using the semi-analytical models such as IMRPhenomXPHM and SEOBNRv4PHM, we find that for more than 20\% of events the NRSur7dq4 model recovers noticeably different measurements of black hole properties like the masses and spins, as well as extrinsic properties like the binary inclination and distance. For instance, GW150914_095045 exhibits noticeable differences in spin precession and spin magnitude measurements. Other notable findings include one event (GW191109_010717) that constrains the effective spin $\chi_{eff}$ to be negative at a 99.3\% credible level and two events (GW191109_010717 and GW200129_065458) with well-constrained kick velocities. Furthermore, compared to the models used in the LVK analyses, NRSur7dq4 recovers a larger signal-to-noise ratio and/or Bayes factors for several events.

Surrogate models of numerical relativity simulations of merging black holes provide the most accurate tools for gravitational-wave data analysis. Neural network-based surrogates promise evaluation speedups, but their accuracy relies on (often obscure) tuning of settings such as the network architecture, hyperparameters, and the size of the training dataset. We propose a systematic optimization strategy that formalizes setting choices and motivates the amount of training data required. We apply this strategy on NRSur7dq4Remnant, an existing surrogate model for the properties of the remnant of generically precessing binary black hole mergers and construct a neural network version, which we label NRSur7dq4Remnant_NN. The systematic optimization strategy results in a new surrogate model with comparable accuracy, and provides insights into the meaning and role of the various network settings and hyperparameters as well as the structure of the physical process. Moreover, NRSur7dq4Remnant_NN results in evaluation speedups of up to 8 times on a single CPU and a further improvement of 2,000 times when evaluated in batches on a GPU. To determine the training set size, we propose an iterative enrichment strategy that efficiently samples the parameter space using much smaller training sets than naive sampling. NRSur7dq4Remnant_NN requires O(10^4) training data, so neural network-based surrogates are ideal for speeding-up models that support such large training datasets, but at the moment cannot directly be applied to numerical relativity catalogs that are O(10^3) in size. The optimization strategy is available through the gwbonsai package.

We investigate the late-time tail behavior in gravitational waves from merging eccentric binary black holes (BBH) using black hole perturbation theory. For simplicity, we focus only on the dominant quadrupolar mode of the radiation. We demonstrate that such tails become more prominent as eccentricity increases. Exploring the phenomenology of the tails in both spinning and nonspinning eccentric binaries, with the spin magnitude varying from = -0.6 to = & thorn;0.6 and eccentricity as high as e = 0.98, we find that these tails can be well approximated by a slowly decaying power law. We study the power law for varying systems and find that the power law exponent lies close to the theoretically expected value -4. Finally, using both plunge geodesic and radiation-reaction-driven orbits, we perform a series of numerical experiments to understand the origin of the tails in BBH simulations. Our results suggest that the late-time tails are strongly excited in eccentric BBH systems when the smaller black hole is in the neighborhood of the apocenter, as opposed to any structure in the strong field of the larger black hole. Our analysis framework is publicly available through the GWTAILS Python package.

gwNRHME is a recently developed framework that seamlessly converts a multimodal (i.e., with several spherical harmonic modes) quasicircular waveform into a multimodal eccentric waveform if the quadrupolar eccentric waveform is known. Here, we employ the gwNRHME framework to combine a multimodal quasicircular numerical relativity surrogate waveform model NRHybSur3dq8 and quadrupolar nonspinning post-Newtonian eccentric waveform model EccentricIMR to construct multimodal nonspinning eccentric model NRHybSur3dq8-gwNRHME. Using a total of 35 eccentric numerical relativity (NR) simulations obtained from the SXS and RIT catalogs, we demonstrate that NRHybSur3dq8-gwNRHME model predictions agree well with NR (with typical relative L2 errors of similar to 0.01 for the dominant quadrupolar mode) for mass ratios 1 <= q <= 4 and eccentricities up to similar to 0.2 measured about 10 cycles before the merger. Our frequency-domain mismatches (calculated assuming advanced LIGO design sensitivity curve) are mostly below 0.01. To demonstrate the modularity of the gwNRHME framework, we further combine EccentricIMR with the BHPTNRSur1dq1e4BHPTNRSur1dq1e4 model and develop a nonspinning eccentric model named BHPTNRSur1dq1e4-gwNRHME. Finally, we develop a different variant of these models by replacing EccentricIMR with EccentricTD. Both the gwNRHME framework and associated models are available through the gwModels package.

LISA, the Laser Interferometer Space Antenna, will usher in a new era in gravitational-wave astronomy. As the first anticipated space-based gravitational-wave detector, it will expand our view to the millihertz gravitational-wave sky, where a spectacular variety of interesting new sources abound: from millions of ultra-compact binaries in our Galaxy, to mergers of massive black holes at cosmological distances; from the beginnings of inspirals that will venture into the ground-based detectors' view to the death spiral of compact objects into massive black holes, and many sources in between. Central to realising LISA's discovery potential are waveform models, the theoretical and phenomenological predictions of the pattern of gravitational waves that these sources emit. This white paper is presented on behalf of the Waveform Working Group for the LISA Consortium. It provides a review of the current state of waveform models for LISA sources, and describes the significant challenges that must yet be overcome.

Highly precise and robust waveform models are required as improvements in detector sensitivity enable us to test general relativity with more precision than ever before. In this work, we introduce a spin-aligned surrogate ringdown model. This ringdown surrogate, NRSur3dq8_RD, is built with numerical waveforms produced using Cauchy-characteristic evolution. In addition, these waveforms are in the superrest frame of the remnant black hole allowing us to do a correct analysis of the ringdown spectrum. The novel prediction of our surrogate model is complex-valued quasinormal mode (QNM) amplitudes, with median relative errors of 10^-2-10^-3 over the parameter space. Like previous remnant surrogates, we also predict the remnant black hole's mass and spin. The QNM mode amplitude errors translate into median errors on ringdown waveforms of 10^-4. The high accuracy and QNM mode content provided by our surrogate will enable high-precision ringdown analyses such as tests of general relativity. Our ringdown model is publicly available through the python package surfinBH.

The anticipated observation of the gravitational microlensing of gravitational waves (GWs) promises to shed light on a host of astrophysical and cosmological questions. However, extracting the parameters of the lens from the modulated GWs requires accurate modeling of the lensing amplification factor, accounting for wave-optics effects. Analytic solutions to the lens equation have not been found to date, except for a handful of simplistic lens models. While numerical solutions to this equation have been developed, the time and computational resources required to evaluate the amplification factor numerically make large-scale parameter estimation of the lens (and source) parameters prohibitive. On the other hand, surrogate modeling of GWs has proven to be a powerful tool to accurately, and rapidly, produce GW templates at arbitrary points in parameter space, interpolating from a finite set of available waveforms at discrete parameter values. In this work, we demonstrate that surrogate modeling can also effectively be applied to the evaluation of the time-domain microlensing amplification factor F(t). We show this by constructing F(t) for two lens models, viz. point-mass lens, and singular isothermal sphere, which notably includes logarithmic divergence behavior. We find both surrogates reproduce the original lens models accurately, with mismatches less than or similar to 5 x 10-4 across a range of plausible microlensed binary black hole sources observed by the Einstein Telescope. This surrogate is between 5 and 103 times faster than the underlying lensing models, and can be evaluated in about 100 ms. The accuracy and efficiency attained by our surrogate models will enable practical parameter estimation analyses of microlensed GWs.

Past studies have empirically demonstrated a surprising agreement between gravitational waveforms computed using adiabatic–driven–inspiral point–particle black hole perturbation theory (ppBHPT) and numerical relativity (NR) following a straightforward calibration step, sometimes referred to as α-β scaling. Specifically focusing on the quadrupole mode, this calibration technique necessitates only two time-independent parameters to scale the overall amplitude and time coordinate. In this article, part of a Special Issue, we investigate this scaling for non-spinning binaries at the equal-mass limit. Even without calibration, NR and ppBHPT waveforms exhibit an unexpected degree of similarity after accounting for different mass scale definitions. Post-calibration, good agreement between ppBHPT and NR waveforms extends nearly up to the point of the merger. We also assess the breakdown of the time-independent assumption of the scaling parameters, shedding light on current limitations and suggesting potential generalizations for the α-β scaling technique.

We present a surrogate model of SEOBNRv4PHM, a fully precessing time-domain effective-one-body waveform model including subdominant modes. We follow an approach similar to that used to build recent numerical relativity surrogate models. Our surrogate is 5000 M in duration, covers mass ratios up to 1:20 and dimensionless spin magnitudes up to 0.8. Validating the surrogate against an independent test set, we find that the median mismatch error is less than 10-3, which is typically smaller than the modeling error of SEOBNRv4PHM itself. At high total mass, a few percent of configurations can exceed mismatches of 10-2 if they are highly precessing and exceed a mass ratio of 1:4. This surrogate is nearly 2 orders of magnitude faster than the underlying time-domain SEOBNRv4PHM model and can be evaluated in similar to 50 ms. Bayesian inference analyses with SEOBNRv4PHM are typically very computationally demanding and can take from weeks to months to complete. The 2 order of magnitude speedup attained by our surrogate model enables practical parameter estimation analyses with this waveform family. This is crucial because Bayesian inference allows us to recover the masses and spins of binary black hole mergers given a model of the emitted gravitational waveform along with a description of the noise.

Gravitational wave signals from extreme mass ratio inspirals are a key target for space-based gravitational wave detectors. These systems are typically modeled as a distributionally-forced Teukolsky equation, where the smaller black hole is treated as a Dirac delta distribution. Time-domain solvers often use regularization approaches that approximate the Dirac distribution that often introduce small length scales and are a source of systematic error, especially near the smaller black hole. We describe a multi-domain discontinuous Galerkin method for solving the distributionally-forced Teukolsky equation that describes scalar fields evolving on a Kerr spacetime. To handle the Dirac delta, we expand the solution in spherical harmonics and recast the sourced Teukolsky equation as a first-order, one-dimensional symmetric hyperbolic system. This allows us to derive the method's numerical flux to correctly account for the Dirac delta. As a result, our method achieves global spectral accuracy even at the source's location. To connect the near field to future null infinity, we use the hyperboloidal layer method, allowing us to supply outer boundary conditions and providing direct access to the far-field waveform. We document several numerical experiments where we test our method, including convergence tests against exact solutions, energy luminosities for circular orbits, the scheme's superconvergence properties at future null infinity, and the late-time tail behavior of the scalar field. We also compare two systems that arise from different choices of the first-order reduction variables, finding that certain choices are numerically problematic in practice. The methods developed here may be beneficial when computing gravitational self-force effects, where the regularization procedure has been developed for the spherical harmonic modes and high accuracy is needed at the Dirac delta's location.

We present BHPTNRSur2dq1e3, a reduced order surrogate model of gravitational waves emitted from binary black hole (BBH) systems in the comparable to large mass ratio regime with aligned spin (chi 1) on the heavier mass (m1). We trained this model on waveform data generated from point particle black hole perturbation theory (ppBHPT) with mass ratios varying from 3 <= q <= 1000 and spins from - 0.8 <= chi 1 <= 0.8. The waveforms are 13,500 m 1 long and include all l <= 4 spin-weighted spherical harmonic modes except the (4,1) and m = 0 modes. We find that, for binaries with chi 1 less than or similar to-0.5, retrograde quasinormal modes are significantly excited, thereby complicating the modeling process. To overcome this issue, we introduce a domain decomposition approach to model the inspiral and merger-ringdown portion of the signal separately. The resulting model can faithfully reproduce ppBHPT waveforms with a median time-domain mismatch error of 8 x 10-5. We then calibrate our model with numerical relativity (NR) data in the comparable mass regime (3 <= q <= 10). By comparing with spin-aligned BBH NR simulations at q = 15, we find that the dominant quadrupolar (subdominant) modes agree to better than approximate to 10 - 3 (approximate to 10-2) when using a time-domain mismatch error, where the largest source of calibration error comes from the transition-to-plunge and ringdown approximations of perturbation theory. Mismatch errors are below approximate to 10 - 2 for systems with mass ratios between 6 <= q <= 15 and typically get smaller at larger mass ratio. Our two models-both the ppBHPT waveform model and the NR-calibrated ppBHPT model-will be publicly available through GWS urrogate and the B lack H ole P erturbation T oolkit packages.

During binary black hole (BBH) mergers, energy and momenta are carried away from the binary system as gravitational radiation. Access to the radiated energy and momenta allows us to predict the properties of the remnant black hole. We develop a Python package gw_remnant to compute the remnant mass, remnant spin, peak luminosity, and the final kick imparted on the remnant black hole from the gravitational radiation. Using this package, we compute the remnant properties of the final black hole in case of non-spinning BBH mergers with mass ratios ranging from q=2.5 to q=1000 using waveform modes generated from BHPTNRSur1dq1e4, a recently developed numerical-relativity-informed surrogate model based on black-hole perturbation theory framework. We validate our results against the remnant properties estimated from numerical relativity (NR) surrogate models in the comparable mass ratio regime and against recently available high-mass-ratio NR simulations at q=[15,32,64]. We find that our remnant property estimates computed from fluxes at future null infinity closely match the estimates obtained from the NR surrogate model of apparent horizon data. Using Gaussian process regression fitting methods, we train a surrogate model, {BHPTNR_Remnant, for the properties of the remnant black hole arising from BBH mergers with mass ratios from q=2.5 to q=1000. Finally, we discuss potential improvements in the BHPTNRSur1dq1e4 waveform model when including remnant information. We make both the gw_remnant and BHPTNR_Remnant packages publicly available.

Numerical relativity simulations provide the most precise templates for the gravitational waves produced by binary black hole mergers. However, many of these simulations use an incomplete waveform extraction technique -- extrapolation -- that fails to capture important physics, such as gravitational memory effects. Cauchy-characteristic evolution (CCE), by contrast, is a much more physically accurate extraction procedure that fully evolves Einstein's equations to future null infinity and accurately captures the expected physics. In this work, we present a new surrogate model, NRHybSur3dq8$\_$CCE, built from CCE waveforms that have been mapped to the post-Newtonian (PN) BMS frame and then hybridized with PN and effective one-body (EOB) waveforms. This model is trained on 102 waveforms with mass ratios $q\leq8$ and aligned spins $\chi_{1z}, \, \chi_{2z} \in \left[-0.8, 0.8\right]$. The model spans the entire LIGO-Virgo-KAGRA (LVK) frequency band (with $f_{\text{low}}=20\text{Hz}$) for total masses $M\gtrsim2.25M_{\odot}$ and includes the $\ell\leq4$ and $(\ell,m)=(5,5)$ spin-weight $-2$ spherical harmonic modes, but not the $(3,1)$, $(4,2)$ or $(4,1)$ modes. We find that NRHybSur3dq8$\_$CCE can accurately reproduce the training waveforms with mismatches $\lesssim2\times10^{-4}$ for total masses $2.25M_{\odot}\leq M\leq300M_{\odot}$ and can, for a modest degree of extrapolation, capably model outside of its training region. Most importantly, unlike previous waveform models, the new surrogate model successfully captures memory effects.

The non-linear gravitational wave (GW) memory effect is a distinct prediction in general relativity. While the effect has been well studied for comparable mass binaries, it has mostly been overlooked for intermediate mass ratio inspirals (IMRIs). We offer a comprehensive analysis of the phenomenology and detectability of memory effects, including contributions from subdominant harmonic modes, in heavy IMRIs consisting of a stellar mass black hole and an intermediate mass black hole. When formed through hierarchical mergers, for example when a GW190521-like remnant captures a stellar mass black hole, IMRI systems have a large total mass, large spin on the primary, and possibly residual eccentricity; features that potentially raise the prospect for memory detection. We compute both the displacement and spin non-linear GW memory from the $m \neq 0$ gravitational waveforms computed within a black hole perturbation theory framework that is partially calibrated to numerical relativity waveforms. We probe the dependence of memory effects on mass ratio, spin, and eccentricity and consider the detectability of a memory signal from IMRIs using current and future GW detectors. We find that (i) while eccentricity introduces additional features in both displacement and spin memory, it does not appreciatively change the prospects of detectability, (ii) including higher modes into the memory computation can increase singal-to-noise (SNR) values by about 7\% in some cases, (iii) the SNR from displacement memory dramatically increases as the spin approaches large, positive values, (iv) spin memory from heavy IMRIs would, however, be difficult to detect with future generation detectors even from highly spinning systems. Our results suggest that hierarchical binary black hole mergers may be a promising source for detecting memory and could favorably impact memory forecasts.

The final black hole left behind after a binary black hole merger can attain a recoil velocity, or a "kick," reaching values up to 5000  km/s. This phenomenon has important implications for gravitational wave astronomy, black hole formation scenarios, testing general relativity, and galaxy evolution. We consider the gravitational wave signal from the binary black hole merger GW200129_065458 (henceforth referred to as GW200129), which has been shown to exhibit strong evidence of orbital precession. Using numerical relativity surrogate models, we constrain the kick velocity of GW200129 to v_{f}∼1542_{-1098}^{+747}  km/s or v_{f}≳698  km/s (one-sided limit), at 90% credibility. This marks the first identification of a large kick velocity for an individual gravitational wave event. Given the kick velocity of GW200129, we estimate that there is a less than 0.48% (7.7%) probability that the remnant black hole after the merger would be retained by globular (nuclear star) clusters. Finally, we show that kick effects are not expected to cause biases in ringdown tests of general relativity for this event, although this may change in the future with improved detectors.

We present a reduced-order surrogate model of gravitational waveforms from non-spinning binary black hole systems with comparable to large mass-ratio configurations. This surrogate model, \texttt{BHPTNRSur1dq1e4}, is trained on waveform data generated by point-particle black hole perturbation theory (ppBHPT) with mass ratios varying from 2.5 to 10,000. \texttt{BHPTNRSur1dq1e4} extends an earlier waveform model, \texttt{EMRISur1dq1e4}, by using an updated transition-to-plunge model, covering longer durations up to 30,500 $m_1$ (where $m_1$ is the mass of the primary black hole), includes several more spherical harmonic modes up to $\ell=10$, and calibrates subdominant modes to numerical relativity (NR) data. In the comparable mass-ratio regime, including mass ratios as low as $2.5$, the gravitational waveforms generated through ppBHPT agree surprisingly well with those from NR after this simple calibration step. We also compare our model to recent SXS and RIT NR simulations at mass ratios ranging from $15$ to $32$, and find the dominant quadrupolar modes agree to better than $\approx 10^{-3}$. We expect our model to be useful to study intermediate-mass-ratio binary systems in current and future gravitational-wave detectors.

Linear wave equations sourced by a dirac delta distribution δ(x) and its derivative(s) can serve as a model for many different phenomena. We describe a discontinuous Galerkin (DG) method to numerically solve such equations with source terms proportional to ∂nδ∕∂xn. Despite the presence of singular source terms, which imply discontinuous or potentially singular solutions, our DG method achieves global spectral accuracy even at the source’s location. Our DG method is developed for the wave equation written in fully first-order form. The first-order reduction is carried out using a distributional auxiliary variable that removes some of the source term’s singular behavior. While this is helpful numerically, it gives rise to a distributional constraint. We show that a time-independent spurious solution can develop if the initial constraint violation is proportional to δ(x). Numerical experiments verify this behavior and our scheme’s convergence properties by comparing against exact solutions.

We develop and use a novel mixed-precision weighted essentially non-oscillatory (WENO) method for solving the Teukolsky equation, which arises when modeling perturbations of Kerr black holes. We show that WENO methods outperform higher-order finite-difference methods, standard in the discretization of the Teukolsky equation, due to the need to add dissipation for stability purposes in the latter. In particular, as the WENO scheme uses no additional dissipation, it is well suited for scenarios requiring long-time evolution such as the study of price tails and gravitational wave emission from extreme mass ratio binaries. In the mixed-precision approach, the expensive computation of the WENO weights is performed in reduced floating-point precision that results in a significant speedup factor of ≈ 3.3 . In addition, we use state-of-the-art Nvidia general-purpose graphics processing units and cluster parallelism to further accelerate the WENO computations. Our optimized WENO solver can be used to quickly generate accurate results of significance in the field of black hole and gravitational wave physics. We apply our solver to study the behavior of the Aretakis charge—a conserved quantity, that if detected by a gravitational wave observatory like LIGO/Virgo would prove the existence of extremal black holes.

We develop new strategies to build numerical relativity surrogate models for eccentric binary black hole systems, which are expected to play an increasingly important role in current and future gravitational-wave detectors. We introduce a new surrogate waveform model, \texttt{NRSur2dq1Ecc}, using 47 nonspinning, equal-mass waveforms with eccentricities up to $0.2$ when measured at a reference time of $5500M$ before merger. This is the first waveform model that is directly trained on eccentric numerical relativity simulations and does not require that the binary circularizes before merger. The model includes the $(2,2)$, $(3,2)$, and $(4,4)$ spin-weighted spherical harmonic modes. We also build a final black hole model, \texttt{NRSur2dq1EccRemnant}, which models the mass, and spin of the remnant black hole. We show that our waveform model can accurately predict numerical relativity waveforms with mismatches $\approx 10^{-3}$, while the remnant model can recover the final mass and dimensionless spin with absolute errors smaller than $\approx 5 \times 10^{-4}M$ and $\approx 2 \times10^{-3}$ respectively. We demonstrate that the waveform model can also recover subtle effects like mode-mixing in the ringdown signal without any special ad-hoc modeling steps. Finally, we show that despite being trained only on equal-mass binaries, \texttt{NRSur2dq1Ecc} can be reasonably extended up to mass ratio $q\approx3$ with mismatches $\simeq 10^{-2}$ for eccentricities smaller than $\sim 0.05$ as measured at a reference time of $2000M$ before merger. The methods developed here should prove useful in the building of future eccentric surrogate models over larger regions of the parameter space.