# Stellar binary black holes in the LISA band: a new class of standard sirens

Walter del Pozzo, Alberto Sesana, Antoine Klein

► To cite this version:

Walter del Pozzo, Alberto Sesana, Antoine Klein. Stellar binary black holes in the LISA band: a new class of standard sirens. Monthly Notices of the Royal Astronomical Society, 2018, 475 (3), pp.3485-3492. 10.1093/mnras/sty057. hal-01768103

HAL Id: hal-01768103

https://hal.science/hal-01768103v1

Submitted on 5 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Stellar binary black holes in the LISA band: a new class of standard sirens

Walter Del Pozzo,<fnref n="1" />,<fnref n="2" />* Alberto Sesana<fnref n="2" /> and Antoine Klein<fnref n="3" />,<fnref n="4" />,<fnref n="5" />

Accepted 2018 January 5. Received 2017 November 16; in original form 2017 March 6

## ABSTRACT

The recent Advanced LIGO detections of coalescing black hole binaries (BHBs) imply a large population of such systems emitting at milli-Hz frequencies, accessible to the *Laser Interferometer Space Antenna* (LISA). We show that these systems provide a new class of cosmological standard sirens. Direct *LISA* luminosity distance – *D*₁–measurements, combined with the inhomogeneous redshift – *z*–distribution of possible host galaxies provide an effective way to populate the *D*₁–*z* diagram at *z* < 0.1, thus allowing a precise local measurement of the Hubble expansion rate. To be effective, the method requires a sufficiently precise *LISA* distance determination and sky localization of a sizeable number of BHBs, which is best achieved for a six-link detector configuration. We find that, for a BHB population consistent with current fiducial LIGO rates, the Hubble constant *H*₀ can be determined at the ~5 per cent and ~2 per cent level (68 per cent confidence), assuming two and five million kilometre arm-length, respectively.

**Key words:** black hole physics – gravitational waves.

## 1 INTRODUCTION

The first direct gravitational wave (GW) observations by Advanced LIGO and Virgo (Abbott et al. 2016b,d, 2017a,b,c) opened a new era in observational astronomy. The detected events provide a formidable laboratory for astrophysics (Abbott et al. 2016e), fundamental physics (Abbott et al. 2016c) and, thanks to the optical counterpart identification for GW170817 (Abbott et al. 2017e), a first ‘distance scale ladder’-free determination of the Hubble constant $H_0$ (Abbott et al. 2017d).

Regarding massive binary black holes (BHBs), as first pointed out by Sesana (2016), the high merger rates inferred by the LIGO detections imply an unexpectedly vast population of fairly massive BHBs (MBHBs) emitting at millihertz frequencies, accessible to space-based interferometers. This is an appealing prospect for the *Laser Interferometer Space Antenna* (LISA; eLISA Consortium et al. 2013) planned by the European Space Agency (ESA), with a possible launch window within 2029–2035.

Besides observing several BHBs few years before they merge in the LIGO band, thus realizing the promise of multiband GW astronomy (Barausse, Yunes & Chamberlain 2016; Sesana 2016; Vitale 2016), many more systems will be observed as slowly inspiralling sources around 5–10 mHz, hundreds of years before their final coalescence (see also Kyutoku & Seto 2016b). This new population

might prove useful in a number of astrophysical and cosmological contexts. For example, eccentricity measurements that are possible in the *LISA* band (Nishizawa et al. 2016b) can be used to discriminate among different astrophysical formation channels (Nishizawa et al. 2016a; Breivik et al. 2016). In this paper, we investigate the prospects of using *LISA* measurements of BHBs to constrain cosmological parameters, and, in particular, to provide an independent local measurement of the Hubble expansion rate $H_0$. This is of paramount importance in light of the recent tension between local and cosmological determination of this parameter. If we write $h = H_0/(100 \text{ km s}^{-1} \text{ Mpc}^{-1})$, a recent analysis of Type Ia supernovae with improved distance calibrators yields $h = 0.732 \pm 0.017$ (1$\sigma$ error) (Riess et al. 2016), largely inconsistent with Planck measurements of $h = 0.678 \pm 0.009$ (Planck Collaboration XIII 2016). Although local measurements might suffer from cosmic variance effects, Riess et al. (2016) showed that this is unlikely to bias the measurement of $h$ at more than 1 per cent. Therefore, a completely independent local measurement, not affected by the same distance calibration systematics, might prove useful in resolving the tension.

The argument of using GW sources as ‘standard sirens’ for cosmology measurements dates back thirty years (Schutz 1986) and goes as follows. Contrary to electromagnetic (EM) probes, GW observations provide a direct measurement of the luminosity distance – *D*₁– to the source. If the source redshift can then be determined through the detection of an associated EM counterpart, one can populate the *D*₁–*z* diagram, thus extracting the cosmological parameters that provide the best fit to the observed points. This idea has been

<footnote n="1">¹Dipartimento di Fisica ‘Enrico Fermi’, Università di Pisa, I-56127 Pisa, Italy</footnote>

<footnote n="2">²School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK</footnote>

<footnote n="3">³Department of Physics and Astronomy, The University of Mississippi, University, MS 38677, USA</footnote>

<footnote n="4">⁴CENTRA, Departamento de Física, Instituto Superior Técnico, Universidade de Lisboa, Avenida Rovisco Pais 1, P-1049 Lisboa, Portugal</footnote>

<footnote n="5">⁵CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis Bd Arago, F-75014 Paris, France</footnote>

<footnote>* E-mail: walter.delpozzo@unipi.it</footnote>
extensively exploited in different contexts both for ground-based (e.g. Nissanke et al. 2010a; Del Pozzo 2012) and space-based detectors (Holz & Hughes 2005; Tamanini et al. 2016). It has been noted, however, that for specific classes of sources, distinctive EM counterparts might be difficult to identify, like in the case of MBHBs at cosmological distances, or might not be expected in first place, like in the case of extreme mass ratio inspirals (EMRIs), i.e. compact objects inspiralling into an MBH. Even in this case, however, one can still take advantage of the clustering properties of galaxies, an idea that was first proposed by MacLeod & Hogan (2008) in the context of EMRIs observations and then applied by Petiteau, Babak & Sesana (2011) to MBHBs. To illustrate the basic idea, let us assume that LISA provides perfect measurements of $D_{i,i}$ on a set of GW sources $i = 1, \dots, n$. If cosmological parameters are unknown within a given prior range, each measured $D_{i,i}$ will translate in a possible $\Delta z_i$ range for the host galaxy. Therefore, each event will have $j = 1, \dots, m_i$ possible hosts at different redshifts within $\Delta z_i$. One can then cross-correlate the N events and search for the unique set of cosmological parameters that aligns one galaxy in each $\Delta z_i$ on a single $D_{i-}z$ relation. LISA measurement errors and the presence of peculiar velocities add significant complications to the method, adding uncertainties in the determination of the parameters, but the basic concept remains unchanged.

Here, we apply the method of MacLeod & Hogan (2008) and Del Pozzo (2012) to the putative population of stellar BHBs observable by LISA. The idea of using these systems as standard sirens has independently been proposed by Kyutoku & Seto (2016a). Instead of using a statistical method, they estimated the number of BHBs that can be localized well enough so that only one galaxy falls in the GW measurement error cube. Based on analytic scaling relations for the average sky location and distance accuracy, they concluded that LISA might measure $H_0$ to a few per cent precision. This result, however, relies on average estimates of quantities that vary by orders of magnitude from one source to another. Here we show, in fact, that association of the BHB with a single galaxy host is generally possible only in a handful of cases, and statistical methods to reconstruct the $D_{i-}z$ relation are crucial.

The paper is organized as follows. In Section 2, we present the analytical framework and introduce the statistical tools needed for the analysis. In Section 3, we construct the BHB population observed by LISA and the associated distribution of galaxy counterparts, which are the necessary input to our analysis. Results on the determination of $h$ are presented in Section 4, and final remarks are left to Section 5.

## 2 THEORETICAL FRAMEWORK

In this section, we summarize the general mathematical treatment given in Del Pozzo (2012) and Del Pozzo, Li & Messenger (2017), and specialize it to the case in which the redshifts of the potential galaxy hosts are given by the sum of the cosmological redshift and the peculiar one given by the galaxy proper motion.

Consider a set of n GW observations $d \equiv d_1, \dots, d_n$ and an all sky galaxy catalogue g. Each galaxy $g_j$ in the galaxy catalogue is defined (for simplicity) in terms of its right ascension and declination coordinates on the celestial sphere $\varphi_j, \theta_j$ and its observed redshift $z_j$ which is given by the sum of the cosmological redshift and the apparent redshift due to the galaxy peculiar velocity. Under the assumption that each GW event is located in one of the galaxies in the Universe, we want to construct the posterior probability distributions for an arbitrary set of cosmological parameters $H \equiv H_0$,

$\Omega_m, \dots$. From Bayes' theorem:

$$p(H|d, I) \propto p(H|I)p(d|H, I), \quad (1)$$

where $p(H|I)$ is the prior probability distribution for the cosmological parameters $H$ and $p(d|H, I)$ is the (quasi-)likelihood function for the observations $d$, we indicate by $I$ all background information that is relevant for the problem at hand. Since each GW observation is independent, the likelihood factorizes in the product of the likelihoods for the individual events:

$$p(d|H, I) = \prod_{i=1}^{n} p(d_i|H, I). \quad (2)$$

The (quasi-)likelihood for a GW event is obtained by marginalization over its celestial coordinates $\theta, \varphi$, luminosity distance $D_l$, and redshift $z$. Define $\mathbf{x} \equiv \{D_l, z, \theta, \varphi\}$:

$$p(d_i|H, I) = \int d\mathbf{x} p(d_i|\mathbf{x}, H, I) p(\mathbf{x}|H, I). \quad (3)$$

We are going to assume that the integral over the celestial coordinates can be factorized and performed analytically.<fnref n="1" /> We are left then with:

$$p(d_i|H, I) = \int dD_l dz \, p(d_i|D_l, z, H, I) p(D_l|z, H, I) p(z|H, I). \quad (4)$$

As in Del Pozzo (2012),

$$p(D_l|z, H, I) = \delta(D_l - D_l(z, H)) \quad (5)$$

with $D_l(z, H)$ given by Hogg (1999). Under the simplifying assumption that the likelihood for the GW event is given by a Gaussian distribution in luminosity distance $D_l$ centred around some mean value $\bar{D}_L$ and with standard deviation $\sigma_D$, after integrating over $D_l$, we get

$$p(d_i|D_l(z, H), z, H, I) \propto \exp \left[ -\frac{1}{2} \left( \frac{D_l(z, H) - \bar{D}_L}{\sigma_D} \right)^2 \right]. \quad (6)$$

Finally, we need to define the prior over the GW event redshift $p(z|H, I)$. Differently from Del Pozzo (2012), we do not simply choose a sum of delta functions over the individual galaxy redshifts, but we account also for the contribution to the redshift due to the proper motion. We choose a prior of the form

$$p(z|H, I) \propto \sum_{j=1}^{K} \exp \left[ -\frac{1}{2} \left( \frac{z_j - z}{\sigma_{pv}} \right)^2 \right] \quad (7)$$

with $\sigma_{pv} = 0.015$. This value is chosen by inspection of the distribution of line of sight velocities in the Millennium simulation. We finally obtain the likelihood for a single GW event and the corresponding potential galactic hosts:

$$p(d_i|H, I) = \int_{z_{\min}}^{z_{\max}} dz \exp \left[ -\frac{1}{2} \left( \frac{D_l(z, H) - \bar{D}_L}{\sigma_D} \right)^2 \right] \times \sum_{j=1}^{K} \exp \left[ -\frac{1}{2} \left( \frac{z_j - z}{\sigma_{pv}} \right)^2 \right]. \quad (8)$$

The limits of integration $z_{\min}, z_{\max}$ are chosen as described in Section 3.3 for the host galaxies selection.

<footnote n="1">¹ The marginalization over sky coordinates corresponds to the selection of only the galaxies within the cone subtended by the measured solid angle. We find that relaxing this assumption and marginalizing over the galaxies sky position numerically does not affect the results presented in Section 4.</footnote>
## 2.1 Numerical implementation

We explore equations (2) and (8) numerically using a parallel implementation of a generic Nested sampling algorithm<fnref n="2" /> that has been specifically adapted by the authors to the problem at hand. For each of the scenarios presented below, we compute posterior probability distributions for the reduced Hubble constant $h = H_0/100 \text{ km s}^{-1} \text{ Mpc}^{-1}$ and for the matter density $\Omega_m$. In addition to these parameters, because of the peculiar velocities of the galaxies in the Millennium simulation, we also compute posterior probability functions for each GW event included in the analysis. We choose uniform prior distributions in [0.1, 1.5] for $h$ and [0.04, 2.0] for $\Omega_m$. The prior for each GW event redshift is chosen to be uniform between [$z_{min}, z_{max}$]. We explore the parameter space using an ensemble of 1000 Live points. At every iteration of the Nested Sampling algorithm, we generate a new Live point evolving the worst likelihood one via a Markov chain Monte Carlo (MCMC) (see Veitch et al. 2015). We limit the maximum number of MCMC step to 5000. The posterior distribution for $h$ and $\Omega_m$ is constructed from the individual events marginal likelihoods in the same way described in Del Pozzo et al. (2017): Once samples from the marginal likelihood are available, we use a Dirichlet Process Gaussian Mixture model (DPGMM) (Ferguson 1973; Hjort et al. 2010) to obtain a density estimate of the likelihood density. The joint distribution is then computed by straightforward multiplication.

# 3 ASTROPHYSICAL MODEL

The practical application of the formalism developed in the previous section requires the numerical evaluation of the likelihood given by equation (8). This requires the specification of several astrophysical inputs:

(i) a population of GW sources defining the vector of observed events $d_1, \dots, d_n$;

(ii) an estimate of the errors in the determination of the source sky location and luminosity distance, which define the redshift support of possible galaxy counterparts in equation (8); and

(iii) the ensemble of galaxies $g_{i,1}, \dots, g_{i,m}$ populating the error box of each GW event $d_i$, and their apparent redshifts.

In the following, we describe how we construct each relevant quantity.

## 3.1 Black hole binary population

To generate a cosmic population of inspiralling BHBs in the milli-Hz window, we follow the procedure described in Sesana (2016). For the given BHB mass function, one can compute the associated chirp mass probability distribution $p(M_r)$, where $M_r = (M_{1,r} M_{2,r})^{3/5}/(M_{1,r} + M_{2,r})^{1/5}$. In the following, the subscript $r$ refers to quantities computed in the BHB rest frame. Those are related to quantities measured on Earth as $M_r = M/(1+z)$, $f_r = f(1+z)$, $t_r = t/(1+z)$. The merger rate density per unit mass can be then written as

$$ \frac{d^2 n}{d M_r dt_r} = \mathcal{R} \times p(M_r), \qquad (9) $$

where $\mathcal{R}$ is the total rate units of $\text{yr}^{-1} \text{ Gpc}^{-3}$. The number of BHBs per unit mass, redshift, and frequency is therefore given by

$$ \frac{d^3 N}{d M_r dz d f_r} = \frac{d^2 n}{d M_r dt_r} \frac{dV}{dz} \frac{dt_r}{d f_r}, \qquad (10) $$

where $dt_r/df_r$ is the binary frequency evolution implied by the standard quadrupole emission formula (we assume circular binaries), and the volume shell $dV/dz$ is given by a fiducial $\Lambda$ cold dark matter (CDM) cosmology with $h = 0.73$, $\Omega_m = 0.25$, $\Omega_\Lambda = 0.75$. Although those parameters are not up to date, this is irrelevant for the purpose of our exercise. The choice has been made to be consistent with the simulated universe employed to select galaxy counterparts, as described in Section 3.3.

A Monte Carlo sampling of the differential distribution given by equation (10) in the relevant frequency ($10^4 \text{Hz} < f < 0.1 \text{Hz}$) and redshift ($z < 0.5$) range provides the BHB population observable in the LISA band. All we need to specify is the intrinsic rate $\mathcal{R}$ and the BHB mass function, defining $p(M_r)$. LIGO observations imply a cosmic BHB merger rate in the range $\mathcal{R} = 9-240 \text{yr}^{-1} \text{ Gpc}^{-3}$ (90 per cent confidence Abbott et al. 2016a,f). With different assumptions in the underlying BHB mass function leading to different rate estimates. For the results shown in Section 4, we considered a log-flat rest-frame BHB mass function in the range $5 < M_{1,2,r} < 100 M_\odot$ and three values of the intrinsic merger rate $\mathcal{R} \approx 12, 34$, and $70 \text{yr}^{-1} \text{ Gpc}^{-3}$. Those roughly correspond to the 10, 50, and 90th percentile of the rate probability distribution for this mass model (see fig. 10 in Abbott et al., 2016a), thus providing a conservative, fiducial, and optimistic scenario. We stress that, as shown in Sesana (2016), the number of BHBs in the LISA band implied by the LIGO results is almost independent on the details of the BHB mass function. For example, a Salpeter mass function, although resulting in an higher estimate of $\mathcal{R}$, would essentially give the same number of detections in the LISA band. Therefore, our results do not strongly depend on the specific assumptions of a log-flat mass function.

## 3.2 LISA observations

We consider two LISA detector configurations, namely N2A2M5L6 and N2A5M5L6, in the nomenclature of Klein et al. (2016). Those correspond to detector baselines with low-frequency noise set to the original LISA requirement (N2), two and five million kilometre arm-length (A2 and A5, respectively), a nominal 5-yr mission lifetime (M5), and six laser links (L6). The first assumption is based on the reported performance of the LISA Pathfinder spacecraft (Armano et al. 2016), whereas the remaining options are set following the recommendation of the Gravitational Observatory Advisory Team to the ESA.<fnref n="3" /> Since the two configurations only differ by their arm-length, in the following we will refer to them simply as A2 and A5, respectively.

Individual signal-to-noise ratios (S/N) and estimates of the relevant parameters are produced using the Fisher information matrix (FIM) approximation, and the analytical sensitivity curves provided by Klein et al. (2016). The GW signals are modelled using 3.5 Post Newtonian precessing spinning waveforms for circular adiabatic inspirals based on the shifting asymptotic method (SUA) described in Klein, Cornish & Yunes (2014). For a given detector configuration and BHB population model, we therefore generate a catalogue of $i = 1, \dots, N$ sources observed with S/N > 8. Each observed system

<footnote n="2">² https://github.com/johnveitch/cpnest</footnote>

<footnote n="3">³ http://www.cosmos.esa.int/web/goat</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Properties of the average BHB population observed by LISA assuming the A5 (green) and the A2 (orange) baselines. Top left-hand panel: distribution of luminosity distance errors for all sources with S/N > 8. Top right-hand panel: redshift distribution of the subset of sources with $\Delta D_1/D_1 < 0.2$. Bottom left-hand panel: sky localization error distribution on the subset of sources at $z < 0.1$. Bottom right-hand panel: distribution of the number of possible galaxy hosts in the error-box defined by LISA measurement errors and the assumed prior range in the cosmological parameters. Numbers always refer to the integral of the distributions plotted in each panel.</figcaption></figure>

is defined by a vector of parameters $\lambda_i$ (including masses, spins, sky location, luminosity distance, etc.) with associated errors $\Delta\lambda_i$ estimated in the FIM approximation.

For the purpose of host identification, errors on distance $\Delta D_1$ and sky localization $\Delta\Omega$ are of paramount importance. The latter is defined as $\Delta\Omega = 2\pi \sin\theta \sqrt{\Delta\theta\Delta\phi - (\Sigma^{\theta\phi})^2}$, where $\theta$ and $\phi$ are the latitude and longitude of the source in ecliptic coordinates and $\Sigma^{\theta\phi}$ in the correlation between those two quantities. For the reconstruction of the $D_1-z$ diagram, only events with a decent $D_1$ estimate are useful. We therefore select systems with $\Delta D_1/D_1 < 0.2$. Furthermore, since we are primarily interested in an independent measurement of $h$, we restrict our analysis to events at a distance corresponding to $z < 0.1$ in the chosen cosmology. We checked, in any case, that only few events at $z > 0.1$ have $\Delta D_1/D_1 < 0.2$ and should not significantly improve on the measurement of $h$.

The cuts in $\Delta D_1/D_1$ and $z$ are shown in Fig. 1, where we show parameter estimation results averaged over 20 Monte Carlo realizations of the BHB population with $\mathcal{R}$ sampled from the posterior distribution derived by the LIGO observations. The top right-hand panel shows the $\Delta D_1/D_1$ distribution. Because of the relatively low S/N and of the inclination-distance degeneracy, distance is essentially undetermined for about 50 per cent of the sources, and only about 11 per cent of them (170/1450 and 22/200 in the A5 and A2 case, respectively) satisfy the requirement $\Delta D_1/D_1 < 0.2$. Most of those sources are at $z < 0.1$ (top left-hand panel), and they can be localized with an accuracy of about 1 deg$^2$ (bottom right-hand panel). The average number of sources that can be exploited for cosmological inference is therefore about 130 for the A5 configuration and only 20 for the A2 one.

### 3.3 Host galaxy selection

For a given cosmology, a LISA measurement $D_1 \pm \Delta D_1$ translates into a specific $z \pm \Delta z$. This conversion is complicated by two facts. First, since our goal is to determine the cosmological parameters, we need to increase our $\Delta z$ estimate to account for the prior range in those parameters. In practice, the minimum and maximum redshifts, $z^-$, $z^+$, compatible with a given distance measurement and some cosmological prior are implicitly given by

$$D_1 - \Delta D_1 = (1 + z^-)D_0 \int_0^{z^-} \frac{dz'}{F^-(z')}, \quad (11)$$

$$D_1 + \Delta D_1 = (1 + z^+)D_0 \int_0^{z^+} \frac{dz'}{F^+(z')}, \quad (12)$$

where $D_0 = 9.26 \times 10^{27}$ cm and $F(z) = hE(z)$, with $E(z)$ defined in Hogg (1999). $F^-(z)$ and $F^+(z)$ are the specific realizations of the function $F(z)$ that minimize and maximize the $D_1-z$ conversion within the assumed cosmological prior for the parameters $h$, $\Omega_m$, $\Omega_\Lambda$. As specified above, in this investigation we assume very conservative prior $h \in [0.6, 0.85]$ and $\Omega_m \in [0.04, 1]$ (fixing $\Omega_m + \Omega_\Lambda = 1$). Therefore, $F^-(z) = 0.6[0.04(1+z)^3 + 0.96]^{1/2}$ and $F^+(z) = 0.85(1+z)^{3/2}$. Secondly, because of peculiar velocities $v_p$, the apparent redshift $z_{\text{obs}}$ of the GW host might be different from the cosmological redshift $z$ determined by the Hubble flow, thus adding a further uncertainty term $\Delta z_{v_p} = (1+z)v_p/c$ to the estimated redshift. To account for this effect, for each event defined by $D_1 \pm \Delta D_1$, we assign a possible host redshift range $[z^- - \Delta z_{v_p}^-, z^+ + \Delta z_{v_p}^+]$. We assume $v_p = 500 \text{ km s}^{-1}$, consistent with the standard deviation of the radial peculiar velocity distribution observed in the Millennium run (Springel et al. 2005). We neglect additional contributions to the apparent redshift of the GW event coming from kicks at the compact binary formation, which can be as high as 360 km s$^{-1}$ (Kelley et al. 2010), thus comparable with our choice of the peculiar velocity.

The support $\Theta$ that enters in the evaluation of the galaxy likelihood function in equation (8) is the error-box defined by $\Delta\Omega \times [z^- - \Delta z_{v_p}^-, z^+ + \Delta z_{v_p}^+]$. To perform the proposed statistical $D_1-z$ reconstruction, we now need to attach a realistic distribution of galaxy to each of those error-boxes. We use the simulated observed light cones of Henriques et al. (2012) constructed from the Millennium run. The simulation assumes a $\Lambda$ CDM cosmology with $h = 0.73$, $\Omega_m = 0.25$, $\Omega_\Lambda = 0.75$ (which justifies the cosmological parameters assumed above to construct the BHB population). Because of the resolution of the simulation, galaxy distributions extracted from it are reliable only down to $M_* = 10^9 M_{\odot}$. We checked that at the magnitude cut imposed by Henriques et al. (2012), their reconstruction of the all-sky galaxy map is complete down to this mass limit out to $z \approx 0.15$. The all-sky galaxy map provides sky location, stellar mass, cosmological redshift $z$, and observed redshift $z_{\text{obs}}$ (i.e. inclusive of peculiar velocities) for each galaxy. To select galaxy counterpart candidates to each GW event, we proceed as follows:

(i) we fix the cosmological parameters to the values of the Millennium run and we select all galaxies with a cosmological redshift (i.e. excluding peculiar velocities) in the volume shell $z \pm \Delta z$ consistent with the LISA measurement $D_1 \pm \Delta D_1$;

(ii) we randomly choose one of those galaxies to be the 'true host' of the GW signal. Note that this amounts to say that the true host can be anywhere in the interval $D_1 \pm \Delta D_1$, i.e. we are not centring the error-box around the true host;
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Example of error-box construction and host selection. The red vertical line is the 'true host' redshift. The dashed vertical green line and the shaded green area are the best LISA estimated redshift and the 1σ confidence interval, when the cosmology is fixed to the Millennium run values (i.e. to the true cosmology, in our experiment). The yellow area represents the full allowed redshift interval once prior uncertainty on the cosmological parameters and peculiar velocities are also taken into account, i.e. $[z^{-} - \Delta z_{vp}^{-}, z^{+} + \Delta z_{vp}^{+}]$ as defined in the text. The black histogram is the redshift distribution of all possible galaxies falling in the error-box defined by $\Delta\Omega \times [z^{-} - \Delta z_{vp}^{-}, z^{+} + \Delta z_{vp}^{+}]$.</figcaption></figure>

(iii) we consider the error-box defined by $\Delta\Omega \times [z^{-} - \Delta z_{vp}^{-}, z^{+} + \Delta z_{vp}^{+}]$, where $\Delta\Omega$ is centred at the sky location of the selected 'true host'; and

(iv) we list all the galaxies with $M_{*} > 10^{9} M_{\odot}$ that fall within the error-box. In doing this, we make sure that their observed redshift $z_{obs}$ lies within the selected range. Each galaxy $z_{obs}$ is the relevant $z_{j}$ that enters in equation (8). In fact, in a real cosmological measurement one can only measure $z_{obs}$ (which is affected by peculiar velocities) and does not have access to the true cosmological redshift.

The outcome of this procedure for a typical GW event is shown in Fig. 2. Note that the procedure (and the statistical analysis described in Section 2) assigns to each galaxy equal probability of hosting the GW event. An alternative reasonable assumption would be that the probability of hosting the event is proportional to the stellar mass of the galaxy. This will likely sharpen the features of the likelihood function, leading to a better estimate of the cosmological parameters. In this respect, our procedure can be considered conservative and we defer the investigation of different assumptions to future work.

In summary, for each LISA configuration and BHB population model, we construct a catalogue of d LISA events with $z < 0.1$, observed with S/N > 8 and a luminosity distance error $\Delta D_{1}/D_{1} < 0.2$. Each BHB is assigned an array of possible galaxy hosts $j = 1, \dots, m$ with $z_{j,obs}$ consistent with the measured LISA uncertainty and the prior on the cosmological parameters. The distribution of number of possible hosts in the associated error-box defined by the GW measurement errors plus the prior in the cosmological parameters is shown in the bottom right-hand panel of Fig. 1. Despite the small redshifts, the poor distance determination and the wide cosmology prior cause several galaxies to fall within the error-box $\Delta\Omega \times [z^{-} - \Delta z_{vp}^{-}, z^{+} + \Delta z_{vp}^{+}]$. Differently from estimates by Kyutoku & Seto (2016a), only few per cent of the GW events (on average 9 for A5 and 2 for A2) have a single possible galaxy host,

<table><caption>**Table 1.** For each of the instrumental configurations and coalescence rates considered, we report the number of sources observed as well as the 68 per cent credible intervals over $h$ averaged over 10 realisations of the galaxy hosts.</caption><thead><tr><th>Run</th><th>LISA design</th><th>R (yr<sup>-1</sup> Gpc<sup>-3</sup>)</th><th>N<sub>BHB</sub></th><th>h (68 per cent)</th></tr></thead><tbody><tr><td>A2<sub>10</sub></td><td>N2A2M5L6</td><td>12</td><td>7</td><td>0.716<sup>+0.052</sup><sub>-0.050</sub></td></tr><tr><td>A2<sub>50</sub></td><td>N2A2M5L6</td><td>34</td><td>22</td><td>0.734<sup>+0.037</sup><sub>-0.033</sub></td></tr><tr><td>A2<sub>90</sub></td><td>N2A2M5L6</td><td>70</td><td>39</td><td>0.726<sup>+0.026</sup><sub>-0.024</sub></td></tr><tr><td>A5<sub>10</sub></td><td>N2A5M5L6</td><td>12</td><td>55</td><td>0.730<sup>+0.021</sup><sub>-0.020</sub></td></tr><tr><td>A5<sub>50</sub></td><td>N2A5M5L6</td><td>34</td><td>143</td><td>0.728<sup>+0.013</sup><sub>-0.012</sub></td></tr><tr><td>A5<sub>90</sub></td><td>N2A5M5L6</td><td>70</td><td>259</td><td>0.731<sup>+0.010</sup><sub>-0.009</sub></td></tr></tbody></table>

the vast majority have between 10 and several hundred possible hosts, making our statistical approach necessary.

Note also that we are considering hosts with $M_{*} > 10^{9} M_{\odot}$ out to $z \approx 0.1$, and the number of galaxies within the error-box is generally less than a thousand. Contrary to the MBHB case at cosmological distances investigated by Petiteau et al. (2011), getting precise spectroscopic redshifts for all the host candidates will be a relatively straightforward task.

## 4 RESULTS AND DISCUSSION

As shown in Fig. 1, on average we expect about 20 and 130 useful BHBs in the A2 and A5 LISA baselines, respectively. We investigate here three fiducial BHB population models corresponding to the 10–50–90-th percentile of the rate distribution inferred from the first LIGO detections. The combination of detector baseline and population models results in six scenarios that we name A2$_{10}$, A2$_{50}$, A2$_{90}$, A5$_{10}$, A5$_{50}$, and A5$_{90}$ and are listed in Table 1. For each scenario, we compute posterior distributions for $h$ and for the matter density $\Omega_{m}$. For each of the cases listed in Table 1, we generated 10 realizations of the set of possible galaxy hosts, thus all the results quoted in this section are averaged over the aforementioned realizations. Call $M$ the number of realizations for each event $d_{i}$, we produce average marginal likelihoods by computing

$$ 
\begin{aligned}
p(< d_i > | h, \Omega_m, I) &= \frac{1}{M} \sum_{j=1}^{M} p(d_i^j | h, \Omega_m, I) \\
&= \frac{1}{M} \int dz \sum_{j=1}^{M} p(d_i^j | h, \Omega_m, z, I)
\end{aligned} 
\quad (13) 
$$

and $p(d_i^j | h, \Omega_m, z, I)$ for each event realization is given in equation 8. The joint posteriors for $h$ and $\Omega_m$ are then computed from the averaged single ones. Regardless of the number of events considered (thus of the rate and instrumental configuration), we find that $\Omega_m$ is never measurable and its posterior is identical to the prior. This is not surprising as all the events we considered have a redshift smaller than 0.1. Hence, we will not report results for $\Omega_m$ as uninformative.

The main result of the study is summarized in Fig. 3: LISA observations of BHBs will allow a determination of the Hubble constant $H_{0}$ with a degree of accuracy at the 68 per cent level that ranges from ~7 per cent in the worst case scenario of a 2Mkm arm-length and low merger rate to ~1.3 per cent for the 5Mkm arm-length and large merger rate. Rather than the increased S/N due to the longer baseline, we find that the inference is driven by the number of detected events. In fact, A2$_{90}$ and A5$_{10}$ yield approximately the same
<figure><img src="image_5.png" /><figcaption>**Figure 3.** 90 per cent (black) and 68 per cent (red) credible intervals for $H_0$/km s$^{-1}$ Mpc$^{-1}$ for each of the configurations considered in this study. The credible regions are the result of the average over the galaxy hosts realizations. The main factor driving the inference of $H_0$ is the number of GW events detected, therefore the best inference is guaranteed by a high merger rate and 5Mkm arm-length for LISA.</figcaption></figure>

number of detections and, consequently, very similar uncertainties on the Hubble constant, 3.5 per cent and 2.8 per cent, respectively. The optimal inference for *h* is given by the $A_{50}$. In this case, the 259 detected events yield a measurement of *h* whose 90 per cent accuracy is ~2 per cent, thus a factor of ~2 worse than the current most accurate determination by Planck Collaboration XIII (2016), but competitive with best distance indicators available today (Riess et al. 2016). Although this is an unfair comparison, since *LISA* is still far in the future and the measurements provided either by the cosmic microwave background or by large SNe Ia monitoring campaigns are most likely going to improve substantially, we would like to stress that a GW-based determination of the cosmological parameters is relying on entirely different assumptions and thus affected by independent systematics. This is, in fact, a key point; as mentioned in the introduction, some tension is emerging between CMB-based and local measurements of $H_0$. Before claiming hints of new physics, the possibility of systematics underestimation has to be ruled out. Local BHB standard sirens provide a direct benchmark comparison to SNe Ia measurements. A putative discrepancy between $H_0$ values measured with the two methods (SNe Ia and BHBs) will unveil an underestimation of the systematics at play in the measurements. On the contrary, a consistent measurement at odds with CMB-based inference will strengthen the need of going beyond the standard $\Lambda$ CDM model.

Moreover, in this paper, we are concerned exclusively with *LISA* observations of GW150914-like systems and not taking into account either the contributions coming from the analysis of *LISA* MBHBs (Petiteau et al. 2011) and EMRs (MacLeod & Hogan 2008) or systems observed by LIGO/Virgo (Nissanke et al. 2010b; Del Pozzo 2012; Nissanke et al. 2013) or *Einstein Telescope* (Sathyaprakash, Schutz & Van Den Broeck 2010; Del Pozzo, Li & Messenger 2017). Since the width of the posterior distribution for *h* scales very close to $N_{GW}^{-1/2}$, it is not unreasonable that by the end of the *LISA* mission the accuracy with which *h* can be measured using GW observations alone will be competitive with the measurements coming from more traditional channels.

Our results are moderately affected by the specific realization of the galactic hosts. For each of the cases here considered, we find

that the variance of the median value of *h* inferred is similar to the average joint posterior variance. The width of the joint posterior is, not unexpectedly, unaffected. The latter depends, in fact, only on the sizes of the distance uncertainties and not on the details of the galactic population. Examples of individual realizations of the observed standard sirens sample and on the associated $H_0$ constrains are shown in Fig. 4 for the $A_{250}$ and $A_{50}$ cases. The top panels show that the width of the joint posteriors of the single realizations is, in fact, very similar to the average values reported in Table 1, whereas the bottom panels highlight the contribution of each individual BHB to the $H_0$ measurement and the $N_{GW}^{-1/2}$ improvement in the joint $H_0$ constrain as systems are added up. The central panels provide a visual representation of how our technique works. Each BHB has an associated ‘likelihood blob’ in the $D_1-z$ plane given by the uncertainty in the GW measurement, weighted with the distribution of galaxies within this uncertainty. As we add up blobs, we can find the cosmological parameters that provide the best $D_1-z$ fit to them. Note that systems at $z < 0.05$ have significantly smaller associated uncertainties, contributing the most to the measurement. However, at such small distances, galaxy peculiar velocities can significantly bias the estimate of the intrinsic redshift of the host galaxies, and it is important to average over several events to avoid a significant bias in the final measurement of $H_0$. Moreover, similarly to what observed for EM measurements from the local Universe (e.g Riess et al. 2016), the actual measurement of $H_0$ from the *LISA* sources using the method presented in this study will have to assess potential systematics caused by cosmic variance effects that can have a non-negligible impact (of the order of 1 per cent) over the inferred value of $H_0$.

Finally, from the class of sources considered here, it is not possible to draw any inference over the value of the matter density $\Omega_m$. This is because the redshift range we considered is limited to $z \le 0.1$, over which the $D_1-z$ relation is very little affected. For an inference of the energy density parameters with *LISA* observations, we would have to turn to MBHBs that can be detected at far larger redshifts and indeed have the potential of probing even the dark energy equation of state (e.g. Petiteau et al. 2011). The lack of sensitivity to $\Omega_m$, however, greatly simplified our analysis. For searches that are sensitive to $\Omega_m$ and whose indicators have large distance uncertainties, selection effects can play a significant role (see discussion in Del Pozzo et al. 2017). Therefore, a simultaneous inference of the coalescence rate (irrelevant for the case of $H_0$) and of the cosmological parameters is necessary.

## 5 CONCLUSIONS

We considered the possibility of using GW150914-like sources as observed by *LISA* as cosmological probes. Building on the original idea of MacLeod & Hogan (2008), we employed a statistical methodology that circumvent the need of direct EM counterparts to the GW events to populate the $D_1-z$ diagram. Essentially, the method associates to each BHB a redshift probability distribution defined by the galaxies falling in an error cube given by the GW measurement uncertainty plus a prior uncertainty range for the cosmological parameters. A measurement of the latter is then possible by stacking several such distributions and searching for the combination of cosmological parameters that maximizes the resulting joint probability of the BHB hosts.

We demonstrated the feasibility of such measurement by performing end-to-end simulations including the following: (i) state of the art rates of BHB mergers as inferred by the first Advanced
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Example of the inference for $A_{2_{50}}$ (left-hand column) and $A_{5_{50}}$ (right-hand column). Top row: joint posterior distributions for $h$. The green dashed lines indicate the fiducial value $h = 0.73$. Centre row: median $D_{1-z}$ regression line (solid) and 68 per cent and 90 per cent credible regions in teal and turquoise, respectively. Each set of contours indicate the 50th, 90th percentiles from the single event posterior distribution for the redshift $z$. The dashed line indicates the $D_{1-z}$ relation from the fiducial values of the cosmological panels, $h = 0.73$ and $\Omega_m = 0.25$. The black dots are the redshifts and luminosity distances – computed using the fiducial cosmology – for the individual galaxies that went in the analysis. Bottom row: evolution of the 90 per cent credible region for $h$ as a function of the number of events analysed (black error bars). In blue, the 90 per cent credible region for $h$ from the analysis is the corresponding individual event. The green dotted line indicates the fiducial value $h = 0.73$.</figcaption></figure>

LIGO run; (ii) sky localization and distance determination of each BHB as determined by LISA observations, employing two different baselines and the full response function of the detector; and (iii) realistic distributions of potential galaxy hosts taken from the light cones derived from the Millennium simulation.

We found that LISA can probe $H_0$ with an accuracy that ranges from ~7 per cent to ~1 per cent (68 per cent confidence) depending on the intrinsic BHB merger rate and detector configuration. Except for the most pessimistic cases, figures are comparable with current local measurement of $H_0$ and can therefore play an important role
in addressing the apparent discrepancy between local and CMB based measurements of the Hubble constant. Although EM-based $H_0$ measurements are expected to improve by the time LISA flies, we stress that GW-based methods rely on independent assumptions and systematics. Moreover, LISA measurement of $H_0$ from BHB can be combined with results from other classes of GW standard sirens (MBHBs, EMRs, etc.) to significantly reduce the associated uncertainties.

While conducting this study, a novel LISA baseline has been proposed by the LISA Consortium to address the L3 mission call by ESA (Amaro-Seoane et al. 2017). The new baseline features a 2.5 Mkm arm-length and slightly different technical specifications. The resulting sensitivity curve is similar to the N2A2M5L6 configuration investigated here. The minimum proposed mission lifetime is 4 yr, with the goal of an extension to up to 10 yr. Sesana (2017) derived BHB observed rates for the new configuration, confirming that they are in line with the N2A2M5L6 baseline. If the mission is extended to 10 yr, however, the number of BHB detections increases by a factor of ≈4, implying a factor of ≈2 better determination of $H_0$. Since the determination of $H_0$ in the N2A2M5L6 case ranges from 7 per cent to 3.5 per cent, a 10-yr mission with the new proposed baseline can push these figures down to 3 per cent–1.5 per cent, depending on the intrinsic BHB merger rate.

Finally, we caution that, although quite sophisticated, our simulations rely on a number of simplifying assumptions. First, GW measurement errors are computed in the FIM approximation. Although this might not be accurate for relatively low S/N events, results from the mock LISA data challenge (Babak et al. 2010) showed that FIM results are quite reliable also at low S/N for persistent sources like EMRs or galactic WD binaries. Secondly, we only consider the 1σ GW measurement error-cube and assign to each galaxy within the cube equal probability to host the BHB. In reality, one will have to weigh the probability of each galaxy according to the shape of the 3D joint posterior distribution of the source distance and sky location. Lastly, our host probability distribution might be refined by applying more physically motivated weights for each galaxy. For example, the probability of hosting a BHB can be taken to be proportional to the galaxy mass. This would likely produce more defined spikes in the host probability distributions, with a consequential improvement on the determination of the cosmological parameters.

## ACKNOWLEDGEMENTS

We thank the referee Will Farr for useful comments that improved the manuscript. The simulations were performed on the Tsunami cluster of the University of Birmingham. AS is supported by a University Research Fellowship of the Royal Society. WDP is supported by the program 'Rientro dei Cervelli Rita Levi Montalcini'. AK is supported by NSF CAREER Grant No. PHY-1055103, and by FCT contract IF/00797/2014/CP1214/CT0012 under the IF2014 Programme. This work was supported by the Centre National d'Études Spatiales. We thank John Veitch for useful comments and discussions. Part of the graphs on this manuscript have been produced using Matplotlib (Hunter 2007).

## REFERENCES

Abbott B. P. et al., 2016a, Phys. Rev., X6, 041015

Abbott B. P. et al., 2016b, Phys. Rev. Lett., 116, 061102

Abbott B. P. et al., 2016c, Phys. Rev. Lett., 116, 221101

Abbott B. P. et al., 2016d, Phys. Rev. Lett., 116, 241103

Abbott B. P. et al., 2016e, ApJ, 818, L22

Abbott B. P. et al., 2016f, ApJ, 833, L1

Abbott B. P. et al., 2017a, Phys. Rev. Lett., 118, 221101

Abbott B. P. et al., 2017b, Phys. Rev. Lett., 119, 141101

Abbott B. P. et al., 2017c, Phys. Rev. Lett., 119, 161101

Abbott B. P. et al., 2017d, Nature, 551, 85

Abbott B. P. et al., 2017e, ApJ, 848, L12

Amaro-Seoane P. et al., 2017, preprint (arXiv:1702.00786)

Armano M. et al., 2016, Phys. Rev. Lett., 116, 231101

Babak S. et al., 2010, Class. Quantum Gravity, 27, 084009

Barausse E., Yunes N., Chamberlain K., 2016, Phys. Rev. Lett., 116, 241104

Breivik K., Rodriguez C. L., Larson S. L., Kalogera V., Rasio F. A., 2016, ApJL, 830, L18

Del Pozzo W., 2012, Phys. Rev. D, 86, 043011

Del Pozzo W., Li T. G. F., Messenger C., 2017, Phys. Rev. D, 95, 043502

eLISA Consortium et al., 2013, preprint (arXiv:1305.5720)

Ferguson T. S., 1973, Ann. Statist., 1, 209

Henriques B. M. B., White S. D. M., Lemson G., Thomas P. A., Guo Q., Marleau G.-D., Overzier R. A., 2012, MNRAS, 421, 2904

Hjort N. L., Holmes C., Müller P., Walker S. G., 2010, Bayesian Nonparametrics. Cambridge Univ. Press, Cambridge

Hogg D. W., 1999, preprint (astro-ph/9905116)

Holz D. E., Hughes S. A., 2005, ApJ, 629, 15

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Kelley L. Z., Ramirez-Ruiz E., Zemp M., Diemand J., Mandel I., 2010, ApJ, 725, L91

Klein A., Cornish N., Yunes N., 2014, Phys. Rev. D, 90, 124029

Klein A. et al., 2016, Phys. Rev. D, 93, 024003

Kyutoku K., Seto N., 2016a, Phys. Rev. D, 95, 083525

Kyutoku K., Seto N., 2016b, MNRAS, 462, 2177

MacLeod C. L., Hogan C. J., 2008, Phys. Rev. D, 77, 043512

Nishizawa A., Sesana A., Berti E., Klein A., 2016a, MNRAS, 465, 4375

Nishizawa A., Berti E., Klein A., Sesana A., 2016b, Phys. Rev. D, 94, 064020

Nissanke S., Holz D. E., Hughes S. A., Dalal N., Sievers J. L., 2010a, ApJ, 725, 496

Nissanke S., Holz D. E., Hughes S. A., Dalal N., Sievers J. L., 2010b, ApJ, 725, 496

Nissanke S., Holz D. E., Dalal N., Hughes S. A., Sievers J. L., Hirata C. M., 2013, preprint (arXiv:1307.2638)

Petiteau A., Babak S., Sesana A., 2011, ApJ, 732, 82

Planck Collaboration XIII, 2016, A&A, 594, A13

Riess A. G. et al., 2016, ApJ, 826, 56

Sathyaprakash B. S., Schutz B. F., Van Den Broeck C., 2010, Class. Quantum Gravity, 27, 215006

Schutz B. F., 1986, Nature, 323, 310

Sesana A., 2016, Phys. Rev. Lett., 116, 231102

Sesana A., 2017, J. Phys., 840, 012018

Springel V. et al., 2005, Nature, 435, 629

Tamanini N., Caprini C., Barausse E., Sesana A., Klein A., Petiteau A., 2016, J. Cosmol. Astropart. Phys., 4, 002

Veitch J. et al., 2015, Phys. Rev., D91, 042003

Vitale S., 2016, Phys. Rev. Lett., 117, 051102

This paper has been typeset from a TeX/LA TeX file prepared by the author.