Introduction: The Hubble-Lemaître constant H 0 is a fundamental cosmological quantity that governs the expansion rate of the Universe. Measurements of H 0 from different astronomical observations are at odds with each other. For instance, H 0 inferred from the fluctuation spectrum of the cosmic microwave background [1] disagrees with the value obtained from the measurement of the luminosity distance and redshift to Type Ia supernovae [2] [3] [4] at 4.0-5.8σ significance [2, 5] . Confirming or ruling out this discrepancy is of paramount importance as it may point to new or missing physics from an epoch in the early Universe just before the recombination era [5] .
Gravitational waves (GWs) facilitate a unique way of determining H 0 , without relying on the cosmic distance ladder [6, 7] . The multi-messenger observations of the binary neutron star merger GW170817 [8] led to the first GW-assisted measurement of H 0 , estimating its value to be 70 +12 −8 km s −1 Mpc −1 [9] . This measurement crucially relied on the coincident detection of an electromagnetic (EM) counterpart to the GW source. The source's redshift came from the counterpart while the luminosity distance was inferred from the GW signal. Observation of EM counterparts to ∼ 50 binary neutron star mergers could nail down the Hubble-Lemaître constant to an accuracy of 2%, sufficient to confirm if there are any systematics in the local measurements of H 0 [10] .
H 0 with Dark Sirens: Binary black hole (BBH) mergers are not expected to have EM counterparts but they too can determine the luminosity distance to their hosts independently of the cosmic distance ladder; for this reason they are sometimes referred to as dark sirens. Even so, it may be possible to identify their potential host galaxies either with the help of a galaxy catalog or by follow-up observations. There is no guarantee that this approach could identify the true host as the GW sky localization with current generation detectors is not precise enough [11] . With multiple potential hosts one has to resort to a statistical approach to determine H 0 as first suggested in [6] (also see [12] ). Dark sirens detected in the first and second observing runs of LIGO and Virgo [13] were used in this way to estimate the value of H 0 to be 68 +14 −7 km s −1 Mpc −1 [14] (also see [15] ).
One source of systematic errors in this method arises from the incompleteness of the available galaxy catalogs. Three galaxy catalogs were used in [14] : The first one is the GLADE catalog [16] , which has an all-sky coverage but the probability of the host galaxy to be in the catalog at z = 0.1, as determined in [14] , is ∼ 60%. The other two catalogs are from the Dark Energy Survey (an ongoing five year survey) [17, 18] and GWENS [19] , which are mostly complete up to z = 0.1, but do not cover the complete sky, with the former covering only an eighth of the sky at the end of its mission. Nonetheless, LIGO and Virgo at their design sensitivity could constrain H 0 with dark sirens to 5% accuracy with ∼ 250 detections [20] .
Motivation for Current Work: Recently, LIGO and Virgo have published two compact binary mergers found during the third observing run: GW190412 [21] and GW190814 [22] , which are exceptional due to their large mass-asymmetry, with mass ratios ∼ 3 and ∼ 9, respectively. These systems have led to the detection of subdominant spherical harmonic modes beyond the quadrupole mode [21] [22] [23] . Such higher modes have been argued to be of importance in the parameter estimation of asymmetric binaries [24] [25] [26] , especially that of the luminosity distance D L and orbital inclination ι [27] . Indeed, among all the dark compact binaries detected till date, GW190814 has the best measured luminosity distance (∼ 18 %) and 90% credible sky area (∼ 19 deg 2 ) [22] , albeit its signal to noise ratio (SNR) is similar to that of GW150914 whose uncertainty in distance is ∼ 35% [13] . The first LIGO-Virgo GW transient catalog GWTC-1 [13] includes ten BBHs, of which some are broadly similar to GW150914. Future observations could potentially improve parameter estimations of such systems when higher modes are present in the observed signal.
We show that the recent discoveries have raised the opportunity to measure H 0 with dark sirens to within 2%, the accuracy required to resolve the Hubble tension, in the next five years. There are two principal reasons for this expectation: (1) planned upgrades to LIGO and Virgo would enhance their sensitive volume by a factor of ∼ 3.4, and (2) higher spherical harmonic GW modes, as we found, should help localize dark sirens in the sky by a factor of ∼ 2 better but, more critically, reduce the uncertainty in distance measurement by a factor as large as 4 [27] [28] [29] . These improvements compound together to localize a 'golden' subset of the dark sirens to a small enough patch in the sky that only a single galaxy would be found within the error region. Galaxy catalogs and EM follow-up campaigns could then determine the host and obtain the source's redshift and hence directly measure H 0 , without relying on the statistical method [30, 31] .
We are primarily interested in the local measurement of H 0 with sources close enough that we can neglect the effect of weak lensing [7] as well as dark matter and dark energy and assume the simplest form of the Hubble-Lemaître law. If sources are too close (say, D L 100 Mpc), H 0 measurements will be flawed due to the systematic bias from peculiar velocities v p of host galaxies. While galaxies in clusters have relatively large peculiar velocities (v p ∼ 2000 km s −1 ), this appears not to be the case (v p ∼ 300 km s −1 ) for the majority of galaxies (∼95%) that are found outside of clusters [32] . Therefore, we will consider, in our study, dark sirens distributed uniformly in co-moving volume up to a redshift of z = 0.1 (D L 460 Mpc). This ascertains that most sources are far enough away (∼ 1% will be closer than D L ∼ 100 Mpc) that the Hubble-Lemaître flow will dominate the peculiar velocities.
Dark Siren Populations for H 0 : We consider three types of dark siren populations in our analysis and compute the precision with which H 0 could be measured with golden binaries among those populations. The first population makes use of the rates and mass distributions inferred from the GWTC-1 BBHs [13, 33] . The second population is a sub-population of these binaries, which is GW150914-like, the loudest BBH in GWTC-1 [34] . The projected merger rate and the high companion-masses make it an interesting source in itself, warranting a separate investigation (see [30, 35] for a treatment of this population, especially the accuracy with which luminosity distance and sky position could be measured, albeit without higher modes). Lastly, we consider a population of GW190814like dark sirens. The interest in this class of sources stems from the expected role higher modes may play in parameter estimation.
The simulated populations differ in the choice of companion masses. In the case of GWTC-1 BBH population, we distribute the larger mass m 1 according to a power-law p(m 1 ) ∝ m −α 1 with exponent α = 1.6 and the smaller mass m 2 uniformly in the range [5M , m 1 ] [33] . The component masses are fixed for the other two populations: m 1 = 36 M and m 2 = 29 M for GW150914-like events [35] and m 1 = 23 M and m 2 = 2.6 M for GW190814-like ones [22] . In all the cases, the companion black holes are assumed to be non-spinning, consistent with GWTC-1. Further, the events are uniformly distributed in comoving volume up to a redshift z = 0.1, as well as over sky positions, and orientation angles. Each simulated population contains 10 4 samples.
We consistently employ the IMRPhenomHM [36] model from lalsimulation [37] for the parameter estimation of all the aforementioned classes of sources. This waveform family includes radiative moments with spherical harmonic indices ( , m) = (2, 2), (3, 3) , (4, 4) , (2, 1), (3, 2) , (4, 3) , which ensures that there are no systematic biases due to the neglect of higher modes and helps us carryout meaningful comparisons between the three populations.
The median merger rate for each population reported by LIGO and Virgo is R GW190814-like = 7 +16 −6 Gpc −3 yr −1 for GW190814-like events [22] and R GWTC-1 = 53 +59 −20 Gpc −3 yr −1 for GWTC-1 BBHs [33] . Since we restrict the populations to a maximum redshift of z = 0.1, we obtain the median detection rates for the three populations to beR GW190814-like = 2.9 yr −1 ,R GTWC-1 = 22 yr −1 , and consequentlyR GW150914-like = 2.6 yr −1 for GW150914-like events, where we assume that R GW150914-like = fR GWTC-1 and f to be the fraction of GWTC-1 BBHs with m 1 , m 2 ≥ 25 M : These rates will be used below to estimate the number of dark sirens that could be localized well enough each year to identify their hosts. Detector Networks: The detector networks considered in this study are combinations of seven geographical locations and three technology generations (essentially, the choice of detector's power spectral density) as summarized in Table I . A+, AdV+, and KAGRA+ [38, 39] are planned upgrades of the Advanced LIGO [40] , Advanced Virgo [41] , and KAGRA [42] detectors, referred to as 2G+. Their targeted strain sensitivity should improve the reach for BBHs by a factor of ∼ 1.5. We consider two 2G+ networks, HLV+ and HLVKI+, with three and five detectors, respectively. The third network, Voy+, is heterogeneous and uses two 2G+ and three '2.5G' detectors; the latter is a proposed upgrade of the LIGO facilities to 'Voyager' technology [43] , which will introduce a further improvement in their reach by a factor of ∼ 3.2 compared to 2G+. The final network, ECC, contains three third generation (3G) observatories, namely one Einstein Telescope [44] in Italy and two Cosmic Explorer detectors [45] at fiducial sites in the United States and Australia. The used power spectral densities, ET-D and CE1 [39] , yield strain sensitivity improvements by factors of ∼ 3 and ∼ 5, respectively, as compared to Voyager.
Measurement Accuracies: For each of the above networks, we compute the error in the estimation of the binary parameters using the Fisher-matrix formalism [46, 47] , which is an excellent approximation for the high-SNR events such as the ones in our three populations. The parameter set includes the binary's intrinsic masses, m 1 and m 2 , two angles describing the position of the binary in the sky (right ascension α and declination δ), two more giving its orientation relative to the detector (inclination ι of the binary's orbital angular momentum relative to the line of sight and the polarization angle ψ), the luminosity distance D L , a fiducial coalescence time, and the phase of the signal at that time. Among these, D L and ι are highly correlated when the observed signal is face-on and contains only the = 2 quadrupole mode but the degeneracy is largely lifted by higher order spherical harmonic modes, > 2 [27, 48] . Gravitational waves from coalescing binaries are dominated by the quadrupole mode, however, higher modes are present in systems with unequal mass companions and more prominent for systems observed with large inclination angles. Figure 1 shows the cumulative density functions of the 90% credible sky area Ω 90 and fractional luminosity distance error ∆D L /D L for each of the three dark siren populations in the four studied networks. The GW150914-like population allows for better constrains on D L and Ω 90 due to the large SNRs such massive mergers would accumulate. The GW190814-like events will be able to determine the luminosity distance to similar accuracies as GWTC-1 BBHs, but fall behind in terms of the sky localization. The improved parameter estimation, for GW190814-like signals with relatively low SNRs, is facilitated by the higher modes that are strongly excited for such highly asymmetric systems.
These findings present the two quantities in an "eventindependent" fashion. We tackle this by applying a sky localization condition to each event. A sky patch of size Ω * 4.4 × 10 −2 deg 2 contains, on average, one L 12 galaxy within z = 0.1 [49] . Thus, from the full set of simulated events, we select the fraction * that is resolved to Ω 90 Ω * . This ensures a unique identification of the dark siren's host galaxy. Table II lists, for this sub-population of events, medians of the SNR ρ and the fractional error in the luminosity distance ∆D L /D L . Furthermore, from the merger rateR and the fraction * , we compute the number of well-localized dark sirens detected each year by the different networks. The listed event rates suggest that we can expect to observe one to several such where N is the expected number of events in two years of observing time for a network, assuming 100% duty cycle. The error estimate for some of the cases are not plotted because we do not expect to see any such event in two years of observing time.
will be rarer at only one event every 31, 6.7, and 2.4 years in HLV+, HLVKI+, and Voy+, respectively. The golden binaries described in Table II should all yield an error ≤ 5% in the luminosity distance, with the most accurate distance estimates to be expected from GW150914-like binaries. Finally, the values in Table II clearly show that dark siren localization, both in terms of luminosity distance measurement and host galaxy identification, will be the norm in the 3G network era (ECC). In fact, the distance will be determined to sub-percent accuracy no matter the source population.
Measurement of H 0 with dark sirens: The luminosity distance-redshift relation in the local Universe is well approximated by the relation D L = cz/H 0 . It follows, then, that the fractional error in H 0 is equal to the fractional error in D L for a single event, if errors in redshift measurements are negligible. In the left panel of Figure 2 , we show the error in the measurement of H 0 by a single, golden event for each of the three dark siren sub-populations that obeys the localization condition Ω 90 < Ω * . We plot the median and variance of a distribution of H 0 errors for 100 realizations of a single random event drawn from each sub-population. We find that HLV+, HLVKI+, and Voy+ detectors would estimate H 0 to a few percent accuracy while the 3G network would measure H 0 to sub-percent precision.
A major factor that contributes to the D L errors is the D L -ι degeneracy [27] . This is most significant in the absence of strongly excited higher modes and, therefore, affects the nearly equal-mass binaries in GWTC-1 and the entire GW150914-like population the most. It is precisely due to the importance of higher modes that golden GW190814-like events stay competitive to much higher-SNR signals from massive GW150914-like binaries. This holds especially true for less sensitive networks which are less effective at overcoming the D L -ι degeneracy due to lower SNRs. In fact, the variances in the distribution of H 0 errors do not favor any of the three dark siren populations in the 2G+ networks. It is not until the 3G era, when the GW150914-like events with high-SNRs will yield significantly tighter bounds on H 0 than the GW190814-like ones.
Further, we note that the variance is larger for the GWTC-1 population compared to the GW150914-like and GW190814like populations. This can be attributed to the wide mass distribution of the GWTC-1 population which ranges from 5M to 100M : lower mass binaries with relatively low SNR result in poorer luminosity distance measurements, but could still fulfill the sky localization condition. However, the GWTC-1 dark sirens stay competitive to the other two populations since they include both massive, high-SNR systems as well as very asymmetric ones.
In the right panel of Figure 2 , we calculate the errors in H 0 by taking account of the number of detections in two years of observing time for each network (assuming 100% duty cycle). The H 0 error estimates are missing for GW150914-like (GW190814-like) dark sirens in the case of HLV+ (HLV+, HLVKI+), since we do not expect to observe well-localized dark sirens of either type in the respective networks within a two year time period. We see that, due to the different rates of detections for different populations, the performance of the GW190814-like population is slightly worse (about a factor of 2) than the other two populations. Note also that the GW150914-like population is still competitive with the GWTC-1 BBHs even though its rate is considerably lower than the latter.
Conclusions: We have demonstrated a tantalizing possibility of measuring the Hubble-Lemaître constant to ∼2%-level precision using dark sirens with the imminent upgrades of the LIGO and Virgo detectors to 2G+ sensitivity. The inclusion of higher spherical harmonic modes is crucial to make such a measurement.
Our conclusions rely on controlling the amplitude calibration of the detectors to below 1%, which can be accomplished with photon calibrators [50] , and two assumptions on a dark siren's host galaxy: its unique identification and a neg-ligible uncertainty in peculiar velocity correction. The sky area could be 'contaminated' with faint galaxies or the host itself could be faint and missing from current catalogs. Further, the host's peculiar velocity correction might not meet the desired accuracy, especially for close-by sources with small Hubble-Lemaître flow. Fortunately, EM follow-up observations of such well-localized sirens should be able to identify the host galaxy, model the velocity flow, and constrain the uncertainty in peculiar velocities to ∼ 100 -150 km/s [51] , which is accurate enough for 99% of the sources considered in this study. Such follow-up surveys will be of interest to the entire astrophysics community since they would not only benefit the Hubble-Lemaître constant measurement, but improve our understanding of the correlations between binary coalescences and their environments.
Given the paucity of binary neutron star mergers with EM counterparts so far, dark sirens offer an alternative to resolve the H 0 tension within the next five years. Beyond the 2G+ era, our results are also very encouraging for a possible synergy between the dark sirens and the bright sirens, wherein the H 0 measurement from low redshift may be used as a prior in the measurement of other cosmological parameters at higher redshifts [52] .