A Detailed Analysis of GW190521 with Phenomenological Waveform Models

We use the same notation and conventions as in our reanalysis of GWTC-1 (Mateu-Lucena et al. 2021). We will report all masses in units of the solar mass M_⊙, and except where otherwise noted, we always refer to masses in the source frame, assuming a standard cosmology (Planck Collaboration et al. 2016; see Appendix B of Abbott et al. 2019). Some figures and tables use an "src" index as a more explicit notation for clarity. Individual component masses are denoted by m_i , and the total mass is M = m₁ + m₂. The chirp mass is ${ \mathcal M }={({m}_{1}\,{m}_{2})}^{3/5}{M}^{-1/5}$. We define mass ratios q = m₂/m₁ ≤ 1 and Q = m₁/m₂ ≥ 1.

We also define effective spin parameters, which are commonly used in waveform modeling and parameter estimation. The parameter χ_eff is defined as

$$\begin{eqnarray}&&{\chi }_{\mathrm{eff}}=\displaystyle \frac{{m}_{1}{\chi }_{1}+{m}_{2}{\chi }_{2}}{{m}_{1}+{m}_{2}},\end{eqnarray} \tag{ 1 }$$

where the χ_i are the projection of the spin vectors of the individual black holes onto the instantaneous direction perpendicular to the orbital plane. The effective spin-precession parameter χ_p (Schmidt et al. 2015) is designed to capture the dominant effect of precession, and corresponds to an approximate average over many precession cycles of the spin in the precessing orbital plane, and is defined in terms of the average spin magnitude S_p (Schmidt et al. 2015),

$$\begin{eqnarray}&&{S}_{{\rm{p}}}=\displaystyle \frac{1}{2}\left({A}_{1}{S}_{1,\perp }+{A}_{2}{S}_{2,\perp }+| {A}_{1}{S}_{1,\perp }-{A}_{2}{S}_{2,\perp }| \right),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&=\ \max \left({A}_{1}{S}_{1,\perp },{A}_{2}{S}_{2,\perp }\right),\end{eqnarray} \tag{ 3 }$$

where A₁ = 2 + 3/(2q), and χ_p is then defined as

$$\begin{eqnarray}&&{\chi }_{{\rm{p}}}=\displaystyle \frac{{S}_{{\rm{p}}}}{{A}_{1}{m}_{1}^{2}}.\end{eqnarray} \tag{ 4 }$$

Both χ_eff and χ_p are dimensionless and thus independent of the frame (source or detector).

We will employ waveforms with several multipoles beyond the quadrupolar contribution, always considering pairs of both positive and negative modes when referring to a particular multipole. Thereby, to refer to the example list of multipoles (l, m) = (2, ±2), (2, ±1), we will use the notation (l, ∣m∣) = (2, 2), (2, 1) or simply (2, 2), (2, 1).

For a complete list of all of the waveform models used in the present paper, see Table 1. The original LVC publications on GW190512 (Abbott et al. 2020b, 2020c) used three families of waveform models, which represent incarnations of three well established approaches to compact binary coalescence (CBC) waveform modeling:

• 1.  
  The SEOBNRv4PHM time-domain model (Ossokine et al. 2020; Babak et al. 2017; Pan et al. 2014): constructed within the effective-one-body (EOB) framework (Akcay et al. 2019; Bohé et al. 2017; Cotesta et al. 2020; Damour & Nagar 2016; Damour et al. 2015, 2013; Nagar et al. 2019, 2018; Nagar & Rettegno 2019; Nagar et al. 2017; Del Pozzo & Nagar 2017; Nagar et al. 2016; Pürrer 2016; Taracchini et al. 2012, 2014; Damour 2001; Cotesta et al. 2018).
• 2.  
  The NRSur7dq4 surrogate time-domain model (Varma et al. 2019), which directly interpolates NR waveforms.
• 3.  
  The IMRPhenomPv3HM frequency-domain model, which corresponds to the third generation of models in the IMRPhenom family (Husa et al. 2016; Khan et al. 2016; Hannam et al. 2014; Bohé et al. 2016; London et al. 2018; Khan et al. 2019, 2020).

Here we employ two further recently developed models that represent upgrades of IMRPhenomPv3HM and constitute a fourth generation of IMRPhenom models:

• 1.  
  The frequency-domain model IMRPhenomXPHM (Pratten et al. 2021), which builds upon an underlying non-precessing model IMRPhenomXHM (Pratten et al. 2020; García-Quirós et al. 2020, 2020b) that features calibration of the subdominant harmonics to NR simulations.
• 2.  
  IMRPhenomTPHM (Estellés et al. 2020b, 2021c), building on the non-precessing model IMRPhenomTHM (Estellés et al. 2020b, 2020c), which essentially applies the same phenomenological techniques at the heart of IMRPhenomXPHM to construct a native time-domain model. Working in the time domain allows several key improvements that we will discuss below.

All of these models include a description of precession effects and subdominant harmonics, but do not include eccentricity, and they have complementary strengths and shortcomings that we will detail below.

Only NRSur7dq4 is calibrated to precessing NR waveforms, but its training data set is restricted to mass ratio Q ≤ 4 and dimensionless component spin magnitudes a_1,2 ≤ 0.8). However the model can also be evaluated in the extrapolation region with Q ≤ 6 and a_1,2 ≤ 1. Furthermore, usage of the model is restricted by the length of the original time-domain NR waveforms, and restrictions get tighter in the frequency domain due to the need of windowing before Fourier transforming the template. The limited length of NRSur7dq4 waveforms leads, in particular, to extra constraints on the minimum frequency and total mass allowed in parameter estimation analyses.

The IMRPhenom models describe precession via an approximate map between signals from non-precessing and precessing systems, which we will refer to as "twisting-up" (Schmidt et al. 2011, 2012; Hannam 2014). This approximation exploits the fact that, at least during the inspiral, the precession timescale is much slower than the orbital timescale, and thus the precessing motion mainly acts as an amplitude modulation. The spin of the remnant of the precessing system is however in general significantly different from the final spin of the non-precessing system. For a recent discussion of the approximations of this approach, see Ramos-Buades et al. (2020).

The SEOBNRv4PHM model (Buonanno et al. 2003; Babak et al. 2017; Pan et al. 2014; Cotesta et al. 2018; Ossokine et al. 2020) numerically integrates in the time domain the EOB BBH dynamics, including the spin-precession equations, using a Hamiltonian and GW flux that are tuned to non-precessing NR simulations. Then, the waveforms in the inertial frame are obtained by applying a time-domain rotation ("twisting-up") to the waveforms in the co-precessing frame (Buonanno et al. 2003; Babak et al. 2017; Pan et al. 2014).

Neither the SEOBNRv4PHM nor IMRPhenom models include the asymmetries between positive and negative m spherical harmonic modes in the co-precessing frame, which are related to the large recoil velocities observed in some NR simulations of precessing binaries, as discussed in Bruegmann et al. (2008).

We now turn to describe relevant aspects of the IMRPhenomX and IMRPhenomT families, which we use for the new results presented in this paper. Frequency-domain waveform models are particularly attractive for GW data analysis, since they are naturally adapted to a matched-filter-type analysis, where the noise is characterized in the frequency domain, and accordingly allows the most computationally efficient Bayesian inference analysis. In order to accelerate the evaluation of precessing waveforms, current frequency-domain IMRPhenom models also use the SPA approximation to compute the transfer functions between the frequency-domain non-precessing waveform and the precessing waveform in an inertial frame (see Marsat & Baker 2018 for more accurate alternatives). The assumptions underlying the SPA fail for merger and ringdown, but the method has been found to work surprisingly well, and has been routinely used in GW data analysis; see, e.g., Abbott et al. (2019) and Abbott et al. (2021a). The approximation does however have to be employed with caution when essentially the only observable part of the signal is the merger and ringdown, as is the case for GW190521. This shortcoming has been one of the main reasons for us to also develop a time-domain phenomenological waveform model, IMRPhenomTPHM, which does not rely on the SPA. One of the goals of this paper is to discuss in detail how to avoid misleading conclusions from analyzing GW190521 with IMRPhenomXPHM, and how improved results can be obtained with IMRPhenomTPHM.

Since neither IMRPhenomXPHM nor IMRPhenomTPHM are calibrated to precessing NR waveforms but rather build on the above approximations to describe precession effects, it is essential to incorporate in the models some functionality to test the robustness of results for challenging events like GW190521. As discussed in Appendix F of Pratten et al. (2021) for IMRPhenomXPHM and in the Appendix of Estellés et al. (2021a) for IMRPhenomTPHM, the LALSuite (LIGO Scientific Collaboration 2020) implementation of our models supports several options regarding the choice of precession prescription and final-spin approximations. These options are selected with parameters that take integer values, which we will refer to as PV for the inspiral precession version and FS for final spin.

Our "twisting-up" procedure is based on time/frequency dependent rotations from the co-precessing frame to an inertial one in which we observe the signal. For the IMRPhenom models, this rotation is implemented through three Euler angles. IMRPhenomPv2 only supports an effective single-spin, orbital-averaged description valid at next-to-next-to-leading (NNLO) post-Newtonian order. IMRPhenomXPHM allows us to use the same prescription, but as a default, it relies on a more recent double-spin description that can be derived within the post-Newtonian framework using multiple-scale-analysis (MSA; Chatziioannou et al. 2017; this description is also used by IMRPhenomPv3HM). IMRPhenomTPHM implements both NNLO and MSA Euler angles, but its default behavior is to numerically integrate evolution equations for the component spins as discussed in Estellés et al. (2021a). We will refer to different precession prescriptions with the acronym PV and follow the same convention enforced in LALSuite, where PV = 223 corresponds to the MSA approximation and PV = 300 to the numerically integrated angles.

Another setting of the models that can be specified by the user is the final-spin approximation, as discussed in Section IV.D of Pratten et al. (2021) for IMRPhenomXPHM and Section II.E of Estellés et al. (2021a) for IMRPhenomTPHM. The default choice for IMRPhenomXPHM is to use a precession-averaged equation inspired by the MSA formalism. This version will be referred to as version FS = 3. Alternative versions attach the in-plane spins to the larger mass, either relying on the usual effective precession spin χ_p (FS = 0, which is adopted by all third-generation IMRPhenom models), or by taking the norm of the in-plane spin vectors at the reference frequency (FS = 2). The default version of IMRPhenomTPHM takes the norm of the in-plane spin vectors at the coalescence time (FS = 4).

There are several improvements in the treatment of precession achieved by the time-domain IMRPhenomTPHM in comparison with the frequency-domain IMRPhenomXPHM. First, in order to obtain explicit expressions for the spherical harmonic modes of the precessing frequency-domain models, IMRPhenomXPHM and previous IMRPhenom models use the SPA to compute approximate Fourier transforms. Second, in the time domain, it is simple to incorporate analytical knowledge about the ringdown frequencies in the ringdown portion of a precessing waveform; see Section II.E of Estellés et al. (2021a). This has not yet been achieved in the frequency domain. This is particularly crucial for GW190521, where a large part of the observed signal-to-noise ratio (S/N) is due to the ringdown portion of the signal. Third, the numerical integration of the equations for the spin dynamics in IMRPhenomTPHM also resolves an inconsistency of the MSA Euler angles with the non-precessing limit, which we discuss in Estellés et al. (2021a). The numerical integration also provides further gains in accuracy, and we find it to decrease the computational cost (Estellés et al. 2021a).

Finally, we note that a particularly challenging region of the BBH parameter space arises at larger mass ratios, where the precession cone angle can become large, and the orbital angular momentum L can become smaller than the (sum of the) component spins. Then both L and the total angular momentum J may end up flipping their direction. The latter situation is also known as transitional precession (Apostolatos et al. 1994; Kidder 1995; Schmidt et al. 2012), as opposed to simple precession, when J at least approximately maintains its direction. Very few NR simulations exist for the cases of large angles between J and L, and for the sign of J flipping, and these situations are related to various caveats in the post-Newtonian and MSA approximations that are often used in waveform modeling, and in particular in the IMRPhenom models. While the systematic errors of precessing waveform models are in general not yet very well understood, this is particularly true when the normal to the orbital plane or the final-spin flip sign with respect to the direction at large separation. This situation indeed arises for the results of Nitz & Capano (2021). In such cases, one should proceed with great caution, and test the robustness of results by comparing different waveform models. Nitz & Capano (2021) compared their IMRPhenomXPHM results with the NRSur7dq4 model, but this analysis was limited by the parameter space coverage of NRSur7dq4 (Varma et al. 2019). In contrast, IMRPhenomTPHM will allow us a more global comparison in this work. In addition, future work will aim to improve the robustness of our models for such situations.

We use public GW strain data collected by the Advanced LIGO detectors (Aasi et al. 2015) and Advanced Virgo detector (Acernese et al. 2015) from the Gravitational Wave Open Science Center (GWOSC; Vallisneri et al. 2015; Abbott et al. 2021b), as well as PSDs and calibration uncertainties included in the GWOSC release (LIGO Scientific Collaboration & Virgo collaboration 2020a). From the available GWOSC strain data sets we have selected the data sampled at 16 kHz, with a sampling rate of 1024 Hz chosen for our analysis, consistent with the choice in Abbott et al. (2020b) and Abbott et al. (2020c). The lower and upper cutoff frequencies for the likelihood integration were taken to be 11.0 Hz and 512 Hz (the Nyquist frequency corresponding to the sampling rate), again consistent with Abbott et al. (2020b) and Abbott et al. (2020c).

We have carried out Bayesian parameter estimation of the signal using two publicly available codes, the Python-based parallel Bilby (pBilby, PB) code (Ashton et al. 2019; Smith et al. 2020), which uses the dynesty (Speagle 2020) variant of the nested sampling algorithm (Skilling 2004), and the LALInference (LI) code (Veitch et al. 2015), which is part of the LALSuite (LIGO Scientific Collaboration 2020) package for GW data analysis, using its implementation of Markov Chain Monte Carlo (MCMC) sampling.

Parallel Bilby provides a highly parallel and flexible implementation of nested sampling, and supports a range of priors and choices of sampling parameters and settings. With pBilby, we sample in mass ratio and chirp mass, which is easier than sampling the component masses in that code. We largely use the default settings of the code apart from the following choices: we fix the minimal (walks) and maximal (maxmcmc) number of MCMC steps to 200 and 15,000, respectively. For our final results, we have set the number of autocorrelation times to use before accepting a point (nact) to a value of 30. We have varied the number of nested sampling live points (nlive) between 1024 and 4096 for selected runs to test that we have obtained (sufficiently) converged results, and always show the results for nlive = 4096. In order to speed up calculations, we use distance marginalization as described in Thrane & Talbot (2019). For each of the pBilby runs we quote results for, we have carried out four independent simulations (independent seeds), and then merged the posteriors to a single posterior with PESummary (Hoy & Raymond 2021).

LALInference samples in mass ratio and chirp mass, re-weighting to a prior that is flat in component masses as described in Veitch et al. (2015). We use essentially standard LALInference settings with eight temperatures, but a large number of independent chains, 120 for our production runs. For our LALInference runs, we do not employ the distance marginalization used for our Bilby runs.

We have previously used pBilby as our primary code for our reanalysis of the GW190412 event (Colleoni et al. 2021; Estellés et al. 2021a; Mateu-Lucena et al. 2021), where we found good agreement with LALInference results as reported in Colleoni et al. (2021). We have however found that the computational cost of comparably well sampled pBilby runs is significantly higher than for LALInference runs due to the high required settings of the nact parameter. Here we use LALInference for our primary results and pBilby for comparisons.

Runs performed with pBilby have been sampled with a prior uniform in "inverse" mass ratio Q = 1/q, following Nitz & Capano (2021). This prior emphasizes unequal masses and improves pBilby convergence in the unequal-mass regime. Unlike Nitz & Capano (2021), we choose to sample in the detector frame, to take advantage of distance marginalization, which would require a nonstandard likelihood in Bilby if sampling in source frame. In most of the results shown, and as indicated when stating results, we have performed a post-processing re-weighting from this prior to a prior flat in component masses with the corresponding functions from the Bilby code, to obtain results matching the same prior as the LALInference runs (flat in component masses).

Additionally we have performed some runs with restricted versions of these priors to improve resolution for more unequal masses. In particular, we report in Section 5 results for LALInference runs in a mass-ratio range q ∈ [0.035, 0.15] and pBilby runs in a range q ∈ [0.035, 0.2]. For studying the possible association of the event with the reported active galactic nucleus (AGN) flare (Graham et al. 2020), we have also performed runs fixing the sky location and luminosity distance to the values reported for the AGN flare.

Finally, we have employed mainly three different sets of mass prior ranges for the runs reported in this work, checking that we avoid any significant railing against prior limits. Note that a small amount of railing against the lower-mass-ratio bound is still present; however, we decided to not attempt to precisely map out the tail at low q as models become less reliable and computational cost increases significantly. We comment on the possibility of further low-q posterior modes in Section 6. For IMRPhenomTPHM and IMRPhenomTHM runs with LALInference, we have employed a prior uniform in component masses (in detector frame) with a range of ${m}_{1,2}^{\det }\in [10,400]{M}_{\odot }$ and default mass-ratio constraints q ∈ [0.035, 1.0]. For IMRPhenomTPHM runs with pBilby, we have employed a prior uniform in inverse mass ratio with range q ∈ [0.035, 1.0] and uniform in total mass with range ${M}_{{\rm{T}}}^{\det }\in [80,550]{M}_{\odot }$, and constraints for component masses (in detector frame) ${m}_{1}^{\det }\in [30,400]{M}_{\odot }$ and ${m}_{2}^{\det }\in [5,400]{M}_{\odot }$. For IMRPhenomXHM and IMRPhenomXPHM runs, we have employed the same ranges in LALInference as with IMRPhenomTPHM runs, but for pBilby runs, the parameter ranges are: q ∈ [0.04, 1.0], ${m}_{1}^{\det }\in [30,300]{M}_{\odot }$, ${m}_{2}^{\det }\in [5,200]{M}_{\odot }$, and ${M}_{{\rm{T}}}^{\det }\in [80,550]{M}_{\odot }$. Differences in the mass ranges are due to problems found with railing posteriors for the IMRPhenomTPHM runs, which in general have support for higher component masses.

In some cases we also re-weight the posterior samples obtained with a certain set of priors to obtain an approximation of what the posterior should be when using a different set of priors, without having to run the full alternative inference. Specifically, we use the bilby implementation of this re-weighting procedure, which converts samples with a given prior in chirp mass and mass ratio to a prior flat in component masses by resampling the posterior with weights defined as the ratio in new-over-old prior values times the Jacobian of the transformation. The procedure generates new posteriors that contain only 25% of the original number of samples.

All runs have maximum component spin magnitudes limited to 0.99, and the luminosity distance prior is chosen as uniform in comoving volume, assuming the Planck15 cosmology (Planck Collaboration et al. 2016) with a range D_L ∈ [0.2, 10] Gpc. The LALInference prior contains an additional factor 1/(1 + z_L ), accounting for time dilation, which is not present in the definition for pBilby; but using the re-weighting procedure on pBilby results, we have tested that this does not have any noticeable influence on estimates of D_L or any other quantities.

The main results of Bayesian parameter estimation are the posterior distributions, and point estimates are usually given as medians with error estimates given by the 90% confidence intervals. Sometimes, it can also be enlightening to consider the maximum-likelihood value ($\max { \mathcal L }$) returned by an analysis, and which point in parameter space it correspond to, as the likelihood directly answers the question of how well the employed model (across the sampled part of parameter space) can fit the data. However, there are several caveats in interpreting the $\max { \mathcal L }$ values over a set of posterior samples, since a Bayesian parameter estimation run, such as those we employ here using pBilby and LALInference, is by construction not an optimal $\max { \mathcal L }$ finding algorithm. The prior has significant influence on how densely which parts of the parameter space are evaluated, and the $\max { \mathcal L }$ reported over the final posterior sample may be far from the actual maximum over all likelihoods evaluated while the sampling chains progressed. It is important to note that the number of posterior samples is typically much smaller than the number of likelihood evaluations, and to achieve a good estimate of $\max { \mathcal L }$, much more expensive sampling settings may be required than in order to get good estimates for source parameter values and error estimates. Nevertheless, comparing the $\max { \mathcal L }$ across runs can be a useful additional diagnostic of the behavior of waveforms and samplers, and we highlight them on all posterior plots in this paper.

Indeed for GW190521, we find that the maximum-likelihood parameters do not appear to be stable across runs, and are influenced by statistical fluctuations, sampler settings, as well as waveform models and priors. For an example, see Section 5.3.

Our main results derive from the posterior distributions we have obtained with the IMRPhenomXPHM and IMRPhenomTPHM models. To test the influence of the harmonic content in the templates, we will in general present results for different harmonic content for IMRPhenomTPHM. In Figures 2 and 3 we show these posteriors for some key parameters: the component masses, mass ratio, effective spins χ_eff and χ_p, luminosity distance d_L , and the angle θ_JN between J and the line of sight. We find consistency between the results obtained with LALInference and pBilby after re-weighting the pBilby results (with the uniform-in-Q prior) to a prior that is flat in component masses, giving us confidence in our sampling of parameter space. Recovered S/Ns and signal-versus-noise Bayes factors for our main runs and several different model options are shown in Table 2. Point estimates for key parameters of the main runs are also summarized in Table 3, again comparing with those from Nitz & Capano (2021), Abbott et al. (2020b), and Abbott et al. (2020c). Complete posterior data sets for our standard LALInference IMRPhenomTPHM run with ${{\ell }}_{\max }=4$ and standard IMRPhenomXPHM run can be found in our Zenodo data release (Estellés et al. 2021b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Network Matched-filter S/Ns with 90% Credible Intervals and log Signal-to-noise Bayes Factors ${ \mathcal B }{ \mathcal F }$ for Runs with Waveform Models in the IMRPhenomX and IMRPhenomT Families, Including Several Different Options of the IMRPhenomTPHM Model

Approx.	$\mathrm{log}{ \mathcal B }{ \mathcal F }$	${\rho }_{\mathrm{mf}}^{{\rm{H}}}$	${\rho }_{\mathrm{mf}}^{{\rm{L}}}$	${\rho }_{\mathrm{mf}}^{{\rm{V}}}$	${\rho }_{\mathrm{mf}}^{{\rm{N}}}$
XHM	80.06 ± 0.15	${8.0}_{-0.3}^{+0.2}$	${11.8}_{-0.3}^{+0.5}$	${2.4}_{-1.2}^{+0.7}$	${14.4}_{-0.3}^{+0.3}$
XPHM PV = 223 FS = 3	80.43 ± 0.21	${7.9}_{-0.3}^{+0.2}$	${11.8}_{-0.3}^{+0.5}$	${2.5}_{-1.2}^{+0.7}$	${14.4}_{-0.3}^{+0.3}$
THM	79.10 ± 0.19	${8.0}_{-0.4}^{+0.3}$	${11.8}_{-0.4}^{+0.4}$	${2.4}_{-1.2}^{+0.7}$	${14.4}_{-0.3}^{+0.3}$
TPHM PV = 300 FS = 4 ℓ ≤ 2	83.47 ± 0.14	${7.8}_{-0.3}^{+0.3}$	${12.2}_{-0.4}^{+0.3}$	${2.7}_{-1.1}^{+0.8}$	${14.7}_{-0.3}^{+0.3}$
TPHM PV = 300 FS = 4 ℓ ≤ 3	83.45 ± 0.19	${8.0}_{-0.3}^{+0.3}$	${12.2}_{-0.4}^{+0.3}$	${2.7}_{-1.2}^{+0.8}$	${14.8}_{-0.3}^{+0.3}$
TPHM PV = 300 FS = 4 ℓ ≤ 4	81.93 ± 0.24	${8.0}_{-0.4}^{+0.3}$	${12.0}_{-0.4}^{+0.3}$	${2.6}_{-1.1}^{+0.7}$	${14.6}_{-0.3}^{+0.3}$
TPHM PV = 300 FS = 4 ℓ ≤ 5	81.90 ± 0.21	${8.0}_{-0.4}^{+0.3}$	${12.0}_{-0.4}^{+0.3}$	${2.6}_{-1.1}^{+0.6}$	${14.6}_{-0.3}^{+0.3}$
TPHM PV = 300 FS = 2	81.88 ± 0.23	${8.0}_{-0.4}^{+0.3}$	${12.0}_{-0.4}^{+0.3}$	${2.6}_{-1.1}^{+0.7}$	${14.6}_{-0.3}^{+0.2}$
TPHM PV = 22311 FS = 3	81.86 ± 0.30	${8.0}_{-0.4}^{+0.3}$	${12.0}_{-0.4}^{+0.3}$	${2.6}_{-1.2}^{+0.8}$	${14.6}_{-0.3}^{+0.3}$

Note. We note that the highest ${ \mathcal B }{ \mathcal F }$ values are recovered by IMRPhenomTPHM with reduced mode content. This is consistent with the slightly negative Bayes factor for dominant versus higher modes reported for the NRSur7dq4 model in Abbott et al. (2020b) and Abbott et al. (2020c), but as discussed in the Appendix, the posteriors become much better resolved once including modes up to ℓ ≤ 4, and seem mostly converged in comparison to adding further modes ℓ ≤ 5; hence, we use the ℓ ≤ 4 run as our main result in this paper.

Download table as:  ASCIITypeset image

Table 3. Source Properties for GW190521, Listed as Median Posterior Values with Error Estimates Given by the 90% Credible Intervals

Waveform Model	NRSur7dq4	Pv3HM	v4PHM	XPHM (NC)	XPHM	TPHM ℓ ≤ 4	TPHM ℓ ≤ 5
Primary BH mass m1	${85}_{-14}^{+21}$	${90}_{-16}^{+23}$	${99}_{-19}^{+42}$	${129}_{-37}^{+46}$	${97}_{-21}^{+34}$	${109}_{-22}^{+80}$	${107}_{-20}^{+68}$
Secondary BH mass m2	${66}_{-18}^{+17}$	${65}_{-18}^{+16}$	${71}_{-28}^{+21}$	${32}_{-17}^{+33}$	${59}_{-25}^{+22}$	${65}_{-34}^{+28}$	${68}_{-33}^{+26}$
Total BBH mass M	${150}_{-17}^{+29}$	${154}_{-16}^{+25}$	${170}_{-23}^{+36}$	${169}_{-20}^{+23}$	${154}_{-16}^{+35}$	${181}_{-27}^{+44}$	${179}_{-25}^{+39}$
Binary chirp mass ${ \mathcal M }$	${64}_{-8}^{+13}$	${65}_{-7}^{+11}$	${71}_{-10}^{+15}$	${55}_{-16}^{+14}$	${64}_{-10}^{+15}$	${71}_{-11}^{+16}$	${72}_{-11}^{+16}$
Mass ratio q = m2/m1	${0.79}_{-0.29}^{+0.19}$	${0.73}_{-0.29}^{+0.24}$	${0.74}_{-0.42}^{+0.23}$	${0.23}_{-0.14}^{+0.46}$	${0.61}_{-0.36}^{+0.32}$	${0.63}_{-0.46}^{+0.32}$	${0.66}_{-0.46}^{+0.29}$
Primary BH spin χ1	${0.69}_{-0.62}^{+0.27}$	${0.65}_{-0.57}^{+0.32}$	${0.80}_{-0.58}^{+0.18}$	${0.84}_{-0.46}^{+0.12}$	${0.67}_{-0.59}^{+0.30}$	${0.86}_{-0.56}^{+0.12}$	${0.84}_{-0.56}^{+0.14}$
Secondary BH spin χ2	${0.73}_{-0.64}^{+0.24}$	${0.53}_{-0.48}^{+0.42}$	${0.54}_{-0.48}^{+0.41}$	${0.57}_{-0.44}^{+0.32}$	${0.55}_{-0.49}^{+0.4}$	${0.56}_{-0.50}^{+0.39}$	${0.56}_{-0.50}^{+0.39}$
Primary BH spin tilt angle ${\theta }_{{{LS}}_{1}}$	${81}_{-53}^{+64}$	${80}_{-49}^{+64}$	${81}_{-45}^{+49}$	${132}_{-54}^{+17}$	${117}_{-81}^{+44}$	${80}_{-32}^{+54}$	${85}_{-37}^{+52}$
Secondary BH spin tilt angle ${\theta }_{{{LS}}_{2}}$	${85}_{-55}^{+57}$	${88}_{-58}^{+63}$	${93}_{-60}^{+61}$	${84}_{-44}^{+48}$	${82}_{-57}^{+68}$	${97}_{-64}^{+57}$	${97}_{-63}^{+57}$
Effective spin parameter χeff	${0.08}_{-0.36}^{+0.27}$	${0.06}_{-0.39}^{+0.31}$	${0.06}_{-0.35}^{+0.34}$	$-{0.46}_{-0.14}^{+0.55}$	$-{0.11}_{-0.47}^{+0.43}$	${0.07}_{-0.44}^{+0.32}$	${0.02}_{-0.41}^{+0.36}$
Precession spin parameter χp	${0.68}_{-0.37}^{+0.25}$	${0.60}_{-0.44}^{+0.33}$	${0.74}_{-0.40}^{+0.21}$	${0.57}_{-0.25}^{+0.19}$	${0.49}_{-0.34}^{+0.34}$	${0.78}_{-0.39}^{+0.17}$	${0.76}_{-0.39}^{+0.19}$
Remnant BH mass Mf (M⊙)	${142}_{-16}^{+28}$	${147}_{-15}^{+23}$	${162}_{-22}^{+35}$	⋯	${148}_{-15}^{+35}$	${173}_{-25}^{+46}$	${171}_{-19}^{+28}$
Remnant BH spin χf	${0.72}_{-0.12}^{+0.09}$	${0.72}_{-0.15}^{+0.11}$	${0.74}_{-0.14}^{+0.12}$	⋯	${0.63}_{-0.24}^{+0.17}$	${0.75}_{-0.18}^{+0.13}$	${0.72}_{-0.14}^{+0.13}$
Radiated energy Erad	${7.6}_{-1.9}^{+2.2}$	${7.2}_{-2.2}^{+2.7}$	${7.8}_{-2.3}^{+2.8}$	⋯	${6.1}_{-3.3}^{+3.4}$	${7.2}_{-3.2}^{+3.1}$	${7.3}_{-3.1}^{+3.0}$
Luminosity distance DL	${5.3}_{-2.6}^{+2.4}$	${4.6}_{-1.6}^{+1.6}$	${4.0}_{-1.8}^{+2.0}$	${2.9}_{-1.4}^{+4.1}$	${3.5}_{-2.0}^{+2.4}$	${3.5}_{-1.7}^{+1.9}$	${3.4}_{-1.6}^{+2.0}$
Source redshift z	${0.82}_{-0.34}^{+0.28}$	${0.73}_{-0.22}^{+0.20}$	${0.64}_{-0.26}^{+0.25}$	${0.33}_{-0.15}^{+0.36}$	${0.59}_{-0.3}^{+0.32}$	${0.58}_{-0.25}^{+0.26}$	${0.56}_{-0.23}^{+0.27}$

Note. The first three results columns correspond to the results reported in Abbott et al. (2020b) and Abbott et al. (2020c), the fourth column summarizes results from Nitz & Capano (2021), and the last three columns are the new results from this paper, taken from our LALInference default runs with the standard versions of the IMRPhenomXPHM and IMRPhenomTPHM waveform models (including two choices of mode content for IMRPhenomTPHM, which yield very similar results).

Download table as:  ASCIITypeset image

To demonstrate the overall behavior of the GW190521 signal and the quality of the match with our waveform models, in Figure 4 we show the $\max { \mathcal L }$ templates from several IMRPhenomTPHM runs with different mode content compared to the whitened detector data at the time of GW190521 and a waveform reconstruction from the unmodeled cWB analysis (Klimenko et al. 2016; Abbott et al. 2020b, 2020c). The detector data is best matched by the ringdown region of the IMRPhenomTPHM model, while the cycles before merger are suppressed once whitened by the instrument's noise amplitude spectral density.

Zoom In Zoom Out Reset image size

Below we will discuss details and consequences of the posterior results shown in Figures 2 and 3, and we will analyze further posterior distributions, starting with different versions of the non-precessing versions of our models in Section 5.2, where we find good agreement between them. This serves as the more solid basis for the more challenging analysis with precessing models. We then use the IMRPhenomXPHM frequency-domain model in Section 5.3 and compare with the analysis of Nitz & Capano (2021), and discuss effects of a code change we have implemented in the default version of IMRPhenomXPHM, tracking the flipping of direction of the total angular momentum J in the same way as for non-default versions. We then investigate the case for multi-modality in the mass parameters reported by Nitz & Capano (2021) in Section 5.4, showing that we recover a multimodal mass posterior both with the IMRPhenomXPHM and IMRPhenomTPHM models, although with modified details compared to Nitz & Capano (2021). Then in Section 5.5 we investigate the support for precession in the source system, comparing results obtained with precessing and non-precessing approximants. Finally we study the implications for component masses in the mass gap and the support for the association with the AGN flare as a possible electromagnetic counterpart in Sections 5.6 and 5.7. Further analysis of the importance of the multimode harmonic content is presented in the Appendix. A crucial part of this analysis is to build confidence in our results by showing consistency between results obtained with different priors and different samplers (nested sampling; Skilling 2004, as implemented in pBilby, and MCMC as available through LALInference).

Before turning our attention to precessing models, we will inspect results obtained with non-precessing waveform approximants. In this simplified context, current waveform models have reached a certain level of maturity, where all state-of-the-art versions have been calibrated to NR simulations, including the subdominant harmonics content, to a varying degree. Therefore, we expect good agreement between different models, at least when the same subdominant mode content is included. We show results for LALInference runs with IMRPhenomTHM and IMRPhenomXHM in Figure 5. One can see that there is consistency between IMRPhenomXHM and IMRPhenomTHM when the same set of modes is included, which implies disabling the (3,2) mode in IMRPhenomXHM and restricting to ${{\ell }}_{\max }=4$ in IMRPhenomTHM. We do observe larger differences when including all of the available modes in each model, with a shift toward slightly lower q and mild multi-modality in the distance and inclination parameters for IMRPhenomXHM, although joint distributions still look broadly consistent with IMRPhenomTHM. We also notice that the recovered mass ratio and effective spins are consistent with the values reported in the LVC publications. The same goes for other key parameters, such as source-frame masses, distance, and inclination (see Figures 1 and 6 in Abbott et al. 2020c). For all of these results, both component masses lie confidently within the PISN mass gap (at 90% credible intervals).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The results reported in Nitz & Capano (2021) for the IMRPhenomXPHM model were obtained with the default version of the model (corresponding to MSA Euler angles and final-spin version FS = 3). The posteriors obtained by Nitz & Capano (2021) have nonzero support in regions of parameter space where the direction of the total angular momentum J flips (see Section 2.2) and would thus require careful cross-checks for robustness, as discussed previously. This is due to the fact that for the default version, we had initially implemented a different behavior as for other options: instead of attempting to track the direction of the total angular momentum J, a warning message was to be printed, alerting the user that the model is less reliable in case of flipped J. After the publication of Nitz & Capano (2021) we however realized that the warning messages had not been printed correctly when the calculation of subdominant harmonics was activated. To avoid confusion, we have more recently implemented a change harmonizing the behavior of the different final-spin versions, and the code now always tracks the direction of J for all parameter settings; this is now also described in the recently updated Section IV D of Pratten et al. (2021).

With this change, all final-spin versions now produce consistent results, as shown in Figure 6, with a much reduced support for the parameter region where the mass ratio is high and the effective spin negative, and thus where J may flip its sign. In particular, we note that, using the latest code version, the support for both masses being outside the mass gap is drastically reduced; see Table 5. Consistent results with the updated default version have also been reported by Capano et al. (2021) and Mehta et al. (2021), but here we present the first direct comparison using multiple final-spin versions. We also find in Figure 6 that when changing the final-spin version, the position of the maximum-likelihood sample changes considerably; this is however not surprising, as discussed in Section 4.4.

We now turn to examining the results obtained with the default settings of IMRPhenomXPHM and IMRPhenomTPHM with LALInference and pBilby, where in both cases the two approximants were run with the same priors and sampler settings. As we have already mentioned, pBilby posteriors have been re-weighted to allow for a direct comparison with LALInference results; see Section 4.3. Results are shown in Figure 2. One can appreciate a remarkable consistency between the two sampling codes. It is also clear that mass-ratio posteriors have a multimodal behavior for both models. The main difference here is that more unequal-mass ratios (q ∼ 0.25) in IMRPhenomXPHM are correlated with large negative χ_eff while the unequal-mass-ratio support for IMRPhenomTPHM is correlated with moderate positive χ_eff. Compared to inference with the aligned-spin models, support for the components to lie within the mass gap is reduced (see left panel); however, we defer an extensive discussion of this point to Section 5.6. In line with Nitz & Capano (2021), we find evidence for at least one high-mass-ratio mode at Q ≈ 5, in addition to the mode with near-equal masses as originally reported by the LVC analysis (right panel).

In Figure 7 we can see a comparison of the highest-S/N values for the default LALInference runs with both IMRPhenomXPHM and IMRPhenomTPHM. We can observe that the q ∼ 0.2 region produces similar S/Ns for both models, but IMRPhenomTPHM has support for higher S/Ns in the close-to-equal mass region. IMRPhenomTPHM also recovers a small strip at q ∼ 0.1 at more or less the same height as the q ∼ 0.2 bulk. Differences in network matched-filter S/N here are only about 0.15 (which corresponds to 3.86 in $\max { \mathcal L }$). For IMRPhenomXPHM, the maximum S/N is located at q = 0.26 while for IMRPhenomTPHM, it is located at q = 0.975, in agreement with the $\max { \mathcal L }$ positions shown in the right panel of Figure 2.

Zoom In Zoom Out Reset image size

In order to better explore the regime of very unequal masses, we have performed IMRPhenomTPHM runs with restricted mass-ratio priors, both with LALInference and with pBilby. In Figure 8 one can see that results are consistent with not finding particular support for another mode below the one at q ∼ 0.2. The full-prior run is poorly sampled in this region, with only 3.6% of samples below q = 0.15, so the small peak at q = 0.06 is probably an artifact from insufficient resolution in this region, and it is not recovered by the restricted runs. The maximum S/N recovered by the restricted runs is also lower than the S/N recovered by the full run near-equal masses.

Zoom In Zoom Out Reset image size

We have also checked that results are robust for different IMRPhenomTPHM versions, as can be seen in Figure 9. As discussed in the Appendix, the bimodality is not recovered when using only the dominant ℓ ≤ 2 spherical harmonic modes, but is robust under inclusion of different subsets of the higher-order modes implemented in IMRPhenomTPHM. Therefore, we overall find clear evidence for a bimodal mass-ratio posterior with a secondary unequal-mass-ratio peak near Q ∼ 5, but while there is also nonvanishing posterior support reaching to even more extreme ratios, we find no clear evidence for a third peak as originally reported in Nitz & Capano (2021). These results are consistent with Capano et al. (2021) and Mehta et al. (2021).

Zoom In Zoom Out Reset image size

In terms of spins, one can see in Figure 3 that while IMRPhenomXPHM recovers a posterior that is approximately symmetric around χ_p = 0.5, IMRPhenomTPHM strongly favors high values of χ_p. For χ_eff, the IMRPhenomXPHM posterior is bimodal, with peaks close to χ_eff = 0 and at large negative values, while the IMRPhenomTPHM posterior is broad but unimodal with the median close to small positive values. This is reflected also in the position of the $\max { \mathcal L }$ points from the posteriors: while for IMRPhenomTPHM results, both with pBilby and LALInference, $\max { \mathcal L }$ lies at low positive χ_eff and high χ_p ≥ 0.8, for IMRPhenomXPHM, the $\max { \mathcal L }$ is located around χ_p ∼ 0.5, and its χ_eff is large and negative.

Another thing to notice in relation to Figure 3 is that adding precession does not seem to significantly affect the distance–θ_JN joint distribution for IMRPhenomXPHM with respect to IMRPhenomXHM, while for the time-domain IMRPhenomT model, family precession indeed helps to better constrain the posteriors. This is supported by the Bayes factors for the precessing hypothesis (computed from the difference in log signal-to-noise Bayes factors between precessing and non-precessing results) shown in Table 4, where the time-domain models show stronger support for precession than the frequency-domain models.

Table 4. Comparison of Bayes Factors between Precessing and Non-precessing Approximants

	TPHM vs THM (${{\ell }}_{\max }=4$)	TPHM vs THM (${{\ell }}_{\max }=5$)	XPHM vs XHM	TPHM vs XHM	NRSur7dq4 (LVC)
Bayes factor	${17.0}_{-4.6}^{+6.2}$	${13.5}_{-3.0}^{+3.9}$	${1.5}_{-0.3}^{+0.4}$	${6.5}_{-1.6}^{+2.2}$	${11.5}_{-1.1}^{+1.1}$

Download table as:  ASCIITypeset image

In Table 5 we show the posterior probabilities for the component objects of the source of GW190521 being inside the PISN gap. The exact boundaries of this gap are not known exactly and depend in a highly nontrivial way on nuclear reaction rates and several aspects of stellar dynamics. For simplicity, we will provide probabilities computed assuming a gap either in the range [50, 120] M_⊙ (Woosley 2017) or, following more recent estimates (Woosley & Heger 2021), in the range [70, 161] M_⊙ (we report results for the latter case in parentheses). Shifting the estimated boundaries toward higher masses increases (decreases) the probability of the heavier (lighter) component lying in the gap; it also drastically decreases the probability of both masses being in the mass gap.

Table 5. Probabilities of Component Objects to Be in the PISN Mass Gap, Assumed to Be Either in the Range [50, 120] M_⊙ or in [70, 161] M_⊙ (in Parentheses)

	$P({m}_{1,\mathrm{src}}^{\mathrm{PISN}})$	$P({m}_{2,\mathrm{src}}^{\mathrm{PISN}})$	$P({m}_{1,\mathrm{src}}^{\mathrm{PISN}}\&{m}_{2,\mathrm{src}}^{\mathrm{PISN}})$	$P({m}_{1,\mathrm{src}}^{\mathrm{PISN}}| | {m}_{2,\mathrm{src}}^{\mathrm{PISN}})$
TPHM PV = 300 FS = 4 ℓ ≤ 3	0.672(0.919)	0.843(0.612)	0.663(0.612)	0.852(0.919)
TPHM PV = 300 FS = 4 ℓ ≤ 4	0.666(0.872)	0.715(0.415)	0.646(0.415)	0.736(0.872)
TPHM PV = 300 FS = 4 ℓ ≤ 5	0.741(0.916)	0.782(0.460)	0.719(0.460)	0.804(0.916)
TPHM PV = 300 FS = 2	0.661(0.871)	0.705(0.414)	0.640(0.414)	0.726(0.871)
TPHM PV = 22311 FS = 3	0.629(0.870)	0.713(0.463)	0.613(0.463)	0.728(0.870)
TPHM PV = 300 FS = 4 fixed sky location	0.648(0.807)	0.644(0.442)	0.625(0.442)	0.667(0.808)
TPHM PV = 300 FS = 4 fixed 3D localization	0.896(0.914)	0.897(0.780)	0.886(0.780)	0.907(0.914)
THM	0.993(0.997)	0.857(0.182)	0.856(0.182)	0.996(0.997)
XPHM N&C	0.470(0.675)	0.372(0.048)	0.365(0.048)	0.477(0.675)
XPHM LI PV = 223 FS = 3	0.910(0.993)	0.785(0.166)	0.750(0.166)	0.946(0.993)
XPHM PB 1/q PV = 223 FS = 3	0.681(0.923)	0.495(0.072)	0.474(0.072)	0.701(0.924)
XPHM PB 1/q PV = 223 FS = 2	0.711(0.970)	0.510(0.069)	0.493(0.069)	0.728(0.970)
XPHM PB 1/q PV = 223 FS = 0	0.658(0.959)	0.479(0.067)	0.464(0.067)	0.674(0.959)
XHM LI	0.994(0.993)	0.847(0.125)	0.845(0.125)	0.997(0.993)

Note. Numbers reported here are computed as the ratio between the posterior samples inside the gap and the total sample size. The runs denoted with 1/q use priors uniform in Q = 1/q = m₁/m₂, while all others use the default prior uniform in q = m₂/m₁; see Section 5.6 for the full discussion.

Download table as:  ASCIITypeset image

We observe that IMRPhenomXPHM with a prior that is flat in component masses (run with LI) has more support for the PISN gap hypothesis, as do the non-precessing models (both IMRPhenomXHM and IMRPhenomTHM). One can see that for the IMRPhenomXPHM LI run, sampled uniform in component masses, the support for at least one component in the mass gap is greater than 90%. The alternative prior uniform in Q = 1/q enhances the high-mass-ratio region of the posterior, which typically makes both components lie outside the gap. In the runs we performed with pBilby and this prior, the support for components in the mass gap drops for IMRPhenomXPHM. With IMRPhenomTPHM, the mass-gap probabilities are generally lower than for IMRPhenomXPHM.

Hence, we conclude that inference with the non-precessing models, or with IMRPhenomXPHM and the uniform-in-q prior, prefer the hypothesis of at least one object in the PISN gap (>9:1 probability ratio), matching the original conclusions of Abbott et al. (2020b) and Abbott et al. (2020c); but IMRPhenomXPHM runs with uniform-in-Q prior and IMRPhenomTPHM inference have only mild preference (ranging from ∼2:1 to ∼3:1 depending on model version and mode content) for this scenario, allowing more readily for the "straddling" hypothesis from Fishbach & Holz (2020) of the heavier BH above and the lighter BH below the PISN gap. As seen in Table 5, results for IMRPhenomTPHM with different mode content (ℓ ≤ 3, 4, 5) are broadly consistent, except for a larger probability with ℓ ≤ 3 of having the smaller mass, or either of the two masses, in the mass gap.

The tentative association of GW190521 with the AGN flare ZTF19abanrhr (Graham et al. 2020) would be hugely impactful as the first detected electromagnetic counterpart to a BBH merger and open up cosmological applications (Chen et al. 2020; Mukherjee et al. 2020; Gayathri et al. 2020b; Finke et al. 2021). This association is however far from certain (Ashton et al. 2020; Palmese et al. 2021). Still, below we report updated association probabilities based on our new parameter estimation results. It is also illustrative to study how the GW190521 posteriors would change when factoring in knowledge of the sky location and distance of this counterpart as priors on the GW parameter estimation. To do so, we consider two scenarios: in the first, we only fix the R.A. and decl. of the source, while in the second, we also fix its luminosity distance, which can be calculated from the flare's redshift assuming a standard cosmology (Planck Collaboration et al. 2016).

First, looking back at the d_L –θ_JN panel of Figure 3 for the runs with standard priors (no constraints on source location), we can see that while none of the models are able to break the hemisphere degeneracy, IMRPhenomTPHM results are able to better constrain a particular inclination in each hemisphere, albeit some bimodality in each hemisphere is also present. Similarly, IMRPhenomTPHM appears to deliver a better constraint on the luminosity distance than IMRPhenomXPHM, for which there are also hints of bimodality in the posterior. We note also that the $\max { \mathcal L }$ sample of the IMRPhenomXPHM posterior has a lower luminosity distance than that of the IMRPhenomTPHM posterior, and that these results are consistent between LALInference and pBilby runs.

Switching to IMRPhenomTPHM runs with counterpart-informed priors, in Figure 10 we present the results for masses and spins. Fixing either the 2D sky location or the full 3D localization to those of the AGN flare, the support at mass ratio Q greater than 5 is enhanced in both cases. However, fixing the 3D localization also constrains the mass posteriors more overall, while fixing only the 2D sky location mostly extends support to higher Q (lower q). The 2D-fixed prior produces a slightly bimodal posterior in χ_eff that is however broadly consistent with that from the standard prior, while the run with fixed 3D localization shifts the posterior mode of χ_eff toward 0. For the precessing spin parameter χ_p, fixing only the sky position has no noticeable effect, while fixing the 3D localization shifts the recovered posterior distribution to milder values. It is also worth noticing that fixing the 3D localization increases the probability for the companions to be inside the PISN mass gap (see Table 5) with respect to the standard run, while fixing only the sky location prior does not seem to have an effect on this. We also report the Bayes factors for these fixed-prior runs over the default runs in Table 6, which are quite high, but require a full analysis of the actual association probability to interpret.

Zoom In Zoom Out Reset image size

Table 6. Comparison of Bayes Factors for the Runs with Either 2D Sky Location or Full 3D Localization (Sky Position Plus Luminosity Distance) Fixed to the AGN Counterpart Candidate (Graham et al. 2020), Against the Default Runs with Unconstrained Localization Priors

	BF Fixed Sky Location	BF Fixed 3D Location
TPHM ${{\ell }}_{\max }=4$	${14}_{-5}^{+7}$	${62}_{-14}^{+18}$
TPHM ${{\ell }}_{\max }=5$	${13}_{-5}^{+7}$	${64}_{-16}^{+21}$
XPHM ${{\ell }}_{\max }=4$	${23}_{-3}^{+3}$	${95}_{-13}^{+15}$

Note. While these are quite high, it is important to take into account a full analysis of the multi-messenger coincidence significance, see the text and Table 7.

Download table as:  ASCIITypeset image

To quantify this probability of a physical association between the GW and electromagnetic (EM) signals, in Table 7 we show the results from performing the same multi-messenger coincidence significance analysis as presented in Ashton et al. (2020; and using their public code), based on the localization overlap. We give results for the posterior overlap integrals for the sky location (${{ \mathcal I }}_{{\rm{\Omega }}}$), for the distance alone (${{ \mathcal I }}_{{D}_{L}}$), and for the combined 3D localization (${{ \mathcal I }}_{{D}_{L}}{{ \mathcal I }}_{{\rm{\Omega }}}$). Assuming the same prior odds as in Ashton et al. (2020), based on the reported number of flares similar to ZTF19abanrhr in the ZTF alert stream, we compute the odds ${{ \mathcal O }}_{{\rm{C}}/{\rm{R}}}$ for the coincident hypothesis. The resulting odds vary by a factor ∼2 between runs with IMRPhenomXPHM and different versions of IMRPhenomTPHM (precession prescriptions, final-spin versions, and mode content) and overall show mild preference for association. However, our results are consistent with the findings of Ashton et al. (2020), with the highest ${{ \mathcal O }}_{{\rm{C}}/{\rm{R}}}$ we find at the same level as that reported for a run with the SEOBNRv4PHM model in Ashton et al. (2020). Essentially, despite the high Bayes factors, the evidence is insufficient to confidently associate the events. For further illustration, Figure 11 shows the position of ZTF19abanrhr compared with the sky location posterior density for GW190521 recovered by our default IMRPhenomTPHM run.

Zoom In Zoom Out Reset image size

Table 7. Posterior Overlap Integrals and Odds for the Association between GW190521 and ZTF19abanrhr (Graham et al. 2020), Following the Method from Ashton et al. (2020)

	${{ \mathcal I }}_{{D}_{L}}$	${{ \mathcal I }}_{{\rm{\Omega }}}$	${{ \mathcal I }}_{{D}_{L}}{{ \mathcal I }}_{{\rm{\Omega }}}$	${{ \mathcal O }}_{{\rm{C}}/{\rm{R}}}$
TPHM PV = 300 FS = 4 ${{\ell }}_{\max }=4$	4.7	17	140	6
TPHM PV = 300 FS = 4 ${{\ell }}_{\max }=5$	5.1	30	140	12
TPHM PV = 300 FS = 2 ${{\ell }}_{\max }=4$	4.6	33	140	12
TPHM PV = 22311 FS = 3 ${{\ell }}_{\max }=4$	4.7	20	110	7
XPHM LI PV = 223 FS = 3	3.8	19	110	6

Download table as:  ASCIITypeset image