# Fixing the BMS Frame of Numerical Relativity Waveforms

Keefe Mitman  0000-0003-0276-3856<EMAIL_ADDRESS>Theoretical Astrophysics
350-17, California Institute of Technology, Pasadena, California 91125, USA
Neev Khera  0000-0003-3515-2859 Institute for Gravitation and the Cosmos &
Physics Department, Penn State, University Park, Pennsylvania 16802, USA
Dante A. B. Iozzo  0000-0002-7244-1900 Cornell Center for Astrophysics and
Planetary Science, Cornell University, Ithaca, New York 14853, USA  Leo C.
Stein  0000-0001-7559-9597 Department of Physics and Astronomy, University of
Mississippi, University, Mississippi 38677, USA  Michael Boyle
0000-0002-5075-5116 Cornell Center for Astrophysics and Planetary Science,
Cornell University, Ithaca, New York 14853, USA  Nils Deppe
0000-0003-4557-4115 Theoretical Astrophysics 350-17, California Institute of
Technology, Pasadena, California 91125, USA  Lawrence E. Kidder
0000-0001-5392-7342 Cornell Center for Astrophysics and Planetary Science,
Cornell University, Ithaca, New York 14853, USA  Jordan Moxon
0000-0001-9891-8677 Theoretical Astrophysics 350-17, California Institute of
Technology, Pasadena, California 91125, USA
Harald P. Pfeiffer  0000-0001-9288-519X Max Planck Institute for Gravitational
Physics (Albert Einstein Institute), Am Mühlenberg 1, Potsdam 14476, Germany
Mark A. Scheel  0000-0001-6656-9134 Theoretical Astrophysics 350-17,
California Institute of Technology, Pasadena, California 91125, USA  Saul A.
Teukolsky  0000-0001-9765-4526 Theoretical Astrophysics 350-17, California
Institute of Technology, Pasadena, California 91125, USA  Cornell Center for
Astrophysics and Planetary Science, Cornell University, Ithaca, New York
14853, USA  William Throwe  0000-0001-5059-4378 Cornell Center for
Astrophysics and Planetary Science, Cornell University, Ithaca, New York
14853, USA

###### Abstract

Understanding the Bondi-Metzner-Sachs (BMS) frame of the gravitational waves
produced by numerical relativity is crucial for ensuring that analyses on such
waveforms are performed properly. It is also important that models are built
from waveforms in the same BMS frame. Up until now, however, the BMS frame of
numerical waveforms has not been thoroughly examined, largely because the
necessary tools have not existed. In this paper, we show how to analyze and
map to a suitable BMS frame for numerical waveforms calculated with the
Spectral Einstein Code (SpEC). However, the methods and tools that we present
are general and can be applied to any numerical waveforms. We present an
extensive study of 13 binary black hole systems that broadly span parameter
space. From these simulations, we extract the strain and also the Weyl scalars
using both SpECTRE’s Cauchy-characteristic extraction module and also the
standard extrapolation procedure with a displacement memory correction applied
during postprocessing. First, we show that the current center-of-mass
correction used to map these waveforms to the center-of-mass frame is not as
effective as previously thought. Consequently, we also develop an improved
correction that utilizes asymptotic Poincaré charges instead of a Newtonian
center-of-mass trajectory. Next, we map our waveforms to the post-Newtonian
(PN) BMS frame using a PN strain waveform. This helps us find the unique BMS
transformation that minimizes the $L^{2}$ norm of the difference between the
numerical and PN strain waveforms during the early inspiral phase. We find
that once the waveforms are mapped to the PN BMS frame, they can be hybridized
with a PN strain waveform much more effectively than if one used any of the
previous alignment schemes, which only utilize the Poincaré transformations.

††preprint: APS/123-QED

###### Contents

  1. I Introduction
    1. I.1 Overview
    2. I.2 Conventions
  2. II Fixing the $\ell<2$ Transformations
  3. III Fixing the $\ell\geq 2$ Transformations
    1. III.1 Mapping to the super rest frame
    2. III.2 Mapping to the PN BMS frame
  4. IV Results
    1. IV.1 $\ell<2$ Results
      1. IV.1.1 A note on transient effects
      2. IV.1.2 Determining the best method for fixing the Poincaré frame
      3. IV.1.3 Comparing the Newtonian trajectory and charge-based frame fixing methods
      4. IV.1.4 Examining improvements to waveforms
    2. IV.2 $\ell\geq 2$ Results
  5. V Conclusion
  6. A Fixing the supertranslation freedom with the Moreschi supermomentum

##  I Introduction

As more and more astrophysical events are observed by gravitational wave
detectors such as LIGO,111The Laser Interferometer Gravitational-Wave
Observatory. Virgo, and KAGRA,222The Kamioka Gravitational Wave Detector.
accurate models of gravitational waves for these systems are indispensable for
conducting precise parameter estimation as well as tests of general
relativity. As of now, the most accurate models of gravitational waves are the
waveforms from numerical relativity (NR) simulations, with the most prevalent
being those that correspond to binary black hole (BBH) mergers [1, 2, 3, 4].

Like any observable system in nature, however, the waveforms produced by
numerical relativity simulations are functions of the frame that they are in.
Currently, the Simulating eXtreme Spacetimes (SXS) Catalog [1]—the largest of
the publicly available waveform catalogs—and the RIT Catalog [5] attempt to
fix the Poincaré frame of their asymptotic waveforms by mapping them to the
unique center-of-mass (CoM) frame through what they call a _center-of-mass
correction_ [1, 2, 6, 5].333This fixing of the Poincaré frame was called a
“correction” in the work of Ref. [6] because the BBH system’s center-of-mass
drift is an unexpected phenomenon in numerical relativity simulations. It
appears to be related to an imperfect boundary condition on the gauge degrees
of freedom at the outer boundary.

This correction—or this fixing of the Poincaré frame—uses the masses and the
trajectories of the two black holes to construct the Newtonian center-of-mass
trajectory

$\displaystyle\vec{x}_{\text{CoM}}\equiv\frac{m_{a}}{M}\vec{x}_{a}+\frac{m_{b}}{M}\vec{x}_{b},$
(1)

where $m_{a}$ and $m_{b}$ are the Christodoulou masses [7] of the primary and
secondary black holes, $M\equiv m_{a}+m_{b}$ is the total initial mass of the
BBH system, and $\vec{x}_{a}$ and $\vec{x}_{b}$ describe the motion of the
black holes’ centers, i.e., the simulation coordinate averages over each one
of the black hole’s apparent horizon. With this trajectory, the translation
and boost that best map the system to the Newtonian center-of-mass frame can
then be found by finding the transformation that minimizes the average square
of the distance between the center-of-mass and the origin of the corrected
frame [6].

While this method for fixing the Poincaré frame has proven useful, there are
certain aspects of this correction that are concerning. For one, when more
post-Newtonian (PN) terms are included in the Newtonian calculation of the
center of mass in Eq. (1), the correction does not improve and, for higher
mass-ratio or precessing systems, becomes noticeably worse [6]. Apart from
this, it is also worrisome that the frame of these waveforms is being fixed
based on information from the bulk of spacetime, rather than information from
the waveforms themselves, especially as it has been shown that the coordinate
velocity of the black holes does not accurately correspond to the asymptotic
velocity for complicated systems [8].

A better way to fix the Poincaré frame is to measure the Poincaré charges for
asymptotic quantities and try to find the Poincaré transformation that changes
these charges in a prescribed way. For example, to map a general relativistic
system to the center-of-mass frame, one could compute the linear momentum
charge and the boost charge from which the system’s boost and translation away
from the origin can be established. Or, as we will do in this work, one could
instead compute just the center-of-mass charge and then determine which
Poincaré transformation minimizes this charge.

Apart from the Poincaré frame, however, there is an additional freedom in
general relativity arising from the extra symmetries of asymptotically flat
spacetimes that extends the usual Poincaré group. This infinite group of
symmetries, which was found by Bondi, van der Burg, Metzner, and Sachs [9,
10], is known as the _BMS group_.444There is also a proposed generalization of
the BMS group, which promotes the Lorentz transformations to be the infinite
group of local diffeomorphisms on $S^{2}$ [11, 12]. In this work, however, we
will focus on just the BMS group and reserve an examination of the generalized
BMS group for future study. Fundamentally, the BMS group is a semidirect
product of the usual Lorentz group with an infinite-dimensional Abelian group
of transformations called _supertranslations_ , which contain the more
familiar spacetime translations as a normal subgroup.

The most straightforward way to understand how a supertranslation affects
coordinates is via the following. First define the Bondi time $u\equiv t-r$.
Under an arbitrary spacetime translation $(\delta t,\delta\vec{x})=(\delta
t,\delta x,\delta y,\delta z)$, we can write the corresponding transformation
of $u$ as

$\displaystyle
u^{\prime}=u-\sum\limits_{\ell<2}\sum\limits_{m\leq|\ell|}\alpha_{\ell
m}Y_{\ell m}(\theta,\phi),$ (2)

where

$\displaystyle\alpha_{0,0}$ $\displaystyle=\sqrt{4\pi}\delta t,$ (3a)
$\displaystyle\alpha_{1,\pm 1}$
$\displaystyle=-\sqrt{\frac{2\pi}{3}}(\mp\delta x+i\delta y),$ (3b)
$\displaystyle\quad\alpha_{1,0}$ $\displaystyle=-\sqrt{\frac{4\pi}{3}}\delta
z.$ (3c)

A proper supertranslation $\alpha(\theta,\phi)$, i.e., a supertranslation that
is not one of the spacetime translations (see Eq. (2)), then acts on $u$ as

$\displaystyle u^{\prime}=u-\alpha(\theta,\phi)$ (4)

for

$\displaystyle\alpha=\sum\limits_{\ell\geq
2}\sum\limits_{m\leq|\ell|}\alpha_{\ell m}Y_{\ell m}(\theta,\phi)$ (5)

with $\alpha_{\ell,m}=(-1)^{m}\overline{\alpha}_{\ell,-m}$ to make sure that
$u^{\prime}$ is real. Consequently, a supertranslation can be understood as a
direction-dependent time translation on the boundary of asymptotically flat
spacetimes, e.g., future null infinity. For example, if there exists a network
of observers on a sphere surrounding a source, then ideally they could combine
their received signals with some understanding of their clocks’
synchronization. At future null infinity, such a synchronization becomes
impossible and we could supply a separate time offset, i.e., a
supertranslation, to each observer without changing the observable physics. An
outline of how supertranslations transform the typical gravitational wave
quantities, such as the strain $h$, the news $\dot{h}$, and the Weyl scalars
$\Psi_{i}$, can be found in [13, 14].

There are really only two reasonable ways to fix this supertranslation
freedom. The first, and simplest, is to find the supertranslation that
minimizes the difference between a NR waveform and a PN waveform. The second,
which is often the most common in the literature [15, 16, 17], can be
understood as follows.

Like fixing the Poincaré frame by making use of the Poincaré charges, a fairly
similar scheme can be executed to fix a BBH system’s proper supertranslation
freedom by utilizing the proper supertranslation charge. We refer to the frame
that fixes this supertranslation freedom as the Bondi frame, while the frame
that corresponds to fixing all of the freedom of our waveforms is called the
BMS frame. Put differently, the BMS frame captures both the Poincaré and the
proper supertranslation freedom, whereas the Bondi frame only involves the
proper supertranslations. Just as the Bondi four-momentum is the charge that
is related to spacetime translations, the _supermomentum_ —an infinite
extension of the usual Bondi four-momentum—is the charge that corresponds to
supertranslations. Thus, the supertranslation freedom can also be uniquely
fixed by finding the supertranslation that minimizes the proper supermomentum
charge. The Bondi frame under this transformation is then related to the
Poincaré frame that corresponds to the transformation that minimizes the
three-momentum, i.e., the rest frame. In this case, though, we call the Bondi
frame with minimal555The reason we say minimal rather than no supermomentum is
because of a subtlety having to do with gravitational memory, which we discuss
in the Appendix A. supermomentum the _nice section_ [17] or the _super rest
frame_.

In this work, we do exactly this. That is, we fix the Poincaré frame of our
asymptotic waveforms by working with particular Poincaré charges to map the
BBH systems to the center-of-mass frame. We then also fix the proper
supertranslation freedom by finding the supertranslation that maps the
waveforms to the Bondi frame that is ideal for current observations. Thus, we
not only improve upon the work of [6] by using relativistic charges instead of
Newtonian trajectories, but we also wholly fix the proper supertranslation
freedom, thereby fixing the complete BMS frame of our numerical waveforms.

Most importantly, we find that the new and improved fixing of the Poincaré
frame is a drastic improvement over the previous method that uses Newtonian
trajectories. Based on our observations, we conclusively find that the
previous center-of-mass correction appears to have only approximately worked
for the equal mass nonspinning, aligned spins, and superkick systems, and
fails for the nonequal mass, anti-aligned spins, or precessing systems.
Further, even in the equal mass nonspinning systems where the previous
correction was nearly the same as the charge-based one, this new method
nonetheless shows obvious improvements, such as reducing the leakage of the
strain $(2,2)$ mode into other, subdominant strain modes. We show an example
of this improvement in Fig. 5.

Apart from improving the fixing of the Poincaré frame, we also make a few
important observations regarding the fixing of the supertranslation freedom.
Even though it is often mentioned that the supertranslation freedom can be
fixed by mapping to the super rest frame [15, 16, 17], we find that the most
practical way to fix the supertranslations is by mapping our BBH systems to
the PN Bondi frame. This is because LIGO expects numerical waveforms to agree
with PN waveforms during the early inspiral phase of a BBH merger. Note,
however, that for conducting any kind of analyses on quasinormal modes, the
preferred BMS frame is actually the super rest frame, since this is the frame
expected by the Teukolsky formalism.666This will be covered in a future work.

Previously, the BMS freedom of numerical waveforms has not been important
because the SXS Collaboration’s waveforms did not exhibit the displacement
memory [18]. However, because waveforms with memory effects can now be easily
produced by numerical relativity [18] or can even have memory effects added to
them via a correction [19], fixing this Bondi frame has become crucial, seeing
as it is absolutely necessary for performing hybridizations between numerical
and PN strain waveforms. We find that by mapping numerical waveforms to the PN
BMS frame, we can significantly improve NR/PN strain hybridizations. Our main
result regarding this is shown in Fig. 7.

Last, we also discover that by completely fixing the BMS frame of our
waveforms, we can perform noticeably better convergence tests of numerical
relativity waveforms, since waveforms from different resolutions can now be
compared while in the same BMS frame.

###  I.1 Overview

We organize our computations and results as follows.
In Sec. II we introduce the four main Poincaré charges that are useful when
examining asymptotic waveforms: the linear momentum, angular momentum, boost,
and center-of-mass charges. Apart from this, we also discuss how the center-
of-mass charge can be used to obtain the system’s velocity and translation
away from the origin. Following this, in Sec. III we discuss the two most
natural Bondi frames and conclude that for practical purposes, the most useful
frame to map to is the PN Bondi frame. Finally, in Secs. IV.1 and IV.2 we
present our results for mapping to the center-of-mass frame as well as the PN
BMS frame and thus illustrate how the BMS frame of numerical waveforms should
be fixed for future analyses and, most importantly, future surrogate models
[20, 21, 22, 23].

###  I.2 Conventions

We set $c=G=1$ and take $\eta_{\mu\nu}$ to be the $(-,+,+,+)$ Minkowski
metric. When working with complex dyads, following the work of [24], we use

$\displaystyle q_{A}=-\frac{1}{\sqrt{2}}(1,i\sin\theta)\text{ and
}q^{A}=-\frac{1}{\sqrt{2}}(1,i\csc\theta),$ (6)

and write the round metric on the two-sphere $S^{2}$ as $q_{AB}$. The complex
dyad obeys the following properties

$\displaystyle q_{A}q^{A}=0,\quad q_{A}\bar{q}^{A}=1,\quad
q_{AB}=q_{A}\bar{q}_{B}+\bar{q}_{A}q_{B}.$ (7)

Note that this convention differs from the related works of [19, 18, 25],
which in contrast do not include the $1/\sqrt{2}$ normalization factor on the
dyads in Eq. (6). We choose this convention because it makes our expressions
for the asymptotic charges in Eq. (15) more uniform. Nonetheless, for
transparency we provide the conversion between our quantities and those of
these previous works in Eq. (13). We build spin-weighted fields with the dyads
as follows. For a tensor field $W_{A\cdots D}$, the function

$\displaystyle W=W_{A\cdots BC\cdots D}q^{A}\cdots
q^{B}\bar{q}^{C}\cdots\bar{q}^{D}$ (8)

with $m$ factors of $q$ and $n$ factors of $\bar{q}$ has a spin weight of
$s=m-n$. When raising and lowering spin weights we use the Geroch-Held-Penrose
differential spin-weight operators $\eth$ and $\bar{\eth}$ [26],

$\displaystyle\eth W$ $\displaystyle=(D_{E}W_{A\cdots BC\cdots D})q^{A}\cdots
q^{B}\bar{q}^{C}\cdots\bar{q}^{D}q^{E},$ (9a) $\displaystyle\bar{\eth}W$
$\displaystyle=(D_{E}W_{A\cdots BC\cdots D})q^{A}\cdots
q^{B}\bar{q}^{C}\dots\bar{q}^{D}\bar{q}^{E}.$ (9b)

Here, $D_{A}$ is the covariant derivative on the two-sphere. The $\eth$ and
$\bar{\eth}$ operators in spherical coordinates are then

$\displaystyle\eth W(\theta,\phi)$
$\displaystyle=-\frac{1}{\sqrt{2}}(\sin\theta)^{+s}(\partial_{\theta}+i\csc\theta\partial_{\phi})$
$\displaystyle\phantom{=.-\frac{1}{\sqrt{2}}(\sin(\theta))^{s}}\left[(\sin\theta)^{-s}W(\theta,\phi)\right],$
(10a) $\displaystyle\bar{\eth}W(\theta,\phi)$
$\displaystyle=-\frac{1}{\sqrt{2}}(\sin\theta)^{-s}(\partial_{\theta}-i\csc\theta\partial_{\phi})$
$\displaystyle\phantom{=.-\frac{1}{\sqrt{2}}(\sin(\theta))^{s}}\left[(\sin\theta)^{+s}W(\theta,\phi)\right].$
(10b)

Thus, when acting on spin-weighted spherical harmonics, these operators
produce

$\displaystyle\eth(\phantom{}_{s}Y_{\ell m})$
$\displaystyle=+\frac{1}{\sqrt{2}}\sqrt{(\ell-s)(\ell+s+1)}_{s+1}Y_{\ell m},$
(11a) $\displaystyle\bar{\eth}(\phantom{}_{s}Y_{\ell m})$
$\displaystyle=-\frac{1}{\sqrt{2}}\sqrt{(\ell+s)(\ell-s+1)}_{s-1}Y_{\ell m}.$
(11b)

We denote the gravitational wave strain777We explicitly define the strain as
described in Appendix C of [2]. by $h$, which we represent in a spin-weight
$-2$ spherical harmonic basis,

$\displaystyle h(u,\theta,\phi)=\sum\limits_{\ell,m}h_{\ell
m}(u)\,{}_{-2}Y_{\ell m}(\theta,\phi),$ (12)

where, again, $u\equiv t-r$ is the Bondi time. We denote the Weyl scalars by
$\Psi_{0-4}$. The conversion from the convention of [25, 18, 19] (denoted
NR888NR because this is the convention that corresponds to the outputs of the
SXS simulations.) to ours (denoted MB999MB because this corresponds to the
Moreschi-Boyle convention used in the works [24, 13, 27] and the code scri
[28, 29, 30, 13].) is

$\displaystyle
h^{\texttt{NR}}=2\bar{\sigma}^{\texttt{MB}}\quad\text{and}\quad\Psi_{i}^{\texttt{NR}}=\frac{1}{2}(-\sqrt{2})^{i}\Psi_{i}^{\texttt{MB}}.$
(13)

Note that we will omit these superscripts and henceforth assume that
everything is in the MB convention.

##  II Fixing the $\ell<2$ Transformations

As discussed in the Introduction, for the past few years the method for fixing
the Poincaré frame of BBH systems in numerical relativity has relied on using
the Newtonian center-of-mass trajectory, i.e., Eq. (1) [6, 2, 5]. While this
has served as a successful first step, we can improve upon this by using
certain Poincaré charges: specifically, the center-of-mass charge.

We first present the main asymptotic Poincaré charges. These charges are the
linear momentum charge $P_{\Psi}$, the angular momentum charge $J_{\Psi}$, the
boost charge $K_{\Psi}$, and the energy moment charge $E_{\Psi}$. Others have
just called $E_{\Psi}$ the center-of-mass charge, but this is misleading
because $E_{\Psi}$ really measures the center of mass scaled by the energy of
the system. Thus, we instead refer to $E_{\Psi}$ as the energy moment charge
and reserve $G_{\Psi}$ to represent the center-of-mass charge. These charges
are computed by integrating the Bondi mass aspect $m$, the Lorentz aspect $N$,
and the energy moment aspect $E$, which are derived in the NR convention in
Appendixes A and B of [18]. In the MB convention these are

$\displaystyle m$
$\displaystyle\equiv-\text{Re}\left[\Psi_{2}+\sigma\dot{\bar{\sigma}}\right],$
(14a) $\displaystyle N$
$\displaystyle\equiv-\left(\Psi_{1}+\sigma\eth\bar{\sigma}+u\eth
m+\frac{1}{2}\eth\left(\sigma\bar{\sigma}\right)\right),$ (14b) $\displaystyle
E$ $\displaystyle\equiv N+u\eth
m=-\left(\Psi_{1}+\sigma\eth\bar{\sigma}+\frac{1}{2}\eth\left(\sigma\bar{\sigma}\right)\right).$
(14c)

Consequently, the main Poincaré charges are

$\displaystyle P_{\Psi}$
$\displaystyle=\frac{1}{4\pi}\int_{S^{2}}\Psi\,m\,d\Omega,$ (15a)
$\displaystyle J_{\Psi}$
$\displaystyle\equiv\frac{1}{4\pi}\int_{S^{2}}\epsilon^{AB}\left(\partial_{A}\Psi\right)\left(\bar{q}_{B}N\right)\,d\Omega,$
$\displaystyle=\frac{1}{4\pi}\int_{S^{2}}\text{Re}\left[\left(\bar{\eth}\Psi\right)\left(-iN\right)\right]\,d\Omega,$
(15b) $\displaystyle K_{\Psi}$
$\displaystyle\equiv\frac{1}{4\pi}\int_{S^{2}}q^{AB}\left(\partial_{A}\Psi\right)\left(\bar{q}_{B}N\right)\,d\Omega,$
$\displaystyle=\frac{1}{4\pi}\int_{S^{2}}\text{Re}\left[\left(\bar{\eth}\Psi\right)N\right]\,d\Omega,$
(15c) $\displaystyle E_{\Psi}$
$\displaystyle\equiv\frac{1}{4\pi}\int_{S^{2}}q^{AB}\left(\partial_{A}\Psi\right)\left(\bar{q}_{B}\left(N+u\eth
m\right)\right)\,d\Omega,$
$\displaystyle=\frac{1}{4\pi}\int_{S^{2}}\text{Re}\left[\left(\bar{\eth}\Psi\right)\left(N+u\eth
m\right)\right]\,d\Omega,$ $\displaystyle=K_{\Psi}+uP_{\Psi}.$ (15d)

where $\epsilon^{AB}$ is the usual Levi-Civita tensor and $\Psi$ is a real
spin-weight $0$ function on the two-sphere. This scalar $\Psi$ is typically
taken to be a unique combination of the $\ell\leq 1$ spherical harmonics so as
to represent one of the four Cartesian coordinates $t,x,y,z$, i.e.,

$\displaystyle t$ $\displaystyle=1$ $\displaystyle=\sqrt{4\pi}Y_{(0,0)},$
(16a) $\displaystyle x$ $\displaystyle=\sin\theta\cos\phi$
$\displaystyle=\sqrt{\frac{4\pi}{3}}\left[\frac{1}{\sqrt{2}}\left(Y_{(1,-1)}-Y_{(1,+1)}\right)\right],$
(16b) $\displaystyle y$ $\displaystyle=\sin\theta\sin\phi$
$\displaystyle=\sqrt{\frac{4\pi}{3}}\left[\frac{i}{\sqrt{2}}\left(Y_{(1,-1)}+Y_{(1,+1)}\right)\right],$
(16c) $\displaystyle z$ $\displaystyle=\cos\theta$
$\displaystyle=\sqrt{\frac{4\pi}{3}}Y_{(1,0)}.$ (16d)

By utilizing some properties of the spherical harmonics, we can create a four-
vector from the projection of a charge along each Cartesian direction:

$\displaystyle A^{t}$ $\displaystyle=\frac{1}{\sqrt{4\pi}}A_{(0,0)},$ (17a)
$\displaystyle A^{x}$
$\displaystyle=\frac{1}{\sqrt{4\pi}}\frac{1}{\sqrt{6}}\text{Re}\left[A_{(1,-1)}-A_{(1,+1)}\right],$
(17b) $\displaystyle A^{y}$
$\displaystyle=\frac{1}{\sqrt{4\pi}}\frac{1}{\sqrt{6}}\text{Im}\left[A_{(1,-1)}+A_{(1,+1)}\right],$
(17c) $\displaystyle A^{z}$
$\displaystyle=\frac{1}{\sqrt{4\pi}}\frac{1}{\sqrt{3}}\text{Re}\left[A_{(1,0)}\right],$
(17d)

where $A_{(\ell,m)}$ is the $(\ell,m)$ mode of the aspect $A$ when $A$ is
written in a spin-weight $0$ or $1$ spherical harmonic basis.101010Note that
$A=m$ is a spin-weight $0$ function and $A=-iN$, $N$, and $N+u\eth m$ are
spin-weight $1$ functions. Note that all of these charges and aspects, which
we have defined in Eq. (15), have been previously examined in earlier works
such as [31, 32, 33, 34, 18].

As mentioned in the introduction of this section, we are primarily interested
in the center-of-mass charge $G_{\Psi}$, which is closely related to the
energy moment charge $E_{\Psi}$ in Eq. (15). By definition, the center-of-mass
charge is the energy moment divided by the energy of the system,

$\displaystyle G_{\Psi}$
$\displaystyle\equiv\frac{E_{\Psi}}{P^{t}}=\frac{E_{\Psi}}{\gamma M_{B}},$
(18)

where $P^{t}$ is computed according to Eq. (17a) with $A=m$, $\gamma$ is the
Lorentz factor

$\displaystyle\gamma=\sqrt{1+\left|\vec{P}/P^{t}\right|^{2}},$ (19)

and $M_{B}$ is the Bondi mass [9]

$\displaystyle M_{B}\equiv\sqrt{-\eta_{\mu\nu}P^{\mu}P^{\nu}}.$ (20)

The charge $G_{\Psi}$ measures the center of mass that evolves linearly as a
function of momentum in the absence of gravitational radiation [34]. The
reason why we are mainly interested in this Poincaré charge is because

$\displaystyle\vec{G}|_{u=0}$ $\displaystyle=\left(\vec{K}/\gamma
M_{B}\right)|_{u=0},$ (21a) $\displaystyle\dot{\vec{G}}$
$\displaystyle=\frac{d}{du}\left[\left(\vec{K}+u\vec{P}\right)/\left(\gamma
M_{B}\right)\right]$
$\displaystyle=\vec{v}+u\dot{\vec{v}}-\frac{\vec{K}}{\gamma
M_{B}}\left[\frac{\dot{M}_{B}}{M_{B}}-\gamma^{2}\left(\vec{v}\cdot\dot{\vec{v}}\right)\right]+\frac{\dot{\vec{K}}}{\gamma
M_{B}}$ $\displaystyle\approx\vec{v}+\text{(oscillations about $\vec{v}$)},$
(21b)

since

$\displaystyle\dot{M}_{B}\approx 0,\quad\dot{\vec{v}}\approx
0,\quad\text{and}\quad\dot{\vec{K}}\approx 0,$ (22)

in the inspiral and ringdown phases of the BBH merger.111111We address issues
with this assumption in Sec. IV.1. Therefore, the intercept of
$\vec{G}\equiv(E^{x},E^{y},E^{z})/P^{t}$ will be the center of mass of the
system at $u=0$ and the slope will be the three-velocity. This means that by
fitting a degree one polynomial to this three-vector, we can obtain the amount
a system is translated and boosted out of the center-of-mass frame and then
apply the opposite Poincaré transformation to map our waveforms to the center-
of-mass frame.121212While we could also use the boost and linear momentum
charges to obtain these transformations, computing both of these charges is
more computationally expensive than if we just compute the center-of-mass
charge, since there are fewer products of waveforms that need to be taken when
finding $G_{\Psi}$. Note that these Poincaré charges in Eq. (15) can also be
used to measure properties of the BBH system or the remnant black hole, as was
illustrated in the recent and related work of [8].

##  III Fixing the $\ell\geq 2$ Transformations

We now move on to a discussion about the proper supertranslation freedom.

###  III.1 Mapping to the super rest frame

While it is often mentioned that a system’s Bondi frame should be fixed by
minimizing the supermomentum [16],
current observatories expect their waveform models to resemble PN expansions.
Therefore, even though fixing a system’s Bondi frame using the supermomentum
is a well-motivated option with unique benefits, we reserve a discussion of
this for Appendix A, since none of the results that we present involve this
technique.

Name | CCE radius |  $q$ |  $\chi_{A}$: | $(\hat{x},\,$ | $\hat{y},\,$ |  $\hat{z})$ |  $\chi_{B}$: | $(\hat{x},\,$ | $\hat{y},\,$ | $\hat{z})$
---|---|---|---|---|---|---|---|---|---|---
q1_nospin | 292 |  $1.0$ |  | $(0,\,$ | $0,\,$ |  $0)$ |  | $(0,\,$ | $0,\,$ | $0)$
q1_aligned_chi0_2 | 261 |  $1.0$ |  | $(0,\,$ | $0,\,$ |  $0.2)$ |  | $(0,\,$ | $0,\,$ | $0.2)$
q1_aligned_chi0_4 | 250 |  $1.0$ |  | $(0,\,$ | $0,\,$ |  $0.4)$ |  | $(0,\,$ | $0,\,$ | $0.4)$
q1_aligned_chi0_6 | 236 |  $1.0$ |  | $(0,\,$ | $0,\,$ |  $0.6)$ |  | $(0,\,$ | $0,\,$ | $0.6)$
q1_antialigned_chi0_2 | 274 |  $1.0$ |  | $(0,\,$ | $0,\,$ |  $0.2)$ |  | $(0,\,$ | $0,\,$ | $-0.2)$
q1_antialigned_chi0_4 | 273 |  $1.0$ |  | $(0,\,$ | $0,\,$ |  $0.4)$ |  | $(0,\,$ | $0,\,$ | $-0.4)$
q1_antialigned_chi0_6 | 270 |  $1.0$ |  | $(0,\,$ | $0,\,$ |  $0.6)$ |  | $(0,\,$ | $0,\,$ | $-0.6)$
q1_precessing | 305 |  $1.0$ |  | $(0.487,\,$ | $0.125,\,$ |  $-0.327)$ |  | $(-0.190,\,$ | $0.051,\,$ | $-0.227)$
q1_superkick | 270 |  $1.0$ |  | $(0.6,\,$ | $0,\,$ |  $0)$ |  | $(-0.6,\,$ | $0,\,$ | $0)$
q4_nospin | 235 |  $4.0$ |  | $(0,\,$ | $0,\,$ |  $0)$ |  | $(0,\,$ | $0,\,$ | $0)$
q4_aligned_chi0_4 | 222 |  $4.0$ |  | $(0,\,$ | $0,\,$ |  $0.4)$ |  | $(0,\,$ | $0,\,$ | $0.4)$
q4_antialigned_chi0_4 | 223 |  $4.0$ |  | $(0,\,$ | $0,\,$ |  $0.4)$ |  | $(0,\,$ | $0,\,$ | $-0.4)$
q4_precessing | 237 |  $4.0$ |  | $(0.487,\,$ | $0.125,\,$ |  $-0.327)$ |  | $(-0.190,\,$ | $0.051,\,$ | $-0.227)$
|  |  |  |  |  |  |  |  |  |
Table 1: Parameters of the BBH mergers used in our results. The mass ratio is
$q=M_{A}/M_{B}$, and the initial dimensionless spins of the two black holes
are $\chi_{A}$ and $\chi_{B}$. These simulations have been made publicly
available at [35].

###  III.2 Mapping to the PN BMS frame

Besides using the supermomentum to map a system to its super rest frame, the
supertranslation freedom of waveforms can also be fixed by mapping them to
their corresponding PN Bondi frame. Unlike PN waveforms, the NR waveforms are
finite in length and do not contain information from the BBH system’s entire
past history. As a result, when numerical strain waveforms that contain memory
are created, either by using Cauchy-characteristic extraction (CCE) [36, 18]
or by correcting extrapolated waveforms [19], their average during the
inspiral phase will tend to be close to zero. By contrast, post-Newtonian
waveforms typically have a memory contribution that is monotonically
increasing with time and only approaches zero as $u\rightarrow-\infty$.
Therefore, if hybridizations of the numerical waveforms are to be made with PN
waveforms, a mapping to the PN Bondi frame is essential to ensure that the
waveforms and their memory contributions can be properly aligned. While we
examined the results of mapping our systems to the super rest frame, we find
that, because of this hybridization concern, mapping to the PN Bondi frame is
the more sensible procedure for fixing the Bondi frame. Consequently, in Sec.
IV we only present the results for mapping our various BBH systems to their
corresponding PN BMS frame, as defined by a PN strain waveform, and reserve a
study of the benefits of mapping to the super rest frame for future work.

##  IV Results

For the following results, we numerically evolved a set of 13 binary black
hole mergers with various mass ratios and spin configurations using the
Spectral Einstein Code (SpEC) [37]. We list the important parameters of these
various BBH systems in Table 1. Each simulation contains roughly 19 orbits
prior to merger and is evolved until the waves from ringdown leave the
computational domain. Unlike the evolutions in the SXS catalog, the full set
of Weyl scalars and the strain have been extracted from these runs and the
waveforms have been computed using the extrapolation technique described in
[24] and the CCE procedure described in [25, 38]. Extrapolation is performed
with the python module scri [28, 29, 13, 30] and CCE is run with SpECTRE’s CCE
module [25, 38, 39].

For the CCE extractions, the four world tubes that are available have radii
that are equally spaced between $2\lambdabar_{0}$ and $21\lambdabar_{0}$,
where $\lambdabar_{0}=1/\omega_{0}$ is the initial reduced gravitational
wavelength as determined by the orbital frequency of the binary from the
initial data. Based on the recent work of [19], however, we choose to use only
the waveforms that correspond to the world tube with the second-smallest
radius, since these waveforms have been shown to minimally violate the BMS
balance laws. For clarity, we provide the world tube radius used for each
system in Table 1. All of these 13 BBH systems’ waveforms have been made
publicly available at [35].

As mentioned above, the asymptotic strain waveforms are computed using two
methods: extrapolation and CCE. The first method utilizes Regge-Wheeler-
Zerilli (RWZ) extraction to compute the strain on a series of concentric
spheres of constant coordinate radius and then proceeds to extrapolate these
values to future null infinity $\mathcal{I}^{+}$ using $1/r$ approximations
[40, 41, 42, 2, 24, 43]. This is the strain that can be found in the SXS
catalog. The other and more faithful extraction method, which is known as CCE,
computes the strain by using the world tube data provided by a Cauchy
evolution as the inner boundary data for a nonlinear evolution of the Einstein
field equations on null hypersurfaces [25, 38]. CCE requires freely specifying
the strain on the initial null hypersurface labeled $u=0$. Like [18, 19], we
choose this field to match the value and the first radial derivative of $h$
from the Cauchy data on the world tube using the ansatz

$\displaystyle
h(u=0,r,\theta^{A})=\frac{A(\theta^{A})}{r}+\frac{B(\theta^{A})}{r^{3}},$ (23)

where the two coefficients $A(\theta^{A})$ and $B(\theta^{A})$ are fixed by
the Cauchy data on the world tube. Unfortunately, constructing a satisfactory
initial null hypersurface for CCE is currently an open issue in numerical
relativity. Consequences of this choice manifest as transient effects arising
at early times [18]. We address these in Sec. IV.1.1.

As for the extrapolated strain waveforms that we use, these have been
postprocessed so that they exhibit the displacement memory effect and are thus
more on par with the waveforms produced by CCE [19].

Last, when performing our analysis, we predominantly use the code scri [28,
29, 30, 13] to compute BMS charges and transform asymptotic waveform
quantities.

###  IV.1 $\ell<2$ Results

Figure 1: Example of what the center-of-mass charge looks like for a system
with large center-of-mass velocity; in this case a $q=4$, precessing system.
The charges plotted are computed directly from Eq. (18), where $\vec{E}$ and
$P^{t}$ are obtained from Eq. (17) by using the energy moment aspect $E\equiv
N+u\eth m$ and mass aspect $m$. We define the peak time to be the peak of the
$L^{2}$ norm of the Geroch supermomentum [see Eq. (32b)], since this quantity
is a supertranslation-invariant quantity [44, 17]. The waveform used is a CCE
waveform.
BBH merger: q4_precessing (see Table 1).

As discussed in Sec. II, a BBH system can be mapped to its center-of-mass
frame by utilizing the center-of-mass charge $\vec{G}$. In Fig. 1, we show
this charge for a $q=4$, precessing BBH system (see Table 1). We compute the
plotted charges from Eq. (18) by using a CCE waveform. As can be seen, in the
$\hat{x}$, $\hat{y}$, and $\hat{z}$ directions the average of the center-of-
mass charge is not constant with respect to the Bondi time. Further, if one
imagines tracing these curves back in time then it can easily be observed that
they begin with a nonzero value. Because of these results, we may assert that,
with time, the BBH system is drifting through space away from a point that is
not the origin. If the system under consideration were not drifting and were
in the center-of-mass frame, then we would expect our charges to have both
zero slope and zero intercept. Fortunately, because of the nature of this
charge, to map to the center-of-mass frame one can simply boost and translate
the system by the negative of the charge’s slope and intercept. In the
subsequent discussion, we first check to see if any of our of 13 binary
systems, either before or after the Newtonian center-of-mass correction [6],
are in the center-of-mass frame. After this, we then proceed to apply our
charge-based center-of-mass correction and evaluate the improvement it has on
our waveforms.

####  IV.1.1 A note on transient effects

Because numerical relativity simulations are evolved from imperfect initial
data [2], the output waveforms contain unphysical effects referred to as _junk
radiation_. In extrapolated waveforms, the junk radiation appears at early
times and typically decays after an orbit or two. For CCE waveforms, however,
the junk radiation tends to persist longer into the waveform [18]. As a
result, because we examine extrapolated and CCE waveforms, we perform our
analysis on the part of our waveforms that is three orbits past the start of
the simulation, since the CCE-specific transient effects have decayed by then.
We refer to this time throughout the results as $u_{1}$.

####  IV.1.2 Determining the best method for fixing the Poincaré frame

Figure 2: The convergence of the boost and translation vectors obtained by
fitting to the center-of-mass charge. Here, $\beta_{i}$ ($\alpha_{i}$)
represents the boost (translation) obtained after $i$ iterations of fitting to
the charge $\vec{G}$, transforming the untransformed waveform with the fit
result, fitting to the new charge, and then composing that fit with the
previous vector. For example, $\beta_{1}$ is obtained from fitting to the
charge once, and $\beta_{2}$ is obtained from fitting to the charge
transformed using $\beta_{1}$ and then composing that fit with $\beta_{1}$.
Note that we show the result for q1_aligned_chi0_6 because this simulation
exhibits the slowest convergence out of the 13 systems that we examined. The
waveform used is a CCE waveform.
BBH merger: q1_aligned_chi0_6 (see Table 1).

Figure 3: Examining the boost and translation vectors as measured by
iteratively fitting to the center-of-mass charge $\vec{G}$, which is computed
by using Eq. (17) to obtain the vector components of Eq. (18). In the top row
we show the magnitude of these vectors, which we label $\vec{v}_{G}$ with
$\vec{v}=\vec{\beta}$ for boosts (left) and $\vec{v}=\vec{\alpha}$ for
translations (right). In the middle row we compute the relative error between
these vectors and the vectors $\vec{v}_{\text{CoM}}$, which are obtained by
the previous center-of-mass correction that relies on the Newtonian
trajectories of the black holes. In the bottom row we plot the angle between
these vectors using Eq. (24).

We now compute the boosts and translations needed to map the 13 BBH systems to
the center-of-mass frame. We first compute the charges according to Eq. (18).
Next, we define an initial time $u_{1}$ to be three orbits past the start of
the simulation and a final time that is three orbits before the peak time
$u_{\text{peak}}$. We define the peak time to be the peak of the $L^{2}$ norm
of the Geroch supermomentum [see Eq. (32b)], since this quantity is a
supertranslation-invariant quantity [44, 17]. We choose this final time to
ensure that we are only working with the inspiral phase of the binary, rather
than the merger phase. Note that one could fix the Poincaré frame using the
remnant BH, i.e., mapping the kick velocity to zero, but this is not as
instinctive as using the inspiral phase to fix the frame, even though it would
matter for fitting quasinormal modes. Equipped with the charges and boundary
times, we then linearly fit to the center-of-mass charge in the $\hat{x}$,
$\hat{y}$, and $\hat{z}$ directions and take the needed boost to be the
negative of the slope and the needed translation to be the negative of the
vertical intercept at $u=0$.

When we first applied this new Poincaré frame fixing, the improvements that we
saw in the extreme systems, e.g., the $q=4$ and fast-spinning systems, were
remarkable. For the simpler systems, however, the improvements were not as
large as we expected them to be. This is because these systems are already
reasonably close to being in the center-of-mass frame. Thus the measured
boosts and translations are more susceptible to errors that are introduced by
oscillations in the charge and also by the failure of the assumptions in Eq.
(22) to hold because of non-Newtonian effects. Fortunately, because this
method for mapping to the center-of-mass frame just involves the computation
of charges, it can be iterated. That is, after the initial transformation is
found, it can be applied, and then the center-of-mass charge can be computed
again from the new asymptotic waveform. With this new charge, a new
transformation can then be found, which may be composed with the previous
transformation to obtain a more accurate mapping to the center-of-mass frame.

In Fig. 2 we plot the convergence of the boosts and translations for various
iterations of this fitting procedure. We only show the results for
bbh_q1_aligned_chi0_6 because this system has the slowest convergence of the
13 binaries examined. In this figure, $\beta_{i}$ ($\alpha_{i}$) represents
the boost (translation) obtained after $i$ iterations of this charge fitting
process. More specifically, this procedure is as follows:

  1. I.

Take $\vec{\beta}_{0},\vec{\alpha}_{0}=0$.

  2. II.

If the iteration number $i$ is not $0$, transform the waveform with the
Poincaré transformation $\vec{\beta}_{i},\vec{\alpha}_{i}$.

  3. III.

Compute $\vec{G}$ from the transformed waveform.

  4. IV.

Obtain $\vec{\beta}$ and $\vec{\alpha}$ by fitting to $\vec{G}$ with a degree
one polynomial ($\vec{\beta}=\text{slope}$, $\vec{\alpha}=\text{vertical
intercept}$).

  5. V.

Compute $\vec{\beta}_{i+1}=\vec{\beta}_{i}+\vec{\beta}$ and
$\vec{\alpha}_{i+1}=\vec{\alpha}_{i}+\vec{\alpha}$.

  6. VI.

Repeat.

As can be seen, even for this slowly converging system, the transformations
converge rather quickly. Thus, based on this simulation, we choose a fixed
number of five iterations for every one of our systems. With this many
iterations, the improvements in the simpler systems then become as large as we
would like to see.

####  IV.1.3 Comparing the Newtonian trajectory and charge-based frame fixing
methods

Figure 4: Comparing the center-of-mass charge from waveforms transformed using
the Newtonian center-of-mass correction (red/dashed curves) to the charge
obtained from those transformed using the charge-based method (black/solid
curves). As a reference, we also plot the charge obtained from the
untransformed waveforms (blue/dotted curves). The charges plotted are computed
directly from Eq. (18), where $\vec{E}$ and $P^{t}$ are obtained from Eq. (17)
by using the energy moment aspect $E\equiv N+u\eth m$ and mass aspect $m$. The
waveform used is a CCE waveform.
BBH merger: q4_precessing (see Table 1).
Figure 5: The strain $(2,1)$ mode of the original waveform (blue/dotted curve,
scaled by $10^{-3}$) compared to that of the Newtonian (red/dashed curve,
scaled by $10^{-1}$) and the charge-based (black/solid curve) center-of-mass
corrections. In the inset plot we provide the strain $(2,2)$ mode of the
original waveform to illustrate that unless the charge-based correction is
used, the $(2,1)$ mode exhibits the same frequency as the $(2,2)$ mode. Note
that the time axis of the inset plot matches up with the main plot’s. The
waveform used is a CCE waveform.
BBH merger: q1_nospin (see Table 1). Figure 6: The modes that produce the
largest mismatches between the charge-based center-of-mass corrected strains
and both the original strains (blue) and the strains corrected using the
Newtonian center-of-mass correction (red). The mismatch in each mode is
computed by using Eq. (IV.1.4). As shown in the plot’s legend, the shape of
each point represents the mode while the color represents the strains used in
the mismatch. Note that the mismatches for the $q=1$ nonspinning and aligned-
spin systems are so large because these systems have a rotation by $\pi$
symmetry, and thus when the charge-based center-of-mass correction is used the
$m=\text{odd}$ modes become much closer to their expected value of zero. The
waveforms used in these computations are CCE waveforms.

We now compare our iterative center-of-mass correction that uses the center-
of-mass charge to the previous version that uses Newtonian trajectories. In
Fig. 3 we show that of the 13 binaries examined, none of them are exactly in
the center-of-mass frame after applying the correction that relies on
Newtonian trajectories since all of the new Poincaré transformations are
nonzero. More specifically, we plot two columns, one for the boosts on the
left and one for the translations on the right. In the top row, we plot the
magnitude of these vectors as obtained from fitting to the center-of-mass
charge, which is a proxy for how much the system fails to be in the center-of-
mass frame. In the middle row, we plot the relative difference between this
charge-based vector $\vec{v}_{\text{G}}$, and the vector obtained from the
Newtonian center-of-mass correction $\vec{v}_{\text{CoM}}$. This serves as a
proxy for how much the system fails to be in the center-of-mass frame, even
after the Newtonian correction. Finally, in the bottom row, we plot the angle
between these two vectors. That is,

$\displaystyle\Delta\theta(\vec{v}_{\text{G}},\vec{v}_{\text{CoM}})\equiv\arccos\left(\frac{\vec{v}_{\text{G}}}{|\vec{v}_{\text{G}}|}\cdot\frac{\vec{v}_{\text{CoM}}}{|\vec{v}_{\text{CoM}}|}\right).$
(24)

As can be seen in the top row, the equal mass nonspinning and aligned systems
are reasonably close to being in the center-of-mass frame, while the other
systems, especially the $q=4$ systems, are not. According to the other rows,
though, because the differences between the vectors are so large, we realize
that the Newtonian method for mapping to the center-of-mass frame is not
nearly as successful as previously thought. Consequently, it is now evident
that this Newtonian method for mapping a BBH system to its center-of-mass
frame does not achieve its objective and the method based on the center-of-
mass charge, which we explore in more detail now, is necessary.

####  IV.1.4 Examining improvements to waveforms

At this point, we examine how the strain waveforms change under the center-of-
mass charge-based mapping to the center-of-mass frame. First, though, in Fig.
4 we show how the center-of-mass charge changes under the Newtonian center-of-
mass correction versus the new charge-based method. What this plot shows is
that while the Newtonian center-of-mass correction only corrects the $\hat{x}$
component of the center-of-mass charges for the q4_precessing system, the new
method produces an average value of exactly zero for every vector component.

In Fig. 5 we show the most important consequence of improving the fixing of
the Poincaré frame by using asymptotic Poincaré charges. For this example, we
show the strain $(2,1)$ mode of q1_nospin as it is, after the Newtonian
center-of-mass transformation, and after our new center-of-mass
transformation. We show this mode because it exhibits the largest mismatch
(see Eq. (IV.1.4)) when comparing strains that have been transformed using
both the previous and the new center-of-mass corrections. Based on PN theory,
during the system’s inspiral phase we expect the frequency of this mode to be
half the strain $(2,2)$ mode’s frequency.131313Really we expect this mode to
be exactly zero because of the symmetry of the system; however, because the
two spins and the eccentricity of the black holes in the bbh_q1_nospin
simulation are not precisely zero, owing to numerical error, there is an
unexpected nonzero contribution to this mode. However, as can be seen by
comparing the original and transformed waveforms, after correctly mapping to
the center-of-mass frame the frequency of the strain $(2,1)$ mode is roughly
half of what it was before. This is because, previously, the system was not
truly in the center-of-mass frame so the strain $(2,2)$ mode was leaking into
the $(2,1)$ mode. Note that in Fig. 5 we have scaled the original waveform and
the waveform transformed with the Newtonian center-of-mass correction by
factors of $10^{-3}$ and $10^{-1}$ to make them more comparable to the
waveform transformed with the charge-based center-of-mass correction.

Last, to show the impact this charge-based method has more broadly, we provide
Fig. 6. This figure shows the mismatch between the newly transformed waveforms
and both the original waveforms and the waveforms that have been transformed
using the previous Newtonian correction. In this plot, we compute the mismatch
between a mode of two strain waveforms via

$\displaystyle\mathcal{M}($ $\displaystyle
h_{(\ell,m)}^{A},h_{(\ell,m)}^{B})\equiv$
$\displaystyle\phantom{=.}1-\frac{\langle
h_{(\ell,m)}^{A},h_{(\ell,m)}^{B}\rangle}{\sqrt{\langle
h_{(\ell,m)}^{A},h_{(\ell,m)}^{A}\rangle\langle
h_{(\ell,m)}^{B},h_{(\ell,m)}^{B}\rangle}},$ (25)

in which the inner product is given by

$\displaystyle\langle
h_{(\ell,m)}^{A},h_{(\ell,m)}^{B}\rangle\equiv\int_{u_{1}}^{+\infty}h_{(\ell,m)}^{A}\overline{h_{(\ell,m)}^{B}}du,$
(26)

where $u=+\infty$ is the final time of the simulation.

###  IV.2 $\ell\geq 2$ Results

Figure 7: Comparing the BMS alignment results between NR and PN strain
waveforms. In red, we show the initial misalignment. In blue, we show the
misalignment if the usual alignment procedure is used, i.e., finding the time
translation and frame rotation that best aligns a NR strain waveform that has
undergone the previous center-of-mass correction to a PN strain waveform.
Finally, in green, we show the misalignment after the new BMS frame alignment
procedure has been used, that is, finding the BMS transformation (up to
$\ell=4$) that minimizes the average $L^{2}$ norm of the difference of the NR
and PN strain waveforms. Note that the measure of this misalignment,
$||f(u,\theta,\phi)||^{2}$, is defined to be
$||f(u,\theta,\phi)||^{2}=\int_{u_{1}}^{u_{1}+\text{four
orbits}}\int_{S^{2}}f\bar{f}\,d\Omega\,du$.

With the $\ell<2$ modes of our waveforms properly fixed using the new center-
of-mass correction, we now explore how to fix the $\ell\geq 2$ modes, i.e.,
choosing a Bondi frame. As discussed in Sec. III, there are really only two
options: mapping the systems to their super rest frame or their PN Bondi
frame. As described earlier, we prefer to map to the PN BMS frame since this
tends to best improve the hybridization between two NR and PN strain
waveforms.

To map our various systems to their PN BMS frame, we begin by creating a 3PN-
order strain waveform from the orbital frequency of the two black holes using
the code GWFrames [28, 29]. To generate this PN waveform, we obtain the
orbital frequency of the system from the horizon information and evolve it
backward in time using the PN evolution equations. We then simultaneously find
the Poincaré transformation and the $2\leq\ell\leq 4$ supertranslation that
minimize the $L^{2}$ norm of the difference between the strain and the PN
waveform. The norm is computed over the time interval starting at $u_{1}$ and
continuing for four orbits. To perform this minimization, we use SciPy’s
minimize function corresponding to the sequential least squares programming
algorithm (SLSQP) [45] and define the following function:

  1. I.

Take, as inputs, a NR strain, a PN strain, and also a center-of-mass
transformation.

  2. II.

Consider the ordered list of transformations

    1. 1.

$\ell=2$ supertranslation,

    2. 2.

time translation,

    3. 3.

frame rotation,

    4. 4.

$\ell=3$ supertranslation,

    5. 5.

$\ell=4$ supertranslation.

  3. III.

Begin with iteration $n=1$.

  4. IV.

For iteration $n$, include all transformations from the above list up to
transformation $n$ as free parameters. Use the findings from iteration $n-1$
as initial guesses for the $n-1$ transformations.

  5. V.

Use SciPy’s SLSQP minimize function to find the collection of $n$
transformations that, when composed with the center-of-mass transformation
(see below), best map the NR strain to the frame of the PN strain, i.e., the
transformations that minimize the $L^{2}$ norm of the difference of the two
waveforms, integrated over the time interval from $u_{1}$ to four orbits past
$u_{1}$.

  6. VI.

Repeat until $n=5$.

We find that it is important to start this procedure with the $\ell=2$
supertranslation, because this tends to be the largest source of error. Beyond
this, however, the order of the transformations is fairly inconsequential and
chosen as such to minimize run time. Note though that it is also important to
allow the previous transformations to be free parameters in the next iteration
because each new transformation tends to influence the previous ones.

With this function, we then run the iterative procedure:

  1. I.

Find the center-of-mass transformation that maps the NR strain waveform to the
center-of-mass frame using the charge-based center-of-mass correction.

  2. II.

Provide this transformation and the NR and PN strain waveforms to the
minimizing function.

  3. III.

Apply the optimized BMS transformation to the raw NR strain141414We apply the
BMS transformation to the raw waveform to ensure that we are not transforming
the same object more than once, which would introduce numerical error. and
find the center-of-mass transformation needed to map this new waveform to the
center-of-mass frame.

  4. IV.

Compose this new transformation with the original center-of-mass
transformation and then repeat steps II. - IV. with the previous BMS
transformation as an initial guess until a desired precision is obtained.

We find that by running this procedure 4 times we can obtain rather impressive
alignments between the input NR and PN strain waveforms for most of our
systems. We choose to run this method, rather than optimizing over all BMS
transformations, because this method not only produces the best alignment but
it also tends to keep the system much closer to its center-of-mass frame.

In Fig. 7, we show the results of mapping the various CCE and extrapolated
waveforms from our 13 systems to their corresponding PN BMS frame. In red, we
show the initial misalignment. In blue, we show the misalignment if the usual
alignment procedure is performed, i.e., finding the time translation and frame
rotation that best aligns a PN strain waveform to a NR strain waveform that
has undergone the previous center-of-mass correction. Finally, in green, we
show the misalignment after using the new BMS frame alignment procedure. As is
clearly illustrated, by capitalizing on the full BMS freedom of NR waveforms
one can perform substantially better alignment between NR and PN strain
waveforms. Apart from this, though, one can also observe the failure of the PN
waveform to accurately model the BBH system, e.g., as the total spin or the
mass ratio of the system increases, or if the system is precessing, the
success of the BMS alignment between the NR and PN strains tends to worsen.
This is expected, however, since the PN waveform is only of 3PN order.
Regardless, this shows that not only can the Bondi frame, and really the whole
BMS frame, be fixed by utilizing a PN strain waveform, but doing so is
critically important for aligning, and thus hybridizing, NR and PN waveforms.

Figure 8: The NR and PN strain hybrid. We present three plots, one for each of
the $(2,2)$, $(2,1)$, and $(2,0)$ modes. In each plot, in the top panel we
show the NR waveform (black/solid curve), the PN waveform (blue/dashed curve),
and also the hybrid waveform (orange/dotted curve), while in the bottom panel
we plot the relative error between the hybrid waveform and both the NR
waveform (black curve) and the PN waveform (blue curve). In each panel, we
also show the hybridization interval in purple. The NR waveform in each of
these plots is a CCE waveform.
BBH merger: q4_aligned_chi0_4 (see Table 1). Figure 9: Illustrating a caveat
about NR and PN hybridizations, as shown through the strain $(2,1)$ mode. In
the top panel we plot the NR waveform (black/solid curve), the PN waveform
(blue/dashed curve), and also the hybrid waveform (orange/dotted curve), while
in the bottom panel we plot the absolute error between the hybrid waveform and
both the NR waveform (black curve) and the PN waveform (blue curve). In each
panel, we also show the hybridization interval in purple. See the text for
more details. The NR waveform in this plot is a CCE waveform.
BBH merger: q1_nospin (see Table 1).

With our NR waveforms now optimally mapped to the PN BMS frame, we perform
strain hybridizations between NR and PN to illustrate the operations this
alignment procedure allows for. To create these hybridizations, we use the
smoothing function

$\displaystyle f(x)=\begin{cases}0&x\leq 0,\\\
\left(1+\exp\left[\frac{1}{x-1}+\frac{1}{x}\right]\right)^{-1}&0<x<1,\\\
1&x\geq 1.\end{cases}$ (27)

so that before the hybridization interval, which is the same as the alignment
interval, the hybrid is equal to the PN waveform and after it is equal to the
NR waveform. Put differently, we build the hybrid waveform $h^{\text{hybrid}}$
via

$\displaystyle
h^{\text{hybrid}}=h^{\text{PN}}+f\left(\frac{u-u_{1}}{u_{2}-u_{1}}\right)\left(h^{\text{NR}}-h^{\text{PN}}\right),$
(28)

where $u_{2}$ is the time that is four orbits past $u_{1}$.

In Fig. 8, we show three plots, one for each of the strain $(2,2)$, $(2,1)$,
and $(2,0)$ modes. In each plot, we compare the hybrid waveform to the NR and
PN strains for simulation q4_aligned_chi0_4. In each top panel, we show the
strain modes, while in each bottom panel, we show the relative error between
the hybrid and the NR and PN strain waveforms. As is expected from our
alignment results in Fig. 7, there is fairly impressive agreement in every
mode. Furthermore, the plot of the strain $(2,0)$ mode shows that the initial
value of the hybrid now agrees with PN, i.e., it exhibits the memory that we
expect to be there due to the emission of radiation throughout the entire past
history of the binary’s inspiral. Apart from these important positive results,
though, there is one minor caveat regarding this hybridization procedure that
is worth mentioning to avoid confusion.

Figure 10: The magnitudes of the boosts and translations needed to map the NR
strains in the PN BMS frame to their center-of-mass frame. The reason these
are not exactly zero is because a few of the BMS transformations that are not
involved in the center-of-mass correction do not commute with the center-of-
mass transformation, causing the system to be pushed slightly away from the
center-of-mass frame to obtain better alignment with the PN strain waveform.
Nonetheless, these transformations are fairly negligible, especially when
compared to those in the top plots of Fig. 3. Figure 11: Comparing the BMS
alignment results between a high resolution and a lower resolution numerical
relativity waveform. In red, we show the initial misalignment, while in green,
we show the misalignment after using the new BMS frame alignment: that is,
finding the BMS transformation (up to $\ell=4$) that minimizes the average
$L^{2}$ norm of the difference of the two waveforms. Note that the measure of
this misalignment, $||f(u,\theta,\phi)||^{2}$, is defined to be
$||f(u,\theta,\phi)||^{2}=\int_{u_{1}}^{u_{1}+\text{four
orbits}}\int_{S^{2}}f\bar{f}\,d\Omega du$.

Recall that when mapping a NR strain waveform to its corresponding PN BMS
frame, we apply certain $2\leq\ell\leq 4$ supertranslations. Because
supertranslations also affect the Bondi time $u$, however, when we
supertranslate the strain we not only have to act on it with the
supertranslation, but we also have to interpolate the waveform on to a new
series of Bondi times. So, when applying a supertranslation with only one
nonzero mode, not only will that mode of the strain change, but so will every
other mode because of the time interpolation. The reason why this is important
is illustrated in Fig. 9.

In Fig. 9, we show the $(2,1)$ mode of the strain hybrid for simulation
q1_nospin. Based purely on symmetry, we would normally expect this mode to be
zero, but because of numerical error, the spins and eccentricity of this
system are not exactly zero. We therefore observe nonzero values in this mode
of the NR waveform, even though they are negligible. Thus, the $(2,1)$ mode of
the hybrid looks a bit strange because of this difference between NR and PN.

Apart from this, however, one may also notice that the zero average value of
the PN strain is not matched by the nonzero average of the NR strain. This
oddity results from the supertranslation’s broader influence on the whole of
the strain waveform because of the needed time interpolation, as mentioned
earlier. In other words, supertranslating the $(2,1)$ mode away from an
average of zero improves the error in other modes more than it worsens the
error in the $(2,1)$ mode. Consequently, while this kind of behavior in the
hybridized waveform is certainly undesirable, it is a natural consequence of
the PN waveform not being a perfect theoretical model for the numerical BBH
system. Furthermore, we note that this behavior only occurs in modes whose
amplitude is rather negligible, i.e., $10^{-6}$ or less.

It is also fairly important to note that, because the center-of-mass
transformation does not commute with the other BMS transformations, this
procedure of mapping to the PN BMS frame does not exactly map our systems to
the center-of-mass frame, even though they are very close. We illustrate this
in Fig. 10, which shows the remaining boosts and translations needed to map
the waveform in the PN BMS frame to its center-of-mass frame. As is shown, the
necessary boosts and translations are nearly zero, especially when compared to
the transformations in the top plots of Fig. 3. Therefore, we consider this
mapping to the PN BMS frame to be remarkably successful, especially since the
improvements in Fig. 7 are so notable.

Finally, we show our last result regarding the benefit of utilizing the whole
BMS freedom of numerical waveforms. When examining physical quantities that
are output by numerical relativity, it is important to run convergence tests
to ensure that conclusions can be made with respect to numerical error. In
Fig. 11 we show how convergence tests can be improved by mapping waveforms
that are from simulations of different numerical resolutions to the same BMS
frame before they are compared with one another. To do this, we perform the
same iterative process as described earlier, but we now optimize over every
BMS transformation, rather than everything but the center-of-mass
transformation. In red, we show the initial misalignment, while in green we
show the misalignment after the BMS frame alignment procedure. As is shown,
the improvements are relatively minor, but could still prove to be important
for numerical simulations run with newer codes, such as SpECTRE [39], which
will be more accurate than the SXS Collaboration’s current code SpEC.

##  V Conclusion

Like any physical system, understanding the frame that a binary black hole
merger is in relative to a family of observers is essential. Understanding
this frame will help us ensure that any analyses on waveform models are
performed properly and no misleading assertions are made. For gravitational-
wave physics, fixing the frame is not as simple as fixing the Poincaré frame,
since the symmetries of asymptotically flat spacetimes are characterized by an
infinite extension of the Poincaré group: the BMS group.

Currently, gravitational-wave physicists who analyze models of gravitational
waves expect those models to be in the center-of-mass frame and the PN Bondi
frame, since this is the BMS frame that analytic models are in. However, the
waveforms that are currently produced by numerical relativity—the supplier of
the most accurate models of gravitational waves—typically are not in such a
frame because of an unexpected center-of-mass drift in numerical simulations
and a lack of initial data that contains information about the entire past
history of the binary black hole’s inspiral. Consequently, they are instead in
some other BMS frame. But with a proper understanding of the BMS group, one
can postprocess these waveforms and map them to the desired BMS frame after
the BBH simulation is complete. As of now, such a postprocessing technique is
used for the Poincaré frame by using the Newtonian center of mass [6], but
there is no such postprocessing for the Bondi frame.

In this paper, by utilizing asymptotic Poincaré charges, i.e., the center-of-
mass charge (see Eq. (18)), we show that this method that relies on Newtonian
trajectories for mapping to the center of mass Frame is not as successful as
previously thought. As a result, we develop an improved procedure for fixing
the Poincaré frame, which shows large benefits in terms of exhibiting the
expected behavior in the asymptotic Poincaré charges and also in minimizing
certain modes of the strain waveform that are expected to be zero because of
axisymmetry.

We also found that we can meaningfully fix the whole BMS frame of our
numerical waveforms by mapping them to their corresponding PN BMS frame using
a 3PN-order strain waveform. With this BMS frame fixing procedure, we observe
that we can produce much more favorable hybridizations between NR and PN
strain waveforms than if one were to use an alignment scheme that only
utilized the Poincaré transformations. Last, we also find that such a BMS
frame alignment will prove important for future numerical relativity codes
that will be able to run simulations at higher resolutions and will need to
properly test the convergence of their waveforms.

With this new method of fully fixing the BMS frame of asymptotic waveforms,
many important improvements to gravitational wave modeling can be made. For
example, by correctly mapping to the center-of-mass Poincaré frame and PN BMS
frame, we can produce much better PN and NR hybridizations. Furthermore,
because we can now ensure that waveforms are in the same BMS frame, surrogate
models built from such waveforms should be more accurate since they are no
longer trying to link waveforms in radically different BMS frames. It would be
very interesting to see how parameter estimation using a NR surrogate changes
depending on whether the waveforms used to build the surrogate have been
mapped to a consistent BMS frame, such as the PN BMS frame.

###### Acknowledgements.

Computations were performed with the High Performance Computing Center and the
Wheeler cluster at Caltech. This work was supported in part by the Sherman
Fairchild Foundation and by NSF Grants No. PHY-2011961, No. PHY-2011968, and
No. OAC-1931266 at Caltech, NSF Grants No. PHY-1912081 and No. OAC-1931280 at
Cornell, and NSF Grant No. PHY-1806356, Grant No. UN2017-92945 from the Urania
Stott Fund of the Pittsburgh Foundation, and the Eberly research funds of Penn
State at Penn State.

##  Appendix A Fixing the supertranslation freedom with the Moreschi
supermomentum

As mentioned earlier, an alternative method for fixing the Bondi frame of NR
waveforms, apart from mapping to the PN BMS frame, is to map to the super rest
frame. Therefore, since this technique could prove rather useful in future
analyses—such as quasinormal mode fitting— we now briefly illustrate how to
perform this procedure. According to Moreschi [15, 16], a reasonable choice
for the supermomentum is

$\displaystyle\Psi^{\text{M}}(u,\theta,\phi)\equiv\sum\limits_{\ell\geq
0}\sum\limits_{m\leq|\ell|}\Psi_{\ell m}^{\text{M}}(u)Y_{\ell
m}(\theta,\phi),$ (29)

where

$\displaystyle\Psi_{\ell
m}^{\text{M}}(u)\equiv-\frac{1}{\sqrt{4\pi}}\int_{S^{2}}Y_{\ell
m}\Psi^{\text{M}}(u)\,d\Omega,$ (30)

and

$\displaystyle\Psi^{\text{M}}(u)\equiv\Psi_{2}+\sigma\dot{\bar{\sigma}}+\eth^{2}\bar{\sigma}.$
(31)

While this is the supermomentum that correctly defines the Bondi frame, it is
important to note that there are many other supermomenta that have been
proposed in the literature, e.g., the Bondi-Sachs (BS), Geroch (G), and
Geroch-Winicour (GW) supermomenta [46, 47, 31, 17]

$\displaystyle\Psi^{\text{BS}}(u)$
$\displaystyle\equiv\Psi_{2}+\sigma\dot{\bar{\sigma}},$ (32a)
$\displaystyle\Psi^{\text{G}}(u)$
$\displaystyle\equiv\Psi_{2}+\sigma\dot{\bar{\sigma}}+\left(\eth^{2}\bar{\sigma}-\bar{\eth}^{2}\sigma\right),$
(32b) $\displaystyle\Psi^{\text{GW}}(u)$
$\displaystyle\equiv\Psi_{2}+\sigma\dot{\bar{\sigma}}-\bar{\eth}^{2}\sigma.$
(32c)

While all of these agree on their $\ell<2$ modes,151515The $\ell<2$ modes of
the supermomenta in Eq. (32) all correspond to the usual Bondi four-momentum.
the Moreschi supermomentum is the supermomentum that can most easily be used
to map to the super rest frame for reasons discussed in [15, 16, 17, 48].

It is important to note, however, that

$\displaystyle\Delta\Psi^{\text{M}}$
$\displaystyle\equiv\int_{u_{1}}^{u_{2}}\dot{\Psi}^{\text{M}}(u)\,du$
$\displaystyle=\int_{u_{1}}^{u_{2}}\left[\left(\dot{\Psi}_{2}+\left[\sigma\ddot{\bar{\sigma}}+\eth^{2}\dot{\bar{\sigma}}\right]\right)+|\dot{\sigma}|^{2}\right]du$
$\displaystyle=\int_{u_{1}}^{u_{2}}|\dot{\sigma}|^{2}\,du$ (33)

since

$\displaystyle\dot{\Psi}_{2}=-\left[\sigma\ddot{\bar{\sigma}}+\eth^{2}\bar{\sigma}\right]$
(34)

by the Bianchi identities. Consequently, since Eq. (A) is proportional to the
radiated energy, this means that the Moreschi supermomentum can never be made
zero so long as there is energy radiated in gravitational waves.161616This
also corresponds to the electric part of the null memory [18].

With $\Psi^{\text{M}}$, a similar procedure as the one presented in Sec. IV.2
for mapping to the PN BMS frame can then be performed. But, instead of
minimizing the $L^{2}$ norm of the difference of NR and PN strain waveforms,
one would simply find the proper supertranslation that minimizes the $L^{2}$
norm of $\Psi^{\text{M}}$ at a certain time or over some finite time interval,
such as the ringdown phase.

## References

  * [1] SXS Gravitational Waveform Database, http://www.black-holes.org/waveforms.
  * Boyle  _et al._ [2019] M. Boyle _et al._ , The SXS Collaboration catalog of binary black hole simulations, Class. Quant. Grav. 36, 195006 (2019), arXiv:1904.04831 [gr-qc] .
  * Jani  _et al._ [2016] K. Jani, J. Healy, J. A. Clark, L. London, P. Laguna, and D. Shoemaker, Georgia Tech Catalog of Gravitational Waveforms, Class. Quant. Grav. 33, 204001 (2016), arXiv:1605.03204 [gr-qc] .
  * Healy  _et al._ [2017] J. Healy, C. O. Lousto, Y. Zlochower, and M. Campanelli, The RIT binary black hole simulations catalog, Class. Quant. Grav. 34, 224001 (2017), arXiv:1703.03423 [gr-qc] .
  * Healy and Lousto [2020] J. Healy and C. O. Lousto, Third RIT binary black hole simulations catalog, Phys. Rev. D 102, 104018 (2020), arXiv:2007.07910 [gr-qc] .
  * Woodford  _et al._ [2019] C. J. Woodford, M. Boyle, and H. P. Pfeiffer, Compact Binary Waveform Center-of-Mass Corrections, Phys. Rev. D 100, 124010 (2019), arXiv:1904.04842 [gr-qc] .
  * Christodoulou [1970] D. Christodoulou, Reversible and irreversible transforations in black hole physics, Phys. Rev. Lett. 25, 1596 (1970).
  * Iozzo  _et al._ [2021a] D. A. B. Iozzo, N. Khera, L. C. Stein, K. Mitman, M. Boyle, N. Deppe, F. Hébert, L. E. Kidder, J. Moxon, H. P. Pfeiffer, M. A. Scheel, W. Throwe, and S. A. Teukolsky, Comparing remnant properties from horizon data and asymptotic data in numerical relativity, In Production (2021a).
  * Bondi  _et al._ [1962] H. Bondi, M. G. J. van der Burg, and A. W. K. Metzner, Gravitational waves in general relativity. 7. Waves from axisymmetric isolated systems, Proc. Roy. Soc. Lond. A 269, 21 (1962).
  * Sachs [1962] R. K. Sachs, Gravitational waves in general relativity. 8. Waves in asymptotically flat space-times, Proc. Roy. Soc. Lond. A 270, 103 (1962).
  * Campiglia and Laddha [2014] M. Campiglia and A. Laddha, Asymptotic symmetries and subleading soft graviton theorem, Phys. Rev. D 90, 124028 (2014), arXiv:1408.2228 [hep-th] .
  * Campiglia and Laddha [2015] M. Campiglia and A. Laddha, New symmetries for the Gravitational S-matrix, JHEP 04, 076, arXiv:1502.02318 [hep-th] .
  * Boyle [2016] M. Boyle, Transformations of asymptotic gravitational-wave data, Phys. Rev. D 93, 084031 (2016), arXiv:1509.00862 [gr-qc] .
  * Gómez López and Quiroga [2017] L. A. Gómez López and G. D. Quiroga, Asymptotic structure of spacetime and the Newman-Penrose formalism: a brief review, Rev. Mex. Fis. 63, 275 (2017), arXiv:1711.11381 [gr-qc] .
  * Moreschi [1988] O. M. Moreschi, Supercenter of Mass System at Future Null Infinity, Class. Quant. Grav. 5, 423 (1988).
  * Moreschi and Dain [1998] O. M. Moreschi and S. Dain, Rest frame system for asymptotically flat space-times, J. Math. Phys. 39, 6631 (1998), arXiv:gr-qc/0203075 .
  * Dain and Moreschi [2000] S. Dain and O. M. Moreschi, General existence proof for rest frame systems in asymptotically flat space-time, Class. Quant. Grav. 17, 3663 (2000), arXiv:gr-qc/0203048 .
  * Mitman  _et al._ [2020] K. Mitman, J. Moxon, M. A. Scheel, S. A. Teukolsky, M. Boyle, N. Deppe, L. E. Kidder, and W. Throwe, Computation of displacement and spin gravitational memory in numerical relativity, Phys. Rev. D 102, 104007 (2020), arXiv:2007.11562 [gr-qc] .
  * Mitman  _et al._ [2021] K. Mitman _et al._ , Adding gravitational memory to waveform catalogs using BMS balance laws, Phys. Rev. D 103, 024031 (2021), arXiv:2011.01309 [gr-qc] .
  * Blackman  _et al._ [2017a] J. Blackman, S. E. Field, M. A. Scheel, C. R. Galley, D. A. Hemberger, P. Schmidt, and R. Smith, A Surrogate Model of Gravitational Waveforms from Numerical Relativity Simulations of Precessing Binary Black Hole Mergers, Phys. Rev. D 95, 104023 (2017a), arXiv:1701.00550 [gr-qc] .
  * Blackman  _et al._ [2017b] J. Blackman, S. E. Field, M. A. Scheel, C. R. Galley, C. D. Ott, M. Boyle, L. E. Kidder, H. P. Pfeiffer, and B. Szilágyi, Numerical relativity waveform surrogate model for generically precessing binary black hole mergers, Phys. Rev. D 96, 024058 (2017b), arXiv:1705.07089 [gr-qc] .
  * Varma  _et al._ [2019] V. Varma, S. E. Field, M. A. Scheel, J. Blackman, L. E. Kidder, and H. P. Pfeiffer, Surrogate model of hybridized numerical relativity binary black hole waveforms, Phys. Rev. D 99, 064045 (2019), arXiv:1812.07865 [gr-qc] .
  * Williams  _et al._ [2020] D. Williams, I. S. Heng, J. Gair, J. A. Clark, and B. Khamesra, Precessing numerical relativity waveform surrogate model for binary black holes: A Gaussian process regression approach, Phys. Rev. D 101, 063011 (2020), arXiv:1903.09204 [gr-qc] .
  * Iozzo  _et al._ [2021b] D. A. B. Iozzo, M. Boyle, N. Deppe, J. Moxon, M. A. Scheel, L. E. Kidder, H. P. Pfeiffer, and S. A. Teukolsky, Extending gravitational wave extraction using Weyl characteristic fields, Phys. Rev. D 103, 024039 (2021b), arXiv:2010.15200 [gr-qc] .
  * Moxon  _et al._ [2020] J. Moxon, M. A. Scheel, and S. A. Teukolsky, Improved Cauchy-characteristic evolution system for high-precision numerical relativity waveforms, Phys. Rev. D 102, 044052 (2020), arXiv:2007.01339 [gr-qc] .
  * Geroch  _et al._ [1973] R. Geroch, A. Held, and R. Penrose, A space-time calculus based on pairs of null directions, J. Math. Phys. 14, 874 (1973).
  * Moreschi [1986] O. M. Moreschi, On angular momentum at future null infinity, Class. Quantum Gravity 3, 503 (1986).
  * Boyle  _et al._ [2020] M. Boyle, D. Iozzo, and L. C. Stein, moble/scri: v1.2 (2020).
  * Boyle [2013] M. Boyle, Angular velocity of gravitational radiation from precessing binaries and the corotating frame, Phys. Rev. D 87, 104006 (2013), arXiv:1302.2919 [gr-qc] .
  * Boyle  _et al._ [2014] M. Boyle, L. E. Kidder, S. Ossokine, and H. P. Pfeiffer, Gravitational-wave modes from precessing black-hole binaries, arXiv:1409.4431 [gr-qc] .
  * Dray and Streubel [1984] T. Dray and M. Streubel, Angular momentum at null infinity, Class. Quant. Grav. 1, 15 (1984).
  * Dray [1985] T. Dray, Momentum flux at null infinity, Class. Quant. Grav. 2, L7 (1985).
  * Flanagan and Nichols [2017] E. E. Flanagan and D. A. Nichols, Conserved charges of the extended Bondi-Metzner-Sachs algebra, Phys. Rev. D 95, 044002 (2017), arXiv:1510.03386 [hep-th] .
  * Compère  _et al._ [2020] G. Compère, R. Oliveri, and A. Seraj, The Poincaré and BMS flux-balance laws with application to binary systems, JHEP 10, 116, arXiv:1912.03164 [gr-qc] .
  * [35] https://data.black-holes.org/waveforms/extcce_catalog.html.
  * Pollney and Reisswig [2011] D. Pollney and C. Reisswig, Gravitational memory in binary black hole mergers, Astrophys. J. Lett. 732, L13 (2011), arXiv:1004.4209 [gr-qc] .
  * [37] https://www.black-holes.org/code/SpEC.html.
  * Moxon  _et al._ [tion] J. Moxon, N. Deppe, L. E. Kidder, G. Lovelace, K. Mitman, M. A. Scheel, S. Teukolsky, and W. Throwe, The spectre cauchy characteristic evolution system for rapid, precise waveform extraction.
  * Deppe  _et al._ [2020] N. Deppe, W. Throwe, L. E. Kidder, N. L. Fischer, C. Armaza, G. S. Bonilla, F. Hébert, P. Kumar, G. Lovelace, J. Moxon, E. O’Shea, H. P. Pfeiffer, M. A. Scheel, S. A. Teukolsky, I. Anantpurkar, M. Boyle, F. Foucart, M. Giesler, D. A. B. Iozzo, I. Legred, D. Li, A. Macedo, D. Melchor, M. Morales, T. Ramirez, H. R. Rüter, J. Sanchez, S. Thomas, and T. Wlodarczyk, SpECTRE (2020).
  * Sarbach and Tiglio [2001] O. Sarbach and M. Tiglio, Gauge invariant perturbations of Schwarzschild black holes in horizon penetrating coordinates, Phys. Rev. D 64, 084016 (2001), arXiv:gr-qc/0104061 .
  * Regge and Wheeler [1957] T. Regge and J. A. Wheeler, Stability of a Schwarzschild singularity, Phys. Rev. 108, 1063 (1957).
  * Zerilli [1970] F. J. Zerilli, Effective potential for even parity Regge-Wheeler gravitational perturbation equations, Phys. Rev. Lett. 24, 737 (1970).
  * Boyle and Mroue [2009] M. Boyle and A. H. Mroue, Extrapolating gravitational-wave data from numerical simulations, Phys. Rev. D 80, 124045 (2009), arXiv:0905.3177 [gr-qc] .
  * Ashtekar and Streubel [1981] A. Ashtekar and M. Streubel, Symplectic geometry of radiative modes and conserved quantities at null infinity, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 376, 585 (1981).
  * Virtanen  _et al._ [2020] P. Virtanen _et al._ , SciPy 1.0–Fundamental Algorithms for Scientific Computing in Python, Nature Meth. 17, 261 (2020), arXiv:1907.10121 [cs.MS] .
  * Geroch [1977] R. Geroch, Asymptotic structure of space-time, in  _Asymptotic Structure of Space-Time_, edited by F. P. Esposito and L. Witten (Springer US, Boston, MA, 1977) pp. 1–105.
  * Geroch and Winicour [1981] R. P. Geroch and J. Winicour, Linkages in general relativity, J. Math. Phys. 22, 803 (1981).
  * Gallo and Moreschi [2014] E. Gallo and O. M. Moreschi, Intrinsic angular momentum for radiating spacetimes which agrees with the Komar integral in the axisymmetric case, Phys. Rev. D 89, 084009 (2014), arXiv:1404.2475 [gr-qc] .