Compact Binary Waveform Center-of-Mass Corrections

TL;DR

This paper analyzes the center-of-mass motion in binary black hole simulations, identifying it as largely a gauge effect that causes mode mixing in gravitational waveforms, and proposes methods to improve c.m. correction techniques.

Contribution

It investigates the physical versus gauge origins of c.m. motion and introduces a new waveform-based method for optimal c.m. correction.

Findings

Large c.m. motion is mainly a gauge effect, not physical.

Mode mixing persists even at large distances from the source.

Waveform-based c.m. correction methods can improve waveform accuracy.

Abstract

We present a detailed study of the center-of-mass (c.m.) motion seen in simulations produced by the Simulating eXtreme Spacetimes (SXS) collaboration. We investigate potential physical sources for the large c.m. motion in binary black hole simulations and find that a significant fraction of the c.m. motion cannot be explained physically, thus concluding that it is largely a gauge effect. These large c.m. displacements cause mode mixing in the gravitational waveform, most easily recognized as amplitude oscillations caused by the dominant (2,$\pm$2) modes mixing into subdominant modes. This mixing does not diminish with increasing distance from the source; it is present even in asymptotic waveforms, regardless of the method of data extraction. We describe the current c.m.-correction method used by the SXS collaboration, which is based on counteracting the motion of the c.m. as measured by the trajectories of the apparent horizons in the simulations, and investigate potential methods to improve that correction to the waveform. We also present a complementary method for computing an optimal c.m. correction or evaluating any other c.m. transformation based solely on the asymptotic waveform data.

Equations

$$h(t)\to$$

$$\approx h(t)+{\partial }_{t}h(t)(\vec{\alpha }+\vec{\beta }t)\cdot \hat{n},$$

$$h(t,\theta ,\varphi ):=h_{+}(t,\theta ,\varphi )-\mathrm{i}{h}_{\times }(t,\theta ,\varphi ).$$

$$h(t,\theta ,\varphi )=\sum _{l,m}{h}^{l,m}(t){}_{-2}{Y}_{l,m}(\theta ,\varphi ),$$

$${\vec{x}}_{\text{c.m.}}=\frac{m_a}{M}\overrightarrow{x_a}+\frac{m_b}{M}\overrightarrow{x_b},$$

$$\Xi (\vec{\alpha },\vec{\beta })={\int }_{t_i}^{t_f}|{\vec{x}}_{c.m.}-(\vec{\alpha }+\vec{\beta }t){|}^{2}dt.$$

$$⟨{\vec{x}}_{\text{c.m.}}⟩=\frac{1}{t_f-t_i}{\int }_{t_i}^{t_f}{\vec{x}}_{\text{c.m.}}(t)dt,$$

$$⟨t{\vec{x}}_{\text{c.m.}}⟩=\frac{1}{t_f-t_i}{\int }_{t_i}^{t_f}t{\vec{x}}_{\text{c.m.}}(t)dt.$$

$$\vec{\alpha }=\frac{4(t_f^2+t_f{t}_i+t_i^2)⟨{\vec{x}}_{\text{c.m.}}⟩-6(t_f+t_i)⟨t{\vec{x}}_{\text{c.m.}}⟩}{(t_f-t_i{)}^2},$$

$$\vec{\beta }=\frac{12⟨t{\vec{x}}_{\text{c.m.}}⟩-6(t_f+t_i)⟨{\vec{x}}_{\text{c.m.}}⟩}{(t_f-t_i{)}^2}.$$

$$\vec{\mu }=\{\begin{array}{ll}\vec{\alpha }& \text{if }|\vec{\alpha }|>|\vec{\alpha }+\vec{\beta }{t}_{\text{merger}}|,\\ \vec{\alpha }+\vec{\beta }{t}_{\text{merger}}& \text{if }|\vec{\alpha }|\le |\vec{\alpha }+\vec{\beta }{t}_{\text{merger}}|.\end{array}$$

$${\chi }_{\text{eff}}=\frac{c}{GM}\left(\frac{{\vec{S}}_a}{m_a}+\frac{{\vec{S}}_b}{m_b}\right)\cdot \frac{\vec{L}}{|\vec{L}|}=\frac{{\vec{\chi }}_a{m}_a+{\vec{\chi }}_b{m}_b}{M}\cdot \frac{\vec{L}}{|\vec{L}|},$$

$${\chi }_p=\frac{c}{B_{1}G{m}_a^2}\max (B_1|{\vec{S}}_{a,\perp }|,B_2|{\vec{S}}_{b,\perp }|).$$

$$⟨L_{(a}{L}_{b)}⟩:={\int }_{S^2}{L}_{(a}{\left\{{\mathcal{T}}_{\vec{\alpha },\vec{\beta }}[h]\right\}}^{\ast }{L}_{b)}\left\{{\mathcal{T}}_{\vec{\alpha },\vec{\beta }}[h]\right\}dA,$$

$${\Phi }_{z^{\prime }}(t_{i,2})-{\Phi }_{z^{\prime }}(t_{i,1})=\phi .$$

$${\Phi }_{z^{\prime }}(t_{i+1,1})-{\Phi }_{z^{\prime }}(t_{i,1})=\pi .$$

$$\begin{gathered} {\mu }_i(t,\theta ,\varphi )=\sum _{l,m}[s_i^{l,m}+{\dot{s}}_i^{l,m}(t-t_{i,c}) \\ +\left(a_i^{l,m}+{\dot{a}}_i^{l,m}(t-t_{i,c})\right){\mathrm{e}}^{i\sigma (m,l){\Phi }_{z^{\prime }}(t)}]{}_{-2}{Y}_{l,m}(\theta ,\varphi ). \end{gathered}$$

$$\sigma (m,l)=\{\begin{array}{ll}1& |m|<l,\\ -1& |m|=l.\end{array}$$

$$s_i^{l,m}=(-1{)}^l{\bar{s}}_i^{l,-m},{\dot{s}}_i^{l,m}=(-1{)}^l{\bar{\dot{s}}}_i^{l,-m},$$

$$a_i^{l,m}=(-1{)}^{l+1}{\bar{a}}_i^{l,-m},{\dot{a}}_i^{l,m}=(-1{)}^{l+1}{\bar{\dot{a}}}_i^{l,-m}.$$

$${\rm Y}\left(\vec{\alpha },\vec{\beta }\right)$$

$$=\underset{\mathbf{s},\dot{\mathbf{s}},\mathbf{a},\dot{\mathbf{a}}}{\min }\sum _i{\int }_{t_{i,1}}^{t_{i,2}}\sum _{l,m}{\left|\mathcal{R}{\left\{{\mathcal{T}}_{\vec{\alpha },\vec{\beta }}\left[h\right]\right\}}^{l,m}-{\mu }_i^{l,m}\right|}^{2}dt.$$

$$\vec{G}=-\frac{48}{5}\frac{G^4}{c^7{r}^4\omega }{m}_a^2{m}_b^2(m_a-m_b)\hat{\lambda },$$

$$r_{\text{recoil}}=|\vec{G}|=\frac{48}{5}\frac{G^4}{c^7{r}^4\omega }{m}_a^2{m}_b^2(m_a-m_b)$$

$$r_{\mathrm{measured}}=\frac{1}{t_f-t_i}{\int }_{t_i}^{t_f}|{\vec{x}}_{\text{c.m.}}-(\vec{\alpha }+\vec{\beta }t)|dt.$$

$${\vec{a}}_{\text{c.m.}}=\frac{d^2}{d{t}^2}\left({\vec{x}}_{\text{c.m.}}-(\vec{\alpha }+\vec{\beta }t)\right).$$

$$\frac{d\vec{p}}{d\Omega dt}=\frac{c^2}{G}\frac{R^2}{16\pi }{\left|\frac{dh}{dt}\right|}^2\hat{r},$$

$$\hat{r}=\sqrt{\frac{2\pi }{3}}(Y_{1,-1}-Y_{1,1},i{Y}_{1,-1}+i{Y}_{1,1},\sqrt{2}{Y}_{1,0}).$$

$$\begin{gathered} {\dot{p}}_j=\frac{c^2}{G}\frac{R^2}{16\pi }\sum _{l,l^{\prime },m,m^{\prime }}{\hat{r}}_j^{1,m^{\prime }-m}{\dot{h}}^{l,m}{\bar{\dot{h}}}^{l^{\prime },m^{\prime }}\sqrt{\frac{3(2l+1)(2l^{\prime }+1)}{4\pi }} \\ \times (-1{)}^{m^{\prime }}\left(\begin{array}{ccc}l& l^{\prime }& 1\\ m& -m^{\prime }& m^{\prime }-m\end{array}\right)\left(\begin{array}{ccc}l& l^{\prime }& 1\\ 2& -2& 0\end{array}\right), \end{gathered}$$

$$\hat{n}$$

$$\hat{k}$$

$$\hat{\lambda }$$

$${\vec{c}}_1={\vec{x}}_{\text{c.m.}}-(\vec{\alpha }+\vec{\beta }t).$$

$${\vec{c}}_2={\vec{x}}_{\text{c.m.}}-{\Delta }_n\hat{n}-{\Delta }_{\lambda }\hat{\lambda }-{\Delta }_k\hat{k},$$

Full text

I Introduction

Binary black hole (BBH) systems have been studied for decades, beginning with analytic work and branching out into numerical relativity. With the introduction of gravitational-wave detectors, particularly LIGO, the pursuit of BBH gravitational waveforms has intensified in an attempt to create and fill vast waveform template banks. Gravitational waveforms created through numerical relativity are generally the most accurate waveforms available, and are used for parameter estimation and to compare and improve semianalytic and analytic models of BBHs, which in turn are used for gravitational-wave detection and parameter estimation Lovelace et al. (2016); Healy et al. (2018); Kumar et al. (2014).

While numerical-relativity waveforms are the most accurate BBH waveforms, there are concerns regarding their validity, accuracy, and reproducibility. There have been numerous discussions on how to measure the accuracy of a numerical relativity simulation and some sources of error in the simulations have been investigated, including numerical truncation errors, error due to extraction at finite radius or imperfect extrapolation to infinite radius, and errors between simulations of different lengths that otherwise have identical parameters Blackman et al. (2017); Chu et al. (2016); Boyle and Mroué (2009).

Though not strictly a source of error like those named above, there are also consequences due to the gauge freedom of general relativity that may be confused for errors if not properly understood, and will effectively become sources of error if ignored Gualtieri et al. (2008); Kelly and Baker (2013); Boyle (2016a). Gauge freedom in general relativity affects all numerical relativity simulations, and thus far no numerical relativity results in the literature have been in a completely specified gauge due to this inherent gauge freedom. While a full accounting of the effects of general gauge freedom is beyond the scope of this work, we will address the translation and boost degrees of freedom. Reference Boyle (2016a) identified these transformations as important for counteracting the observed motion of the center of mass in simulations produced by the Simulating eXtreme Spacetimes (SXS) collaboration. Here, we expand on that analysis, using the recently updated catalog of 2,018 SXS simulations SXS ; Boyle and others (SXS Collaboration), and investigating possible improvements to the correction method.

A translation $\vec{\alpha}$ and a boost $\vec{\beta}$ will transform the waveform $h$, measured at some point distant from the source, as

[TABLE]

where $\hat{n}$ is the direction to the observer from the source Boyle (2016a). Note that this is independent of the distance to the source; even the asymptotic waveform will exhibit this dependence regardless of any extrapolation, Cauchy-characteristic extraction, or similar techniques that may be applied to the data. We can understand this intuitively by thinking about a sphere surrounding the source. If we displace the source away from the center of the sphere, an emitted signal will arrive at the part of the sphere closest to the source before it will arrive at the opposite side of the sphere. The difference in arrival times is independent of the radius of the sphere; it only depends on the size of the displacement. The additional term in Eq. (1b) introduces an angular dependence that is not generally included in waveform models.

Figure 1 demonstrates the effects of c.m. motion for two systems from the SXS catalog. The most striking example is the upper pair of panels, which show data from a nonprecessing system with mass ratio 1.23. On physical grounds, there is nothing to suggest modulations in the mode amplitudes on the orbital timescale; this is a relatively symmetric system with very low eccentricity. The dominant physical behavior on the orbital timescale is simply rotation, which should have no effect on the amplitudes of the modes, along with a secular increase toward merger. Nonetheless, the raw waveform data from the simulation (thin dark lines in the upper left panel) shows very clear amplitude modulations of the subdominant modes on the orbital timescale. These modulations—like the c.m. trajectories seen in the upper right panel—show no signs of convergence with increasing numerical resolution in the simulation, even though the initial data for each resolution is created from identical high-resolution initial data.

As we discuss below, the c.m. motions found in the SXS catalog are effectively random and apparently independent of any physical parameters of the systems. Therefore, they comprise an essentially random source of unmodeled and unphysical contributions to waveforms from numerical relativity. In particular, they are not systematic; the modulations found in waveforms for one set of physical parameters will be uncorrelated with the modulations in waveforms even for nearly identical physical parameters. Clearly, expecting waveform models such as effective-one-body (EOB)  Buonanno and Damour (1999); Pan et al. (2014); Taracchini et al. (2014); Babak et al. (2017); Bohé et al. (2017); Cotesta et al. (2018), phenomenological Ajith et al. (2007); Khan et al. (2019); Nagar et al. (2018), and surrogate models Varma et al. (2019) to accurately represent these features across a range of physical parameters is tantamount to expecting them to fit large, discontinuous, random signals.

However, by simply compensating for the inertial part of the measured c.m. motion, the modulations can be almost completely eliminated (thick transparent lines in the upper left panel). It is notable that the c.m. only drifts by roughly $0.1M$ during almost the entire inspiral for the system shown in the upper panels of Fig. 1, but still has such a drastic effect on the waveform’s modes. Even though this is only a gauge choice—which we have been trained to consider irrelevant in principle—in practice, gauge choices must be made consistently and systematically for the waveforms to be really useful. In this sense, we might suggest that the c.m.-centered gauge is really an optimal choice.

A more difficult comparison is for the precessing system shown in the lower panels of Fig. 1. The precession already mixes the modes drastically, leading to a complicated waveform with pronounced amplitude modulations, even after c.m. correction. Clearly the c.m. correction changed the data, but it is not obvious that we can say it was a change for the better—at least from looking at this plot alone. To make the comparison more quantitative, we introduce a new measure of a waveform’s “simplicity” in Sec. IV. Essentially, this quantity measures the residual when the waveform is modeled by simple linear-in-time amplitudes in the corotating frame Boyle (2013). The value of this residual is 117 times smaller for the c.m.-corrected data than for the raw data in this precessing system, showing that the corrected waveform is clearly and objectively better in this sense at least.

To address the miscalculation of the c.m. and its correction, this paper is organized as follows:

1. i

In Sec. II, we discuss the current definition of the c.m. and the consequences this definition and its use have on SXS gravitational waveforms. 2. ii

In Sec. III, we discuss the current method for correcting waveform data and selecting an optimal gauge. Any correlations found between simulation parameters and the c.m. correction factors are discussed. 3. iii

In Sec. IV, we discuss a quantitative method for evaluating the “correctness” of the gauge a waveform is currently in. We compare the waveform data in its original, unoptimized gauge to the c.m.-corrected gauge described in Sec. III. 4. iv

In Sec. V, we discuss how we may improve the definition of the c.m. to find a better correction, potentially leading to a further optimized choice of gauge. This section also investigates alternative definitions of the c.m. with a focus on potential physical causes of c.m. motion like that seen in Fig. 1, including post-Newtonian definitions and considerations of linear-momentum recoil. 5. v

Finally, we present our findings and results in Sec. VI.

II The need for c.m. corrections

One of the primary concerns with BBH simulations with regards to gravitational-wave astronomy is the validity of their gravitational waveforms. Above all, the output from a BBH simulation should result in a reliable, reproducible waveform that can then be released for public usage. In the case of the SXS collaboration, many of the waveforms produced are also compressed into a catalog that is released to LIGO for data analysis and waveform comparisons with their gravitational wave detector data.

Gravitational waveforms in the SXS catalog are given in terms of the gravitational-wave strain $h$, or the Weyl component $\Psi_{4}$. In regards to detecting gravitational waves, $h$ and $\Psi_{4}$ contain the same information, and the analysis and corrections applied in this work may be applied to either with the same results. For simplicity, we will focus on $h$.

Waveforms from SXS are represented by mode weights, or amplitudes, for spin-weighted spherical harmonics (SWSHs). The gravitational-wave strain may be represented by the transverse-traceless projection of the metric perturbation caused by the gravitational waves at time $t$ and location $(\theta,\phi)$ relative to the binary, and can be combined into a single complex quantity, given by

[TABLE]

For each slice in time, the combined perturbation $h(t,\theta,\phi)$ is measured on the coordinate sphere. The angular dependence of this measurement can then be expanded in SWSHs. The quantity $h(t,\theta,\phi)$ has a spin weight of $-2$ Newman and Penrose (1966), and may be represented as

[TABLE]

where the complex quantities $h^{l,m}(t)$ are referred to as modes or mode weights, and are much more convenient when analyzing BBH than the total perturbation in any particular direction Boyle et al. ; Ajith et al. (2007). Spin-weighted Spherical harmonics are further discussed in Appendix A.

The expansion in Eq. (3) depends on orientation of the spherical coordinates $\theta,\phi$ and their origin. The customary choice places the binary at the origin with the binary’s initial orbital plane coinciding with the equatorial plane $\theta=\pi/2$. For comparable mass, nonprecessing binaries, the quadrupolar $(l,m)=(2,\pm 2)$ modes then dominate the waveform. While the $h^{2,\pm 2}$ modes are dominant, it is important to consider the behavior of the other modes present in the waveform. The other modes may not be used directly for BBH detection currently, as they are much smaller in magnitude compared to the $h^{2,\pm 2}$ modes for most systems, but can be important for parameter estimation Varma and Ajith (2017); Lange et al. (2017); Chatziioannou et al. (2019). Additionally, there are proposals for using higher-order modes in BBH searches Pekowsky et al. (2013); Calderón Bustillo et al. (2017); Harry et al. (2018). Higher-order modes are also useful for verifying the reliability and potentially the accuracy of the waveform. If the shape, variability, magnitude, or any other characteristic of the higher order, or subdominant, modes are found to not suitably match with theory, then this could indicate a possible flaw in the simulation.

One clear issue is the coordinate system, or gauge choice, for the simulation, as spherical harmonics and hence SWSHs depend on the coordinates. The center chosen for the simulation is the c.m. of the system, as calculated and set in the initial data. It is expected that the c.m. will move slightly throughout the simulation; however, large movements are not expected and suggest a flaw in the choice of gauge. If the c.m. moves significantly, there is mode mixing Boyle (2016a). The dominant effect Boyle (2016a) is leaking of the $h^{2,\pm 2}$ modes of BBH waveforms into the higher modes, and this leakage can be at least partially removed through c.m. drift corrections, as described in Sec. III.

Mode mixing is manifested in the waveforms as oscillating amplitudes, which can clearly be seen in the left column panels of Fig. 1, especially for the $(2,1)$, $(3,1)$, and $(3,3)$ modes in the top panel for SXS:BBH:0314. Precessing simulations, like SXS:BBH:0622, are expected to have some amplitude modulations purely due to the orientation of the system. The worse the c.m. calculation is for a simulation, the more altered the SWSH representations are, and the worse the mode mixing becomes.

It is easily seen in Fig. 1 that the applied c.m. correction removes what we will find to be unphysical waveform amplitude modulations for nonprecessing systems, however for precessing systems it is not so clear or obvious. Therefore, especially for precessing systems, a quantitative method is required for evaluating the “correctness” of the c.m. and gauge. This is discussed in Sec. IV.

The current definition used during BBH simulations for the c.m. is the usual Newtonian definition:

[TABLE]

where $M=m_{a}+m_{b}$ is the total mass of the system, $m_{a}$ and $m_{b}$ are the Christodoulou masses Christodoulou (1970) of the primary and secondary black holes respectively, and $\vec{x}_{a}$ and $\vec{x}_{b}$ are the coordinates of the centers for the primary and secondary black hole respectively. This is a Newtonian expression for the c.m., and from output of the simulations like in Fig. 1, we know is not a perfect description of the optimal c.m.. The tracking of the c.m. throughout the simulations can be seen for SXS:BBH:0314 and SXS:BBH:0622 in the right column panels of Fig. 1.

The c.m. motion is an effect of the initial data. One aspect of the initial data construction method proposed in Ref. Ossokine et al. (2015) is the elimination of Arnowitt-Deser-Misner (ADM), or spatial, linear momentum in the initial data for precessing systems, namely enforcing $\vec{P}_{\mathrm{ADM}}=0$. The work done in Ref. Ossokine et al. (2015) proposed a new, and now adopted, method for calculating and constructing the initial data for BBH simulations. The improved method for calculating initial data has far-reaching effects in Spectral Einstein (SpEC) simulations and most of the simulations in the SXS simulation catalog were completed using this relatively new method. This had the effect of reducing specific components of mode mixing as seen in the gravitational waveforms, however as showcased in Fig. 1, significant mode mixing is still present.

As is further discussed in Sec. V.2, linear-momentum recoil is an expected physical contribution to the motion of the c.m.. However, unphysical contributions to the linear momentum in the initial data of simulations introduce unphysical motion in the c.m., essentially imparting spurious linear-momentum kicks. By controlling the linear momentum and removing it, this effect from the initial data is removed. However, even for simulations with initial data constructed using the method described in Ref. Ossokine et al. (2015), significant translations and boosts, and the resulting mode mixing, are still present in the gravitational waveforms. This warrants further investigation into the c.m. motion and the application of a c.m. correction.

It had been suggested in Refs. Boyle (2016a); Ossokine et al. (2015) that much of the c.m. motion depicted in the right column panels of Fig. 1, and seen in all SXS simulations, was largely unphysical and could be removed from the data. The c.m. correction used to remedy the unphysical c.m. motion is discussed in the following section. Additionally, there are alternative definitions of the c.m. and physical effects that are expected to cause the c.m. to move, or imply that the c.m. is not moving at all. The more obvious of these physical effects are post-Newtonian (PN) corrections for the c.m. which may include effects explaining the c.m. motion, and linear-momentum recoil from the system. PN and linear momentum contributions are examined in Secs. V.1 and V.2.

III Centre-of-mass correction method

Previous work Boyle (2016a); Ossokine et al. (2015) suggests that the c.m. motion is largely a result of gauge choice. Therefore, understanding the c.m. correction begins with understanding the permissible gauge transformations. More specifically, we are interested in the gauge transformations that will affect the waveform measured by distant observers. Because a gravitational-wave detector will typically be very far from the source, only the asymptotic behavior of the waves is generally considered relevant—specifically at future null infinity, $\mathscr{I}^{+}$.

While the asymptotic gauge of waveforms from numerical relativity has not been extensively investigated, it is certainly fair to say that no results in the literature thus far have been in a completely specified gauge. Even the strongest claims of “gauge-invariant” asymptotic waveforms Handmer et al. (2015) are only invariant modulo the infinite-dimensional Bondi-Metzner-Sachs (BMS) gauge group Bondi et al. (1962); Sachs (1962), which is the asymptotic symmetry group corresponding to the Bondi gauge condition. An important feature of Bondi gauge is that the gravitational waves measured by any distant inertial observer (at least over a duration short compared to the distance to the source) are approximately given by the asymptotic metric perturbation at fixed spatial coordinates as a function of retarded time in some member of this gauge class—and conversely, any such function corresponds to a signal that could be measured by some distant inertial observer Strominger (2018). Essentially, we might think of Bondi gauge to be “as simple as possible, but not simpler” for the purposes of gravitational-wave detection. Because the BMS group alters the waveform while preserving Bondi gauge, we consider it to be the fundamental symmetry group relevant to gravitational-wave modeling.111Other possible gauge choices exist. For example, Newman-Unti gauge Newman and Unti (1962) is closely related to Bondi gauge, and is invariant under the same asymptotic symmetry algebra Barnich and Lambert (2012). More generally, it is not even clear that waveforms from numerical relativity are actually expressed in either of these well-defined gauge classes, in which case more general gauge transformations may be of interest. Ultimately, the gauge freedoms relevant to counteracting c.m. motion are simply space translations and boosts. As long as these transformations are allowed, this discussion of c.m. motion remains relevant. Previous work Taylor et al. (2013) suggests that SXS waveforms are consistent with waveforms in Bondi gauge, though further research is warranted.

Because BMS transformations preserve the inertial property of observers, we cannot expect to counteract all of the c.m. motion seen in Fig. 1—particularly the cyclical behavior. However, in addition to the cyclical behavior, these coordinate tracks begin with some overall displacement from the origin, and then drift away from that initial location over the entire course of the inspiral. Thus, we expect that a space translation and a boost are needed to negate the effects of some of the c.m. motion. In particular, we will choose the translation $\vec{\alpha}$ and boost $\vec{\beta}$ to minimize the average of the square of the distance between the measured c.m. and the origin.222This measure will be invariant under time translation and rotation, which are generally dealt with separately during gravitational-wave analysis, so we simply ignore those degrees of gauge freedom. Furthermore, it is not at all clear how a higher-order supertranslation should affect the coordinates close to the center of a simulation, and so we leave discussion of more general supertranslations to future work.

III.1 Choosing the translation and boost

We follow Appendix E of Ref. Boyle (2016a) in choosing the translation $\vec{\alpha}$ and boost $\vec{\beta}$ to minimize the average square of the distance between the c.m. measured in the raw data and the origin of the corrected frame. That is, we choose $\vec{\alpha}$ and $\vec{\beta}$ to minimize the function

[TABLE]

It is not hard to find the minimum of this quantity analytically. We define two moments of the c.m. position

[TABLE]

Then, the minimum of Eq. (5) is achieved with

[TABLE]

We then apply this transformation to the asymptotic waveform using the method described in Ref. Boyle (2016a). Note that this rests on implicit assumptions about how directly comparable the coordinates of the apparent-horizon data and the asymptotic coordinates are. For example, this assumes that the time coordinate of the apparent-horizon data and the asymptotic retarded-time coordinate are equal. While there is no rigorous motivation for this assumption, the results of Sec. IV bear out its approximate validity.

Using this minimization method, the c.m. offsets for every public waveform in the SXS catalog have been corrected in the waveform data. The first instance of the c.m. corrections to waveforms in the SXS public waveform catalog was in January of 2017. Center-of-mass corrected waveform data is recommended over non corrected data in all cases, and corresponding files are listed in the SXS public waveform catalog as files ending in CoM.h5. An overview of the c.m. correction values is shown in Fig. 2.

It is clear from the upper-left panel of Fig. 1 that c.m. removal is “helpful” in the sense that it reduces the amplitude oscillations, which are not expected on physical grounds. Unfortunately, this by-eye analysis is not quantitative, and it is not clear how it would apply to a precessing system, as seen in the lower-left panel of the same figure. We discuss a better measure of how the waveform quality is impacted by c.m. corrections in Sec. IV.

III.2 Choosing the integration region

The determination of $\vec{\alpha}$ and $\vec{\beta}$ is made over a subset of the total simulation time, from $t_{i}$ to $t_{f}$ [see Eq. (5)]. Choosing different values of $t_{i}$ and $t_{f}$ may affect the resulting $\vec{\alpha}$ and $\vec{\beta}$ values. For the corrections performed on the SXS catalog, a standard subset of the simulation time was chosen. All waveforms had their c.m.-correction values calculated from $t_{i}=t_{\text{relax}}$, the “relaxation” time after which the initial transients have dissipated, to 10% before the end of the inspiral: $t_{f}=0.9t_{\text{merger}}$. These time bounds were chosen to avoid including periods of junk radiation as well as the merger and ringdown stages.

However, changing $t_{i}$, $t_{f}$ by small amounts could change the c.m. correction values. As there is epicyclical motion of the c.m. (as seen in Fig. 1, for example), changing the beginning or ending time may cause the resulting $\vec{\alpha}$ and $\vec{\beta}$ to change, depending on where $t_{i}$ and $t_{f}$ fall on an epicycle. For example, if $t_{i}$ and $t_{f}$ are separated by an integer number of epicycles, then we might expect any effect from the epicycles on the calculation of $\vec{\alpha}$ and $\vec{\beta}$ to cancel out. However, if $t_{i}$ and $t_{f}$ are separated by a noninteger number of epicycles, especially by a half-integer number of epicycles, the epicycles may induce significant bias in $\vec{\alpha}$ and $\vec{\beta}$. The overall number of epicycles included in the calculation of $\vec{\alpha}$ and $\vec{\beta}$ may also affect how sensitive they are to this bias.

Here, we compare the size of the c.m. correction using the standard prescription to the size of the correction when $t_{i}$ and/or $t_{f}$ are changed by half an orbit. This will give us some idea of the stability of the c.m.-correction procedure. However, it must be noted that, at a larger scale, the choices of $t_{i}$ and $t_{f}$ are quite arbitrary. For some purposes, it may be preferable to choose those values to range over only the first half of the inspiral, or even just the ringdown stage. The values used in the SXS catalog were chosen for robustness and easy reproducibility.

To simplify the comparison, we describe the c.m. motion using the quantity $\vec{\mu}$, which gives the most distant position of the corrected origin of coordinates throughout the inspiral, relative to the origin used in the simulation. Specifically, we can define $\vec{\mu}$ according to

[TABLE]

For 96% of the simulations in the SXS catalog, we find that $\vec{\mu}=\vec{\alpha}+\vec{\beta}t_{\text{merger}}$. The 4% of simulations with $\vec{\mu}=\vec{\alpha}$ have no apparent correlations with system parameters, and are effectively random. We also introduce subscripts, so that $\vec{\mu}_{00}$ is the result of this calculation when using the original values of $t_{i}$ and $t_{f}$; $\vec{\mu}_{10}$ is the result when moving $t_{i}$ later by half an orbit; $\vec{\mu}_{01}$ is the result when moving $t_{f}$ later by half an orbit; and $\vec{\mu}_{11}$ is the result when moving both $t_{i}$ and $t_{f}$ later by half an orbit.

In Fig. 3, we examine $\max_{jk}|\vec{\mu}_{00}-\vec{\mu}_{jk}|$ as a measure of how robust the c.m. corrections are with respect to these small adjustments in the choices of $t_{i}$ and $t_{f}$. In the great majority of systems the c.m. changes by less than $10^{-2}M$. This is, for example, just one tenth the size of the displacements seen in the upper panels of Fig. 1. The median change is $3.1\times 10^{-3}M$, and in all cases is smaller than the median value of $|\vec{\mu}_{00}|$ itself, which is $6.9\times 10^{-2}M$. The systems with the largest c.m. corrections in the SXS catalog change by fractions of a percent, suggesting that the results are certainly stable in the cases where applying a c.m. correction is most important. There are several cases where the fractional change is greater than 100%, though these are systems with relatively small values of $\vec{\mu}_{00}$. The median fractional change is $4.3\%$. It is also notable that the data points separate roughly into three groups. The group in the lower left corner of Fig. 3 is comprised exclusively of equal-mass nonspinning simulations with various eccentricities, though several of these are also found in the central group. The central group is where all equal-mass simulations with equal but nonzero spins are found, which includes ten systems with significant precession. Every other type of system is in or near the largest group, on the upper right.

III.3 Correlations between c.m. correction values and physical

parameters

Along with having c.m. corrected the waveforms, we have also performed an analysis of the values of the boosts and translations needed by each simulation in the SXS public waveform catalog.

No obvious correlations can be seen in Fig. 2 between spin aligned and precessing systems. We also show a more in-depth correlation plot in Fig. 4, taking more of the simulation parameters into consideration. It can be seen that typically precessing simulations may have larger overall c.m. displacement, $\vec{\delta}=\vec{\alpha}+\vec{\beta}t_{\text{merger}}$, and that larger boost values $\vec{\beta}$ correspond with larger overall displacement values $\vec{\delta}$ for both spin aligned and precessing systems.

Outside of the correlations between the boost $\vec{\beta}$ and total displacement $\vec{\delta}$ of the c.m., there does not appear to be any other strong correlations present for the current SXS simulation catalog. It was expected that precessing, high mass ratio, and eccentric systems should have vastly different c.m. correction values than spin aligned, low mass ratio, and more circular systems, however no such correlations are present with this data set.

For Fig. 4, we use the eccentricity $e$, number of orbits, and mass ratio $q$ reported by SpEC at the end of the simulation. We calculate the effective spin Ajith et al. (2011); Santamaria et al. (2010); Abbott et al. (2016)

[TABLE]

and an effective precession parameter Schmidt et al. (2015); Abbott et al. (2016)

[TABLE]

Here, $M=m_{a}+m_{b}$ is the total mass of the system, $\vec{S_{i}}=G/c\,\vec{\chi}_{i}m_{i}^{2}$ is the angular momentum of the $i$-th black hole and $\vec{\chi}_{i}$ its dimensionless spin, $B_{1}=2+3m_{b}/2m_{a}$, $B_{2}=2+3m_{a}/2m_{b}$, and the subscript $\perp$ indicates the quantity perpendicular to the orbital angular momentum $\vec{L}$, e.g., $\vec{S}_{a,\perp}=\vec{S}_{a}-(\vec{S}_{a}\cdot\hat{L})\hat{L}$. Note that $\chi_{p}$ gives a measure of how much a system is precessing during a simulation.

IV Quantifying c.m. Correction Using Waveforms Alone

Any discussion of c.m. based on the positions of the individual black holes will suffer from the same fundamental ambiguity: reliance on coordinates—specifically in the highly dynamical region between the two black holes—that are subject to unknown gauge ambiguities. The only region of the spacetime where the gauge freedom is limited in any useful sense is the asymptotic region, in which we assume the only freedom is given by the BMS group (described in Sec. III). While there are many suggestions in the literature Newman and Penrose (1966); Adamo and Newman (2010); Adamo et al. (2012); Deadman and Stewart (2009); Kozameh et al. (2008); Szabados (2004); Helfer (2007); Dray and Streubel (1984); Ashtekar and Winicour (1982); Flanagan et al. (2016) for using asymptotic information to specify the asymptotic gauge more narrowly, they all require more information than is available from most catalogs of numerical-relativity waveforms—such as additional Newman-Penrose quantities or more precise characterization of the asymptotic behavior of the various fields.

Here, we present a simplistic but effective measure of c.m. effects that can be applied exclusively to asymptotic waveform data $h$ or $\Psi_{4}$. The basic idea is that we expect to be able to decompose a waveform measured in c.m.-centered coordinates into modes that are, at least for small portions of the inspiral, given by a slowly changing complex amplitude times a complex phase that varies proportionally with the orbital phase. When the waveform is decomposed in off-center coordinates, those well-behaved modes mix, so that the amplitude and phase do not behave as expected. Therefore, we will attempt to model a given waveform in a sort of piecewise fashion that assumes the expected behavior, and simply measure the residual between the model and the waveform itself. For a given transformation applied to the waveform, we will minimize the residual by adjusting the parameters to the model while keeping the waveform fixed. The smaller the minimized residual, the more accurately the waveform with that transformation can be modeled in this simple way, and the more nearly we expect that the waveform is decomposed in c.m.-centered coordinates. Roughly speaking, we can think of this as a measure of the “simplicity” of the waveform, which is not only in line with our basic expectations for waveforms in the appropriate coordinates, but also a useful measure of how accurately simple waveform models (EOB, surrogate, etc.) will be able to capture features in the numerical waveforms. This criterion is obviously totally distinct from any criteria involving the BH positions, but is important precisely because it provides a complementary way of looking at the data. Finding agreement between the results of this method and another will lend support to the idea that the other method is suitable.

IV.1 Defining the method

We now describe this method more precisely. The initial inputs are some translation $\vec{\alpha}$ and boost $\vec{\beta}$ that we wish to evaluate. We transform the waveform by those inputs and denote the result $\mathcal{T}_{\vec{\alpha},\vec{\beta}}[h]$. We then transform to a “corotating frame”, which is a time-dependent frame chosen so that the waveform in that frame is varying as slowly as possible Boyle (2013). Only the angular velocity of this frame, $\vec{\Omega}$, is determined by the condition that the waveform vary slowly; it is integrated in time to obtain one such frame Boyle (2017), but the result is only unique up to an overall rotation. We choose that overall rotation so that the $\vec{z}\,^{\prime}$ axis of the final corotating frame is aligned as nearly as possible throughout the inspiral portion of the waveform with the dominant eigenvector O’Shaughnessy et al. (2011); Boyle et al. (2011) of the matrix

[TABLE]

where $L_{a}$ is the usual angular-momentum operator. This still leaves the frame defined only up to an overall rotation about $\vec{z}\,^{\prime}$, but such a rotation will have no effect on our results. The transformed waveform in this corotating frame will be denoted $\mathcal{R}\left\{\mathcal{T}_{\vec{\alpha},\vec{\beta}}\left[h\right]\right\}(t,\theta,\phi)$, though we will usually suppress the parameters, and may decompose the angular dependence in terms of SWSH mode weights as usual. This is the quantity we will be attempting to model.

For the model itself, we first break the inspiral up into smaller spans of time; we will be modeling the waveforms using simple linear-in-time approximations, so we cannot expect to accurately reproduce the nonlinear evolution over a very long portion of the inspiral using just one such model. The relevant measurement of the waveform’s dynamical behavior is the angular velocity of the corotating frame. More specifically, we define $\Omega_{z^{\prime}}=\vec{\Omega}\cdot\vec{z}\,^{\prime}$, and use that to determine a phase333This phase is loosely related to the orbital phase of the binary. The angular velocity $\vec{\Omega}$, however, is defined solely with respect to the waveform at $\mathscr{I}^{+}$, and entirely without reference to any quantities at finite distance in the system. Nonetheless, for reasonably well-behaved coordinate systems, we would expect it to agree roughly with the orbital phase deduced from the trajectories of the black holes, especially during the early inspiral regime. $\Phi_{z^{\prime}}=\int\Omega_{z^{\prime}}\,dt$. An obvious span of time would be a single cycle of this phase, which would include enough data so that the fit would actually reflect the behavior of the waveform, but not so much that we would expect a poor fit due to evolution on the inspiral timescale. However, we will essentially be fitting oscillatory terms with linear models. In the simple case of fitting a line to a basic sine function, it is not hard to see that the optimal line has the expected slope of zero—independent of the phase of the sine function—when the fit region is such that the argument of the sine function goes through a phase change of $\varphi\approx 8.9868$ [or other solutions of $\varphi=2\tan(\varphi/2)$]. Therefore, we select each span of time so that it extends over a phase $\Phi_{z^{\prime}}$ of approximately $\varphi$, thereby determining the difference in time between $t_{i,1}$ and $t_{i,2}$ so that they satisfy

[TABLE]

We find that this choice does drastically reduce the oscillations in the optimal fit parameters as we shift the fitting window. While the individual time spans extend over this range, we find that remaining effects from oscillation are minimized by selecting successive time spans to be separated by half of a period—so that $\Phi_{z^{\prime}}$ changes by exactly $\pi$ between $t_{i,1}$ and $t_{i+1,1}$:

[TABLE]

So that the model may be reasonably accurate, without encountering excessive numerical noise or excessively dynamical behavior at merger, we limit the region over which we choose these time spans to be the central 80% of time between the “relaxation time” listed in the waveform metadata and the time of maximum signal power in the waveform. This establishes $t_{0,1}$, and all successive times can be computed from that using Eqs. (12) and (13).

Now, we model the waveform “piecewise” on these spans of time, though the pieces are overlapping. The advantage of transforming the waveform as described above is that each mode separates Boyle et al. into two parts that are symmetric and antisymmetric under reflection along the $z$ axis. The symmetric part varies on an inspiral timescale because the primary rotational behavior has been factored out by transforming to the corotating frame; the antisymmetric part is mostly due to spin-orbit coupling and therefore varies most rapidly by a complex phase with frequency equal to the rotational frequency of the frame itself, though possibly with opposite sign. We model these two parts separately as simple linear-in-time complex quantities, with an additional phase-evolution term for the antisymmetric parts. For each time span $i$, we write

[TABLE]

Here, each of the coefficients $s_{i}^{l,m}$, $\dot{s}_{i}^{l,m}$, $s_{i}^{l,m}$, and $\dot{a}_{i}^{l,m}$ is a complex constant, we use $t_{i,c}=(t_{i,1}+t_{i,2})/2$ to mitigate degeneracy between the constant-in-time and linear-in-time terms, and the function $\sigma$ is given by

[TABLE]

These signs are chosen because they represent the dominant behavior of the corresponding terms in the data. Note that the symmetry properties imply that once the quantities $s_{i}^{l,m}$, etc., are chosen for positive $m$, they are automatically known for negative $m$ from the relations

[TABLE]

Because the $m=0$ modes of the SXS waveforms we use are not considered reliable Boyle and Mroué (2009); Boyle and others (SXS Collaboration), we simply ignore those modes in both the model and the data. That is, the sum in Eq. (14) does not include any $m=0$ modes. If the sum over modes extends from $l=2$ to some maximum $l=L$, the total number of (real) degrees of freedom in this model is $4L(L+1)-8$ for each span of time. While the data we use contains up to $l=8$, the highest-order modes contribute little to the result, and drastically increase the number of degrees of freedom in the problem (and therefore the time taken to optimize the model). Therefore, we use only up to $l=6$ in constructing the model and evaluating the residual, reducing the degrees of freedom from 280 to 160 per time span. Finally, because this is still such a large number of degrees of freedom, we limit the evaluation to only the first two and last two time spans; we find that including the rest has no significant effect on the result, but vastly increases the amount of processing time required. This leaves us with a manageable 640 degrees of freedom in this model.

Now, using this model, we define the objective function

[TABLE]

We will use this function in two ways: first, to simply evaluate $\Upsilon$ for given values of $(\vec{\alpha},\vec{\beta})$, where those values are obtained from the methods described in other sections; second, to minimize $\Upsilon$ over possible values of $(\vec{\alpha},\vec{\beta})$ to find the optimum c.m. correction.

IV.2 Results for the standard c.m.-correction method

We can now compare the value of $\Upsilon$ defined in Eq. (17) for all the waveforms discussed in the previous sections. First, we compare its value $\Upsilon_{\mathrm{raw}}$ in the raw data to its value $\Upsilon_{\mathrm{0PN}}$ using $\vec{\alpha}$ and $\vec{\beta}$ as given by Eqs. (7), where $\vec{x}_{\mathrm{c.m.}}$ is given by the Newtonian (0PN) formula. The latter corresponds to the technique actually used in the current SXS data, for waveforms found in the SXS simulation catalog with file names ending in CoM. The results of this comparison are shown in Fig. 5. One unusually short simulation (SXS:BBH:1145 Lovelace (2019b)) in the SXS catalog did not have enough GW cycles to evaluate $\Upsilon$ properly, leaving a total of 2,017 systems shown in these figures. The vast majority of systems improve significantly by this measure. The notable feature is that even though the naive 0PN c.m. trajectory is so fundamentally different from $\Upsilon$, this plot suggests that they agree in the sense that the 0PN correction improves the waveforms for all but three systems—and even for those three the change is very small.444These three systems are unusual, in that they are quite short (having 13 to 15 orbits before merger, compared to an average of 22), and have eccentricities (0.215 and higher) that place them among the 12 most eccentric in the SXS catalog. Furthermore, the magnitude of the change in $\Upsilon$ for each of them is very small—in the lowest percentile for the entire catalog—which suggests that the negative results may be consistent with numerical error, and in any case are not cause for much concern.

We can also actively optimize $\Upsilon$ over the values of $\vec{\alpha}$ and $\vec{\beta}$. The results are shown in Fig. 6. Naturally, $\Upsilon$ improves in every case because it is specifically being optimized. In Fig. 6, we see the pattern that the vast majority of systems are changing by small fractions. In this case, there are just three systems in which $\Upsilon$ changes at the percent level. These are some of the same systems that changed the most in going from the raw data to the 0PN-corrected data. These particular systems also happen to be extremely long, with significant overall accelerations during the inspiral. This suggests that the corrections will be sensitive to the precise span of times over which the corrections are being made, which may explain why they continue to change so much by optimization. However, as discussed in Sec. III.2, changing the beginning and ending fractions does not significantly change $\vec{\alpha}$ and $\vec{\beta}$. Nevertheless, the overall scale of the changes seen in this plot suggests once again that the naive 0PN c.m. correction is achieving near-optimal results in the vast majority of cases.

V Improving the c.m. correction

V.1 Post-Newtonian c.m. definition

To characterize the motion seen in the c.m. in the raw simulation data, an obvious first step to finding a more accurate definition of the c.m. during the simulation is to try low orders of post-Newtonian (PN) corrections. Note that the c.m. should be near the origin of the coordinate system, with minimal motion around the origin from linear-momentum recoil, as discussed in Sec. V.2. 1PN order corrections are analytically trivial, and are zero for circular systems. However, SXS simulations are not perfectly circular and so we investigate 1PN and 2PN order corrections to the c.m.. We implement the PN corrections given by Eq (4.5) in de Andrade et al. (2001). This formalism goes up to 3.5 PN for the c.m. vector in time. The 2PN version of the expression that we used can be found in Appendix B.

Using this formalism, the effects of the correction on the c.m. are small but measurable at 1PN and significantly different at 2PN for many systems. We are not concerned with the coordinates of the c.m. itself changing, but of the c.m. correction values changing, as discussed in Sec. III. In general, we assume that the c.m. drifts in a linear fashion and can be corrected with a translation $\vec{\alpha}$ and boost $\vec{\beta}$, Our results from the 1PN and 2PN analysis are shown in Fig. 7, which shows the relative difference between $\vec{\alpha}$, $\vec{\beta}$, and total c.m. displacement $\vec{\delta}=\vec{\alpha}+\vec{\beta}t_{\text{merger}}$ for the 1PN (top panels) correction to the c.m. and the 2PN (bottom panels) correction to the c.m.. The 1PN corrections show small changes for most simulations. However, the 2PN correction shows more sizable changes in the c.m. correction values, which may indicate that including at least up to 2PN corrections to the c.m. will give better accuracy either to the c.m. itself or to the correction factors.

To see any potential correlations with large 1PN or 2PN corrections and simulation parameters, we compared the relative difference in the c.m. corrections to the eccentricity $e$ of the system. As shown in Fig. 8, no correlations between the magnitude of the 1PN and 2PN corrections to the c.m. and the eccentricity are apparent, despite the definition of the 1PN contribution to the c.m. being dependent on $e$.

Using the method described in Sec. IV, comparisons of $\Upsilon$, as defined in Eq. (17), between 1PN, 2PN, and the original c.m. correction method (dubbed 0PN), can be found in Fig. 9. The striking feature of these plots is that a significant number of systems actually have larger values of $\Upsilon$ when including either of these corrections. While it is reassuring that the majority of systems in each case only exhibit quite small changes—changes of order $10^{-4}$ or less in a quantity that already improved significantly from the raw data—the 1PN and 2PN corrections plots include a large group of systems that change at the percent level. These systems also happen to be the same systems that changed most drastically in going from the raw data to the 0PN-corrected data (found near the upper-right corner of the plots), and are particularly biased towards increasing values of $\Upsilon$. That is to say, it appears that the 1PN and 2PN corrections do worst for the most extreme systems. This should not come as a great surprise, since those systems tend to be the ones with the most extreme mass ratios and precession, so that post-Newtonian analysis is also expected to be at its least accurate.

V.2 Linear-momentum recoil

Any binary with asymmetric components will emit net linear momentum in the form of gravitational waves, which will cause a recoil of the binary itself. As the system rotates, the direction of recoil will also rotate, pushing the c.m. roughly in a circle Blanchet and Faye (2019); Fitchett (1983). In principle, this effect could cause the epicyclical motion apparent in the c.m. trajectories, which is further characterized in Sec. V.4. To see if recoil is responsible, we use methods described in Refs. Blanchet and Faye (2019); Fitchett (1983); Boyle et al. ; Gerosa et al. (2018) to investigate the size of the linear-momentum recoil implied by the gravitational-wave emission in these systems.

As shown in the right-column panels of Figure 1, the c.m. motion follows an overall linear track with additional epicyclical motion. The linear motion of the c.m. is well understood, and discussed in Sec. III. For this analysis, we assume that the linear part of the c.m. motion may be removed from the data without loss of information, leaving the epicyclical behavior about the coordinate origin.

Blanchet and Faye Blanchet and Faye (2019) consider the motion of the c.m. from linear momentum flux and the flux of the c.m. itself up to 3.5PN order to calculate the instantaneous c.m. motion induced by these effects, finding the position of the c.m. relative to its average location over an orbit to be

[TABLE]

[cf. Eq. (6.9) in Ref. Blanchet and Faye (2019)] leading to a circular motion of the c.m. with radius

[TABLE]

for a system comprised of nonspinning BHs in a circular orbit. Here, $r=|\vec{x}_{a}-\vec{x}_{b}|$ is the distance between the two black holes, $\omega=\sqrt{GM/r^{3}}$ is the Newtonian orbital frequency, and $\hat{\lambda}$ is the unit vector in the direction of motion of $m_{a}$ [cf. Eq. (25)]. An earlier analysis using a simplified model and lower-order approximations can be found in Ref. Fitchett (1983).

Subtracting the motion described by Eq. (19) out of the c.m. motion and comparing the radius of the measured motion to that of Eq. (19) immediately shows that the measured c.m. motion—specifically the epicyclical motion—is significantly larger than what can be explained by linear-momentum and c.m. reactions for all simulations. Figure 10 shows a comparison across the SXS public waveform catalog for the ratio between the measured c.m. radius about the coordinate origin and the estimated c.m. radius given by Eq. (19). The measured c.m. radius is calculated by averaging the distance of the c.m. at time $t$ away from the line $\vec{\alpha}+\vec{\beta}t$ between times $t_{i}$ and $t_{f}$:

[TABLE]

The results show that the measured radius is typically hundreds of times larger than the radius implied by Blanchet and Faye’s analysis.

It is also possible to go beyond the analysis in Ref. Blanchet and Faye (2019) by using the measured gravitational waves to compute the linear-momentum flux, and compare the acceleration that would cause to the measured acceleration of the c.m. in the simulation. The acceleration of the c.m., $\vec{a}_{\text{c.m.}}$ was calculated by taking two time derivatives of the coordinate position of the c.m. after removing the linear motion:

[TABLE]

The linear-momentum flux may be calculated from the gravitational radiation as Ruiz et al. (2008)

[TABLE]

where $\hat{r}$ is the direction from the source to the point in question, $R$ is the distance from the source to the observation sphere, and $\Omega$ represents all angles on the sphere. Integrating over all angles to find the total linear momentum flux $d\vec{p}/dt=\dot{\vec{p}}$, we need to expand $h$ into SWSHs as done in Eq. (3). We can also decompose $\hat{r}$ accordingly as

[TABLE]

Then integrating over all angles gives

[TABLE]

where the last two factors are Wigner 3-$\!j$ symbols. The sum over $m^{\prime}$ and most terms in the sum over $l^{\prime}$ can be eliminated using properties of the 3-$\!j$ symbols. Explicit expressions for the calculation of the linear momentum flux are given in Appendix C. Our calculation of the linear momentum flux can be found in the open-source package spherical_functions Sph .

Of course, there are multiple methods for calculating the linear momentum for a BBH system. We compared our method with the method proposed in Ref. Gerosa et al. (2018) and originally stated in Ref. Ruiz et al. (2008), and found that the two derivations give the same results to within numerical accuracy. This is unsurprising given that the methods are identical up to the choice of notation. In particular, the core definition of the total linear momentum flux given in Eq. (22) is the same as that given in Eq.(2.11) of Ref. Ruiz et al. (2008).

We compared the linear-momentum flux divided by the total mass of the system with the measured acceleration of the c.m. and find that the values for most of the SXS public waveform catalog do not agree. We further found that the c.m. acceleration for nonequal mass systems is consistently larger than the acceleration found through the linear momentum flux, confirming that most of the c.m. motion is not due to linear-momentum recoil.

An overview of the ratios between $\vec{a}_{\text{c.m.}}$ and $\dot{\vec{p}}/M$ can be seen in Figure 10. Note that linear-momentum recoil in equal mass or near-equal mass BBH systems is expected to be very small, and so these systems are isolated as a separate case in the middle panel of Figure 10.

V.3 Causes of unphysical c.m. motion

If the c.m. motion seen in the SXS BBH simulations cannot be explained by physical processes such as inclusion of PN terms or linear-momentum recoil, then why is it there? There are two potential causes for the appearance of unphysical or erroneous c.m. motion: (i) the presence of uncontrolled residual linear momentum in the initial data, and (ii) the emission of unresolved junk radiation at the beginning of the BBH simulation causing effectively random and unphysical coordinate kicks to the system. The presence of uncontrolled residual linear momentum was addressed and partially rectified in Ref. Ossokine et al. (2015), leading to a new method of creating initial data for the BBH simulations. However, this method was not used for all systems in the SXS catalog, and does not completely resolve the issues of spurious translations and boosts even when it is applied.

Another factor that seems to cause c.m. motion starting from early in the simulations appears to be junk radiation, which is an inherent part of BBH simulations. It is the radiation emitted when a BBH relaxes from its initial-data “snapshot”, which is only an approximation to the true state of a long binary inspiral at the time the simulation starts. Junk radiation is physical, in the sense that if the entire spacetime were actually in the configuration given by the initial data, it would indeed emit this radiation as the system evolved. However, it is not astrophysical, in the sense that no real system in the universe is expected to contain this type of radiation. Junk radiation contains high frequencies that are largely unresolved in BBH simulations because of limits of computational power and time. As the resolution increases in the simulations, more of the junk radiation is accurately treated. This is a potential cause for the difference in initial kicks of the c.m., as seen in Fig 2. Even in our higher-resolution simulations, not all of the junk radiation is accounted for. This failure to resolve all of the junk radiation possibly contributes to the c.m. kicks observed in SpEC BBH simulations.

Fortunately, kicks from the emission of unresolved junk radiation can be corrected using the gauge transformations discussed in the previous section. The large epicyclical motion in the c.m. cannot be accounted for using a BMS gauge transformation like the linear motion of the c.m.. The cause for such large, seemingly unphysical epicyclic motion is unknown, and is left for future work.

V.4 Epicycle quantification

As seen in the right column panels of Figure 1, the c.m. motion has both a linear and epicyclic component. The linear component of the c.m. motion has been discussed and is already considered in the current c.m. correction technique.

The size of the epicyclic motion in the c.m. cannot be solely explained by linear-momentum recoil, as discussed in Sec. V.2. The leftmost plot in Figure 10 illustrates that the expected radius of the epicycles from linear-momentum recoil, even on an approximate basis, is orders of magnitude smaller than what is actually seen, given by $r_{\mathrm{measured}}$ defined in Eq. (20). The actual size of $r_{\mathrm{measured}}$ is fairly consistent across the SXS simulation catalog regardless of simulation parameters or initial data construction, and tends to be between 0.01 and 0.1 with an average value of 0.026 in simulation units. There is not a significant change in the distribution or magnitude of $r_{\mathrm{measured}}$ between the 0PN, 1PN, or 2PN representations of the c.m..

We assume that the epicycle motion seen in the c.m. after the translation and boost are applied is from the calculated c.m., $\vec{x}_{\text{c.m.}}$, being displaced from the optimal c.m. by a small amount. “Displaced” here means that we assume $\vec{x}_{\text{c.m.}}$ is displaced from the optimal c.m. along the separation vector $\vec{r}_{ab}=\vec{x}_{a}-\vec{x}_{b}$, and not in any other direction.

Regardless of the origin of these epicycles, removing them to calculate the linear c.m. correction should improve the quality of the waveforms. Their removal should reduce the error associated with the averaging done to calculate the translation and boost correction values. As discussed regarding the choice of beginning and ending times in Sec. III.2, the presence of large epicycles has the potential to affect the reproducibility and reliability of the calculation of $\vec{\alpha}$ and $\vec{\beta}$. Removing the epicycles, or at the very least minimizing them, to calculate a more accurate c.m. correction should further reduce mode mixing. Of course, epicycle motion cannot be completely subtracted from the waveform itself as this would require an acceleration correction, which is not an allowed BMS transformation.

To accurately describe the epicycles, we need to define the corotating coordinate frame. For our simulations, we have three unit vectors that describe the rotating coordinate frame:

[TABLE]

where $\hat{n}$ points along the separation vector $\vec{r}_{ab}$, $\hat{k}$ points along the orbital angular velocity $\vec{\omega}$, and $\hat{\lambda}$ points along the direction of rotation.

This leads us to a potential method for epicycle removal. First, we calculate the estimated spatial translation $\vec{\alpha}$ and boost $\vec{\beta}$ from the original c.m. $\vec{x}_{\rm c.m.}$ using the current averaging method

[TABLE]

Note that $\langle|\vec{c}_{1}|\rangle=r_{\rm measured}$ as defined in Eq. (20), using the angle-bracket notation to denote averaging over time, as for the moments of the c.m. position in Eq. (6). We then calculate the corresponding corotating coordinate frame unit vectors given in Eq. (25). These unit vectors are then applied to the original c.m. $\vec{x}_{\text{c.m.}}$ as

[TABLE]

where $\Delta_{n}=\langle\vec{c}_{1}(t)\cdot\hat{n}(t)\rangle$, $\Delta_{\lambda}=\langle\vec{c}_{1}(t)\cdot\hat{\lambda}(t)\rangle$, $\Delta_{k}=\langle\vec{c}_{1}(t)\cdot\hat{k}(t)\rangle$, are the time averaged projections of the linearly corrected c.m. onto the rotational coordinate system unit vectors. The epicycle-corrected c.m. $\vec{c}_{2}$ is then fed back into the averaging method to get the final spatial-translation $\vec{\alpha}_{\mathrm{epi}}$ and boost $\vec{\beta}_{\mathrm{epi}}$ values, which can then be applied to the waveform data.

A visual representation of our epicycle-correction method is shown in Figure 11, which uses the spin aligned system SXS:BBH:0314 and the precessing system SXS:BBH:0622 as sample cases. As seen in the upper two panels on the right, the removal of the epicyclic motion as calculated by the time-averaged values $\Delta_{n}$, $\Delta_{\lambda}$, and $\Delta_{k}$ greatly diminishes the large size of the epicycles and allows for potentially better optimization of the c.m.-correction values $\vec{\alpha}$ and $\vec{\beta}$. In these same panels, it can also be seen that not all of the epicycle motion is removed by our method, and in particular there are larger deviations towards the beginning and end of the simulation data, which are a result of our time-averaged method capturing most but not all of the epicycle motion. Specifically, the epicycle radius tends to grow with time in most (81% of) SXS simulations.

Figure 12 shows $\vec{c}_{1}$ [Eq. (26)] and the projections of $\vec{c}_{1}$ onto the $\hat{n}$, $\hat{\lambda}$ unit vectors. These panels show very similar behavior between resolutions of the same simulation, implying that the cause of the size of the epicycles is not random, from initial conditions, or the junk radiation phase, and that the randomness of the initial kick has been completely removed by the correction applied in Eq. (26).

Applying this method to the BBH simulations in the SXS public simulation catalog, we can calculate $\vec{\alpha}_{\mathrm{epi}}$ and $\vec{\beta}_{\mathrm{epi}}$ c.m. correction values. Figure 13 compares the usual “0PN” c.m. correction with the epicycle-removed values. The values plotted involve $\vec{\mu}$, defined in Eq. (8), which is the largest displacement between the origin of coordinates in the simulation and the corrected origin. We see that epicycle removal changes the c.m. correction values at a scale comparable to the changes caused by varying the end points of integration used to determine the c.m., as seen in Figure 3. The changes due to epicycle removal are generally somewhat smaller than changes due to varying end points. If we remove the epicycles before applying those variations, the changes seen in Figure 3 are reduced by a typical factor of two—though there is no apparent effect on roughly 10% of systems. Systems that change by more than 10% in this figure are also typically changing by more than the post-Newtonian changes shown in Figure 7.

Using the method outlined in Sec. IV, we also find that we cannot reliably conclude that the epicycle correction actually makes a significant improvement in the waveforms. The results of the $\Upsilon$ comparisons between the original c.m. correction method outlined in Sec. LABEL:sec:SOMCorrectionMethod and the epicycle removal method in this section can be found in Figure 14. This plot shows that approximately $30\%$ of simulations improve using the epicycle method, and approximately $70\%$ get worse. This is not enough of a benefit to warrant the use of the epicycle removal step in all simulations, and implies that the epicycle removal, at this stage of BBH simulations, is an unnecessary step in calculating and applying the c.m. correction, a somewhat disappointing conclusion.

As mentioned in Sec. III.2, epicycles are also a potential source of instability regarding choice of beginning and ending times $t_{i}$, $t_{f}$. Initial investigation into how the epicycle removal method affects changes in the c.m. correction values due to differing $t_{i}$, $t_{f}$ implies that the epicycle removal method does not diminish changes in the c.m. correction values. This is also an unintuitive and disappointing result, and may imply that other methods are required to calculate the c.m. correction values $\vec{\alpha}$ and $\vec{\beta}$ after epicycle removal or a different method for epicycle removal entirely. We leave such an investigation to future work.

V.5 Position of the c.m.

During the analysis of the epicycles present in the c.m., we investigated the position of the c.m. relative to the two black holes. We still assumed, that $\vec{x}_{\text{c.m.}}$ would lie along the separation vector between the two black holes or completely in the rotation direction as given in Eq. (18). As seen in Figure 15, this is not the case. The c.m. deviates significantly between the rotation vector and the separation vector between the two black holes but typically lies in the negative rotation direction $-\hat{\lambda}$, as predicted by Eq. (18) and the analysis in Ref. Blanchet and Faye (2019). The projection in the $\pm\hat{\lambda}$ direction, $\Delta_{\lambda}$, averages at $-0.44r_{\rm measured}$ when considering all simulations in the SXS Catalog. On average, the projection of $\vec{c}_{1}$ into the $\pm\hat{n}$ direction, $\Delta_{n}$, is smaller than $\Delta_{\lambda}$, with an average ratio $\Delta_{\lambda}/\Delta_{n}$ for spin aligned systems of $-1.43$ and $-2.48$ for precessing systems. The projection in the $\pm\hat{k}$ direction, $\Delta_{k}$, is significantly smaller than $\Delta_{\lambda}$, with an average ratio $\Delta_{\lambda}/\Delta_{k}$ of $-3.22\times 10^{3}$ for spin aligned systems and $-2.78$ for precessing systems.

Having typically most but not all of the corrected c.m. vector $\vec{c}_{1}$ in the direction of $-\hat{\lambda}$ does not have an obvious cause. This behavior could possibly indicate unaccounted spin-orbit effects on the c.m., unknown effects from unresolved junk radiation, or additional gauge effects that cannot be compensated for with BMS transformations. Attributes of $\vec{c}_{1}$ warrant further investigation, and are left to future work.

Correlations between $\Delta_{n}$, $\Delta_{\lambda}$, and $\Delta_{k}$ with pertinent simulation parameters are shown in Figure 15. This plot shows the correlations between $\Delta_{n}$, $\Delta_{\lambda}$, $\Delta_{k}$, $\chi_{\mathrm{eff}}$, $\chi_{p}$, and $q$. A few notable correlations are apparent, the most obvious being the correlation between $\Delta_{n}$ and $\Delta_{\lambda}$. Most simulations tend to have negative $\Delta_{\lambda}$ values that grow larger in magnitude with increasing $\Delta_{n}$, however there is also a cluster of aligned-spin simulations with $\Delta_{n}$, $\Delta_{\lambda}$ values close to zero. Additionally, $\Delta_{\lambda}$ becomes more negative with increasing $q$, for $q<5$, and there are some weak correlations between both c.m. position offsets $\Delta_{n}$, $\Delta_{\lambda}$ with $\chi_{\mathrm{eff}}$, but not with $\chi_{p}$. $\Delta_{k}$ does not appear to have any strong correlations. It is apparent that spin aligned simulations tend to have $\Delta_{k}$ values which are much smaller than for precessing simulations. One possible explanation for the c.m. to move out of the orbital plane is momentum flow between the gravitational fields and black holes Keppel et al. (2009), however we leave analysis regarding this mechanism to future work. It can also be seen in Figure 15 that larger $\Delta_{k}$ values cluster around $\Delta_{n}=0$ and $\Delta_{\lambda}<0$, which may only be due to $\Delta_{n}$ values being symmetric around [math] and $\Delta_{\lambda}$ values being mostly negative. No apparent correlations are present for other simulation parameters.

VI Conclusions

In this work, we have investigated the effects of c.m. motion waveforms, removed unphysical c.m. motion through allowed gauge transformations to the waveforms, and have investigated methods for improving the c.m. correction. Having unphysical motion in the c.m. causes mode mixing in the gravitational waveforms, and thus a power loss from the dominant $(2,\pm 2)$ modes to the less-dominant, higher-order modes—which is typically visible as amplitude modulation in the higher-order modes.

We found that the c.m. motion observed in the SXS simulations cannot be entirely accounted for by PN corrections or linear-momentum recoil. We also found that the motion of the c.m. does not lie along any one basis vector describing the rotating coordinate frame as defined in Eq. (25), and is offset from the estimated c.m. within and out of the orbital plane—which is not expected on physical grounds.

The current method for correcting the c.m. motion uses allowed BMS transformations, namely a spatial translation and boost that counteracts the linear motion from the c.m. and removes a large amount of the mode mixing from the waveform. The translation and boost are calculated for all simulations at all resolutions, as the c.m. motion is not consistent between different resolutions of the same system.

We attempted to improve the c.m. correction by developing a method to remove the large epicycles from the c.m. motion before calculating the BMS translation and boost. We found that the resulting changes to the translation and boost values were not significant and did not improve the waveforms compared to the originally calculated values.

Last, we introduced a complementary method to quantify the effect of the c.m. correction on the waveforms. We used this method to determine that PN corrections and the epicycle removal technique did not improve the c.m. correction transformations, and thus would not further improve the waveforms or accurately describe the c.m. physically.

Future work includes investigating spin-orbit effects on the c.m. and the peculiarity of the c.m. position. Further investigation is required specifically on the unaccounted for size of the radius of the epicycles seen in the c.m. motion, which may be due to unknown spin-orbit or unresolved junk radiation effects, and may be corrected with additional gauge transformations that minimize the epicycles.

Acknowledgements.

It is our pleasure to thank Saul Teukolsky for useful conversations and guidance on this project. We also thank Leo Stein for helpful insight and suggestions regarding notation. This project was supported in part by the Sherman Fairchild Foundation, by NSF Grant No. PHY-1606654 at Cornell, and by NSERC of Canada Grant No. PGSD3-504366-2017 as well as the Canada Research Chairs Program. The computations described in this paper were performed on the Wheeler cluster at Caltech, which is supported by the Sherman Fairchild Foundation and by Caltech, and on the GPC and Gravity clusters at the SciNet HPC Consortium, funded by the Canada Foundation for Innovation under the auspices of Compute Canada, the Government of Ontario, Ontario Research Fund–Research Excellence, and the University of Toronto.

Appendix A Spin-Weighted Spherical Harmonics

Spin-weighted spherical harmonics (SWSHs) are typically used to generalize the well-known standard spherical harmonics. Specifically, SWSHs provide a decomposition of general spin-weighted spherical functions (SWSFs) into a sum of SWSHs. Spin-weighted spherical functions themselves provide a vital way to study waves radiating from bounded regions, and so have an obvious and important application in gravitational-wave astronomy, which is the focus of this work. Spin-weighed spherical functions play two key roles in this field: (1) describing the magnitude of the wave given any direction of emission or observation, and (2) providing polarization information. There are a number of subtleties in defining SWSFs and hence SWSHs, including dependencies on the chosen coordinate system and the explicit definition of SWSFs.

The spin weight of a function is defined by how it transforms under rotation of the spacelike vectors $\Re(m^{\mu})$, $\Im(m^{\mu})$ where $m^{\mu}$ is a complex null vector tangent to the coordinate sphere $S^{2}$. The rotation of these spacelike vectors is given by

[TABLE]

A function $\eta$ is then said to have a spin weight $s$ if it transforms as

[TABLE]

In the case of gravitational waves, the metric perturbation $h$ has a spin weight of $-2$ Newman and Penrose (1966); Ajith et al. (2007) and this decomposition has been used in numerical relativity extensively.

The classic definition of SWSHs Newman and Penrose (1966) writes the functions in terms of spherical coordinates for $S^{2}$, giving them as explicit formulas using polar and azimuthal angles $(\theta,\phi)$ and using two integer variables that define the order of spherical harmonic to be used. The convention used here and in SXS is $Y_{l,m}$.

Spin-weighted spherical harmonics are thus classically defined as

[TABLE]

where $\eth$ is effectively a covariant differentiation operator in the surface of the sphere. $\eth$ is defined Newman and Penrose (1966) as

[TABLE]

when operating on some function $\eta$ that has a spin weight $s$.

The above classic method inherits an unfortunate dependency on the chosen coordinates. In particular, SWSHs cannot be written as functions on the sphere $S^{2}$; at best they can only be written as functions on coordinates of $S^{2}$. As such, SWSHs as defined in Eq. (30) do not transform among themselves under rotation of the sphere (or, equivalently, rotation of coordinates of the sphere). That is, a SWSH in a given coordinate system cannot generally be expressed as a linear combination of SWSHs in another coordinate system. A more correct method for defining SWSFs that does not inherit these coordinate-system dependencies is to represent them as functions from $\mathrm{Spin}(3)\approx\mathrm{SU}(2)$, which in turn maps onto $S^{2}$ Boyle (2016b). By forming a representation of $\mathrm{Spin}(3)$, SWSHs defined in this way do transform among themselves and still agree with the classic definition. SWSHs may then be defined as

[TABLE]

where $\mathfrak{D}$ is a Wigner matrix, which are representations of the spin group, and R is the $\mathrm{Spin}(3)$ argument. Taking R to be in the unit-quaternion representation of $\mathrm{Spin}(3)$, the $\mathfrak{D}$ matrices may be expressed as

[TABLE]

where $\rho_{1}={\max(0,m^{\prime}-m)}$, $\rho_{2}={\min(l+m^{\prime},l-m)}$, and $R_{s}$ and $R_{a}$ are the geometric projections of the quaternion into symmetric and antisymmetric parts under reflection along the $z$ axis, which are essentially complex combinations of components of the quaternion:

[TABLE]

The definition given in Eq. (32) is consistent with the definition of SWSHs typically used within the SXS collaboration and is the assumed formulation for this work. For further information on SWSFs and SWSHs, a comprehensive in-depth discussion of the history, details, and additional formulations of SWSHs can be found in Boyle (2016b).

Appendix B Post-Newtonian Correction to the c.m.

As discussed in Sec. V.1, we used the 1PN and 2PN corrections to the c.m. as outlined in Ref. de Andrade et al. (2001). The c.m. up to 2PN order is given in Eq. (4.5) of Ref. de Andrade et al. (2001) as

[TABLE]

where the superscript $i$ designates the vector component being considered; subscripts $a,b$ designate which object is being considered; $\vec{y}$ is the position of the body being considered; $r_{ab}=|\vec{y}_{a}-\vec{y}_{b}|$ is the distance between body $a$ and $b$; $\vec{v}$ is the velocity of the body being considered and likewise $v$ is the magnitude of the velocity; and $\vec{n}_{ab}=\vec{r}_{ab}/r_{ab}$. Parentheses here represent the scalar product of the interior values, e.g., $(n_{ab}v_{b})=\vec{n}_{ab}\cdot\vec{v}_{b}$.

Note that this representation of the c.m. position does not include an overall division by the total mass of the system, and so our calculations deviate from Eq. (35) only by including an overall denominator $M=m_{a}+m_{b}$.

Appendix C Linear Momentum Flux from $h^{l,m}$ modes

As mentioned in Sec. V.2, our calculation of the linear momentum flux from the simulations is based on the formalism outlined in Boyle et al. , which uses the SWSH structure of the gravitational strain $h$. Starting from Eq. (24), which gives the general form of the linear momentum flux in $h^{l,m}$ modes, we can evaluate the components of the linear momentum flux:

[TABLE]

where $x,y,z$ refer to the simulation coordinates, with the orbit typically lying in the $x-y$ plane, $R$ is the distance to the observation sphere, and $\dot{h}$, $\bar{\dot{h}}$ are the time derivative and its conjugate of the mode amplitudes. The matrices throughout the summations are Wigner 3-$\!j$ symbols.