Identifying Strong Gravitational-wave Lensing during the Second Observing Run of Advanced LIGO and Advanced Virgo

When gravitational waves (GWs) propagate near massive galaxies or galaxy clusters, similar to light, the GWs can be strongly lensed. If the massive galaxies or galaxy clusters are along the line of sight of the GW source, GW observatories are expected to see multiple images with a time delay of hours to weeks (Haris et al. 2018) from the same astrophysical source as long as both images are above the GW detection threshold. Based on predictions of the number of expected GW sources and the distribution of lenses in the universe, Li et al. (2018) and Oguri (2018) suggested that around one in 1000 events observed by Advanced LIGO and Advanced Virgo (Harry 2010; Aasi et al. 2015; Acernese et al. 2015; Abbott et al. 2018) at design sensitivity will be lensed. The rate computations typically assume that a single image is detected, and the majority of the lenses are galaxy lenses. The lensing rate is expected to be lower at O2 sensitivity and even lower when considering double images. Galaxy cluster lenses have been investigated by Smith et al. (2017, 2018, 2019), Dai et al. (2020), and Robertson et al. (2020), who found that the rate of galaxy cluster lensing is around 10⁻⁵ yr⁻¹ at O1 sensitivity, and given the small sensitivity improvement in O2, we do not expect that the rate has notably increased. Lensing event rates can also be inferred with the measured amplitude of the binary black hole (BBH) background as shown in Buscicchio et al. (2020) and Mukherjee et al. (2021).

Under the presence of a lens, the corresponding strongly lensed GW signal is magnified such that ${\rho }_{l}=\sqrt{\mu }\rho$ (Wang et al. 1996), where ρ_l and ρ are the signal-to-noise ratios (S/Ns) under the lensed and unlensed models, respectively, and μ is the relative magnification factor (Narayan & Bartelmann 1996; Hannuksela et al. 2019). Since the GW frequency evolution is not affected by strong lensing, the lensing magnification is equivalent to a scaling of the source luminosity distance by a factor of $1/\sqrt{\mu }$ (Wang et al. 1996; Dai et al. 2017). Thus, a loud and nearby GW source could potentially be lensed and thus appear to be more distant than it seems. Also, inferring the source frame chirp mass ${{ \mathcal M }}_{c}^{\mathrm{source}}$ relies on the detector frame chirp mass ${{ \mathcal M }}_{c}$ and redshift z: ${{ \mathcal M }}_{c}=(1+z){{ \mathcal M }}_{c}^{\mathrm{source}}$. In the geometrical optics regime, lensing does not affect frequencies; hence, the detector frame chirp mass is not affected. Rather, the inferred amplitude is larger (smaller) if the event is magnified (demagnified), and if we assume the event is an unlensed event in the background cosmology, the relation between distance (through the amplitude) and the redshift would lead us to infer a larger (smaller) source frame mass.

In addition, according to Takahashi & Nakamura (2003) and Dai & Venumadhav (2017), lensing shifts the original phase of the waveform by Δϕ in such a way that the shift is absorbed into the phase of the coalescence Δϕ_c in the case of GWs with the relation Δϕ = 2Δϕ_c (except if precession, eccentricity, or higher modes are present; Ezquiaga et al. 2020). The shift depends on the type of lensed image: type I induces no phase shift, type II induces a +π/2 phase shift, and type III induces a +π phase shift. In electromagnetic observations, type III images are typically suppressed and rarely seen (save for some rare exceptions; Dahle et al. 2013; Collett et al. 2017) because they occur near the maximum of the lens potential, where there is usually a lens galaxy or galaxy cluster; hence, it is hard to pick out faint images on top of a luminous lens. In GW observations, this type of image is expected to be demagnified; therefore, one would typically expect lensed GWs to consist of type I or type II images.

During the second observing run of Advanced LIGO and Advanced Virgo, seven BBHs (Abbott et al. 2017a, 2017b, 2017c, 2019) and one binary neutron star (Abbott et al. 2017d) were detected. A search for GW lensing signatures on the GWTC-1 catalog (Abbott et al. 2019) was performed in Hannuksela et al. (2019), but no good evidence of strong lensing was found. Note that the highest Bayes factor event pair in the analysis was the GW170104–GW170814 pair, but this was disfavored due to (1) the large time delay between the events and (2) the prior probability of lensing being low, around ∼10⁻⁵ yr⁻¹ at O1 sensitivity (Smith et al. 2018) for galaxy cluster lensing, and a relative lensing rate of ≲10⁻² for galaxy lensing (Ng et al. 2018; Hannuksela et al. 2019). If one or more of the observed images is of type III, the rate is understood to be significantly lower.

The same event pair was studied in more detail in Dai et al. (2020), which appeared at the time of writing of this paper. A third, subthreshold image consistent with the lensing hypothesis was found, GWC 170620, which was first discovered in the PyCBC subthreshold search (Nitz et al. 2019). The authors further analyzed the image configurations required for the lensing hypothesis, finding that it would consist of either one or two type III images and require a galaxy cluster lens due to the long time delay. Neglecting the a priori probability of lensing, the false-alarm probability for the double (triplet) was estimated at ∼10⁻⁴–10⁻², depending on specific O2 GW events (Dai et al. 2020). However, when accounting for the prior probability of lensing and the fact that the observed images would require a very peculiar image configuration, the lensing hypothesis is disfavored; the authors concluded that there is not sufficient evidence to conclude that the event pair is lensed. The authors excluded the mass and other binary parameters that rely on the knowledge of the BBH population parameters to determine the false-alarm rate of the pair due to lensing.

If the double/triplet events were lensed, then it would likely imply that the existing estimates of the lensing statistics are likely incorrect in predicting the relative fraction of galaxy cluster lenses and the total rate of lensed events (and hence the merger rate density of BBHs at high redshift). Another likely implication is a population of lenses that can form type III images more frequently than observed in the electromagnetic spectrum. To reconcile the discrepancy, one would likely require all of the following: (1) the merger rate density of BBHs rising at a higher rate than existing estimates from the usual formation channels, (2) galaxy cluster lenses making up a significant portion of the lensing optical depth, and (3) prominence of lensing configurations that can form heavily magnified type III images in GW channels but not electromagnetic channels.

In this paper, we present a Bayesian model selection method similar to Haris et al. (2018), but instead of computing the lensing model evidence using kernel density estimation from independent event posterior samples, we calculate the lensing evidence directly with parameter estimation by jointly fitting both images. We explicitly test the expected phase shifts and use an astrophysically motivated prior for the relative magnification factor. Moreover, we calculate the Bayes factors between the lensed and unlensed hypotheses, and, from the measured time delays, we determine the prior odds for any two events to be likely images of each other to produce an odds ratio that we can use to test the lensed and unlensed hypotheses.

For a GW signal at luminosity distance D_L , the amplitudes of the corresponding lensed images are magnified by a factor of $\sqrt{{\mu }_{i}}$, where i labels the corresponding absolute magnification factor for each image so that the observed luminosity distances will be

$$\begin{eqnarray}&&{D}_{L}^{(i)}={D}_{L}/\sqrt{{\mu }_{i}}.\end{eqnarray} \tag{ 1 }$$

Since the magnification factors and luminosity distances are degenerate, the individual magnification factors are difficult to constrain. Thus, we instead use the relative magnification factor μ,

$$\begin{eqnarray}&&\mu ={\left(\displaystyle \frac{{D}_{L}^{(1)}}{{D}_{L}^{(2)}}\right)}^{2}=\displaystyle \frac{{\mu }_{2}}{{\mu }_{1}},\end{eqnarray} \tag{ 2 }$$

where ${D}_{L}^{(1)}$ and ${D}_{L}^{(2)}$ are the observed luminosity distances of the first and second images, respectively, and μ₁ and μ₂ are the corresponding absolute magnification factors.

For strong lensing, the probability distribution for the individual magnifications is well known in the high-magnification limit and given by $p({\mu }_{i})\propto {\mu }_{i}^{-3}$ (Blandford & Narayan 1986). To estimate the prior distribution of the relative magnification factor μ, we follow the simulation in (Haris et al. 2018) and obtain p(μ) ∝ μ for μ < 1 and p(μ) ∝ μ⁻³ for μ ≥ 1.

Given two observed detector strains, d₁(t) and d₂(t), with confirmed GW detections, we want to determine whether these two signals are lensed or not. The lensed hypothesis ${{ \mathcal H }}_{L}$ states that the two signals come from the same astrophysical source and are thus lensed. Meanwhile, the unlensed model ${{ \mathcal H }}_{U}$ assumes that the two signals are from independent astrophysical sources. Under the lensing hypothesis, we first introduce a set of common parameters for the two events, η = {m₁, m₂, a₁, a₂, ι, α, δ, ψ}, where m₁ and m₂ are the component source frame masses, a₁ and a₂ are the component spins, ι is the inclination angle of the binary, α and δ are the R.A. and decl., and ψ is the polarization angle. We also introduce lensing-dependent parameters, ζ = {D_L , ϕ_c , t_c }, where D_L is the luminosity distance to the source, ϕ_c is the coalescence phase, and t_c is the time at coalescence. Hence, for the lensed hypothesis, we expect the common parameters η to be the same for the two events and only the lensing-dependent parameters to differ. Thus, the likelihood under the lensed hypothesis, given GW strain data d₁ and d₂, can be written as

$$\begin{eqnarray}&&P({d}_{1},{d}_{2}| {{\boldsymbol{\theta }}}_{1},{{ \mathcal H }}_{L})=P({d}_{1}| \eta ,{\zeta }_{1},{{ \mathcal H }}_{L})P({d}_{2}| \eta ,{\zeta }_{2},{{ \mathcal H }}_{L}),\end{eqnarray} \tag{ 3 }$$

where d₁, ζ₁, and d₂, ζ₂ are the data and independent parameters for the first and second images, respectively. Note that for this model, we sample the magnification factor μ instead of ${D}_{L}^{(2)}$ in the independent parameters ζ₂.

For the unlensed hypothesis, the parameters of the two events are sampled independently. The likelihood in the unlensed hypothesis ${{ \mathcal H }}_{U}$ is simply the product of the likelihoods of the two events, since they are independent of each other,

$$\begin{eqnarray}&&P({d}_{1},{d}_{2}| {{\boldsymbol{\theta }}}_{2},{{ \mathcal H }}_{U})=P({d}_{1}| {\eta }_{1},{\zeta }_{1},{{ \mathcal H }}_{U})P({d}_{2}| {\eta }_{2},{\zeta }_{2},{{ \mathcal H }}_{U}),\end{eqnarray} \tag{ 4 }$$

where η₁, ζ₁, and η₂, ζ₂ are the parameters for the first and second GW events, respectively.

To compare the two models, we compute the ratio of the evidence $P({d}_{1},{d}_{2}| {n}_{j},{{ \mathcal H }}_{L})$ and $P({d}_{1},{d}_{2}| {{ \mathcal H }}_{U})$, also known as the Bayes factor,

$$\begin{eqnarray}\begin{array}{rcl}{B}_{U}^{L} & = & \displaystyle \frac{P({d}_{1},{d}_{2}| {n}_{j},{{ \mathcal H }}_{L})}{P({d}_{1},{d}_{2}| {{ \mathcal H }}_{U})}\\ & = & \displaystyle \frac{\displaystyle \int P({d}_{1}| \eta ,{\zeta }_{1},{{ \mathcal H }}_{L})P({d}_{2}| \eta ,{\zeta }_{2},{{ \mathcal H }}_{L})P(\eta ,{\zeta }_{1},{\zeta }_{2}| {{ \mathcal H }}_{L},{n}_{j})d\eta d{\zeta }_{1}d{\zeta }_{2}}{\displaystyle \int P({d}_{1}| {\eta }_{1},{\zeta }_{1},{{ \mathcal H }}_{U})P({d}_{2}| {\eta }_{2},{\zeta }_{2},{{ \mathcal H }}_{U})P({\eta }_{1},{\zeta }_{1},{\eta }_{2},{\zeta }_{2}| {{ \mathcal H }}_{U})d{\eta }_{1}d{\zeta }_{1}d{\eta }_{2}d{\zeta }_{2}},\end{array}\end{eqnarray} \tag{ 5 }$$

where n_j is the Morse index, which determines the type of lensing image (type I/II/III) and thus the expected phase difference for the pair, while $P(\eta ,{\zeta }_{1},{\zeta }_{2}| {{ \mathcal H }}_{L},{n}_{j})$ and $P({\eta }_{1},{\zeta }_{1},{\eta }_{2},{\zeta }_{2}| {{ \mathcal H }}_{U})$ are the priors for the lensed and unlensed hypothesis, respectively. We calculate the lensing evidence with nested sampling (Skilling 2006; Veitch et al. 2015) using lalinference_nest (LIGO Scientific Collaboration 2018).

The time delay between any two events can also be used to compute a corresponding timing Bayes factor (Haris et al. 2018; Hannuksela et al. 2019),

$$\begin{eqnarray}&&{B}_{t}=\displaystyle \frac{P({\rm{\Delta }}t| {{ \mathcal H }}_{L})}{P({\rm{\Delta }}t| {{ \mathcal H }}_{U})}.\end{eqnarray} \tag{ 6 }$$

We estimate the probability distribution $P({\rm{\Delta }}t| {{ \mathcal H }}_{L})$ through simulation following the methodology of Haris et al. (2018). We compute $P({\rm{\Delta }}t| {{ \mathcal H }}_{U})$ by assuming that independent (unlensed) events are Poisson distributed.

To obtain the odds ratio for the lensed and unlensed hypotheses, we compute

$$\begin{eqnarray}&&{O}_{U}^{L}=\displaystyle \frac{P({d}_{1},{d}_{2}| {n}_{j},{{ \mathcal H }}_{L})}{P({d}_{1},{d}_{2}| {{ \mathcal H }}_{U})}\displaystyle \frac{P({\rm{\Delta }}t| {{ \mathcal H }}_{L})}{P({\rm{\Delta }}t| {{ \mathcal H }}_{U})}\displaystyle \frac{P({{ \mathcal H }}_{L})}{P({{ \mathcal H }}_{U})},\end{eqnarray} \tag{ 7 }$$

where the ratio $P({{ \mathcal H }}_{L})/P({{ \mathcal H }}_{U})$ is the prior odds for lensing compared to the unlensed event model. However, type III images are very rare; hence, it is likely that the image shift corresponding to Δϕ = ± π is heavily disfavored: $p({n}_{j}=1| {{ \mathcal H }}_{L})\ll p({n}_{j}=\{0,1/2\}| {{ \mathcal H }}_{L})$.

The prior odds reflect our belief in the probability of lensing for any two events and can be estimated through simulations, as well as from electromagnetic observations. As such, we compute this via the ratio of expected lensed event to independent event rate. The relative lensed event rate has been estimated for galaxy lenses to be around $p({{ \mathcal H }}_{L})/p({{ \mathcal H }}_{U})\sim {10}^{-3}$ at design sensitivity (Li et al. 2018; Oguri 2018), while Smith et al. (2018) found the relative rate of galaxy cluster lensing to be $p({{ \mathcal H }}_{L})\sim {10}^{-5}\,{\mathrm{yr}}^{-1}$ at O1 sensitivity; we expect this to be somewhat larger at O2 sensitivity.

When comparing two models, the Bayesian evidence penalizes a more complex model. If we compare two models, the one with smaller prior volume or fewer parameters would be favored. This penalty is known as the Occam factor (Thrane & Talbot 2019), which is automatically achieved by Bayesian inference. In our analysis, the lensing model has fewer parameters due to the parameter sharing. Indeed, when a signal is consistent with the lensed hypothesis, the magnitude of the Bayes factor is entirely set by the prior volume; a larger prior can increase the Bayes factor by several orders of magnitude, and vice versa. In order to reduce the prior volume difference between the two models, we impose a uniform in $\mathrm{log}({m}_{1})$ and $\mathrm{log}({m}_{2})$ priors, instead of the typical uniform priors on m₁ and m₂ (both in the detector frame). We impose a prior within the mass range of 1–100 M_⊙. The difference between the prior volumes can be reduced by a factor of 10²–10³ when using the uniform in log space prior. Also, the posteriors are almost unaffected by the change from one prior to another because masses can be well constrained from data. In other words, this choice of prior has a much greater effect on the prior volume than the posteriors. We note that an astrophysically motivated mass prior could be used instead, such as the power-law model used by the LVC (Abbott et al. 2019); however, for the reasons stated above, we decided to use the uniform in log prior instead.

We also take care of selection effects in the joint parameter estimation, since GW detectors are not sensitive enough to detect all the binaries in the prior parameter space. Thus, we incorporate a selection function in the parameter estimation directly that keeps a sample if it is above the detection S/N threshold and otherwise rejects it. We also note that taking into account selection effects is important when one or both events are below the detection threshold (Li et al. 2019; McIsaac et al. 2020), as will be the case when performing the joint parameter estimation using the potential third image, GWC 170620 (Dai et al. 2020).

We analyze potential pairs of lensed events from the second observing run of Advanced LIGO and Virgo. Due to the high computational cost of the parameter estimation, we select pairs of events that have similar sky localizations. We then run lalinference_nest (LIGO Scientific Collaboration 2018) to obtain the lensed and unlensed model evidence. We apply the selection function implemented in lalinference_nest to the parameter estimation and set the network S/N threshold to 10 for single events and 14 for joint events, except for the subthreshold event GWC 170620. We sample uniformly in $\mathrm{log}(m1)$ and $\mathrm{log}(m2)$ for both lensed and unlensed models in order to mitigate the prior volume difference between models. The waveform used in our analysis is the IMRPhenomD approximant (Husa et al. 2016; Khan et al. 2016), a nonprecessing and spin-aligned (22-mode only) frequency domain BBH waveform that enables us to test the different coalescence phase shifts due to different image types. Since there is no evidence that precession has been observed in any of the events detected in O1 and O2 (Abbott et al. 2019), we expect that the IMRPhenomD model is sufficiently accurate for this analysis.

In Figure 1, we show the Bayes factors B^L _U and B_t for pairs of events in O2. The lensed models for each pair of events are evaluated with the four different possible coalescence phase shifts. For the unlensed model, we sample the phases for each event independently. The GW170104–GW170814 pair with a π/2 coalescence phase shift has the largest Bayes factor, ${B}_{U}^{L}\sim 1.98\times {10}^{4}$, favoring the lensed hypothesis in the absence of prior probability. The results in which the prior probability is taken into account are shown in Figure 2, which will be discussed later. The event could still be an unlensed event, however, as the two events could be from independent sources that have similar parameters. Therefore, we note that a high Bayes factor is not necessarily indicative of lensing. Nevertheless, it is intriguing that the event favors the lensing hypothesis even when including all of the binary parameters. The GW170809–GW170814 pair, which was suggested as a lensed event by Broadhurst et al. (2019), is clearly disfavored by the model selection.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The sky localization posterior for the GW170104–GW170814 pair is shown in Figure 3. The sky localization posterior inferred under the lensing hypothesis (joint parameter estimation) is better constrained and lies within the overlap region of the GW170104 and GW170814 independent parameter estimation runs. The 90% confidence region is better constrained because the joint run has a higher S/N than each individual run and benefits from the "extra detectors" (more baselines for localization) due to the different times of arrival of each image with respect to the rotation of the Earth. We also show the posterior distributions over the parameters that we expect to be unchanged due to lensing in Figure 4. Similarly, the parameters inferred under the lensing hypothesis are better constrained, compared to those inferred independently.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Assuming an observation time of 9 months for O2 and the simulations in Haris et al. (2018), we obtain an estimate for P(Δt∣L). For the unlensed case, we assume that the detected event rate follows a Poisson distribution, that is, P(Δt∣U) =2(T − Δt)/T², where T is the observation time. In Figure 1, we show the Bayes factor B_t of the lensing model with four different phase shifts compared to the unlensed model. Due to the ∼7 month time delay between GW170104 and GW170814, the Bayes factor B_t is ∼ 8.7 × 10⁻². The Bayes factor B_t can be large if the time delay is only a few days, such as for the GW170809–GW170814 and GW170818–GW170823 pairs, which results in B_t ∼ 6.

The two Bayes factors B^L _U and B_t , together with the prior odds for lensing, can be combined to compute the odds ratio (Equation (7)). The lensed event rate is estimated to be ∼0.1 yr⁻¹ for O2 (Ng et al. 2018; Hannuksela et al. 2019). During O2, there were eight events detected in 9 months of observing time; therefore, the prior odds are ∼0.009. In Figure 2, we show the odds ratio for events detected in O2. The odds ratio for GW170104 and GW170814 with a π/2 coalescence phase shift is ${O}_{U}^{L}\sim 20$. It is the only pair of events that moderately prefers the lensing hypothesis, even with prior information folded in. For other pairs detected in O2, we do not see any lensed evidence, as the odds ratios are much less than 1.

Our B_t and prior odds estimations are based on galaxy lensing. However, the long time delay may point to lensing by a galaxy cluster. Galaxy cluster lensing is expected to be rare at O2 sensitivity. The prior probability of galaxy cluster lensing is ∼10⁻⁵ yr⁻¹ (Smith et al. 2018). Meanwhile, because the phase shift corresponds to π or a type III image as pointed out in Dai et al. (2020), the probability of lensing is further disfavored; the probability of observing type III images should be very low, such that $p({n}_{j}| {{ \mathcal H }}_{L})\ll 1$. Indeed, type III images are rarely observed in the electromagnetic band.

To summarize, we used Bayesian model selection to identify lensed events in the second observing run of Advanced LIGO and Virgo. We test the lensing model with 0, π/4, π/2, and 3π/4 coalescence phase shifts and take care of selection effects. The most significant event pair is found to be GW170104–GW170814 at a very high Bayes factor; however, the two signals may come from two independent sources and have similar parameters, and therefore it is not clear if this is indicative of lensing. Moreover, the event is disfavored as a lensed candidate based on the understanding of lens configurations and the BBH and lens populations by an overwhelming amount. Indeed, the probability of observing lensing configurations with these types of time delays and image configurations is low for both galaxies and galaxy clusters. The estimates of lensed rates and the relative contribution of galaxies and galaxy clusters may vary to some degree, but no current estimate predicts that galaxy cluster lensing should become prominent at O2 sensitivity.

Let us then entertain the possibility that the event was lensed. If this were the case, then unless we were to accept that we were simply incredibly lucky, it would likely imply all of the following.

• 1.  
  The relative contribution of galaxy cluster lensing is more important than previously believed, which would explain the high time delay between the events. Smith et al. (2018, 2017, 2019) and Robertson et al. (2020) studied galaxy cluster lensing and argued that highly magnified events have historically been observed more prominently lensed by galaxy cluster scale lenses.
• 2.  
  The merger rate density of BBHs must rise at a significantly higher rate than previously predicted. Indeed, the current lensing rate estimates rely on black holes tracing the star formation rate density. For example, the Belczynski distribution is often used to model the merger rate density (Belczynski et al. 2016). Out of the studied models of black hole formation, none predict high enough merger rates that galaxy cluster lensing would become observable at O2 sensitivity (Oguri 2018).
• 3.  
  Type III images are more prominent for GW sources than they are for electromagnetic sources.

Let us therefore state that extraordinary claims require extraordinary evidence. Based on the prior probability of lensing by these types of systems, we advise the reader to be very careful in interpreting the results. Indeed, in the absence of clear-cut evidence to the contrary, we must conclude that there is insufficient evidence to claim that the event is lensed, in agreement with Dai et al. (2020) and Hannuksela et al. (2019).

However, two pieces of evidence could together possibly determine if the events were lensed. First, it is vital to perform injection campaigns to determine the probability of a nonlensed event. A similar study was conducted in Dai et al. (2020) with false-alarm probability between 10⁻⁴ and 10⁻² for O2 events. Moreover, another intriguing possibility for cross-verification is through searches in the electromagnetic channels, as pointed out by Dai et al. (2020); if the events are lensed, then their host galaxy must also be lensed (Hannuksela & Collett 2020). If the third event proposed as a lensed candidate in Dai et al. (2020) and the pair are lensed images, then the triplets would likely afford us three time delays, which would allow for a unique opportunity to localize the host galaxy and the galaxy cluster that lensed it in an electromagnetic follow-up. Hannuksela & Collett (2020) and Yu et al. (2020) demonstrated that such a search is possible for galaxies. We note that due to the rarity of galaxy clusters, the search is expected to be even more powerful for galaxy clusters. In the case of doubly lensed events, such as the GW170104–GW170814 pair, the single time-delay estimate may be quite degenerate with the lens parameters and source alignment.

The authors would like to thank Otto Hannuksela, K. Haris, John Veitch, Rico K. L. Lo, Ajit Mehta, Jose Ezquiaga, Daniel Holz, David Keitel, Suvodip Mukherjee, Ken K. Y. Ng, Leo Singer, and Will Farr for useful discussions. I.M.H. is supported by the NSF Graduate Research Fellowship Program under grant DGE-17247915. This work was supported by NSF awards PHY-1607585 and PHY-1912649. The authors are grateful for computational resources provided by the LIGO Laboratory and supported by National Science Foundation grants PHY-0757058 and PHY-0823459 and those provided by the Leonard E. Parker Center for Gravitation, Cosmology and Astrophysics at the University of Wisconsin–Milwaukee. We thank the LIGO and Virgo Collaboration for providing the data for this work. This article has been assigned LIGO document number LIGO-P1900017.