EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 646 (February 2021) A&A, 646 (2021) A65 Full HTML
Open Access
Issue
A&A
Volume 646, February 2021
Article Number A65
Number of page(s) 11
Section Cosmology (including clusters of galaxies)
DOI https://doi.org/10.1051/0004-6361/201936724
Published online 08 February 2021

© S. Mukherjee et al. 2021

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Efforts to take an accurate measurement of the expansion rate of the Universe at the present epoch, known as the Hubble constant (H0), have been ongoing since the discovery of the expanding Universe1 by Lemaître (1927, 1931) and Hubble (1929). Several complementary approaches have measured its value with a high level of precision (Hinshaw et al. 2013; Bennett et al. 2013; Planck Collaboration XVI 2014; Anderson et al. 2014; Planck Collaboration XIII 2016; Cuesta et al. 2016; Planck Collaboration VI 2020; Riess et al. 2019). However, current measurements of H0 obtained using standard rulers anchored in the early Universe, such as cosmic microwave background (CMB), baryon acoustic oscillations (BAO, Planck Collaboration XVI 2014; Anderson et al. 2014; Aubourg et al. 2015; Planck Collaboration XIII 2016; Macaulay et al. 2019), and Big Bang nucleosynthesis (BBN, Addison et al. 2018; Abbott et al. 2018a), differ from late-Universe probes using standard candles, such as supernovae (SN type-Ia, Reid et al. 2009; Riess et al. 2019) and strong lensing from the H0LiCOW project (Wong et al. 2019), along with the use of the angular diameter distance between the lensed images as a calibrator (Jee et al. 2019). A recent measurement of H0 from the Carnegie-Chicago Hubble program, based on the use of the tip of the red giant branch (TRGB) as a calibrator for the SN type-Ia, has reduced the tension (Freedman et al. 2019). However, a more recent analysis by Yuan et al. (2019) has claimed inaccuracies in the calibration of Freedman et al. (2019), which once again aggravates the tension.

Taken at face value, this tension, statistically significant by more than 4σ, would necessitate a revision of the flat Lambda cold dark matter (ΛCDM) model of cosmology (Verde et al. 2013, 2019; Bernal et al. 2016; Di Valentino et al. 2017; Kreisch et al. 2020; Poulin et al. 2019; Lin et al. 2019; Agrawal et al. 2019; Knox & Millea 2020). Whether this discrepancy is associated with systematic or calibration errors in either of the data sets or whether it indicates new physics is currently a subject of intense debate. In this context, the spotlight has turned to standard sirens (Schutz 1986; Abbott et al. 2017a), with binary neutron star mergers as an independent probe with the potential to reach the percent-level precision needed to validate the low redshift (z) determination of H0 (Dalal et al. 2006; Nissanke et al. 2010, 2013a; Chen et al. 2018; Feeney et al. 2019; Seto & Kyutoku 2018; Mortlock et al. 2019). This potential crucially depends on whether the contamination from the peculiar velocity can be corrected at the required level of accuracy. The estimation of the cosmic velocity field was established in several earlier studies (Kaiser et al. 1991; Shaya et al. 1992; Hudson 1994; Davis et al. 1996, 2011; Branchini et al. 1999, 2001; Saunders et al. 2000; Nusser et al. 2001; Hudson et al. 2004; Radburn-Smith et al. 2004; Pike & Hudson 2005; Erdoǧdu et al. 2006; Lavaux & Hudson 2011; Ma et al. 2012; Turnbull et al. 2012).

In this paper, we propose a new framework to obtain an unbiased and accurate measurement of H0 from gravitational wave (GW) observations. Our method relies on the BORG framework, named for the Bayesian Origin Reconstruction from Galaxies (Jasche & Wandelt 2013; Jasche et al. 2015; Lavaux & Jasche 2016), to reconstruct the cosmic velocity field using the galaxy field observed in the redshift space. This framework is quite different from the methods used by Abbott et al. (2017a), which depend on the linear velocity estimates from Carrick et al. (2015). Along with the complete Bayesian posterior distribution of the velocity field available from the BORG framework, it is also useful in capturing the non-linear velocity component (as discussed in Sect. 4), which goes beyond the framework from Carrick et al. (2015). Our method is an alternative way to incorporate the peculiar velocity corrections to the future GW sources in the sky patch of the 2M++ and SDSS-III samples, as discussed in Sect. 4. To make a reliable estimation of H0 from GW sources, it is essential to correct for peculiar velocity using multiple independent approaches to minimize the chance of any systematic bias. The proposal made in this work can be included along with other methods, such as those of Carrick et al. (2015), and Springob et al. (2014) for future GW sources.

The paper is organized as follows. In Sect. 2, we discuss the low redshift probes to H0 using standard candles and standard sirens. In Sects. 3 and 4, we discuss the effect of peculiar velocity on luminosity distance and redshift and we discuss a BORG framework in the context of estimating it from the cosmological observations. In Sect. 5, we discuss the algorithm for incorporating a peculiar velocity correction to the host of standard sirens. Finally, in Sects. 6 and 7, respectively, we obtain the revised value of H0 from the event GW170817 (Abbott et al. 2019a) and discuss the applicability of our method to the future GW events.

2. Low redshift probes of Hubble constant from standard sirens

All direct low-z measurements of H0 depend on measuring the luminosity distance to the source, which is given by

d L = c ( 1 + z ) H 0 ∫ 0 z d z E ( z ) , $$ \begin{aligned} d_{\rm L}= \frac{c(1+z)}{H_0}\int _0^z \frac{\mathrm{d}z}{E(z)}, \end{aligned} $$(1)

in a homogeneous Friedmann-Lemaître-Robertson-Walker (FLRW) metric when assuming a stationary source and observer. Here, c is the speed of light, and E(z)≡H(z)/H0= Ω m ( 1 + z ) 3 + ( 1 − Ω m ) $ \sqrt{\Omega_{\mathrm{m}} (1+z)^3 + (1-\Omega_{\mathrm{m}})} $, within the framework of flat ΛCDM (Hinshaw et al. 2013; Bennett et al. 2013; Planck Collaboration XVI 2014; Anderson et al. 2014; Planck Collaboration XIII 2016; Cuesta et al. 2016; Planck Collaboration VI 2020), with Ωm= ρm/ρc (matter density ρm today divided by the critical density ρ c =3 H 0 2 /8πG $ \rho_{\rm c} = 3H_0^2/8\pi G $). At very low redshift, z the Hubble parameter H(z) is nearly constant and this relation is simplified greatly while becoming independent of the background cosmological model:

d L = cz H 0 · $$ \begin{aligned} d_{\rm L}= \frac{cz}{H_0}\cdot \end{aligned} $$(2)

GW sources provide a new avenue wherby we can measure the luminosity distance as they provide remarkable standard sirens that assume only that the general theory of relativity serves as a valid description; that is, they do not need to be “standardized” using other astrophysical sources. Indeed, as pointed out by Schutz (1986), the distance to GW sources can be measured without a degeneracy using its redshifted chirp mass ℳz, which is related to the physical chirp mass in the source frame by the relation ℳz = (1 + z)ℳ. Source frame chirp mass is related to the mass of each of the compact objects, m1 and m2, via the relation ℳ = (m1m2)3/5/(m1 + m2)1/5. The GW strain in the frequency domain can be expressed for the two polarization states, h+,  h×, via the relation (Hawking & Israel 1987; Cutler & Flanagan 1994; Poisson & Will 1995; Maggiore 2008)

h + ( f z ) = 5 96 G 5 / 6 M z 2 ( f z M z ) − 7 / 6 c 3 / 2 π 2 / 3 d L ( 1 + cos 2 ( i ) ) e i ϕ z , h × ( f z ) = 5 96 G 5 / 6 M z 2 ( f z M z ) − 7 / 6 c 3 / 2 π 2 / 3 d L cos ( i ) e i ϕ z + i π / 2 , $$ \begin{aligned}&h_+ (f_z)= \sqrt{\frac{5}{96}}\frac{G^{5/6}\mathcal{M} _z^2 (f_z\mathcal{M} _z)^{-7/6}}{c^{3/2}\pi ^{2/3}d_{\rm L}} \left(1+ \cos ^2(i)\right) e^{i\phi _z},\nonumber\\&h_\times (f_z)= \sqrt{\frac{5}{96}}\frac{G^{5/6}\mathcal{M} _z^2 (f_z\mathcal{M} _z)^{-7/6}}{c^{3/2}\pi ^{2/3}d_{\rm L}} \cos (i) e^{i\phi _z+ i\pi /2}, \end{aligned} $$(3)

where G is the gravitational constant, c is the speed of light, fz = f/(1 + z) redshifted GW frequency, ϕz is the phase of the GW signal, and i ≡ arccos ( L ̂ . n ̂ ) $ i\equiv \arccos(\hat L.\hat n) $ is the inclination angle of the GW source which is defined as the angle between the unit vector of angular momentum, L ̂ $ \hat L $, and the line of sight of the source position, n ̂ $ \hat n $. The above expression indicates the presence of a degeneracy between dL and cos(i) for individual polarization states. However, the polarization states of GW can be measured as they have different dependency on cos(i) (Holz & Hughes 2005; Dalal et al. 2006; Nissanke et al. 2010, 2013a). The measurement of both polarization states of the GW signal requires multiple detectors, with different relative orientation between the arms of the detectors.

Along with the luminosity distance measurement to the standard sirens, a measurement of the Hubble constant also requires an independent measurement of the redshift of the source. The binary neutron star, black hole-neutron star, and supermassive binary black hole mergers are all expected to have electromagnetic counterparts. This can lead to the identification of the host galaxy of the GW source and the redshift to the GW source can be estimated from the electromagnetic spectra of the host using spectroscopic (or photometric) measurement by the relation 1 + z = λo/λe2. As a result, by using dL from the GW signal, and z from the electromagnetic spectra, GW sources provide an excellent avenue for measuring the Hubble constant using Eq. (2).

3. Luminosity distance and redshift in the presence of large-scale structure

The luminosity distance to a source and its observed redshift in a homogeneous FLRW Universe is different from its counterpart in the presence of cosmic perturbations3. The presence of perturbations in the matter density leads to temporal and spatial fluctuations in the metric perturbations (through terms involving Φ ˙ , Φ ¨ $ \dot \Phi,\,\ddot \Phi $, ∇Φ, ∇2Φ), which are related to effects including (but not limited to) the peculiar velocity of the source and of the observer, the gravitational redshift, the integrated Sachs-Wolfe, and gravitational lensing (Sasaki 1987; Kolb et al. 2005; Barausse et al. 2005; Bonvin et al. 2006).

The observed redshift to the source also differs from the cosmological redshift due to the contributions from the difference in the peculiar velocity between the observer, vo, and the source, vs, and also to the gravitational redshift. At low redshift, the most important contributions arise from the difference in the velocity of the observer and source. The observed redshift zobs can be written in terms of the peculiar velocity, v p = ( v s − v o ) . u ̂ $ \mathit{v}_{\mathrm{p}}= ({\boldsymbol{v}}_{\mathrm{s}}{-}{\boldsymbol{v}}_{\mathrm{o}}).\hat u $4, via the relation ( 1 + z obs ) = ( 1 + z ) ( 1 + v p c ) $ (1+z_{\mathrm{obs}})= (1+z)\left(1+ \frac{\mathit{v}_{\mathrm{p}}}{c}\right) $. The corresponding modification in Eq. (2) is

d L = c z + v p H 0 · $$ \begin{aligned} d_{\rm L}= \frac{c z+ { v}_{\rm p}}{H_0}\cdot \end{aligned} $$(4)

This implies that the contribution from the peculiar velocity leads to a bias in the inferred value of H0, if it is not accounted for. For multiple sources with uncorrelated velocities, the effect of ignoring the peculiar velocity component in the average is an excess variance in the measurement of H0. If we assume the distribution of the peculiar velocity field to be Gaussian with a variance σ v 2 $ \sigma_{\it v}^2 $, then the corresponding excess variance in the H0 measurement for a source at distance dL becomes σ v 2 / d L 2 $ \sigma^2_{\it v}/d^2_{\rm L} $. The current framework for obtaining H0 from the standard sirens use a Gaussian prior on the peculiar velocity (Abbott et al. 2017a; Chen et al. 2018; Feeney et al. 2019). We propose an algorithm to correct for the peculiar velocity contribution to standard sirens.

The peculiar velocity of a host galaxy, vp = vh + vvir, can arise from two components, namely, (i) the motion of the halo due to the spatial gradient in the gravitational potential, vh, and (ii) the virial velocity component, vvir, of the host galaxy inside the halo. The non-rotational velocity component of vl can be obtained from linear perturbation theory as

v l ( x , z ) = − 2 3 g ( z ) a H ( z ) Ω m ( z ) ∇ x Φ , $$ \begin{aligned} {\boldsymbol{v}}_{\rm l}({\boldsymbol{x}},z)= - \frac{2}{3}\frac{g(z)}{aH(z)\,\Omega _{\rm m}(z)}\, {\boldsymbol{\nabla }}_{{x}} \Phi , \end{aligned} $$(5)

where Φ is the gravitational potential, g is the growth rate related to the linear growth factor D by the relation g = dlnD/dlna, and a = 1/(1 + z) is the scale factor. At small scales, non-linear effects can be quite important and the previous relation becomes inaccurate. The virial velocity component, vvir, can be related to the mass of the halo, Mh, and the distance to halo center, r, via the relation v vir 2 ∝ M h /r $ {\it v}^2_{\rm vir} \propto M_{\rm h}/r $, which holds for a virialized system. For a system with a constant mass density, the halo mass and r are related by Mh ∝ r3. As a result, the virial velocity vvir is only related to the mass of the halo by the relation v vir 2 ∝ M h 2 / 3 $ \mathit{v}^2_{\mathrm{vir}}\propto M^{2/3}_{\mathrm{h}} $. This relation indicates that a galaxy with a heavier halo has a larger velocity dispersion than if it resides in a smaller halo. We use the following form, fitted to simulations, to estimate the velocity dispersion of the non-linear velocity component (Sheth & Diaferio 2001):

σ vir = 476 g v ( Δ nl ( z ) E ( z ) 2 ) 1 / 6 ( M h / 10 15 M ⊙ h − 1 ) 1 / 3 , $$ \begin{aligned} \sigma _{\rm vir}= 476\, g_{ v}\, (\Delta _{\rm nl}(z)E(z)^2)^{1/6}\,(M_{\rm h}/ 10^{15}\,M_\odot \,\mathrm{h}^{-1})^{1/3}, \end{aligned} $$(6)

where gv = 0.9, Δnl(z) = 18π2 + 60x − 32x2, and x = Ωm(1 + z)3/E2(z)−1.

4. Estimation of the velocity field using BORG

The velocity field that we use in this work is produced by the Bayesian Origin Reconstruction from Galaxies (BORG) probabilistic model of large-scale structure, as currently applied to the 2M++ compilation (Lavaux & Hudson 2011; Jasche & Lavaux 2019) and the SDSS-III/BOSS survey5 (Eisenstein et al. 2011; Dawson et al. 2013; Lavaux et al. 2019). The BORG algorithm is aimed at inferring a fully probabilistic and physically plausible model of the three-dimensional cosmic matter distribution from observed galaxies in cosmological surveys (see e.g., Jasche & Wandelt 2013; Jasche et al. 2015; Lavaux & Jasche 2016). To that effect, the method solves a large-scale Bayesian inverse problem by fitting a dynamical structure formation model to data and thereby inferring the primordial initial conditions from which presently observed structures formed. The BORG forward modeling approach marginalizes automatically over unknown galaxy bias and accounts for selection and evolutionary effects. The variant of BORG that we use for the 2M++ catalog models the galaxy to dark matter bias using a three-parameter function motivated by the analysis of N-body simulations (Neyrinck et al. 2014). We note that this bias relation is a local non-linear relation for which the conventional understanding of bias on large scale does not apply. It is more akin to a halo occupation distribution for high luminosity. The data specifically selects that behavior as can be seen in Sect. 5.1 of Jasche & Lavaux (2019). We should approach with caution any naive numerical comparison of BORG bias parameters to linear bias values found in the literature. Instead of the linear perturbation result in the previous section, BORG uses a full particle-mesh N-body solver to evolve the initial conditions to a dark matter density distribution at z ≈ 0 (Jasche & Lavaux 2019) through direct integration of the Hamiltonian dynamics equation.

The BORG model provides a set of points in the parameter space (dimensionality ≃2563, for initial conditions, and a few more for additional bias parameters) that provides a numerical approximation of the posterior distribution of those parameters given the 2M++ observed galaxy distribution. Once initial and final positions of dark matter particles are known, the velocity field is estimated using the Simplex-in-Cell estimator (SIC, Hahn et al. 2015; Leclercq et al. 2017). The SIC estimator relies on a phase-space interpolation technique, which provides an accurate estimate of velocity fields without the need for commonly used but arbitrary kernel smoothing procedures, such as cloud-in-cell (CIC). It provides a physical picture for velocity fields even in low-density regions, which are poorly-sampled by dark matter particles, and takes into account the multi-valued nature of the velocity field in shell-crossed regions by proper stream averaging.

4.1. Validation of the BORG reconstruction algorithm

In this section, we discuss several validation tests of the BORG velocity field inference. In Sect. 4.1.1, we assess the validity on the sole basis of the statistics of the reconstructed samples. In Sect. 4.1.2, we run a BORG analysis on a set of mock tracers to check the unbiasedness and the typical reconstruction errors.

4.1.1. Comparisons to performance of linear perturbation theory

The model of Carrick et al. (2015) is commonly used for correcting the estimated cosmological redshift of observed objects in the nearby Universe. This model relies on computing the inverse Laplacian operator on a luminosity-weighted distribution of the 2M++ galaxies. This model has been fitted to observed distances of galaxies derived from the SFI++ (Springob et al. 2007) and the First Amendment A1 SNIa data (Turnbull et al. 2012). The BORG model goes beyond this approach. It is based on deriving a physically meaningful velocity field, beyond linear perturbation theory, using several pieces of information, including the same 2M++ data-set and cosmology.

The performance as well as the need to go beyond linear perturbations are two aspects that have been extensively discussed in the past. The differences between the density and velocity power spectra have already been computed analytically on large scales (Scoccimarro 2004) and observed in simulations (Chodorowski & Ciecielag 2002; Scoccimarro 2004; Jennings et al. 2011; Jennings 2012), but it is a challenging task to derive from observational data. This additionally motivates the BORG approach. Leclercq et al. (2015, Fig. 1) showed that BORG’s constrained simulations of the density field are statistically consistent with independent random realizations. We further emphasize this point here by analyzing the capability of going beyond the linear level for the velocity field. In that respect, we show the cross-correlation between the density and θ = −∇ ⋅ v in Fig. 1 for constrained realizations and for random simulations generated using the same particle mesh implementation and cosmological setup. There are no visible differences between the cross-correlations between the inferred 2M++ fields and the ones for random realizations. This provides evidence that the density de-correlates from the velocity divergence in the way that is expected from the non-linear gravitational clustering physics of dark matter. We also show in Fig. 2 the contribution of the different scales to the total variance of the velocity field. We note that pure Eulerian linear perturbation theory (solid blue line) is systematically above the curve obtained from constrained 2M++ realizations (shades of red) and the random velocity field power spectrum (shades of green) for the range of scales of interest. That behavior certainly changes on smaller scales for which random motions become significant, as is hinted by the crossing of the two curves at k ≃ 2 h Mpc−1. For scales of k <  0.8 h Mpc−1, the 2M++ constrained realization matches well with the pure Eulerian linear perturbation theory within the expected variance (denoted in the red shaded region). The drop of power at that scale for the mock and constrained realizations is related to the observation of non-linearities in the power spectrum as pointed out in Jennings et al. (2011).

thumbnail Fig. 1.

Display of the correlation coefficient between the non-linear density and the negative divergence of the non-linear velocity field, θ ≡ −∇ ⋅ v. The red line is the ensemble average over the posterior distribution of the cross-correlation and the shaded dark red area is the one standard deviation according to the mean of that cross-correlation derived from the same distribution obtained by BORG. The green and shaded green region is the result of computing the same quantities on 80 random realizations. Apart from very minor differences, the two curves are on top of each other, which can be expected as the physical model is the same in both cases. The vertical grey band on the right side indicate the resolution limit of the BORG reconstruction. The correlation reaches asymptotically one for the largest scales (small k values). The red line is reaching a value of 99% for the largest probed mode. As expected, this cross-correlation goes to zero at small scales (high k value).

thumbnail Fig. 2.

Ensemble average power-spectrum of the non-linear velocity field scaled, alongside with one standard deviation contour in the shaded area. The line and shaded region in red are computed from the ensemble posterior distribution built from the 2M++ data. The line and shaded region in green are measured from 80 random realizations of a ΛCDM universe of the same size. Those random realizations were computed using the same N-body solver, but starting from a random sample of a Gaussian distribution with a power-spectrum provided by the ΛCDM model. We show in thick blue line the power spectrum assuming linear theory. We note that large scales are unaffected by non-linearities. Intermediate scales (k ≳ 0.1 h Mpc−1, corresponding to ∼60 h Mpc−1) get non-negligible contributions from non-linear dynamics, as highlighted by the area under the curve. The vertical grey band on the right side indicate the resolution limit of the BORG reconstruction (∼2.6 h Mpc−1).

The two above points show the limits of using the density field as a predictor for the velocity field on intermediate scales. The linear model may likely have reached its maximum capacity of prediction. The BORG-derived non-linear model of the velocity field is physically more realistic, but the simpler model in Carrick et al. (2015) has the virtue of having been previously compared to distance observations. We intend to do the same in a forthcoming publication for BORG. In addition, BORG self-consistently models redshift-space distortions, which brings additional information with which to disentangle galaxy bias from dark matter clustering.

4.1.2. Tests with an N-body simulation

Although a full mock catalog analysis is beyond the scope of this article, we ran some tests on N-body simulations to showcase the correctness of the BORG reconstruction algorithm. Using Gadget-2 (Springel 2005), we generated an N-body simulation covering a cubic volume with a side length of 200 h Mpc−1 and sampled with 2563 particles. We used the same cosmology as that which was used for the BORG 2M++ reconstruction in Jasche & Lavaux (2019). To simplify the test case, we generated a mock catalog of 8457 objects directly from the particles of the simulation. The corresponding number density is ∼8 − 10 times smaller than in the observational data. We then ran a BORG-PM reconstruction with 20 time-steps from z = 50 to z = 0 linearly distributed according to the scale factor, with a grid of 323 elements to represent initial conditions. This corresponds to a spatial resolution of 6.25 h Mpc−1, which is two to three times lower than the one used for the run on observations.

We show in Fig. 3, a slice of the density field of the N-body simulation used for the test (top-left panel), as well as the result of the nearest-grid point assignment of the objects of the mock catalog in that same slice (top-right panel). We also show the density contrasts (bottom-left panel) computed from a single BORG sample, after the burn-in phase. We conclude with a difference plot highlighting the difference between the density field shown in the top-left panel and the bottom-left panel. We note the qualitative agreement between that sample and the density contrasts of the full simulation. In Fig. 4, we give a quantitative comparison derived from the entire posterior distribution. The assessment is done by computing the cross-correlation between the BORG density field and the simulation density field and the variance between those two fields. They show that the BORG density field is unbiased and tightly correlated on large scales.

thumbnail Fig. 3.

Tests of the BORG algorithm on a N-body simulation. Top-left panel: center slice of the density contrast computed at a resolution of 323 from the N-body cosmological simulation used in our test. The simulation itself is obtained with 2563 dark matter particles with ΛCDM cosmology as used in Jasche & Lavaux (2019). Top-right panel: central slice of grid obtained by assigning objects of the mock catalog to a 323 grid with the nearest grid point filter. Bottom-left panel: density slice computed from one posterior sample obtained by BORG. We note the clear spatial correlations between the density field of the fiducial simulation and the inferred sample, despite the very low sampling rate of the mock catalog. Bottom-right panel: difference between the BORG reconstruction and the simulation truth divided by the standard deviation per voxel estimated from the posterior distribution. The BORG posterior covers the truth and gives a conservative estimate of the uncertainty.

thumbnail Fig. 4.

Correlation (blue color) and error residual (orange) for the density field inferred by BORG compared to the simulation density field. We compute the mean and one standard deviation using the posterior distribution sampled by the BORG chain. Both curves are normalized by the power spectrum of each field. We see that indeed when correlation is weak, we get the maximum relative variance.

Figure 5 illustrates the distribution of the differences between the velocity field of the simulation and the one provided by BORG inference. The comparison is done for a single sample obtained by BORG and for linear theory directly from the density contrast obtained from the mock catalog (linear theory). The velocity fields are computed using the Simplex-In-Cell algorithm at high resolution and downgraded at 323 using local averages. In both cases, the reconstructions are unbiased but linear theory exhibits heavier tails for the residuals, leading to a typical width of 200 km s−1. On the other hand, BORG has residuals with a typical width of 139 km s−1. To further check the unbiased-ness, this time using the mock tracers themselves, we show in Fig. 6 how the reconstructed velocities, which are obtained through the BORG velocity field, compares to the original mock tracer velocities directly taken from the N-body simulations. The first three panels (top-left, top-right, lower-left) give the scatter plots for each component of the velocities for each mock tracer. The lower-left panel shows the histogram of the residual velocities for the three components together.

thumbnail Fig. 5.

Tests of the BORG-PM algorithm on a N-body simulation and result for the velocity field. We compute the velocity field using the Simplex-In-Cell algorithm using the entire set of particles of the N-body simulation and a resimulation from a single BORG sample. We compute the linear theory prediction for the velocity field based on the mock catalog as well. We show here the residual between each of the reconstructed velocity field and the one of the simulation. The difference is computed for each mesh element. We note that the residual distribution is narrower and have less heavy tails for the BORG reconstruction than for the linear theory one. In both cases, the velocity field is unbiased.

thumbnail Fig. 6.

Components (xtop-left, ytop-right, zbottom left) of the velocities of the mock tracers as reconstructed from a single sample from BORG and compared to the ones provided by the original N-body simulation. The points are color-coded according to their local density in the scatter plot to allow better identification of the central part and the tails of the distribution. The diagonal red line is added for reference and to indicate unbiased velocity reconstruction. Lower-right panel: histogram of the difference between the BORG reconstructed velocities and the original velocities of the mock tracers for the three components.

We conclude from this test that the peculiar velocities reconstructed with the BORG algorithm are unbiased. They also have typically smaller errors compared to velocities derived using a linear perturbation theory approach.

4.2. Implementation of BORG on 2M++

In Fig. 7 (left panel), we show the velocity field in supergalactic coordinates for 2M++ along with the starred spatial position of NGC 4993. The middle panel in Fig. 7 indicates the estimate of the velocity field from 2M++ in the Supergalactic coordinate system. Additionally, we further show in the right panel the velocity field as inferred from the SDSS-III/BOSS survey. The latter uses a simpler dynamical model based on first-order Lagrangian perturbation theory (Lavaux et al. 2019). The contribution of the non-linear velocity component captured by BORG is shown in Fig. 8. This non-linear velocity component is the one captured by the particle mesh solver of BORG, although it is not yet sufficient to resolve virial motions in large scale structures. We denote this regime by “NL” for non-linear. It corresponds to the intermediate “gray” regime when linear modeling of the velocity field is not sufficient and still not describable without a full non-perturbative description of the non-linear dynamics. To assess its importance, we considered two residuals. The left panel shows the relative contribution of the pure NL velocity from BORG with respect to the total velocity contribution by evaluating the expression |(vr, BORG − vr, linear)/vr, BORG|. It is derived from the difference between the total contribution and the velocity field derived from the gradient of the gravitational potential of the matter density field. We note that the difference between the two fields may be due to at least two causes: either the linear theory approximation underestimates the amplitude of flow in infalling regions or it unphysically overestimates the velocity flow close to the peak of the matter density contrast. The right panel shows the estimated standard deviation of the non-linear velocity field, both the total velocity and the residue, obtained by taking the difference between the total non-linear and the linear component of the velocity field. This quantity corresponds to the second moment of the probability density function, conditioned on δm. We note that this moment does not correspond to the full marginal distribution and, thus, it is possible that the conditional standard deviation is reduced compared to the marginal standard deviation. We compute this standard deviation from each bin of the matter density field. We note that the ratio of these two standard deviations of the total velocity field and the non-linear component of the velocity field can go up to ∼(20−30)% at δm ∼ 4, indicating that the non-linear component is far from negligible.

thumbnail Fig. 7.

Global properties of the velocity field that is available from the BORG inference framework applied on the 2M++ compilation and the SDSS-III/BOSS sample of galaxies. The black cross gives the position of the observer. Left panel: velocity field derived from BORG in the Supergalactic plane. The black star indicates the spatial position of NGC 4993 which is the host of GW170817. Middle panel: spatial distribution of the available model of the velocity field in the Supergalactic plane B = 0° for the volume of 2M++. Right panel: same as a middle panel but for the SDSS-III/BOSS survey. The part at the center is the region from 2M++.

thumbnail Fig. 8.

Non-linearities (NL) modeled by the BORG dynamical model. Left panel: relative difference between the velocity field derived from the gravity field assuming linearity and the non-linear velocity field produced by BORG with respect to the total for the ensemble mean fields. Right panel: standard deviation of the velocity field for each density bin (orange curve) and standard deviation of the non-linear velocity, (vr, BORG − vr, linear), for each density bin (blue curve). We see here that the non-virial, non-linear component of the velocity field can represent a substantial portion of the cosmic velocity field: up to ∼25% of the standard deviation of the total velocity field when the local matter density exceeds δm ≳ 4.

We point that the velocity and density field inference is largely independent of the Hubble constant. Our working coordinate system uses the arbitrary value of H0 = 100 km s−1 Mpc−1. In this coordinate system, the value of H0 is irrelevant for computing the dynamics and the observable in the observer coordinate frame. However, what matters is the redshift evolution of E(z), which only depends on relative expansion history. For the portion of Universe that is considered here, the only relevant parameters are Ωm = 0.315 and ΩΛ = 0.685 (Planck Collaboration VI 2020). The only place where the absolute value of H0 enters is in the prior power spectrum of primordial fluctuations that was used to run the analysis. This is, however, at least a second-order effect as this only governs the way BORG has to fill the missing information in the data space, but not the actual mass distribution that is inferred on large scales.

5. Algorithm for the peculiar velocity correction for standard sirens

The observed GW signal from a network of GW detectors is a probe of the distance to the source, as we discussed earlier on. The GW signal parameters comprise the source parameters, such as the chirp mass, individual masses of the compact objects, spin of the individual compact objects, and the object’s tidal deformability, and the parameters that relate to the observer of the GW source, such as the luminosity distance dL, sky location n ̂ $ \hat n $, polarization angle ψ, and the inclination angle i. It is this latter set that is useful to estimate the value of H0 (Dalal et al. 2006; Nissanke et al. 2010; Abbott et al. 2017a). The GW strain, h+ and h×, provides estimates of the luminosity distance and the inclination angle of the GW source. The sky location of the GW source is obtained using the arrival times of the GW signal at different detectors, as well as their antenna function (Fairhurst 2009, 2011; Nissanke et al. 2011, 2013b; Schutz 2011; Veitch et al. 2012). The redshift of the GW source may be measured by identifying the host galaxy using the electromagnetic counterpart from the GW sources.

We propose a way of estimating the peculiar velocity of the host galaxy based on the estimation of the large-scale flow, vh, from BORG, and the small scale motion as a stochastic velocity dispersion using the fitting form given in Eq. (6). We note that vh is itself composed of the velocity field at the linear order in perturbation, with an additional (at least) 30% correction from non-linear evolution provided by BORG for density δm >  0. Using the location of the identified host galaxy (from the electromagnetic counterpart), we estimate the corresponding posterior probability density function (PDF) of the velocity field component, vh, from 118 realizations of BORG for this galaxy. The virial velocity component is estimated using the mass of the galaxy halo in Eq. (6). The PDF of the non-linear velocity field vvir is assumed to have undergone a Gaussian distribution. The combined posterior of the peculiar velocity, vp = vh + vvir, can be obtained by convolving the BORG posterior of vh with the PDF of vvir.

By combining the measurements from GW data, dGW, the redshift of the host galaxy, z ̂ $ \hat z $, and the peculiar velocity estimate, vp, and luminosity distance, dL, we can obtain the posterior of H0 for n GW sources by using the Bayesian framework:

P ( H 0 | { d GW } , { z ̂ } ) ∝ ∏ i = 1 n ∫ d d L i d v p i L ( d L i | H 0 , v p i , z i , u ̂ i , d GW i ) × P ( z i | u ̂ i ) P ( v p i | M , u ̂ i ) Π ( H 0 ) , $$ \begin{aligned} P(H_0|\, \{d_{\rm GW}\},\, \{\hat{z}\}) \propto &\prod _{i=1}^{\boldsymbol{n}}\int \mathrm{d}d^i_{\rm L}\, \mathrm{d}{{ v}^i_{\rm p}} \ \mathcal{L} (d^i_{\rm L}| H_0, { v}^i_{\rm p}, z_i, {\boldsymbol{\hat{u}}}_i, d^i_{\rm GW})\nonumber \\&\times \mathcal{P} (z_i|{\boldsymbol{\hat{u}}}_i)\, \mathcal{P} ({ v}^i_{\rm p}| M,{\boldsymbol{\hat{u}}}_i)\ \Pi (H_0), \end{aligned} $$(7)

where u ̂ i ( RA i , Dec i ) $ {\boldsymbol{\hat{u}}}_i(\mathrm{RA}_i,\,\mathrm{Dec}_i) $ is the sky position6, ℒ denotes the likelihood which is assumed to be Gaussian, P ( z i | u ̂ i ) $ \mathcal{P}(z_i|{\boldsymbol{\hat{u}}}_i) $ is the posterior of the redshift estimate at the GW source, u ̂ i $ \hat{u}_i $, P ( v p i | M , u ̂ i ) $ \mathcal{P}(\mathit{v}^i_{\mathrm{p}}| M,\,{\boldsymbol{\hat{u}}}_i) $ is the peculiar velocity estimate using the cosmological method, and Π(H0) is the prior on the value of H0.

6. Applying our method to GW170817

6.1. Peculiar velocity estimate for NGC 4993

NGC 4993 is the host galaxy of the event GW170817 (the merger of two neutron stars) which was discovered by the LIGO Scientific Collaboration and Virgo Collaboration (LVC Abbott et al. 2017b,a). The distance to the GW170817 inferred from the GW observation is d L = 43 . 8 − 6.9 + 2.9 $ d_{\mathrm{L}}=43.8_{-6.9}^{+2.9} $ Mpc (Abbott et al. 2017b,a) and the recession velocity of its host NGC 4993, with respect to the CMB rest frame, is 3327 ± 72 km s−1 (Abbott et al. 2017a); Crook:2006sw The velocity estimate for NGC 4993 has two components, vh and vvir. We estimate the vh component from BORG, and the non-linear part of the velocity is obtained using the form mentioned in Eq. (6). NGC 4993 has a halo of mass about Mh ∼ 1012 M⊙ (Pan et al. 2017; Ebrová & Bílek 2020). The corresponding estimate of σvir is about 100 km s−1. The excess velocity component is included as an additional velocity dispersion assuming a Gaussian PDF of variance σ vir 2 $ \sigma^2_{\rm vir} $ with a zero mean.

The posterior of the peculiar velocity from BORG is shown in blue in Fig. 9 for the Planck-2015 best-fit cosmological parameters (Planck Collaboration XIII 2016). The PDF of the velocity field is non-Gaussian. The mean value of the velocity field is v ¯ p = 373 $ \mathit{{\bar v}}_{\mathrm{p}}= 373 $ km s−1 and the maximum a posterior value of the velocity field is v p MP = 385 $ \mathit{v}^{\mathrm{MP}}_{\mathrm{p}}= 385 $ km s−1. The standard deviation of the velocity pdf is σh ∼ 82 km s−1. The total standard deviation of the velocity due to both vh and vvir is σ t = σ h 2 + σ vir 2 ∼ 130 $ \sigma_{\mathrm{t}}= \sqrt{\sigma_{\mathrm{h}}^2+ \sigma_{\mathrm{vir}}^2} \sim 130 $ km s−1. The corresponding PDF for vp is shown by the curve in black color in Fig. 9. The inferred value of the peculiar velocity of NGC 4993 differs from the value assumed by LVC (Abbott et al. 2017a). Abbott et al. (2017a) considered a velocity distribution that was Gaussian with a mean and standard deviation given by 310 km s−1 and 150 km s−1, respectively. The estimates of the velocity from our method predict a 24% higher mean value of the velocity and about 13% less standard deviation. The comparison of the distribution of the peculiar velocity between the LVC and our method is shown in red and black, respectively, in Fig. 9.

thumbnail Fig. 9.

Posterior of the peculiar velocity of NGC 4993. The blue curve displays the large scale flow, vp, inferred from BORG while the black curve gives the required total velocity including the virial component within the halo. The posterior of the velocity used by the LIGO Scientific Collaboration and Virgo Collaboration (LVC) is shown with red dashes.

6.2. Revised H0 from GW170817

Using our new, more precise estimate of the peculiar velocity for the NGC 4993, we obtain a revised value of H0 using the Bayesian framework given in Eq. (7) from the data of luminosity distance and inclination angle, assuming either a high-spin or a low-spin prior on each compact object (Abbott et al. 2019a). Figure 10 shows the corresponding posterior of H0. The posterior of H0 peaks at H0 = 68.3 km s−1 Mpc−1, instead of 70 − 8 + 12 $ 70_{-8}^{+12} $ km s−1 Mpc−1 from the results of LVC (Abbott et al. 2017a) with the same 68.3% credible interval as obtained by LVC. The marginal difference in the shape of the posterior between our method (black line) and LVC (red dashed line) arises only as a result of the difference in the velocity correction between our method and LVC. In addition to calibration noise and peculiar velocity measurements, the main source of the error in the measurement of H0 is due to the degeneracy between the inclination angle, i, and the distance, dL (Abbott et al. 2017a). This acts as a limiting factor for the measurement of H0 from GW170817 if there is no independent measurement of inclination angle. Critically, the addition of Virgo (Acernese et al. 2015) to the nearly co-aligned two LIGO detectors (Abramovici et al. 1992; Martynov et al. 2016) allowed for the improvement in the sky localization of the GW source and offered the possibility of improving the measurement of the inclination angle.

thumbnail Fig. 10.

Posterior of H0 for GW170817 is obtained using the peculiar velocity correction by our method and is shown in black when we use the GW data and in blue when using the GW and VLBI data. We also show the H0 measurement from the LIGO Scientific Collaboration and Virgo Collaboration (LVC) in dashed red line (Abbott et al. 2017a) and the magenta dashes presents the result from Hotokezaka et al. (2019) from hydrodynamics simulation jet using a GW+VLBI observation.

A further improvement in the measurement of H0 is possible if the uncertainty in the inclination angle, i, can be reduced. A recent study has used the observed data of the electromagnetic counterpart to GW170817, the superluminal motion measured by the Very Large Baseline Interferometry (VLBI) observations (Mooley et al. 2018a) and the afterglow light curve data (e.g., Mooley et al. 2018b) to constrain the inclination angle of GW170817. Using the constraints on the inclination angle, Hotokezaka et al. (2019) obtained a revised value of H0. We implement our velocity-field correction method for the combined measurement of GW+VLBI, by using a flat prior on the value of inclination angle between 0.25 rad ≤ i ( d L 41 Mpc ) ≤ 0.45 $ i\,(\frac{d_{\mathrm{L}}}{41\,\mathrm{Mpc}}) \leq 0.45 $ rad. The revised H0 by our method is shown in Fig. 10 in blue and is compared with the H0 value from Hotokezaka et al. (2019) (shown by the magenta dashed line). The error in the measurement of H0 from GW+VLBI (Hotokezaka et al. 2019) arises from several sources, such as the GW data, the shape of the light curve, flux centroid motion, and the peculiar velocity. Our new estimate of the peculiar velocity improves the precision of the H0 measurement with GW+VLBI from 70 . 3 − 5.0 + 5.3 $ 70.3^{+5.3}_{-5.0} $ km s−1 Mpc−1 to 68 . 3 − 4.5 + 4.6 $ 68.3^{+ 4.6}_{-4.5} $ km s−1 Mpc−1.

7. Conclusions and future prospects

Cosmology with gravitational waves is a newly emerging field that holds the enormous potential to explore new aspects of the Universe and fundamental physics (Saltas et al. 2014; Lombriser & Taylor 2016; Lombriser & Lima 2017; Sakstein & Jain 2017; Baker et al. 2017; Nishizawa 2018; Belgacem et al. 2018a,b, 2019; Pardo et al. 2018; Abbott et al. 2019b; Mukherjee et al. 2020a,b). One such primary science goal is the measurement of the Hubble constant from the GW sources when the redshift to the source can be identified from the electromagnetic counterpart or the cross-correlation with galaxy surveys (Oguri 2016; Mukherjee & Wandelt 2018). Using the current GW detectors, such as Advanced LIGO (Martynov et al. 2016) and Advanced Virgo (Acernese et al. 2015), recent forecasts predict the ability to measure H0 with an accuracy of less than 2% with about 50 binary neutron star (BNS, Chen et al. 2018; Feeney et al. 2019). To achieve such a precise estimation of H0, it is essential to accurately correct for the peculiar velocity of the hosts of GW sirens. This is especially the case for very low redshift sources, which will have a high signal-to-noise ratio (S/N) and, hence, significantly contribute to the joint posteriors for ensembles of events. In addition, peculiar velocity corrections will be important for events with a more favorable S/N in the Virgo or additional other detectors, as well as neutron star-black hole mergers (Nissanke et al. 2010; Vitale & Chen 2018). The absence of such a correction will affect both the accuracy and precision of the measurement of the value of H0. For a typical value of the peculiar velocity of about ∼300 km s−1 for a GW host at redshift z = 0.01, the contribution from peculiar velocity is comparable to the term related to the Hubble flow. The contribution becomes less (or more) severe at higher (or lower) redshift. As a result, we need to estimate the value of the peculiar velocity with sufficient accuracy to avoid any systematic bias and additional variance in the measurement of H0. Although averaging over large GW samples can lead to an unbiased estimate of H0, the absence of a peculiar velocity correction will increase the error budget on H0 due to the additional variance from the peculiar velocity contribution. In such a case, it’s necessary to have a larger number of GW samples of Ngw to beat the variance as N gw − 1 / 2 $ N_{\mathrm{gw}}^{-1/2} $ (Nissanke et al. 2010, 2013a; Chen et al. 2018; Feeney et al. 2019; Mortlock et al. 2019). Incorporating an accurate correction for the peculiar velocity of the host of GW sources, we can achieve faster and more economically accurate, and also precise, measurements of H0.

In this paper, we use a statistical reconstruction method to correct for the peculiar velocity of the host of the GW sources. The peculiar velocity for the host galaxy arises from the gravitational potential of the cosmic density field. We estimate the posterior of the peculiar velocity for both the linear and the non-linear component. The large-scale velocity flow is estimated using the Bayesian framework called BORG. The stochastic velocity dispersion of the source is obtained using a numerical fitting-form given in Eq. (6), by using the mass of the halo of the host. By combining the results from peculiar velocity estimate for each source, their redshifts, and the inferred luminosity distance from the GW data, we obtain a Bayesian inference of the value of H0 according to the framework discussed around Eq. (7).

We implemented our method on GW170817 (Abbott et al. 2017b), which is in the host NGC 4993 (Pan et al. 2017). The corresponding posterior distribution with the results from LVC (Abbott et al. 2017a) and GW+VLBI (Hotokezaka et al. 2019) are shown in Fig. 10 by the solid line in black and in blue, respectively. While our correction marginally reduces the maximum a posteriori value of H0 to 68 . 3 − 4.5 + 4.6 $ 68.3^{+4.6}_{-4.5} $ km s−1 Mpc−1 for GW+VLBI, it slightly disfavors very low values of H0 ≲ 60 km s−1 Mpc−1 compared to the results from LVC. As we were preparing this manuscript for submission, an analysis was published (Howlett & Davis 2020) that was based on the implementation of an alternative velocity correction approach to the host of GW170817 and recovered a value of H 0 = 64 . 8 − 7.2 + 7.3 $ H_0= 64.8^{+7.3}_{-7.2} $ km s−1 Mpc−1. The mean value of H0 differs from the values obtained in this work ( H 0 = 68 . 3 − 4.5 + 4.6 $ H_0= 68.3^{+ 4.6}_{-4.5} $ km s−1 Mpc−1) and also from the value obtained by Hotokezaka et al. (2019) ( H 0 = 70 . 3 − 5.0 + 5.3 $ H_0= 70.3^{+5.3}_{-5.0} $ km s−1 Mpc−1). However, the value by Howlett & Davis (2020) is consistent within the error-bars of both the estimations. Also, another recent work (Nicolaou et al. 2020) discussed the possible impact of peculiar velocity on the standard siren measurements and obtained a revised value of H 0 = 68 . 6 − 8.5 + 14 $ H_0= 68.6_{-8.5}^{+14} $ km s−1 Mpc−1 for the event GW170817. Furthermore, a recent work by Boruah et al. (2020) also discusses the impact of the line of sight marginalization for models derived using linear theory applied to spectroscopic data calibrated with supernovae data and the velocity field derived from the distance data obtained from Tully–Fisher or Fundamental Plane methods.

The velocity correction is readily available for GW sources with electromagnetic counterparts (such as binary neutron stars; black hole neutron stars) in the cosmic volumes covered by the 2M++ (Lavaux & Hudson 2011; Jasche & Lavaux 2019) and the SDSS-III/BOSS surveys (Eisenstein et al. 2011; Dawson et al. 2013; Lavaux et al. 2019 see Fig. 7). The 2M++ volume covers galactic latitudes |b| > 10° up to a redshift of z ∼ 0.05; the SDSS-III/BOSS survey spans redshifts of z ∼ 0.2 − 0.7 for the sky areas (0° < Dec <  60° and 120° < RA <  240°) or (0° < Dec <  30° and |RA| < 30°). We expect that within a year our method will be available for the SDSS-IV/eBOSS survey (Dawson et al. 2016). In the long term, with the availability of the nearly full-sky data sets jointly from the upcoming cosmology missions such as Dark Energy Spectroscopic Instrument, (DESI, Aghamousa et al. 2016) and the Large Synoptic Survey Telescope (LSST, LSST Science Collaboration 2009), our algorithm will be available for most of the low-redshift GW sources that are to be observed by the currently planned network of ground-based GW detectors (Schutz 2011; Abbott et al. 2018b).

1

The International Astronomical Union has recently renamed the Hubble law as the Hubble-Lemaître law, in recognition of the pioneering contribution of Lemaître (https://www.iau.org/news/pressreleases/detail/iau1812/).

2

λo and λe are the measured and emitted wavelength of the light, respectively.

3

The perturbation in the homogeneous FLRW metric is caused by the presence of structures in the Universe.

4

The line of sight vector u ̂ $ \hat u $ is directed from the observer to the source, i.e., v . u ̂ ≥ 0 $ {\boldsymbol{v}}.\hat u \geq 0 $, implying that the source is moving away from the observer.

5

SDSS-III/BOSS is the Baryon Oscillation Spectroscopic Survey.

6

RA and Dec respectively denote the right accession and declination.

Acknowledgments

SM would like to thank Emanuele Berti, Rahul Biswas, Stephen Feeney, Daniel Mortlock, Daniel Sclonic, and Joseph Silk for fruitful discussions. SM and SMN would like to thank Hiranya Peiris and Will Farr for useful comments on the draft of the paper. The Flatiron Institute is supported by the Simons Foundation. SM is supported by the Lagrange postdoctoral Fellowship at Institut d’Astrophysique de Paris. The work of SM and BDW is also supported by the Labex ILP (reference ANR-10-LABX-63) part of the Idex SUPER, and received financial state aid managed by the Agence Nationale de la Recherche, as part of the programme Investissements d’avenir under the reference ANR-11-IDEX-0004-02. BDW thanks the CCPP at New York University for its hospitality while this work was completed. This work was supported by the ANR BIG4 project, grant ANR-16-CE23-0002 of the French Agence Nationale de la Recherche. SM and SMN are grateful for support from NWO VIDI and Projectruimte Grants of the Innovational Research Incentives Scheme (Vernieuwingsimpuls) financed by the Netherlands Organization for Scientific Research (NWO). FL acknowledges funding from the Imperial College London Research Fellowship Scheme. KH is supported by Lyman Spitzer Jr. Fellowship at Department of Astrophysical Sciences, Princeton University. This work was granted access to the HPC resources of CINES (Centre Informatique National de l’Enseignement Supérieur) under the allocation A0020410153 made by GENCI and has made use of the Horizon cluster hosted by the Institut d’Astrophysique de Paris in which the cosmological simulations were post-processed. We thank Stéphane Rouberol for running smoothly this cluster for us. This work is done within the Aquila. Consortium (https://www.aquila-consortium.org/). We acknowledge the use of following packages in this analysis: Astropy (Astropy Collaboration 2013, 2018), emcee: MCMC Hammer (Foreman-Mackey et al. 2013), IPython (Pérez & Granger 2007), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), and SciPy (Jones et al. 2001). The authors would like to thank the LIGO scientific collaboration and Virgo collaboration for providing the data for the event GW170817.

References

  1. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, Nature, 551, 85 [Google Scholar]
  2. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, Phys. Rev. Lett., 119, 161101 [Google Scholar]
  3. Abbott, T. M. C., Abdalla, F. B., Annis, J., et al. 2018a, MNRAS, 480, 3879 [Google Scholar]
  4. Abbott, B. P., Abbott, R., Abbott, T. D. et al. 2018b, Liv. Rev. Rel., 21, 3 [Google Scholar]
  5. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019a, Phys. Rev. X, 9, 011001 [Google Scholar]
  6. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019b, Phys. Rev. Lett., 123, 011102 [Google Scholar]
  7. Abramovici, A., Althouse, W. E., Drever, R. W. P., et al. 1992, Science, 256, 325 [Google Scholar]
  8. Acernese, F., Agathos, M., Agatsuma, K., et al. 2015, Class. Quant. Grav., 32, 024001 [Google Scholar]
  9. Addison, G. E., Watts, D. J., Bennett, C. L., et al. 2018, ApJ, 853, 119 [Google Scholar]
  10. Aghamousa, A., Aguilar, J., Ahlen, S., et al. 2016, ArXiv e-prints [arXiv:1611.00036] [Google Scholar]
  11. Agrawal, P., Cyr-Racine, F.-Y., Pinner, D., & Randall, L. 2019, ArXiv e-prints [arXiv:1904.01016] [Google Scholar]
  12. Anderson, L., Aubourg, É., Bailey, S., et al. 2014, MNRAS, 441, 24 [Google Scholar]
  13. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  15. Aubourg, É., Bailey, S., Bautista, J. E., et al. 2015, Phys. Rev. D, 92, 123516 [Google Scholar]
  16. Baker, T., Bellini, E., Ferreira, P. G., et al. 2017, Phys. Rev. Lett., 119, 251301 [Google Scholar]
  17. Barausse, E., Matarrese, S., & Riotto, A. 2005, Phys. Rev. D, 71, 063537 [Google Scholar]
  18. Belgacem, E., Dirian, Y., Foffa, S., & Maggiore, M. 2018a, Phys. Rev. D, 97, 104066 [Google Scholar]
  19. Belgacem, E., Dirian, Y., Foffa, S., & Maggiore, M. 2018b, Phys. Rev. D, 98, 023510 [Google Scholar]
  20. Belgacem, E., Calcagnib, G., Crisostomi, M., et al. 2019, JCAP, 1907, 024 [Google Scholar]
  21. Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, A20 [Google Scholar]
  22. Bernal, J. L., Verde, L., & Riess, A. G. 2016, JCAP, 1610, 019 [Google Scholar]
  23. Bonvin, C., Durrer, R., & Gasparini, M. A. 2006, Phys. Rev. D, 73, 023523; erratum: (2012) Phys. Rev. D, 85, 029901 [Google Scholar]
  24. Boruah, S. S., Hudson, M. J., Lavaux, G., et al. 2020, ArXiv e-prints [arXiv:2010.01119] [Google Scholar]
  25. Branchini, E., Teodoro, L., Frenk, C. S., et al. 1999, MNRAS, 308, 1 [Google Scholar]
  26. Branchini, E., Freudling, W., Da Costa, L. N., et al. 2001, MNRAS, 326, 1191 [Google Scholar]
  27. Carrick, J., Turnbull, S. J., Lavaux, G., & Hudson, M. J. 2015, MNRAS, 450, 317 [Google Scholar]
  28. Chen, H.-Y., Fishbach, M., & Holz, D. E. 2018, Nature, 562, 545 [Google Scholar]
  29. Chodorowski, M. J., & Ciecielag, P. 2002, MNRAS, 331, 133 [Google Scholar]
  30. Crook, A. C., Huchra, J. P., Martimbeau, N., et al. 2007, ApJ, 655, 790 [Google Scholar]
  31. Cuesta, A. J., Vargas-Magaña, M., Beutler, F., et al. 2016, MNRAS, 457, 1770 [Google Scholar]
  32. Cutler, C., & Flanagan, E. E. 1994, Phys. Rev. D, 49, 2658 [Google Scholar]
  33. Dalal, N., Holz, D. E., Hughes, S. A., & Jain, B. 2006, Phys. Rev. D, 74, 063006 [Google Scholar]
  34. Davis, M., Nusser, A., & Willick, J. A. 1996, ApJ, 473, 22 [Google Scholar]
  35. Davis, M., Nusser, A., Masters, K. L., et al. 2011, MNRAS, 413, 2906 [Google Scholar]
  36. Dawson, K. S., Schlegel, D. J., Ahn, C. P., et al. 2013, AJ, 145, 10 [Google Scholar]
  37. Dawson, K. S., Kneib, J.-P., Percival, W. J., et al. 2016, AJ, 151, 44 [Google Scholar]
  38. Di Valentino, E., Melchiorri, A., & Mena, O. 2017, Phys. Rev. D, 96, 043503 [Google Scholar]
  39. Ebrová, I., & Bílek, M. 2020, A&A, 634, A73 [CrossRef] [EDP Sciences] [Google Scholar]
  40. Eisenstein, D. J., Weinberg, D. H., Agol, E., et al. 2011, AJ, 142, 72 [Google Scholar]
  41. Erdoǧdu, P., Lahav, O., Huchra, J. P., et al. 2006, MNRAS, 373, 45 [Google Scholar]
  42. Fairhurst, S. 2009, New J. Phys., 11, 123006 [Google Scholar]
  43. Fairhurst, S. 2011, Class. Quant. Grav., 28, 105021 [Google Scholar]
  44. Feeney, S. M., Peiris, H. V., Williamson, A. R., et al. 2019, Phys. Rev. Lett., 122, 061105 [Google Scholar]
  45. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
  46. Freedman, W. L., Madore, B. F., Hatt, D., et al. 2019, ApJ, 882, 34 [Google Scholar]
  47. Hahn, O., Angulo, R. E., & Abel, T. 2015, MNRAS, 454, 3920 [Google Scholar]
  48. Hawking, S. W., & Israel, W. 1987, Three Hundred Years of Gravitation (Cambridge: Cambridge University Press) [Google Scholar]
  49. Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19 [Google Scholar]
  50. Holz, D. E., & Hughes, S. A. 2005, ApJ, 629, 15 [Google Scholar]
  51. Hotokezaka, K., Nakar, E., Gottlieb, O., et al. 2019, ArXiv e-prints [arXiv:1806.10596] [Google Scholar]
  52. Howlett, C., & Davis, T. M. 2020, MNRAS, 492, 3803 [Google Scholar]
  53. Hubble, E. 1929, Proc. Natl. Acad. Sci., 15, 168 [Google Scholar]
  54. Hudson, M. J. 1994, MNRAS, 266, 475 [Google Scholar]
  55. Hudson, M. J., Smith, R. J., Lucey, J. R., & Branchini, E. 2004, MNRAS, 352, 61 [Google Scholar]
  56. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  57. Jasche, J., & Lavaux, G. 2019, A&A, 625, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Jasche, J., & Wandelt, B. D. 2013, MNRAS, 432, 894 [Google Scholar]
  59. Jasche, J., Leclercq, F., & Wandelt, B. D. 2015, JCAP, 2015, 036 [Google Scholar]
  60. Jee, I., Suyu, S., Komatsu, E., et al. 2019, Science, 365, 1134 [Google Scholar]
  61. Jennings, E. 2012, MNRAS, 427, L25 [NASA ADS] [Google Scholar]
  62. Jennings, E., Baugh, C. M., & Pascoli, S. 2011, MNRAS, 410, 2081 [NASA ADS] [Google Scholar]
  63. Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open Source Scientific Tools for Python [Google Scholar]
  64. Kaiser, N., Efstathiou, G., Saunders, W., et al. 1991, MNRAS, 252, 1 [Google Scholar]
  65. Knox, L., & Millea, M. 2020, Phys. Rev. D, 101, 043533 [Google Scholar]
  66. Kolb, E. W., Matarrese, S., Notari, A., & Riotto, A. 2005, Phys. Rev. D, 71, 023524 [Google Scholar]
  67. Kreisch, C. D., Cyr-Racine, F.-Y., & Doré, O. 2020, Phys. Rev. D, 101, 123505 [Google Scholar]
  68. Lavaux, G., & Hudson, M. J. 2011, MNRAS, 416, 2840 [Google Scholar]
  69. Lavaux, G., & Jasche, J. 2016, MNRAS, 455, 3169 [Google Scholar]
  70. Lavaux, G., Jasche, J., & Leclercq, F. 2019, MNRAS, submitted [Google Scholar]
  71. Leclercq, F., Jasche, J., & Wandelt, B. 2015, JCAP, 2015, 15 [Google Scholar]
  72. Leclercq, F., Jasche, J., Lavaux, G., Wandelt, B., & Percival, W. 2017, JCAP, 2017, 049 [Google Scholar]
  73. Lemaître, G. 1927, Ann. Soc. Sci. Bruxelles, 47, 49 [Google Scholar]
  74. Lemaître, G. 1931, MNRAS, 91, 483 [NASA ADS] [Google Scholar]
  75. Lin, M.-X., Benevento, G., Hu, W., & Raveri, M. 2019, Phys. Rev. D, 100, 063542 [Google Scholar]
  76. Lombriser, L., & Lima, N. A. 2017, Phys. Lett. B, 765, 382 [Google Scholar]
  77. Lombriser, L., & Taylor, A. 2016, JCAP, 1603, 031 [Google Scholar]
  78. LSST Science Collaboration (Abell, P. A., et al.) 2009, ArXiv e-prints [arXiv:0912.0201] [Google Scholar]
  79. Ma, Y.-Z., Branchini, E., & Scott, D. 2012, MNRAS, 425, 2880 [Google Scholar]
  80. Macaulay, E., Nichol, R. C., Bacon, D., et al. 2019, MNRAS, 486, 2184 [Google Scholar]
  81. Maggiore, M. 2008, Gravitational Waves: Volume 1: Theory and Experiments, Gravitational Waves (Oxford: Oxford University Press) [Google Scholar]
  82. Martynov, D. V., Hall, E. D., Abbott, B. P., et al. 2016, Phys. Rev. D, 93, 112004; erratum: (2018) Phys. Rev. D, 97, 059901 [Google Scholar]
  83. Mooley, K. P., Deller, A. T., Gottlieb, O., et al. 2018a, Nature, 561, 355 [Google Scholar]
  84. Mooley, K. P., Frail, D. A., Dobie, D., et al. 2018b, ApJ, 868, L11 [Google Scholar]
  85. Mortlock, D. J., Feeney, S. M., Peiris, H. V., Williamson, A. R., & Nissanke, S. M. 2019, Phys. Rev. D, 100, 103523 [Google Scholar]
  86. Mukherjee, S., & Wandelt, B. D. 2018, ArXiv e-prints [arXiv:1808.06615] [Google Scholar]
  87. Mukherjee, S., Wandelt, B. D., & Silk, J. 2020a, Phys. Rev. D, 101, 103509 [Google Scholar]
  88. Mukherjee, S., Wandelt, B. D., & Silk, J. 2020b, MNRAS, 494, 1956 [Google Scholar]
  89. Neyrinck, M. C., Aragón-Calvo, M. A., Jeong, D., & Wang, X. 2014, MNRAS, 441, 646 [Google Scholar]
  90. Nicolaou, C., Lahav, O., Lemos, P., Hartley, W., & Braden, J. 2020, MNRAS, 495, 90 [Google Scholar]
  91. Nishizawa, A. 2018, Phys. Rev. D, 97, 104037 [Google Scholar]
  92. Nissanke, S., Holz, D. E., Hughes, S. A., Dalal, N., & Sievers, J. L. 2010, ApJ, 725, 496 [Google Scholar]
  93. Nissanke, S., Sievers, J., Dalal, N., & Holz, D. 2011, ApJ, 739, 99 [Google Scholar]
  94. Nissanke, S., Holz, D. E., Dalal, N., et al. 2013a, ArXiv e-prints [arXiv:1307.2638] [Google Scholar]
  95. Nissanke, S., Kasliwal, M., & Georgieva, A. 2013b, ApJ, 767, 124 [Google Scholar]
  96. Nusser, A., Da Costa, L. N., Branchini, E., et al. 2001, MNRAS, 320, L21 [Google Scholar]
  97. Oguri, M. 2016, Phys. Rev. D, 93, 083511 [Google Scholar]
  98. Pan, Y. C., Kilpatrick, C. D., Simon, J. D., et al. 2017, ApJ, 848, L30 [Google Scholar]
  99. Pardo, K., Fishbach, M., Holz, D. E., & Spergel, D. N. 2018, JCAP, 1807, 048 [Google Scholar]
  100. Pérez, F., & Granger, B. E. 2007, Comput. Sci. Eng., 9, 21 [Google Scholar]
  101. Pike, R. W., & Hudson, M. J. 2005, ApJ, 635, 11 [Google Scholar]
  102. Planck Collaboration XVI. 2014, A&A, 571, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  103. Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  104. Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  105. Poisson, E., & Will, C. M. 1995, Phys. Rev. D, 52, 848 [Google Scholar]
  106. Poulin, V., Smith, T. L., Karwal, T., & Kamionkowski, M. 2019, Phys. Rev. Lett., 122, 221301 [Google Scholar]
  107. Radburn-Smith, D. J., Lucey, J. R., & Hudson, M. J. 2004, MNRAS, 355, 1378 [Google Scholar]
  108. Reid, M. J., Braatz, J. A., Condon, J. J., et al. 2009, ApJ, 695, 287 [Google Scholar]
  109. Riess, A. G., Casertano, S., Yuan, W., Macri, L. M., & Scolnic, D. 2019, ApJ, 876, 85 [Google Scholar]
  110. Sakstein, J., & Jain, B. 2017, Phys. Rev. Lett., 119, 251303 [Google Scholar]
  111. Saltas, I. D., Sawicki, I., Amendola, L., & Kunz, M. 2014, Phys. Rev. Lett., 113, 191101 [Google Scholar]
  112. Sasaki, M. 1987, MNRAS, 228, 653 [Google Scholar]
  113. Saunders, W., Oliver, S., Keeble, O., et al. 2000, in Cosmic Flows Workshop, eds. S. Courteau, & J. Willick, ASP Conf. Ser., 201, 223 [Google Scholar]
  114. Schutz, B. F. 1986, Nature, 323, 310 [Google Scholar]
  115. Schutz, B. F. 2011, Class. Quant. Grav., 28, 125023 [Google Scholar]
  116. Scoccimarro, R. 2004, Phys. Rev. D, 70, 083007 [Google Scholar]
  117. Seto, N., & Kyutoku, K. 2018, MNRAS, 475, 4133 [Google Scholar]
  118. Shaya, E. J., Tully, R. B., & Pierce, M. J. 1992, ApJ, 391, 16 [Google Scholar]
  119. Sheth, R. K., & Diaferio, A. 2001, MNRAS, 322, 901 [Google Scholar]
  120. Springel, V. 2005, MNRAS, 364, 1105 [Google Scholar]
  121. Springob, C. M., Masters, K. L., Haynes, M. P., Giovanelli, R., & Marinoni, C. 2007, ApJS, 172, 599 [Google Scholar]
  122. Springob, C. M., Magoulas, C., Colless, M., et al. 2014, MNRAS, 445, 2677 [Google Scholar]
  123. Turnbull, S. J., Hudson, M. J., Feldman, H. A., et al. 2012, MNRAS, 420, 447 [Google Scholar]
  124. van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22 [Google Scholar]
  125. Veitch, J., Mandel, I., Aylott, B., et al. 2012, Phys. Rev. D, 85, 104045 [Google Scholar]
  126. Verde, L., Protopapas, P., & Jimenez, R. 2013, Phys. Dark Univ., 2, 166 [Google Scholar]
  127. Verde, L., Treu, T., & Riess, A. G. 2019, Nat. Astron., 3, 891 [Google Scholar]
  128. Vitale, S., & Chen, H.-Y. 2018, Phys. Rev. Lett., 121, 021303 [Google Scholar]
  129. Wong, K. C., Suyu, S. H., Chen, G. C.-F., et al. 2019, MNRAS, 498, 1420 [Google Scholar]
  130. Yuan, W., Riess, A. G., Macri, L. M., Casertano, S., & Scolnic, D. 2019, ApJ, 886, 61 [Google Scholar]

All Figures

thumbnail Fig. 1.

Display of the correlation coefficient between the non-linear density and the negative divergence of the non-linear velocity field, θ ≡ −∇ ⋅ v. The red line is the ensemble average over the posterior distribution of the cross-correlation and the shaded dark red area is the one standard deviation according to the mean of that cross-correlation derived from the same distribution obtained by BORG. The green and shaded green region is the result of computing the same quantities on 80 random realizations. Apart from very minor differences, the two curves are on top of each other, which can be expected as the physical model is the same in both cases. The vertical grey band on the right side indicate the resolution limit of the BORG reconstruction. The correlation reaches asymptotically one for the largest scales (small k values). The red line is reaching a value of 99% for the largest probed mode. As expected, this cross-correlation goes to zero at small scales (high k value).
In the text
thumbnail Fig. 2.

Ensemble average power-spectrum of the non-linear velocity field scaled, alongside with one standard deviation contour in the shaded area. The line and shaded region in red are computed from the ensemble posterior distribution built from the 2M++ data. The line and shaded region in green are measured from 80 random realizations of a ΛCDM universe of the same size. Those random realizations were computed using the same N-body solver, but starting from a random sample of a Gaussian distribution with a power-spectrum provided by the ΛCDM model. We show in thick blue line the power spectrum assuming linear theory. We note that large scales are unaffected by non-linearities. Intermediate scales (k ≳ 0.1 h Mpc−1, corresponding to ∼60 h Mpc−1) get non-negligible contributions from non-linear dynamics, as highlighted by the area under the curve. The vertical grey band on the right side indicate the resolution limit of the BORG reconstruction (∼2.6 h Mpc−1).
In the text
thumbnail Fig. 3.

Tests of the BORG algorithm on a N-body simulation. Top-left panel: center slice of the density contrast computed at a resolution of 323 from the N-body cosmological simulation used in our test. The simulation itself is obtained with 2563 dark matter particles with ΛCDM cosmology as used in Jasche & Lavaux (2019). Top-right panel: central slice of grid obtained by assigning objects of the mock catalog to a 323 grid with the nearest grid point filter. Bottom-left panel: density slice computed from one posterior sample obtained by BORG. We note the clear spatial correlations between the density field of the fiducial simulation and the inferred sample, despite the very low sampling rate of the mock catalog. Bottom-right panel: difference between the BORG reconstruction and the simulation truth divided by the standard deviation per voxel estimated from the posterior distribution. The BORG posterior covers the truth and gives a conservative estimate of the uncertainty.
In the text
thumbnail Fig. 4.

Correlation (blue color) and error residual (orange) for the density field inferred by BORG compared to the simulation density field. We compute the mean and one standard deviation using the posterior distribution sampled by the BORG chain. Both curves are normalized by the power spectrum of each field. We see that indeed when correlation is weak, we get the maximum relative variance.
In the text
thumbnail Fig. 5.

Tests of the BORG-PM algorithm on a N-body simulation and result for the velocity field. We compute the velocity field using the Simplex-In-Cell algorithm using the entire set of particles of the N-body simulation and a resimulation from a single BORG sample. We compute the linear theory prediction for the velocity field based on the mock catalog as well. We show here the residual between each of the reconstructed velocity field and the one of the simulation. The difference is computed for each mesh element. We note that the residual distribution is narrower and have less heavy tails for the BORG reconstruction than for the linear theory one. In both cases, the velocity field is unbiased.
In the text
thumbnail Fig. 6.

Components (xtop-left, ytop-right, zbottom left) of the velocities of the mock tracers as reconstructed from a single sample from BORG and compared to the ones provided by the original N-body simulation. The points are color-coded according to their local density in the scatter plot to allow better identification of the central part and the tails of the distribution. The diagonal red line is added for reference and to indicate unbiased velocity reconstruction. Lower-right panel: histogram of the difference between the BORG reconstructed velocities and the original velocities of the mock tracers for the three components.
In the text
thumbnail Fig. 7.

Global properties of the velocity field that is available from the BORG inference framework applied on the 2M++ compilation and the SDSS-III/BOSS sample of galaxies. The black cross gives the position of the observer. Left panel: velocity field derived from BORG in the Supergalactic plane. The black star indicates the spatial position of NGC 4993 which is the host of GW170817. Middle panel: spatial distribution of the available model of the velocity field in the Supergalactic plane B = 0° for the volume of 2M++. Right panel: same as a middle panel but for the SDSS-III/BOSS survey. The part at the center is the region from 2M++.
In the text
thumbnail Fig. 8.

Non-linearities (NL) modeled by the BORG dynamical model. Left panel: relative difference between the velocity field derived from the gravity field assuming linearity and the non-linear velocity field produced by BORG with respect to the total for the ensemble mean fields. Right panel: standard deviation of the velocity field for each density bin (orange curve) and standard deviation of the non-linear velocity, (vr, BORG − vr, linear), for each density bin (blue curve). We see here that the non-virial, non-linear component of the velocity field can represent a substantial portion of the cosmic velocity field: up to ∼25% of the standard deviation of the total velocity field when the local matter density exceeds δm ≳ 4.
In the text
thumbnail Fig. 9.

Posterior of the peculiar velocity of NGC 4993. The blue curve displays the large scale flow, vp, inferred from BORG while the black curve gives the required total velocity including the virial component within the halo. The posterior of the velocity used by the LIGO Scientific Collaboration and Virgo Collaboration (LVC) is shown with red dashes.
In the text
thumbnail Fig. 10.

Posterior of H0 for GW170817 is obtained using the peculiar velocity correction by our method and is shown in black when we use the GW data and in blue when using the GW and VLBI data. We also show the H0 measurement from the LIGO Scientific Collaboration and Virgo Collaboration (LVC) in dashed red line (Abbott et al. 2017a) and the magenta dashes presents the result from Hotokezaka et al. (2019) from hydrodynamics simulation jet using a GW+VLBI observation.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.