Prospects for Measuring the Hubble Constant with Neutron-Star-Black-Hole Mergers
Stephen M. Feeney, Hiranya V. Peiris, Samaya M. Nissanke, Daniel J. Mortlock
PProspects for Measuring the Hubble Constant with Neutron-Star-Black-Hole Mergers
Stephen M. Feeney, Hiranya V. Peiris,
1, 2
Samaya M. Nissanke,
3, 4 and Daniel J. Mortlock
5, 6, 7 Department of Physics & Astronomy, University College London, Gower Street, London WC1E 6BT, UK Oskar Klein Centre for Cosmoparticle Physics, Department of Physics,Stockholm University, AlbaNova, Stockholm SE-106 91, Sweden GRAPPA, Anton Pannekoek Institute for Astronomy and Institute of High-Energy Physics,University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands Nikhef, Science Park 105, 1098 XG Amsterdam, The Netherlands Astrophysics Group, Imperial College London, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK Department of Mathematics, Imperial College London, London SW7 2AZ, UK Department of Astronomy, Stockholm University, AlbaNova, SE-10691 Stockholm, Sweden (Dated: January 26, 2021)Gravitational wave (GW) and electromagnetic (EM) observations of neutron-star-black-hole(NSBH) mergers can provide precise local measurements of the Hubble constant ( H ), ideal forresolving the current H tension. We perform end-to-end analyses of realistic populations of simu-lated NSBHs, incorporating both GW and EM selection for the first time. We show that NSBHscould achieve unbiased 1.5–2.4% precision H estimates by 2030. The achievable precision is stronglyaffected by the details of spin precession and tidal disruption, highlighting the need for improvedmodeling of NSBH mergers. Introduction.
The current expansion rate of theUniverse – the Hubble constant, H – is at the heartof a significant cosmological controversy. Direct mea-surements in the local Universe by the SH0ES team’sCepheid-supernova distance ladder [1] find H = 74 . ± .
42 km s − Mpc − . This is discrepant at the 4.4- σ levelfrom the 67 . ± .
54 km s − Mpc − value inferred fromthe Planck satellite’s observations of the cosmic mi-crowave background (CMB) anisotropies, assuming thestandard flat cosmological model [2].There are two potential explanations for this discrep-ancy, the most exciting of which derives from the model-dependence of the CMB constraint: could the discrep-ancy be due to physics beyond the standard model? De-spite extensive effort (e.g., Refs. [3, 4]), consensus on acompelling theoretical explanation has not been reached.The more prosaic explanation posits undiagnosed sys-tematic errors or underestimated uncertainties; however,despite multiple investigations of both the distance lad-der [5] and CMB [6] datasets, no study has found incon-trovertible evidence warranting a change of conclusions.In the absence of conclusive evidence of systematic er-rors or consensus on an extended model, independent ver-ifications of the two central measurements offer a promis-ing route to resolving the tension. Independent verifi-cation of the CMB anisotropy constraints comes fromrecent inverse distance ladder datasets [7]. Local verifi-cation has, however, proven more challenging, with somealternative analyses supporting the SH0ES team’s find-ings [8] and others providing contradictory conclusionsof varying significance [9], in some cases using the samedata.A direct, completely independent local measurementwith percent-level precision is therefore needed to resolvethe H tension. Combined gravitational-wave (GW) andelectromagnetic (EM) observations of nearby compact-object mergers are ideal candidates to provide that mea- surement, yielding H estimates that depend on gen-eral relativity alone [10–33]. Thanks to their accompa-nying EM emission, the utility of binary neutron star(BNS) mergers is well established [12–25], but less atten-tion has been paid to the potential contribution of as-yet undiscovered neutron-star-black-hole (NSBH) merg-ers with EM counterparts [13, 16, 34]. Using idealized,fixed-signal-to-noise simulations at indicative parametervalues, Vitale and Chen [34] recently showed that cata-logs of GW-selected NSBH observations may constrain H as well as BNSs, depending on the relative mergerrates and BH spins. In particular, they showed explicitlythat luminosity distance estimates could improve, as mis-aligned BH spins induce spin precession, helping breakthe degeneracy between the luminosity distance and in-clination angle for some NSBH systems.Here, we determine the H constraints realistic NSBHsamples will achieve, by performing end-to-end analysesof simulated NSBH samples incorporating fully specifiedparent populations, combined GW and EM selection, anda complete noise treatment. We use state of the art GWwaveforms [35, 36] and EM outflow models [37], bothcalibrated to a suite of numerical relativity simulations,for our GW and EM signals, and focus on the “A+”era of the mid-to-late 2020s, assuming an expanded GWnetwork including LIGO India and KAGRA. Simulations.
In this work we simulate the results ofa circa-2025 GW detector network, consisting of LIGOA+, Virgo AdV+, KAGRA and LIGO India [38, 39] ob-serving for t obs = 5 yr with duty cycle of ∆ obs = 0 . fid = 610 yr − Gpc − (corresponding to the 90% upperlimit of Ref. [40]), and cosmological parameters matchingRef. [2], with H = 67 .
36 km s − Mpc − and q = − . λ = ∆ obs t obs V Γ fid , where V is the redshifted volume a r X i v : . [ a s t r o - ph . C O ] J a n d e n s i t y ( S E O B N R ) prior sampleGW selectedGW+EM selected − . . .
00 1 2 30 1000 2000 distance , d [Mpc] d e n s i t y ( I M R P h e n o m ) mass ratio , q − BH z spin , a z BH . . . BH spin magnitude , | a BH | inclination , ι . . . . . . . . . . . σ H / H FIG. 1. Distributions of a subset of parameters from our
SEOBNR (top) and
IMRPhenom (bottom) samples, as drawn from theprior (dotted), selected by GW SNR (dashed) and selected by GW and EM emission (colored histograms). The bins are coloredby the fractional H uncertainty the mergers within the bin achieve: the yellowest/lightest bins are most informative. (Eq. A3) calculated using a third-order cosmographicco-moving volume element (Eq. A2). To reduce com-putation time, the volume integral can be truncated atsome redshift z max where there is negligible probabilityof even the loudest merger being detected: we find that z max = 0 .
44 suffices for our setting. For our fiducial pa-rameter set, the mean number of mergers λ = 25160;our particular Poisson draw yields a pre-selection totalof 25241.For each merger, we draw a cosmological redshift, z ,assuming a constant source-frame rate (see, e.g., Eq. 28of Ref. [24]), along with an isotropically-distributed an-gular sky position, inclination angle, phase and polar-ization angle. We draw uniformly distributed BH massesfrom P( m BH ) = U(2 . M (cid:12) , M (cid:12) ), taking the upper limitfrom low-metallicity binary population synthesis simu-lations [41] and extending to low masses to reflect the de-tection of objects in the purported NS/BH mass gap [42].NS masses are drawn from P( m NS ) = U(1 M (cid:12) , . M (cid:12) ),with an upper limit chosen to match that of the DD2equation of state (EOS) [43]. Dimensionless BH andNS spin magnitudes are drawn from the uniform distri-butions P( a BH ) = U(0 , .
99) and P( a NS ) = U(0 , . Using solar metallicities would reduce this upper limit to 12 M (cid:12) . to zero [45]. Finally, we generate a peculiar velocity v foreach merger from a zero-mean normal with a standarddeviation of 500 km s − .With the NSBH parameters in hand, we gener-ate mock data for each merger and apply our selec-tion criteria. To determine the impact of differentphysical effects on our results, we simulate two pop-ulations using different waveform approximants: theBNS-calibrated IMRPhenomPv2 NRTidal [35] and NSBH-specific
SEOBNRv4 ROM NRTidalv2 NSBH [36] (hereafter
IMRPhenom and
SEOBNR ). We refer the reader to the Sup-plemental Material for a complete discussion of the dif-ferences between the two waveforms. The
SEOBNR wave-form requires aligned or anti-aligned spins, so we set thetransverse NS and BH spins to zero after sampling themisotropically (mimicking, in some sense, spins becomingaligned over time). For each merger, we generate a 32-second segment of noisy GW data ˆ x per waveform us-ing the same random seed and a frequency range of 20–2048 Hz, considering it detected if the network signal-to-noise ratio (SNR) is at least ρ ∗ = 12. We assumethat the GW detectors operate in concert with an EMfollowup program capable of detecting all mergers withejecta mass greater than m ∗ ej = 0 . M (cid:12) , modeling for thefirst time a hybrid GW-EM selection function for NSBHs.Finally, we generate noisy measured redshifts and pecu-liar velocities by drawing from P(ˆ z | z ) = N( z, . v | v ) = N( v,
200 km s − ), respectively. Of the 25241simulated mergers, 2477 (2954) are detected in GWs us-ing the IMRPhenom ( SEOBNR ) waveform, 99 (75) of whichhave sufficient ejecta to be detected in EM; 62 appear in Using spectra from Ref. [39]. both samples. The SNRs for
SEOBNR waveforms are, onaverage, 5.9% larger than their
IMRPhenom counterparts,resulting in the GW-detected
SEOBNR sample containing ∼
500 more objects. Setting the transverse spins to zerofor use with the
SEOBNR waveform, however, has the side-effect of reducing the typical ejecta mass [37] and hencethe final GW+EM-detected sample.The impact of our selection function is illustrated inFig. 1, in which we plot histograms of our full popu-lation (dotted lines), GW-selected events (dashed lines)and GW+EM-selected mergers (colored bars) for a sub-set of our parameters. The prior curves are identical inboth cases apart from the BH spin magnitudes, wherezeroing the transverse spins has made the
SEOBNR popu-lation’s distribution non-uniform. The primary impact ofthe GW SNR threshold is, as expected, to select nearbymergers; it also imparts very slight preferences for lowmass ratios and prograde BH z spins [46]. It is interest-ing to note that the GW-selected SEOBNR distance distri-bution is broader than that of
IMRPhenom and peaked atslightly higher distances: this is a direct consequence ofthe
SEOBNR injections’ systematically higher SNRs. Fur-ther, the presence of spin precession permits the detec-tion of more edge-on
IMRPhenom waveforms.The ejecta-mass threshold (i.e., EM selection) stronglyimpacts the observed distributions. The GW+EM-detected distributions are shifted to even smaller dis-tances, particularly for the
SEOBNR waveform, as the low-mass-ratio systems which produce significant ejecta masscan only be detected nearby. There is a very strongpreference for mass ratios under 10 (again, especially sofor
SEOBNR ) and large spins [37, 44], and the preferencefor positive z spins is much more pronounced. The dif-ferences between the two waveforms’ GW+EM-selecteddistributions are slightly obfuscated by the small sam-ple sizes, but the SEOBNR sample is shifted towards lowerdistances and mass ratios. As the
SEOBNR mergers’ BHspin magnitudes are smaller than their
IMRPhenom coun-terparts’, they require smaller mass ratios to produce sig-nificant ejecta [37], and the resulting systems are harderto detect at distance. Our implementation of EM se-lection captures the current best understanding of thedependence on ejecta mass, but we note that a fully self-consistent model of EM selection does not yet exist. Thisselection does not, for example, incorporate any viewing-angle dependence (see Ref. [47] for a treatment in BNSmergers) or EM survey selection effects [e.g., 48–51].
Methods.
The probabilistic inference of the Hubbleconstant from catalogs of compact object mergers hasbeen described in detail in the literature [10, 12–14, 16, We hypothesise that this is due to the effects of generic spin pre-cession on the
IMRPhenom population. Differences in the lengthsof GW signals within the detector frequency bands due to thetwo waveforms’ different merger frequencies would tend to boostthe
IMRPhenom
SNRs. Populations with ∼ solar metallicities produce NSBHs with lowerBH masses [41] and hence more GW+EM-detectable mergers.
19, 21–25, 52, 53]. In the following, we adopt a slightvariant of the formalism set out in Ref. [24], whose Fig.9 depicts a network diagram for the model we use todescribe the data . The posterior we evaluate is definedin Eq. A5.We infer the parameters of this model in two parts,using two sampling methods. First, we process eachmerger individually in order to obtain the GW likeli-hoods marginalized over all parameters θ i other than theluminosity distance d to the merger ( θ i here comprisingthe i th merger’s component masses, spin magnitudes andorientations where used, inclination, polarization angle,NS tidal deformability, and time and phase at coales-cence). We adopt priors identical to the distributionsused in the generative model for all parameters otherthan the masses. Convergence is greatly improved bysampling chirp masses and mass ratios instead of com-ponent masses, and we therefore sample using interimpriors that are uniform in these parameters (over theranges permitted by our component-mass extrema), be-fore importance-sampling the outputs to reinstate ourdesired component-mass priors. The marginal GW like-lihoods are sampled with the pypolychord nested sam-pler [54], wrapped by bilby [55], using 1000 live pointsand bilby ’s marginalize phase , time and distance settings. Each 15 (11)-dimensional IMRPhenom ( SEOBNR )sampling run takes 6-14 (4-6) days to complete on oneIntel Xeon 2.7 GHz CPU.Given the marginal GW likelihoods, we use No-U-TurnSampling as implemented by the pystan package [56]to infer the cosmological and population parameters.To connect to the cosmological parameters, we adopt athird-order distance-redshift relation matching our vol-ume element [57], with H and q allowed to vary but thejerk set to one. We assume a broad Gaussian prior on H , P( H ) = N(70 ,
20) km s − Mpc − , a truncated Gaus-sian prior on q , P( q ) = Θ( q + 2)Θ(1 − q )N( − . , . ∝ / Γ. To use pystan , we must be able to sample all parameters fromanalytic distributions. We therefore perform a GaussianMixture Model fit to each merger’s marginal distancelikelihood using pomegranate [58]. We fit each likeli-hood with an integer grid of 2–10 mixture components,repeating 10 times at each grid point and selecting thebest fit using the Akaike Information Criterion [59].Finally, we must evaluate the expected number of de-tected mergers ¯ N at each sampled value of the cosmo-logical and population parameters. We do so by re-simulating the catalogs 100 times at each point of a 5 × { H , q } assuming our fiducial rate, Γ fid , and in-terpolating the results using a 2D fourth-order interpola-tion. The dependence on the sampled rate is captured bymultiplying the interpolation coefficients by Γ / Γ fid . The The only addition required to the network diagram of Ref. [24]is the dependence of the selection, S , on an intrinsic parameter:the merger’s ejecta mass.
200 300 400 d L (Mpc) . . . . . . . ι ρ IMRP = 53 . ρ IMRP , al = 48 . ρ SEOB = 55 . d L (Mpc) ρ IMRP = 30 . ρ IMRP , al = 34 . ρ SEOB = 36 . d L (Mpc) ρ IMRP = 12 . ρ IMRP , al = 12 . ρ SEOB = 15 . . . . fractional uncertainty f r e q u e n c y IMRP d L IMRP H SEOB d L SEOB H FIG. 2. Left three panels: distance and inclination posteriors for a selection of mergers, simulated and sampled using the
IMRPhenom waveform with precessing (grey filled) and aligned (dark red dashed) spins, and using
SEOBNR with aligned spins(red). The selection includes the highest-SNR merger common to both catalogs (left) and the
IMRPhenom merger whose BH spinis closest to being aligned (second from right). Right: distributions of fractional uncertainties on luminosity distance (dotted)and H (solid) from individual mergers from our IMRPhenom (grey) and
SEOBNR (red) NSBH catalogs. resulting 153 (193)-dimensional pystan inference runstake less than a minute to generate 20,000 well convergedsamples on a 3.1 GHz Intel Core i7 CPU. The set of trueredshifts and peculiar velocities are uninteresting for thepurposes of cosmology inference, and we marginalize overthese parameters when quoting the results below.
Results.
Processing the simulated
SEOBNR and
IMRPhenom catalogs through our two-stage inferencepipeline produces the cosmology and population pa-rameter posteriors shown in the Supplemental Material.In both cases, the recovered H , q and rate poste-riors are completely consistent with the input values,indicating, as expected, that the selection effects arecorrectly accounted for [24]. The 68% credible inter-vals on the near-Gaussian H marginal posteriors are68 . ± . − Mpc − for the SEOBNR sample and 66 . ± . − Mpc − for IMRPhenom . As the
IMRPhenom sample contains 99 objects to the
SEOBNR sample’s 75,we should therefore expect that the
IMRPhenom sam-ple’s H posterior be roughly 13% narrower than thatof the SEOBNR sample. The remaining reduction, there-fore, reflects the ability for precessing spins to break thedistance-inclination degeneracy [e.g., 21, 60–62]. This ad-ditional constraining power is equivalent to an approxi-mate doubling of the catalog size.The H uncertainties we find for both waveforms arecomparable to the current Cepheid-SN distance ladderprecision [1]. NSBH populations – should they produceEM counterparts and occur at rates roughly matchingour assumptions – will therefore strongly inform the out-come of the current H tension, particularly when com-bined with accompanying BNS populations , likely of We note that comparable H precision had been hoped for by2023 from catalogs of BNS mergers [21]. That timescale, how-ever, now appears optimistic due to the lack of EM counterpartdetections following GW170817. comparable size [21, 23, 24, 34]. The mergers are alsoinformative about the deceleration parameter, q , shrink-ing its uncertainty from 0.5 to 0.32 or 0.27, depend-ing on the waveform. This further implies that NSBHcatalogs will be able to begin constraining parameterssuch as the matter density and dark energy EOS (in thecontext of ΛCDM and extended models), complemen-tary to BBH results from higher redshifts [31, 33, 63].The merger rates are recovered with roughly 10% preci-sion [e.g., 40, 42, 64, 65].To obtain a picture of the parameter combinations thatare most important for the H constraints, we return toFig. 1. Here, the colors of the histogram bins indicate thefractional H uncertainty the mergers within each bin at-tain, with the most-constraining bins colored yellow andthe least-constraining blue. For both waveforms, the bulkof the H constraining power comes from mergers out toroughly 700 Mpc, not just the very nearest ( ∼
200 Mpc),loudest events. For the
IMRPhenom mergers, all mass ratiobins less than ∼ M (cid:12) contribute equally, despite the fre-quency dropping rapidly with mass ratio. For the SEOBNR mergers, the constraints are instead driven by the lowest-mass ratio events [e.g., 60, 61]. From the
IMRPhenom spinpanels, it is clear that that highest-spin events constrain H most strongly, with the full H constraint comingalmost entirely from the highest-spin (and most popu-lated) bin. The SEOBNR constraint, on the other hand, issourced by events with a broader range of spins, thoughthis is likely driven by the prior.We further highlight the importance of precession inbreaking the distance-inclination degeneracy in Fig. 2.In the first three panels we plot distance and inclina-tion constraints for a selection of mergers when using the
SEOBNR (red, filled) and
IMRPhenom (grey, filled) wave-forms. For the first two mergers, detected with high SNR,the long degeneracies present in their
SEOBNR posteriorsare almost completely broken when using the
IMRPhenom waveform, in which the spins precess. We illustrate thispoint further by re-running the
IMRPhenom case assum-ing aligned spins, having set the transverse spins to zero.These results are overlaid as dashed dark red contours.Both distance-inclination degeneracies blow up, increas-ing the distance uncertainties by a factor of over three,with commensurate consequences for the mergers’ abil-ity to constrain H . In the third panel, we show equiv-alent posteriors for the IMRPhenom merger whose BHspin is closest to being aligned: the effect of switchingwaveforms here is markedly reduced. The impact onthe population level is clear from the right-hand panelof Fig. 2, in which we plot the distributions of individ-ual mergers’ fractional errors on distance (dashed) and H (solid) when using the SEOBNR (red) and
IMRPhenom (grey) waveforms. The
SEOBNR distributions are shiftedto significantly higher errors than their
IMRPhenom coun-terparts, despite the
SEOBNR mergers typically havinghigher SNRs. The smallest percentage error for any indi-vidual merger is 2.8% for the
IMRPhenom case and 6.1%for
SEOBNR ; the medians are 13.2% and 17.3%, respec-tively. The H constraints imparted by both “golden”and normal events are therefore stronger when spins pre-cess significantly. Finally, we note from the lower limitsof the dashed curves that peculiar velocity and redshift uncertainties strongly suppress the constraining power ofthe nearest and loudest events. Conclusions.
In this Letter, we present the resultsof the first end-to-end inference of H from realistic sim-ulated catalogs of NSBH mergers incorporating GW andEM selection effects. The precision we should expectfrom such catalogs is very promising for resolving thecurrent H tension, with five years of A+ era observa-tions yielding H uncertainties of 1.5–2.4%. We find,however, that the detailed physics of the NSBH wave-forms strongly impacts the achievable precision. Usingthe SEOBNRv4 ROM NRTidalv2 NSBH waveform with non-precessing BH spins results in boosted SNRs, and an in-crease of ∼
500 GW-detected NSBHs. However, includ-ing precessing spins using the
IMRPhenomPv2 NRTidal waveform markedly increases the typical ejecta massand hence the number of combined GW+EM detec-tions. Critically, precessing spins also break the distance-inclination degeneracy in the resulting GW parameterposteriors, yielding a significant improvement ( ∼
40% af-ter accounting for differing catalog sizes) in the resulting H constraint. Our results strongly highlight the need forimproved modeling of NSBH signals in both gravitationalwaves (see e.g., [66]) and the electromagnetic spectrum. Software.
The Python simulation and inference soft-ware developed for this analysis are publicly available at https://github.com/sfeeney/nsbh . Acknowledgments.
We thank Sukanta Bose for pro-viding details on LIGO India, Nikhil Sarin and GregAshton for help with bilby , Will Handley for help with pypolychord , and Tanja Hinderer, Andrew Williamson,Francois Foucart, and Bastien DuBoeuf for useful discus-sions. HVP’s work was partially supported the researchenvironment grant “Gravitational Radiation and Elec- tromagnetic Astrophysical Transients (GREAT)” fundedby the Swedish Research Council (VR) under Dnr 2016-06012 and the research project grant “Gravity MeetsLight” funded by the Knut and Alice Wallenberg Foun-dation Dnr KAW 2019.0112. SMN is grateful for fi-nancial support from the Nederlandse Organisatie voorWetenschappelijk Onderzoek (NWO) through the VIDIand Projectruimte grants. HVP and DJM acknowledgethe hospitality of the Aspen Center for Physics, which issupported by National Science Foundation grant PHY-1607611. The participation of HVP and DJM at theAspen Center for Physics was supported by the SimonsFoundation. This work used computing facilities pro-vided by the UCL Cosmoparticle Initiative; and we thankthe HPC systems manager Edd Edmondson for his dedi-cated support. We are deeply grateful to the
IMRPhenom and
SEOBNR waveform modelers for making these wave-forms public, without which this work would not havebeen possible.
Author contributions.
SMF : conceptualization,methodology, software, investigation, validation, writ-ing (original draft preparation).
HVP : conceptualiza-tion, methodology, validation, writing (review & editing),funding acquisition.
SMN : conceptualization, method-ology, validation, writing (review & editing).
DJM : con-ceptualization, methodology, writing (review & editing).
Appendix A: Supplemental Material
Waveform Approximants.
The waveform approx-imants used in this work are based on the two mainmethods that model fully the inspiral-merger-ringdownphases of binary black holes. Both models extend theanalytically–derived inspiral to incorporate the mergerand ringdown by calibrating the full waveforms with asuite of numerical relativity simulations. With thesebaseline binary black hole models in hand, the two ap-proaches treat the finite-size effects of neutron stars byadding tidal deformation, using BNS numerical relativ-ity simulations, and spin-quadrupole corrections to theGW phase evolution in the waveforms. Specifically,
IMRPhenom uses a phenomenological approach when ex-tending the post-Newtonian approximation of the inspi-ral using numerical relativity simulations, which encap-sulates the dominant effects of generic precessing binaryblack holes [e.g., 67]. Importantly,
IMRPhenom mod-els the merger frequency of the system with a taperfor BNS systems only. In contrast,
SEOBNR is basedon the Effective-One-Body (EOB) formalism to modelthe two body problem in General Relativity [e.g., 68].
SEOBNR specifically applies to NSBH systems where theBH spin is aligned or anti-aligned to the orbital angularmomentum plane. In addition to the finite-size-inducedGW phase correction,
SEOBNR automatically describesthe merger frequency for mass ratios of 1–5 using nu-merical relativity simulations of NSBHs. An updated,NSBH-specific waveform,
IMRPhenomNSBH , was releasedin Ref. [69]; however, we do not use it here as it is onlyvalid for BH spin magnitudes up to 0.5.
Cosmography.
We assume a cosmographic distance-redshift relation based on Taylor expanding the relevantquantities to third order in expansion parameters, i.e.,the Hubble constant, H , the deceleration parameter, q ,and the jerk, j (which we assume to be one throughout).The luminosity distance is then related to the expansionparameters by [57] d (cid:39) czH (cid:20) − q ) z + 16 (cid:0) − q − j + 3 q (cid:1) z (cid:21) , (A1)and the comoving volume element byd V d z (cid:39) π c z H × (A2) (cid:20) − q ) z + 512 (cid:0) q − j + 9 q (cid:1) z (cid:21) . The redshifted volume is defined to be V = 4 π (cid:90) z max d V d z d z z . (A3) Posterior.
Our two-stage Bayesian inference pipeline requires us to first evaluate, for each merger in turn, theGW likelihood marginalized over all parameters θ i otherthan the luminosity distance d to the merger:P( ˆ x i | d [ z i , H , q ]) = (cid:90) d θ i P( θ i )P( ˆ x i | d [ z i , H , q ] , θ i ) . (A4)ˆ x i here are the the i th merger’s GW strain observations,and θ i comprises the merger’s component masses, spinmagnitudes and orientations (if using precessing spins),inclination, polarization angle, NS tidal deformability,and time and phase at coalescence. With these in hand,the joint posterior of the cosmological and NSBH param-eters is given byP (cid:0) H , q , Γ , { z, v }| N, { ˆ x , ˆ z, ˆ v } , ρ ∗ , m ∗ ej (cid:1) ∝ (A5)P ( H , q , Γ) exp (cid:0) − ¯ N (cid:2) H , q , Γ , ρ ∗ , m ∗ ej (cid:3)(cid:1) × N (cid:89) i =1 Γ1 + z i d V d z [ H , q ] P( ˆ x i | d [ z i , H , q ])P( v i )P(ˆ v i | v i )P(ˆ z i | z i ) , where N and ¯ N are the actual and expected number ofmergers detected, respectively, curly brackets denote setsof quantities and bold denotes per-merger vectors. Thespecific posteriors we obtain after processing our simu-lated catalogs are shown in Fig. 3. [1] A. G. Riess, S. Casertano, W. Yuan, L. M. Macri, andD. Scolnic, Large Magellanic Cloud Cepheid StandardsProvide a 1% Foundation for the Determination of theHubble Constant and Stronger Evidence for Physics be-yond ΛCDM, Astrophys. J. , 85 (2019).[2] Planck Collaboration, Planck 2018 results. VI. Cosmolog-ical parameters, arXiv e-prints (2018), arXiv:1807.06209[astro-ph.CO].[3] L. Knox and M. Millea, Hubble constant hunter’s guide,Phys. Rev. D , 043533 (2020), arXiv:1908.03663[astro-ph.CO].[4] S. Vagnozzi, New physics in light of the H tension:An alternative view, Phys. Rev. D , 023518 (2020),arXiv:1907.07569 [astro-ph.CO].[5] G. Efstathiou, H revisited, MNRAS , 1138 (2014);M. Rigault et al. , Confirmation of a Star Formation Biasin Type Ia Supernova Distances and its Effect on theMeasurement of the Hubble Constant, Astrophys. J. ,20 (2015); D. O. Jones, A. G. Riess, and D. M. Scol-nic, Reconsidering the Effects of Local Star Formationon Type Ia Supernova Cosmology, Astrophys. J. , 31(2015); W. Cardona, M. Kunz, and V. Pettorino, Deter-mining H with Bayesian hyper-parameters, JCAP , 056(2017); B. R. Zhang, M. J. Childress, T. M. Davis, N. V.Karpenka, C. Lidman, B. P. Schmidt, and M. Smith,A blinded determination of H from low-redshift TypeIa supernovae, calibrated by Cepheid variables, MNRAS , 2254 (2017); B. Follin and L. Knox, Insensitiv-ity of the distance ladder Hubble constant determinationto Cepheid calibration modelling choices, MNRAS ,4534 (2018); S. M. Feeney, D. J. Mortlock, and N. Dal- masso, Clarifying the Hubble constant tension with aBayesian hierarchical model of the local distance ladder,MNRAS , 3861 (2018); H.-Y. Wu and D. Huterer,Sample variance in the local measurements of the Hub-ble constant, MNRAS , 4946 (2017); S. Dhawan,S. W. Jha, and B. Leibundgut, Measuring the Hubbleconstant with Type Ia supernovae as near-infrared stan-dard candles, A&A , A72 (2018); C. A. P. Ben-galy, U. Andrade, and J. S. Alcaniz, How does an incom-plete sky coverage affect the Hubble Constant variance?,ArXiv e-prints (2018), arXiv:1810.04966; M. Rigault et al. , Strong Dependence of Type Ia Supernova Stan-dardization on the Local Specific Star Formation Rate, ibid . (2018), arXiv:1806.03849; D. O. Jones, A. G. Riess,D. M. Scolnic, Y. C. Pan, E. Johnson, D. A. Coul-ter, K. G. Dettman, M. M. Foley, R. J. Foley, M. E.Huber, S. W. Jha, C. D. Kilpatrick, R. P. Kirshner,A. Rest, A. S. B. Schultz, and M. R. Siebert, ShouldType Ia Supernova Distances Be Corrected for TheirLocal Environments?, Astrophys. J. , 108 (2018),arXiv:1805.05911 [astro-ph.CO]; A. G. Riess, W. Yuan,S. Casertano, L. M. Macri, and D. Scolnic, The Ac-curacy of the Hubble Constant Measurement Verifiedthrough Cepheid Amplitudes, ApJL , L43 (2020),arXiv:2005.02445 [astro-ph.CO]; G. Efstathiou, A Lock-down Perspective on the Hubble Tension (with commentsfrom the SH0ES team), arXiv e-prints , arXiv:2007.10716(2020), arXiv:2007.10716 [astro-ph.CO].[6] D. N. Spergel, R. Flauger, and R. Hloˇzek, Planck datareconsidered, Phys. Rev. D , 023518 (2015); G. E. Ad-dison, Y. Huang, D. J. Watts, C. L. Bennett, M. Halpern,
66 70 74 H (cid:0) km s − Mpc − (cid:1) Γ (cid:0) y r − G p c − (cid:1) − q H = 68 . ± . − q q = − . ± .
32 600 800 Γ (cid:0) yr − Gpc − (cid:1) Γ = 619 +65 −
64 66 68 H (cid:0) km s − Mpc − (cid:1) Γ (cid:0) y r − G p c − (cid:1) − q H = 66 . ± . − q q = − . ± .
27 600 800 Γ (cid:0) yr − Gpc − (cid:1) Γ = 607 +56 − FIG. 3. Cosmological and population parameter posteriors inferred for the mock
SEOBNR (left) and
IMRPhenom (right) catalogs.G. Hinshaw, and J. L. Weiland, Quantifying Discordancein the 2015 Planck CMB Spectrum, Astrophys. J. ,132 (2016); G. Obied, C. Dvorkin, C. Heinrich, W. Hu,and V. Miranda, Inflationary features and shifts in cos-mological parameters from Planck 2015 data, Phys. Rev.D , 083526 (2017), arXiv:1706.09412 [astro-ph.CO];E. Calabrese, R. A. Hloˇzek, J. R. Bond, M. J. De-vlin, J. Dunkley, M. Halpern, A. D. Hincks, K. D. Ir-win, A. Kosowsky, K. Moodley, L. B. Newburgh, M. D.Niemack, L. A. Page, B. D. Sherwin, J. L. Sievers,D. N. Spergel, S. T. Staggs, and E. J. Wollack, Cos-mological parameters from pre-Planck CMB measure-ments: A 2017 update, Phys. Rev. D , 063525 (2017),arXiv:1702.03272 [astro-ph.CO]; G. Efstathiou andS. Gratton, A Detailed Description of the CamSpec Like-lihood Pipeline and a Reanalysis of the Planck High Fre-quency Maps, arXiv e-prints , arXiv:1910.00483 (2019),arXiv:1910.00483 [astro-ph.CO]; P. Motloch and W. Hu,Lensing-like tensions in the Planck legacy release, Phys.Rev. D , 083515 (2020), arXiv:1912.06601 [astro-ph.CO]; S. Aiola et al. , The Atacama Cosmology Tele-scope: DR4 Maps and Cosmological Parameters, arXive-prints , arXiv:2007.07288 (2020), arXiv:2007.07288[astro-ph.CO].[7] G. E. Addison, D. J. Watts, C. L. Bennett, M. Halpern,G. Hinshaw, and J. L. Weiland, Elucidating ΛCDM: Im-pact of Baryon Acoustic Oscillation Measurements on theHubble Constant Discrepancy, ArXiv e-prints (2017),arXiv:1707.06547; DES Collaboration, Dark Energy Sur-vey Year 1 Results: A Precise H0 Measurement from DESY1, BAO, and D/H Data, ibid . (2017), arXiv:1711.00403;O. H. E. Philcox, M. M. Ivanov, M. Simonovi´c, andM. Zaldarriaga, Combining full-shape and BAO analysesof galaxy power spectra: a 1.6% CMB-independent con-straint on H , JCAP , 032 (2020), arXiv:2002.04035[astro-ph.CO].[8] W. Yuan, A. G. Riess, L. M. Macri, S. Casertano, and D. M. Scolnic, Consistent Calibration of the Tip of theRed Giant Branch in the Large Magellanic Cloud on theHubble Space Telescope Photometric System and a Rede-termination of the Hubble Constant, Astrophys. J. ,61 (2019), arXiv:1908.00993 [astro-ph.GA]; C. D. Huang,A. G. Riess, W. Yuan, L. M. Macri, N. L. Zakamska,S. Casertano, P. A. Whitelock, S. L. Hoffmann, A. V.Filippenko, and D. Scolnic, Hubble Space Telescope Ob-servations of Mira Variables in the SN Ia Host NGC1559: An Alternative Candle to Measure the HubbleConstant, Astrophys. J. , 5 (2020), arXiv:1908.10883[astro-ph.CO]; K. C. Wong, S. H. Suyu, G. C. F. Chen,C. E. Rusu, M. Millon, D. Sluse, V. Bonvin, C. D. Fass-nacht, S. Taubenberger, M. W. Auger, S. Birrer, J. H. H.Chan, F. Courbin, S. Hilbert, O. Tihhonova, T. Treu,A. Agnello, X. Ding, I. Jee, E. Komatsu, A. J. Shajib,A. Sonnenfeld, R. D. Blandford, L. V. E. Koopmans, P. J.Marshall, and G. Meylan, H0LiCOW - XIII. A 2.4 percent measurement of H from lensed quasars: 5.3 σ ten-sion between early- and late-Universe p, MNRAS ,1420 (2020), arXiv:1907.04869 [astro-ph.CO]; M. Millon,A. Galan, F. Courbin, T. Treu, S. H. Suyu, X. Ding,S. Birrer, G. C. F. Chen, A. J. Shajib, D. Sluse, K. C.Wong, A. Agnello, M. W. Auger, E. J. Buckley-Geer,J. H. H. Chan, T. Collett, C. D. Fassnacht, S. Hilbert,L. V. E. Koopmans, V. Motta, S. Mukherjee, C. E.Rusu, A. Sonnenfeld, C. Spiniello, and L. Van de Vyvere,TDCOSMO. I. An exploration of systematic uncertain-ties in the inference of H from time-delay cosmography,A&A , A101 (2020), arXiv:1912.08027 [astro-ph.CO];D. W. Pesce, J. A. Braatz, M. J. Reid, A. G. Riess,D. Scolnic, J. J. Condon, F. Gao, C. Henkel, C. M. V. Im-pellizzeri, C. Y. Kuo, and K. Y. Lo, The Megamaser Cos-mology Project. XIII. Combined Hubble Constant Con-straints, ApJL , L1 (2020), arXiv:2001.09213 [astro-ph.CO].[9] W. L. Freedman, B. F. Madore, D. Hatt, T. J. Hoyt, I. S. Jang, R. L. Beaton, C. R. Burns, M. G. Lee, A. J. Mon-son, J. R. Neeley, M. M. Phillips, J. A. Rich, and M. Seib-ert, The Carnegie-Chicago Hubble Program. VIII. An In-dependent Determination of the Hubble Constant Basedon the Tip of the Red Giant Branch, Astrophys. J. ,34 (2019), arXiv:1907.05922 [astro-ph.CO]; W. L. Freed-man, B. F. Madore, T. Hoyt, I. S. Jang, R. Beaton,M. G. Lee, A. Monson, J. Neeley, and J. Rich, Cal-ibration of the Tip of the Red Giant Branch, Astro-phys. J. , 57 (2020), arXiv:2002.01550 [astro-ph.GA];S. Birrer, A. J. Shajib, A. Galan, M. Millon, T. Treu,A. Agnello, M. Auger, G. C. F. Chen, L. Christensen,T. Collett, F. Courbin, C. D. Fassnacht, L. V. E. Koop-mans, P. J. Marshall, J. W. Park, C. E. Rusu, D. Sluse,C. Spiniello, S. H. Suyu, S. Wagner-Carena, K. C. Wong,M. Barnab`e, A. S. Bolton, O. Czoske, X. Ding, J. A.Frieman, and L. Van de Vyvere, TDCOSMO IV: Hi-erarchical time-delay cosmography – joint inference ofthe Hubble constant and galaxy density profiles, arXive-prints , arXiv:2007.02941 (2020), arXiv:2007.02941[astro-ph.CO]; S. S. Boruah, M. J. Hudson, andG. Lavaux, Peculiar velocities in the local Universe:comparison of different models and the implications for H and dark matter, arXiv e-prints , arXiv:2010.01119(2020), arXiv:2010.01119 [astro-ph.CO].[10] B. F. Schutz, Determining the Hubble constant fromgravitational wave observations, Nature (London) ,310 (1986).[11] D. E. Holz and S. A. Hughes, Using Gravitational-Wave Standard Sirens, Astrophys. J. , 15 (2005),arXiv:astro-ph/0504616 [astro-ph].[12] N. Dalal, D. E. Holz, S. A. Hughes, and B. Jain, ShortGRB and binary black hole standard sirens as a probe ofdark energy, Phys. Rev. D , 063006 (2006).[13] S. Nissanke, D. E. Holz, S. A. Hughes, N. Dalal, andJ. L. Sievers, Exploring Short Gamma-ray Bursts asGravitational-wave Standard Sirens, Astrophys. J. ,496 (2010).[14] S. R. Taylor, J. R. Gair, and I. Mandel, Cosmology usingadvanced gravitational-wave detectors alone, Phys. Rev.D , 023535 (2012).[15] C. Messenger and J. Read, Measuring a CosmologicalDistance-Redshift Relationship Using Only GravitationalWave Observations of Binary Neutron Star Coalescences,Phys. Rev. Lett. , 091101 (2012).[16] S. Nissanke, D. E. Holz, N. Dalal, S. A. Hughes,J. L. Sievers, and C. M. Hirata, Determining theHubble constant from gravitational wave observationsof merging compact binaries, ArXiv e-prints (2013),arXiv:1307.2638.[17] M. Oguri, Measuring the distance-redshift relation withthe cross-correlation of gravitational wave standardsirens and galaxies, Phys. Rev. D , 083511 (2016).[18] W. Del Pozzo, T. G. F. Li, and C. Messenger, Cosmolog-ical inference using only gravitational wave observationsof binary neutron stars, Phys. Rev. D , 043502 (2017).[19] B. P. Abbott et al. , A gravitational-wave standard sirenmeasurement of the Hubble constant, Nature (London) , 85 (2017).[20] N. Seto and K. Kyutoku, Prospects of the local Hubbleparameter measurement using gravitational waves fromdouble neutron stars, MNRAS , 4133 (2018).[21] H.-Y. Chen, M. Fishbach, and D. Holz, A 2 per cent Hub-ble constant measurement from standard sirens within 5 years, Nature (London) , 545 (2018).[22] M. Fishbach, R. Gray, I. Maga˜na Hernandez, H. Qi,A. Sur, F. Acernese, L. Aiello, A. Allocca, M. A. Aloy,and A. Amato, A Standard Siren Measurement of theHubble Constant from GW170817 without the Electro-magnetic Counterpart, ApJL , L13 (2019).[23] S. M. Feeney, H. V. Peiris, A. R. Williamson, S. M.Nissanke, D. J. Mortlock, J. Alsing, and D. Scolnic,Prospects for Resolving the Hubble Constant Tensionwith Standard Sirens, Phys. Rev. Lett. , 061105(2019).[24] D. J. Mortlock, S. M. Feeney, H. V. Peiris, A. R.Williamson, and S. M. Nissanke, Unbiased Hubble con-stant estimation from binary neutron star mergers, Phys.Rev. D , 103523 (2019), arXiv:1811.11723 [astro-ph.CO].[25] R. Gray, I. Maga˜na Hernandez, H. Qi, A. Sur, P. R.Brady, H.-Y. Chen, W. M. Farr, M. Fishbach, J. R. Gair,A. Ghosh, D. E. Holz, S. Mastrogiovanni, C. Messenger,D. A. Steer, and J. Veitch, Cosmological Inference us-ing Gravitational Wave Standard Sirens: A Mock DataChallenge, arXiv e-prints (2019), arXiv:1908.06050 [gr-qc].[26] M. Soares-Santos et al. , First Measurement of the Hub-ble Constant from a Dark Standard Siren using the DarkEnergy Survey Galaxies and the LIGO/Virgo Binary-Black-hole Merger GW170814, ApJL , L7 (2019),arXiv:1901.01540 [astro-ph.CO].[27] A. Palmese et al. , A Statistical Standard Siren Mea-surement of the Hubble Constant from the LIGO/VirgoGravitational Wave Compact Object Merger GW190814and Dark Energy Survey Galaxies, ApJL , L33(2020), arXiv:2006.14961 [astro-ph.CO].[28] S. Mukherjee and B. D. Wandelt, Beyond the classi-cal distance-redshift test: cross-correlating redshift-freestandard candles and sirens with redshift surveys, arXive-prints , arXiv:1808.06615 (2018), arXiv:1808.06615[astro-ph.CO].[29] S. Mukherjee, B. D. Wandelt, S. M. Nissanke, andA. Silvestri, Accurate and precision Cosmology withredshift unknown gravitational wave sources, arXive-prints , arXiv:2007.02943 (2020), arXiv:2007.02943[astro-ph.CO].[30] S. S. Vasylyev and A. V. Filippenko, A Measurementof the Hubble Constant Using Gravitational Waves fromthe Binary Merger GW190814, Astrophys. J. , 149(2020), arXiv:2007.11148 [astro-ph.CO].[31] H.-Y. Chen, C.-J. Haster, S. Vitale, W. M. Farr, andM. Isi, A Standard Siren Cosmological Measurementfrom the Potential GW190521 Electromagnetic Counter-part ZTF19abanrhr, arXiv e-prints , arXiv:2009.14057(2020), arXiv:2009.14057 [astro-ph.CO].[32] V. Gayathri, J. Healy, J. Lange, B. O’Brien, M. Szczep-anczyk, I. Bartos, M. Campanelli, S. Klimenko,C. Lousto, and R. O’Shaughnessy, Hubble ConstantMeasurement with GW190521 as an Eccentric BlackHole Merger, arXiv e-prints , arXiv:2009.14247 (2020),arXiv:2009.14247 [astro-ph.HE].[33] S. Mukherjee, A. Ghosh, M. J. Graham, C. Karathanasis,M. M. Kasliwal, I. Maga˜na Hernandez, S. M. Nissanke,A. Silvestri, and B. D. Wandelt, First measurement ofthe Hubble parameter from bright binary black holeGW190521, arXiv e-prints , arXiv:2009.14199 (2020),arXiv:2009.14199 [astro-ph.CO]. [34] S. Vitale and H.-Y. Chen, Measuring the Hubble Con-stant with Neutron Star Black Hole Mergers, Phys. Rev.Lett. , 021303 (2018).[35] T. Dietrich, S. Khan, R. Dudi, S. J. Kapadia, P. Ku-mar, A. Nagar, F. Ohme, F. Pannarale, A. Samaj-dar, S. Bernuzzi, G. Carullo, W. Del Pozzo, M. Haney,C. Markakis, M. P¨urrer, G. Riemenschneider, Y. E.Setyawati, K. W. Tsang, and C. Van Den Broeck, Mat-ter imprints in waveform models for neutron star bina-ries: Tidal and self-spin effects, Phys. Rev. D , 024029(2019), arXiv:1804.02235 [gr-qc].[36] A. Matas, T. Dietrich, A. Buonanno, T. Hinderer,M. P¨urrer, F. Foucart, M. Boyle, M. D. Duez,L. E. Kidder, H. P. Pfeiffer, and M. A. Scheel,Aligned-spin neutron-star-black-hole waveform modelbased on the effective-one-body approach and numerical-relativity simulations, Phys. Rev. D , 043023 (2020),arXiv:2004.10001 [gr-qc].[37] F. Foucart, T. Hinderer, and S. Nissanke, Remnantbaryon mass in neutron star-black hole mergers: Pre-dictions for binary neutron star mimickers and rapidlyspinning black holes, Phys. Rev. D , 081501 (2018),arXiv:1807.00011 [astro-ph.HE].[38] B. P. Abbott et al. , Prospects for observing and localizinggravitational-wave transients with Advanced LIGO, Ad-vanced Virgo and KAGRA, Living Reviews in Relativity , 3 (2018), arXiv:1304.0670 [gr-qc].[39] The LIGO, Virgo and KAGRA Collaborations,https://dcc.ligo.org/ligo-t2000012/public.[40] The LIGO Scientific Collaboration, the Virgo Collabora-tion, et al. , GWTC-1: A Gravitational-Wave TransientCatalog of Compact Binary Mergers Observed by LIGOand Virgo during the First and Second Observing Runs,arXiv e-prints (2018), arXiv:1811.12907 [astro-ph.HE].[41] M. U. Kruckow, T. M. Tauris, N. Langer, M. Kramer,and R. G. Izzard, Progenitors of gravitational wave merg-ers: binary evolution with the stellar grid-based codeCOMBINE, MNRAS , 1908 (2018), arXiv:1801.05433[astro-ph.SR].[42] LVC, GWTC-2: Compact Binary Coalescences Observedby LIGO and Virgo During the First Half of the ThirdObserving Run, arXiv e-prints , arXiv:2010.14527 (2020),arXiv:2010.14527 [gr-qc].[43] S. Typel, G. R¨opke, T. Kl¨ahn, D. Blaschke, and H. H.Wolter, Composition and thermodynamics of nuclearmatter with light clusters, Phys. Rev. C , 015803(2010), arXiv:0908.2344 [nucl-th]; M. Hempel andJ. Schaffner-Bielich, A statistical model for a com-plete supernova equation of state, Nucl. Phys. A ,210 (2010), arXiv:0911.4073 [nucl-th]; T. Fischer,M. Hempel, I. Sagert, Y. Suwa, and J. Schaffner-Bielich,Symmetry energy impact in simulations of core-collapsesupernovae, European Physical Journal A , 46 (2014),arXiv:1307.6190 [astro-ph.HE].[44] F. Foucart, Black-hole-neutron-star mergers: Diskmass predictions, Phys. Rev. D , 124007 (2012),arXiv:1207.6304 [astro-ph.HE].[45] T. Binnington and E. Poisson, Relativistic theoryof tidal Love numbers, Phys. Rev. D , 084018(2009), arXiv:0906.1366 [gr-qc]; P. Pani, L. Gualtieri,A. Maselli, and V. Ferrari, Tidal deformations of aspinning compact object, Phys. Rev. D , 024010(2015), arXiv:1503.07365 [gr-qc]; V. Cardoso, E. Franzin,A. Maselli, P. Pani, and G. Raposo, Testing strong- field gravity with tidal Love numbers, Phys. Rev. D , 084014 (2017), arXiv:1701.01116 [gr-qc]; H. S.Chia, Tidal Deformation and Dissipation of RotatingBlack Holes, arXiv e-prints , arXiv:2010.07300 (2020),arXiv:2010.07300 [gr-qc].[46] K. K. Y. Ng, S. Vitale, A. Zimmerman, K. Chatziioan-nou, D. Gerosa, and C.-J. Haster, Gravitational-waveastrophysics with effective-spin measurements: Asym-metries and selection biases, Phys. Rev. D , 083007(2018), arXiv:1805.03046 [gr-qc].[47] H.-Y. Chen, Systematic Uncertainty of Standard Sirensfrom the Viewing Angle of Binary Neutron Star Inspirals,Phys. Rev. Lett. , 201301 (2020), arXiv:2006.02779[astro-ph.HE].[48] S. Rosswog, U. Feindt, O. Korobkin, M. R. Wu, J. Soller-man, A. Goobar, and G. Martinez-Pinedo, Detectabil-ity of compact binary merger macronovae, Classical andQuantum Gravity , 104001 (2017), arXiv:1611.09822[astro-ph.HE].[49] D. Scolnic, R. Kessler, D. Brout, P. S. Cowperthwaite,M. Soares-Santos, J. Annis, K. Herner, H. Y. Chen,M. Sako, Z. Doctor, R. E. Butler, A. Palmese, H. T.Diehl, J. Frieman, D. E. Holz, E. Berger, R. Chornock,V. A. Villar, M. Nicholl, R. Biswas, R. Hounsell, R. J.Foley, J. Metzger, A. Rest, J. Garc´ıa-Bellido, A. M¨oller,P. Nugent, T. M. C. Abbott, F. B. Abdalla, S. Al-lam, K. Bechtol, A. Benoit-L´evy, E. Bertin, D. Brooks,E. Buckley-Geer, A. Carnero Rosell, M. Carrasco Kind,J. Carretero, F. J. Castander, C. E. Cunha, C. B.D’Andrea, L. N. da Costa, C. Davis, P. Doel, A. Drlica-Wagner, T. F. Eifler, B. Flaugher, P. Fosalba, E. Gaz-tanaga, D. W. Gerdes, D. Gruen, R. A. Gruendl,J. Gschwend, G. Gutierrez, W. G. Hartley, K. Hon-scheid, D. J. James, M. W. G. Johnson, M. D. John-son, E. Krause, K. Kuehn, S. Kuhlmann, O. Lahav,T. S. Li, M. Lima, M. A. G. Maia, M. March, J. L.Marshall, F. Menanteau, R. Miquel, E. Neilsen, A. A.Plazas, E. Sanchez, V. Scarpine, M. Schubnell, I. Sevilla-Noarbe, M. Smith, R. C. Smith, F. Sobreira, E. Suchyta,M. E. C. Swanson, G. Tarle, R. C. Thomas, D. L. Tucker,A. R. Walker, and DES Collaboration, How Many Kilo-novae Can Be Found in Past, Present, and Future Sur-vey Data Sets?, ApJL , L3 (2018), arXiv:1710.05845[astro-ph.IM].[50] P. S. Cowperthwaite, V. A. Villar, D. M. Scolnic, andE. Berger, LSST Target-of-opportunity Observations ofGravitational-wave Events: Essential and Efficient, As-trophys. J. , 88 (2019), arXiv:1811.03098 [astro-ph.HE].[51] C. N. Setzer, R. Biswas, H. V. Peiris, S. Rosswog, O. Ko-robkin, R. T. Wollaeger, and LSST Dark Energy ScienceCollaboration, Serendipitous discoveries of kilonovae inthe LSST main survey: maximizing detections of sub-threshold gravitational wave events, MNRAS , 4260(2019), arXiv:1812.10492 [astro-ph.IM].[52] I. Mandel, W. M. Farr, and J. R. Gair, Extracting distri-bution parameters from multiple uncertain observationswith selection biases, MNRAS , 1086 (2019).[53] S. Vitale, D. Gerosa, W. M. Farr, and S. R. Tay-lor, Inferring the properties of a population of com-pact binaries in presence of selection effects, arXive-prints , arXiv:2007.05579 (2020), arXiv:2007.05579[astro-ph.IM].[54] W. J. Handley, M. P. Hobson, and A. N. Lasenby, poly- chord: nested sampling for cosmology., MNRAS , L61(2015), arXiv:1502.01856 [astro-ph.CO]; POLYCHORD:next-generation nested sampling, MNRAS , 4384(2015), arXiv:1506.00171 [astro-ph.IM].[55] G. Ashton, M. H¨ubner, P. D. Lasky, C. Talbot, K. Ack-ley, S. Biscoveanu, Q. Chu, A. Divakarla, P. J. Easter,B. Goncharov, F. Hernandez Vivanco, J. Harms, M. E.Lower, G. D. Meadors, D. Melchor, E. Payne, M. D.Pitkin, J. Powell, N. Sarin, R. J. E. Smith, andE. Thrane, BILBY: A User-friendly Bayesian InferenceLibrary for Gravitational-wave Astronomy, ApJS ,27 (2019).[56] Stan Development Team, PyStan: the python interfaceto stan, Version 2.19.1.1 (2019).[57] M. Visser, Jerk, snap and the cosmological equation ofstate, Classical and Quantum Gravity , 2603 (2004),gr-qc/0309109.[58] J. Schreiber, Pomegranate: fast and flexible probabilisticmodeling in python, arXiv e-prints , arXiv:1711.00137(2017), arXiv:1711.00137 [cs.AI].[59] H. Akaike, A new look at the statistical model identifica-tion, IEEE Transactions on Automatic Control , 716(1974).[60] T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S.Thorne, Spin-induced orbital precession and its modula-tion of the gravitational waveforms from merging bina-ries, Phys. Rev. D , 6274 (1994).[61] C. Cutler and ´E. E. Flanagan, Gravitational wavesfrom merging compact binaries: How accurately canone extract the binary’s parameters from the inspiralwaveform \ ?, Phys. Rev. D , 2658 (1994).[62] A. Vecchio, LISA observations of rapidly spinning mas-sive black hole binary systems, Phys. Rev. D , 042001(2004), arXiv:astro-ph/0304051 [astro-ph].[63] W. M. Farr, M. Fishbach, J. Ye, and D. E. Holz, A Fu-ture Percent-level Measurement of the Hubble Expansionat Redshift 0.8 with Advanced LIGO, ApJL , L42(2019), arXiv:1908.09084 [astro-ph.CO].[64] The LIGO Scientific Collaboration, the Virgo Collab-oration, et al. , Binary Black Hole Population Proper-ties Inferred from the First and Second Observing Runsof Advanced LIGO and Advanced Virgo, arXiv e-prints(2018), arXiv:1811.12940 [astro-ph.HE].[65] The LIGO Scientific Collaboration and the Virgo Col-laboration, Population Properties of Compact Objectsfrom the Second LIGO-Virgo Gravitational-Wave Tran-sient Catalog, arXiv e-prints , arXiv:2010.14533 (2020),arXiv:2010.14533 [astro-ph.HE]. [66] Y. Huang, C.-J. Haster, S. Vitale, V. Varma, F. Foucart,and S. Biscoveanu, Statistical and systematic uncertain-ties in extracting the source properties of neutron star- black hole binaries with gravitational waves, arXiv e-prints , arXiv:2005.11850 (2020), arXiv:2005.11850 [gr-qc].[67] P. Ajith, S. Babak, Y. Chen, M. Hewitson, B. Krish-nan, J. T. Whelan, B. Br¨ugmann, P. Diener, J. Gon-zalez, M. Hannam, S. Husa, M. Koppitz, D. Pollney,L. Rezzolla, L. Santamar´ıa, A. M. Sintes, U. Sperhake,and J. Thornburg, A phenomenological template fam-ily for black-hole coalescence waveforms, Classical andQuantum Gravity , S689 (2007), arXiv:0704.3764 [gr-qc]; P. Ajith, S. Babak, Y. Chen, M. Hewitson, B. Kr-ishnan, A. M. Sintes, J. T. Whelan, B. Br¨ugmann,P. Diener, N. Dorband, J. Gonzalez, M. Hannam,S. Husa, D. Pollney, L. Rezzolla, L. Santamar´ıa, U. Sper-hake, and J. Thornburg, Template bank for gravita-tional waveforms from coalescing binary black holes:Nonspinning binaries, Phys. Rev. D , 104017 (2008),arXiv:0710.2335 [gr-qc]; M. Hannam, P. Schmidt,A. Boh´e, L. Haegel, S. Husa, F. Ohme, G. Pratten, andM. P¨urrer, Simple Model of Complete Precessing Black-Hole-Binary Gravitational Waveforms, Phys. Rev. Lett. , 151101 (2014), arXiv:1308.3271 [gr-qc].[68] A. Buonanno and T. Damour, Effective one-body ap-proach to general relativistic two-body dynamics, Phys.Rev. D , 084006 (1999), arXiv:gr-qc/9811091 [gr-qc];Transition from inspiral to plunge in binary black holecoalescences, Phys. Rev. D , 064015 (2000), arXiv:gr-qc/0001013 [gr-qc]; A. Boh´e, L. Shao, A. Taracchini,A. Buonanno, S. Babak, I. W. Harry, I. Hinder, S. Os-sokine, M. P¨urrer, V. Raymond, T. Chu, H. Fong, P. Ku-mar, H. P. Pfeiffer, M. Boyle, D. A. Hemberger, L. E.Kidder, G. Lovelace, M. A. Scheel, and B. Szil´agyi, Im-proved effective-one-body model of spinning, nonprecess-ing binary black holes for the era of gravitational-waveastrophysics with advanced detectors, Phys. Rev. D ,044028 (2017), arXiv:1611.03703 [gr-qc]; E. Barausse andA. Buonanno, Improved effective-one-body Hamiltonianfor spinning black-hole binaries, Phys. Rev. D , 084024(2010), arXiv:0912.3517 [gr-qc].[69] J. E. Thompson, E. Fauchon-Jones, S. Khan, E. Nitoglia,F. Pannarale, T. Dietrich, and M. Hannam, Modelingthe gravitational wave signature of neutron star blackhole coalescences, Phys. Rev. D101