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/00046361/201936724  
Published online  08 February 2021 
Velocity correction for Hubble constant measurements from standard sirens
^{1}
Institut d’Astrophysique de Paris (IAP), CNRS & Sorbonne Université, UMR 7095, 98bis bd Arago, 75014 Paris, France
email: s.mukherjee@uva.nl, mukherje@iap.fr
^{2}
Institut Lagrange de Paris (ILP), Sorbonne Universités, 98bis Boulevard Arago, 75014 Paris, France
^{3}
Gravitation Astroparticle Physics Amsterdam (GRAPPA), Anton Pannekoek Institute for Astronomy and Institute for HighEnergy Physics, University of Amsterdam, Science Park 904, 1090 GL Amsterdam, The Netherlands
^{4}
The Oskar Klein Centre, Department of Physics, Stockholm University, AlbaNova University Centre, 106 91 Stockholm, Sweden
^{5}
Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA
^{6}
Imperial Centre for Inference and Cosmology (ICIC) & Astrophysics Group, Imperial College London, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK
^{7}
Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08540, USA
^{8}
Research Center for the Early Universe, Graduate School of Science, University of Tokyo, Bunkyoku, Tokyo 1130033, Japan
Received:
18
September
2019
Accepted:
24
November
2020
Gravitational wave (GW) sources are an excellent probe of the luminosity distance and offer a novel measure of the Hubble constant, H_{0}. This estimation of H_{0} from standard sirens requires an accurate estimation of the cosmological redshift of the host galaxy of the GW source after correcting for its peculiar velocity. The absence of an accurate peculiar velocity correction affects both the precision and accuracy of the measurement of H_{0}, particularly for nearby sources. Here, we propose a framework to incorporate such a peculiar velocity correction for GW sources. A first implementation of our method to the event GW170817, combined with observations taken with Very Large Baseline Interferometry (VLBI), leads to a revised value of H_{0} = 68.3_{−4.5}^{+4.6} km s^{−1} Mpc^{−1}. While this revision is minor, it demonstrates that our method makes it possible to obtain unbiased and accurate measurements of H_{0} at the precision required for the standard siren cosmology.
Key words: gravitational waves / cosmological parameters / galaxies: peculiar / cosmology: observations
© S. Mukherjee et al. 2021
Open 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 (H_{0}), have been ongoing since the discovery of the expanding Universe^{1} 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 H_{0} 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 lateUniverse probes using standard candles, such as supernovae (SN typeIa, 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 H_{0} from the CarnegieChicago Hubble program, based on the use of the tip of the red giant branch (TRGB) as a calibrator for the SN typeIa, 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 percentlevel precision needed to validate the low redshift (z) determination of H_{0} (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; RadburnSmith 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 H_{0} 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 nonlinear 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 SDSSIII samples, as discussed in Sect. 4. To make a reliable estimation of H_{0} 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 H_{0} 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 H_{0} 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 lowz measurements of H_{0} depend on measuring the luminosity distance to the source, which is given by
$$\begin{array}{c}\hfill {d}_{\mathrm{L}}=\frac{c(1+z)}{{H}_{0}}{\displaystyle {\int}_{0}^{z}\frac{\mathrm{d}z}{E(z)},}\end{array}$$(1)
in a homogeneous FriedmannLemaîtreRobertsonWalker (FLRW) metric when assuming a stationary source and observer. Here, c is the speed of light, and E(z)≡H(z)/H_{0}=$\sqrt{{\mathrm{\Omega}}_{\mathrm{m}}{(1+z)}^{3}+(1{\mathrm{\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 ${\rho}_{\text{c}}=3{H}_{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:
$$\begin{array}{c}\hfill {d}_{\mathrm{L}}=\frac{\mathit{cz}}{{H}_{0}}\xb7\end{array}$$(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, m_{1} and m_{2}, via the relation ℳ = (m_{1}m_{2})^{3/5}/(m_{1} + m_{2})^{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)
$$\begin{array}{cc}& {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}_{\mathrm{L}}}(1+{cos}^{2}(i)){e}^{i{\varphi}_{z}},\hfill \\ & {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}_{\mathrm{L}}}cos(i){e}^{i{\varphi}_{z}+i\pi /2},\hfill \end{array}$$(3)
where G is the gravitational constant, c is the speed of light, f_{z} = f/(1 + z) redshifted GW frequency, ϕ_{z} is the phase of the GW signal, and $i\equiv arccos(\widehat{L}.\widehat{n})$ is the inclination angle of the GW source which is defined as the angle between the unit vector of angular momentum, $\widehat{L}$, and the line of sight of the source position, $\widehat{n}$. The above expression indicates the presence of a degeneracy between d_{L} 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 holeneutron 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}/λ_{e}^{2}. As a result, by using d_{L} 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 largescale 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 perturbations^{3}. The presence of perturbations in the matter density leads to temporal and spatial fluctuations in the metric perturbations (through terms involving $\dot{\mathrm{\Phi}},\phantom{\rule{0.166667em}{0ex}}\ddot{\mathrm{\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 SachsWolfe, 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, v_{o}, and the source, v_{s}, 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 z_{obs} can be written in terms of the peculiar velocity, ${\mathit{v}}_{\mathrm{p}}=({\mathit{v}}_{\mathrm{s}}{\mathit{v}}_{\mathrm{o}}).\widehat{u}$^{4}, via the relation $(1+{z}_{\mathrm{obs}})=(1+z)(1+\frac{{\mathit{v}}_{\mathrm{p}}}{c})$. The corresponding modification in Eq. (2) is
$$\begin{array}{c}\hfill {d}_{\mathrm{L}}=\frac{cz+{v}_{\mathrm{p}}}{{H}_{0}}\xb7\end{array}$$(4)
This implies that the contribution from the peculiar velocity leads to a bias in the inferred value of H_{0}, 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 H_{0}. If we assume the distribution of the peculiar velocity field to be Gaussian with a variance ${\sigma}_{v}^{2}$, then the corresponding excess variance in the H_{0} measurement for a source at distance d_{L} becomes ${\sigma}_{v}^{2}/{d}_{\text{L}}^{2}$. The current framework for obtaining H_{0} 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, v_{p} = v_{h} + v_{vir}, can arise from two components, namely, (i) the motion of the halo due to the spatial gradient in the gravitational potential, v_{h}, and (ii) the virial velocity component, v_{vir}, of the host galaxy inside the halo. The nonrotational velocity component of v_{l} can be obtained from linear perturbation theory as
$$\begin{array}{c}\hfill {\mathit{v}}_{\mathrm{l}}(\mathit{x},z)=\frac{2}{3}\frac{g(z)}{aH(z)\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{\mathrm{m}}(z)}\phantom{\rule{0.166667em}{0ex}}{\mathbf{\nabla}}_{x}\mathrm{\Phi},\end{array}$$(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, nonlinear effects can be quite important and the previous relation becomes inaccurate. The virial velocity component, v_{vir}, can be related to the mass of the halo, M_{h}, and the distance to halo center, r, via the relation ${v}_{\text{vir}}^{2}\propto {M}_{\text{h}}/r$, which holds for a virialized system. For a system with a constant mass density, the halo mass and r are related by M_{h} ∝ r^{3}. As a result, the virial velocity v_{vir} is only related to the mass of the halo by the relation ${\mathit{v}}_{\mathrm{vir}}^{2}\propto {M}_{\mathrm{h}}^{2/3}$. 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 nonlinear velocity component (Sheth & Diaferio 2001):
$$\begin{array}{c}\hfill {\sigma}_{\mathrm{vir}}=476\phantom{\rule{0.166667em}{0ex}}{g}_{v}\phantom{\rule{0.166667em}{0ex}}{({\mathrm{\Delta}}_{\mathrm{nl}}(z)E{(z)}^{2})}^{1/6}\phantom{\rule{0.166667em}{0ex}}{({M}_{\mathrm{h}}/{10}^{15}\phantom{\rule{0.166667em}{0ex}}{M}_{\odot}\phantom{\rule{0.166667em}{0ex}}{\mathrm{h}}^{1})}^{1/3},\end{array}$$(6)
where g_{v} = 0.9, Δ_{nl}(z) = 18π^{2} + 60x − 32x^{2}, and x = Ω_{m}(1 + z)^{3}/E^{2}(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 largescale structure, as currently applied to the 2M++ compilation (Lavaux & Hudson 2011; Jasche & Lavaux 2019) and the SDSSIII/BOSS survey^{5} (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 threedimensional 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 largescale 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 threeparameter function motivated by the analysis of Nbody simulations (Neyrinck et al. 2014). We note that this bias relation is a local nonlinear 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 particlemesh Nbody 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 ≃256^{3}, 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 SimplexinCell estimator (SIC, Hahn et al. 2015; Leclercq et al. 2017). The SIC estimator relies on a phasespace interpolation technique, which provides an accurate estimate of velocity fields without the need for commonly used but arbitrary kernel smoothing procedures, such as cloudincell (CIC). It provides a physical picture for velocity fields even in lowdensity regions, which are poorlysampled by dark matter particles, and takes into account the multivalued nature of the velocity field in shellcrossed 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 luminosityweighted 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++ dataset 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 crosscorrelation 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 crosscorrelations between the inferred 2M++ fields and the ones for random realizations. This provides evidence that the density decorrelates from the velocity divergence in the way that is expected from the nonlinear 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 nonlinearities in the power spectrum as pointed out in Jennings et al. (2011).
Fig. 1. Display of the correlation coefficient between the nonlinear density and the negative divergence of the nonlinear velocity field, θ ≡ −∇ ⋅ v. The red line is the ensemble average over the posterior distribution of the crosscorrelation and the shaded dark red area is the one standard deviation according to the mean of that crosscorrelation 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 crosscorrelation goes to zero at small scales (high k value). 
Fig. 2. Ensemble average powerspectrum of the nonlinear 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 Nbody solver, but starting from a random sample of a Gaussian distribution with a powerspectrum 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 nonlinearities. Intermediate scales (k ≳ 0.1 h Mpc^{−1}, corresponding to ∼60 h Mpc^{−1}) get nonnegligible contributions from nonlinear 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 BORGderived nonlinear 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 selfconsistently models redshiftspace distortions, which brings additional information with which to disentangle galaxy bias from dark matter clustering.
4.1.2. Tests with an Nbody simulation
Although a full mock catalog analysis is beyond the scope of this article, we ran some tests on Nbody simulations to showcase the correctness of the BORG reconstruction algorithm. Using Gadget2 (Springel 2005), we generated an Nbody simulation covering a cubic volume with a side length of 200 h Mpc^{−1} and sampled with 256^{3} 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 BORGPM reconstruction with 20 timesteps from z = 50 to z = 0 linearly distributed according to the scale factor, with a grid of 32^{3} 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 Nbody simulation used for the test (topleft panel), as well as the result of the nearestgrid point assignment of the objects of the mock catalog in that same slice (topright panel). We also show the density contrasts (bottomleft panel) computed from a single BORG sample, after the burnin phase. We conclude with a difference plot highlighting the difference between the density field shown in the topleft panel and the bottomleft 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 crosscorrelation 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.
Fig. 3. Tests of the BORG algorithm on a Nbody simulation. Topleft panel: center slice of the density contrast computed at a resolution of 32^{3} from the Nbody cosmological simulation used in our test. The simulation itself is obtained with 256^{3} dark matter particles with ΛCDM cosmology as used in Jasche & Lavaux (2019). Topright panel: central slice of grid obtained by assigning objects of the mock catalog to a 32^{3} grid with the nearest grid point filter. Bottomleft 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. Bottomright 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. 
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 SimplexInCell algorithm at high resolution and downgraded at 32^{3} 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 unbiasedness, 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 Nbody simulations. The first three panels (topleft, topright, lowerleft) give the scatter plots for each component of the velocities for each mock tracer. The lowerleft panel shows the histogram of the residual velocities for the three components together.
Fig. 5. Tests of the BORGPM algorithm on a Nbody simulation and result for the velocity field. We compute the velocity field using the SimplexInCell algorithm using the entire set of particles of the Nbody 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. 
Fig. 6. Components (xtopleft, ytopright, 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 Nbody simulation. The points are colorcoded 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. Lowerright 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 SDSSIII/BOSS survey. The latter uses a simpler dynamical model based on firstorder Lagrangian perturbation theory (Lavaux et al. 2019). The contribution of the nonlinear velocity component captured by BORG is shown in Fig. 8. This nonlinear 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 nonlinear. It corresponds to the intermediate “gray” regime when linear modeling of the velocity field is not sufficient and still not describable without a full nonperturbative description of the nonlinear 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 (v_{r, BORG} − v_{r, linear})/v_{r, 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 nonlinear velocity field, both the total velocity and the residue, obtained by taking the difference between the total nonlinear 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 nonlinear component of the velocity field can go up to ∼(20−30)% at δ_{m} ∼ 4, indicating that the nonlinear component is far from negligible.
Fig. 7. Global properties of the velocity field that is available from the BORG inference framework applied on the 2M++ compilation and the SDSSIII/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 SDSSIII/BOSS survey. The part at the center is the region from 2M++. 
Fig. 8. Nonlinearities (NL) modeled by the BORG dynamical model. Left panel: relative difference between the velocity field derived from the gravity field assuming linearity and the nonlinear 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 nonlinear velocity, (v_{r, BORG} − v_{r, linear}), for each density bin (blue curve). We see here that the nonvirial, nonlinear 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 H_{0} = 100 km s^{−1} Mpc^{−1}. In this coordinate system, the value of H_{0} 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 H_{0} enters is in the prior power spectrum of primordial fluctuations that was used to run the analysis. This is, however, at least a secondorder 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 d_{L}, sky location $\widehat{n}$, polarization angle ψ, and the inclination angle i. It is this latter set that is useful to estimate the value of H_{0} (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 largescale flow, v_{h}, from BORG, and the small scale motion as a stochastic velocity dispersion using the fitting form given in Eq. (6). We note that v_{h} is itself composed of the velocity field at the linear order in perturbation, with an additional (at least) 30% correction from nonlinear 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, v_{h}, 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 nonlinear velocity field v_{vir} is assumed to have undergone a Gaussian distribution. The combined posterior of the peculiar velocity, v_{p} = v_{h} + v_{vir}, can be obtained by convolving the BORG posterior of v_{h} with the PDF of v_{vir}.
By combining the measurements from GW data, d_{GW}, the redshift of the host galaxy, $\widehat{z}$, and the peculiar velocity estimate, v_{p}, and luminosity distance, d_{L}, we can obtain the posterior of H_{0} for n GW sources by using the Bayesian framework:
$$\begin{array}{cc}\hfill P({H}_{0}\phantom{\rule{0.166667em}{0ex}}\{{d}_{\mathrm{GW}}\},\phantom{\rule{0.166667em}{0ex}}\{\widehat{z}\})\propto & {\displaystyle \prod _{i=1}^{\mathit{n}}\int \mathrm{d}{d}_{\mathrm{L}}^{i}\phantom{\rule{0.166667em}{0ex}}\mathrm{d}{v}_{\mathrm{p}}^{i}\phantom{\rule{4pt}{0ex}}\mathcal{L}({d}_{\mathrm{L}}^{i}{H}_{0},{v}_{\mathrm{p}}^{i},{z}_{i},{\widehat{\mathit{u}}}_{i},{d}_{\mathrm{GW}}^{i})}\hfill \\ \hfill & \times \mathcal{P}({z}_{i}{\widehat{\mathit{u}}}_{i})\phantom{\rule{0.166667em}{0ex}}\mathcal{P}({v}_{\mathrm{p}}^{i}M,{\widehat{\mathit{u}}}_{i})\phantom{\rule{4pt}{0ex}}\mathrm{\Pi}({H}_{0}),\hfill \end{array}$$(7)
where ${\widehat{\mathit{u}}}_{i}({\mathrm{RA}}_{i},\phantom{\rule{0.166667em}{0ex}}{\mathrm{Dec}}_{i})$ is the sky position^{6}, ℒ denotes the likelihood which is assumed to be Gaussian, $\mathcal{P}({z}_{i}{\widehat{\mathit{u}}}_{i})$ is the posterior of the redshift estimate at the GW source, ${\widehat{u}}_{i}$, $\mathcal{P}({\mathit{v}}_{\mathrm{p}}^{i}M,\phantom{\rule{0.166667em}{0ex}}{\widehat{\mathit{u}}}_{i})$ is the peculiar velocity estimate using the cosmological method, and Π(H_{0}) is the prior on the value of H_{0}.
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}_{\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, v_{h} and v_{vir}. We estimate the v_{h} component from BORG, and the nonlinear part of the velocity is obtained using the form mentioned in Eq. (6). NGC 4993 has a halo of mass about M_{h} ∼ 10^{12} 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 ${\sigma}_{\text{vir}}^{2}$ with a zero mean.
The posterior of the peculiar velocity from BORG is shown in blue in Fig. 9 for the Planck2015 bestfit cosmological parameters (Planck Collaboration XIII 2016). The PDF of the velocity field is nonGaussian. The mean value of the velocity field is ${\overline{\mathit{v}}}_{\mathrm{p}}=373$ km s^{−1} and the maximum a posterior value of the velocity field is ${\mathit{v}}_{\mathrm{p}}^{\mathrm{MP}}=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 v_{h} and v_{vir} is ${\sigma}_{\mathrm{t}}=\sqrt{{\sigma}_{\mathrm{h}}^{2}+{\sigma}_{\mathrm{vir}}^{2}}\sim 130$ km s^{−1}. The corresponding PDF for v_{p} 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.
Fig. 9. Posterior of the peculiar velocity of NGC 4993. The blue curve displays the large scale flow, v_{p}, 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 H_{0} from GW170817
Using our new, more precise estimate of the peculiar velocity for the NGC 4993, we obtain a revised value of H_{0} using the Bayesian framework given in Eq. (7) from the data of luminosity distance and inclination angle, assuming either a highspin or a lowspin prior on each compact object (Abbott et al. 2019a). Figure 10 shows the corresponding posterior of H_{0}. The posterior of H_{0} peaks at H_{0} = 68.3 km s^{−1} Mpc^{−1}, instead of ${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 H_{0} is due to the degeneracy between the inclination angle, i, and the distance, d_{L} (Abbott et al. 2017a). This acts as a limiting factor for the measurement of H_{0} from GW170817 if there is no independent measurement of inclination angle. Critically, the addition of Virgo (Acernese et al. 2015) to the nearly coaligned 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.
Fig. 10. Posterior of H_{0} 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 H_{0} 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 H_{0} 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 H_{0}. We implement our velocityfield 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\phantom{\rule{0.166667em}{0ex}}(\frac{{d}_{\mathrm{L}}}{41\phantom{\rule{0.166667em}{0ex}}\mathrm{Mpc}})\le 0.45$ rad. The revised H_{0} by our method is shown in Fig. 10 in blue and is compared with the H_{0} value from Hotokezaka et al. (2019) (shown by the magenta dashed line). The error in the measurement of H_{0} 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 H_{0} measurement with GW+VLBI from $70.{3}_{5.0}^{+5.3}$ km s^{−1} Mpc^{−1} to $68.{3}_{4.5}^{+4.6}$ 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 crosscorrelation 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 H_{0} 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 H_{0}, 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 signaltonoise 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 starblack 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 H_{0}. 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 H_{0}. Although averaging over large GW samples can lead to an unbiased estimate of H_{0}, the absence of a peculiar velocity correction will increase the error budget on H_{0} 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 N_{gw} to beat the variance as ${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 H_{0}.
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 nonlinear component. The largescale velocity flow is estimated using the Bayesian framework called BORG. The stochastic velocity dispersion of the source is obtained using a numerical fittingform 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 H_{0} 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 H_{0} to $68.{3}_{4.5}^{+4.6}$ km s^{−1} Mpc^{−1} for GW+VLBI, it slightly disfavors very low values of H_{0} ≲ 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}$ km s^{−1} Mpc^{−1}. The mean value of H_{0} differs from the values obtained in this work (${H}_{0}=68.{3}_{4.5}^{+4.6}$ km s^{−1} Mpc^{−1}) and also from the value obtained by Hotokezaka et al. (2019) (${H}_{0}=70.{3}_{5.0}^{+5.3}$ km s^{−1} Mpc^{−1}). However, the value by Howlett & Davis (2020) is consistent within the errorbars 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}$ 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 SDSSIII/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 SDSSIII/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 SDSSIV/eBOSS survey (Dawson et al. 2016). In the long term, with the availability of the nearly fullsky 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 lowredshift GW sources that are to be observed by the currently planned network of groundbased GW detectors (Schutz 2011; Abbott et al. 2018b).
The International Astronomical Union has recently renamed the Hubble law as the HubbleLemaître law, in recognition of the pioneering contribution of Lemaître (https://www.iau.org/news/pressreleases/detail/iau1812/).
SDSSIII/BOSS is the Baryon Oscillation Spectroscopic Survey.
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 ANR10LABX63) 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 ANR11IDEX000402. 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 ANR16CE230002 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 postprocessed. We thank Stéphane Rouberol for running smoothly this cluster for us. This work is done within the Aquila. Consortium (https://www.aquilaconsortium.org/). We acknowledge the use of following packages in this analysis: Astropy (Astropy Collaboration 2013, 2018), emcee: MCMC Hammer (ForemanMackey 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
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, Nature, 551, 85 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, Phys. Rev. Lett., 119, 161101 [Google Scholar]
 Abbott, T. M. C., Abdalla, F. B., Annis, J., et al. 2018a, MNRAS, 480, 3879 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D. et al. 2018b, Liv. Rev. Rel., 21, 3 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019a, Phys. Rev. X, 9, 011001 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019b, Phys. Rev. Lett., 123, 011102 [Google Scholar]
 Abramovici, A., Althouse, W. E., Drever, R. W. P., et al. 1992, Science, 256, 325 [Google Scholar]
 Acernese, F., Agathos, M., Agatsuma, K., et al. 2015, Class. Quant. Grav., 32, 024001 [Google Scholar]
 Addison, G. E., Watts, D. J., Bennett, C. L., et al. 2018, ApJ, 853, 119 [Google Scholar]
 Aghamousa, A., Aguilar, J., Ahlen, S., et al. 2016, ArXiv eprints [arXiv:1611.00036] [Google Scholar]
 Agrawal, P., CyrRacine, F.Y., Pinner, D., & Randall, L. 2019, ArXiv eprints [arXiv:1904.01016] [Google Scholar]
 Anderson, L., Aubourg, É., Bailey, S., et al. 2014, MNRAS, 441, 24 [Google Scholar]
 Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Astropy Collaboration (PriceWhelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
 Aubourg, É., Bailey, S., Bautista, J. E., et al. 2015, Phys. Rev. D, 92, 123516 [Google Scholar]
 Baker, T., Bellini, E., Ferreira, P. G., et al. 2017, Phys. Rev. Lett., 119, 251301 [Google Scholar]
 Barausse, E., Matarrese, S., & Riotto, A. 2005, Phys. Rev. D, 71, 063537 [Google Scholar]
 Belgacem, E., Dirian, Y., Foffa, S., & Maggiore, M. 2018a, Phys. Rev. D, 97, 104066 [Google Scholar]
 Belgacem, E., Dirian, Y., Foffa, S., & Maggiore, M. 2018b, Phys. Rev. D, 98, 023510 [Google Scholar]
 Belgacem, E., Calcagnib, G., Crisostomi, M., et al. 2019, JCAP, 1907, 024 [Google Scholar]
 Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, A20 [Google Scholar]
 Bernal, J. L., Verde, L., & Riess, A. G. 2016, JCAP, 1610, 019 [Google Scholar]
 Bonvin, C., Durrer, R., & Gasparini, M. A. 2006, Phys. Rev. D, 73, 023523; erratum: (2012) Phys. Rev. D, 85, 029901 [Google Scholar]
 Boruah, S. S., Hudson, M. J., Lavaux, G., et al. 2020, ArXiv eprints [arXiv:2010.01119] [Google Scholar]
 Branchini, E., Teodoro, L., Frenk, C. S., et al. 1999, MNRAS, 308, 1 [Google Scholar]
 Branchini, E., Freudling, W., Da Costa, L. N., et al. 2001, MNRAS, 326, 1191 [Google Scholar]
 Carrick, J., Turnbull, S. J., Lavaux, G., & Hudson, M. J. 2015, MNRAS, 450, 317 [Google Scholar]
 Chen, H.Y., Fishbach, M., & Holz, D. E. 2018, Nature, 562, 545 [Google Scholar]
 Chodorowski, M. J., & Ciecielag, P. 2002, MNRAS, 331, 133 [Google Scholar]
 Crook, A. C., Huchra, J. P., Martimbeau, N., et al. 2007, ApJ, 655, 790 [Google Scholar]
 Cuesta, A. J., VargasMagaña, M., Beutler, F., et al. 2016, MNRAS, 457, 1770 [Google Scholar]
 Cutler, C., & Flanagan, E. E. 1994, Phys. Rev. D, 49, 2658 [Google Scholar]
 Dalal, N., Holz, D. E., Hughes, S. A., & Jain, B. 2006, Phys. Rev. D, 74, 063006 [Google Scholar]
 Davis, M., Nusser, A., & Willick, J. A. 1996, ApJ, 473, 22 [Google Scholar]
 Davis, M., Nusser, A., Masters, K. L., et al. 2011, MNRAS, 413, 2906 [Google Scholar]
 Dawson, K. S., Schlegel, D. J., Ahn, C. P., et al. 2013, AJ, 145, 10 [Google Scholar]
 Dawson, K. S., Kneib, J.P., Percival, W. J., et al. 2016, AJ, 151, 44 [Google Scholar]
 Di Valentino, E., Melchiorri, A., & Mena, O. 2017, Phys. Rev. D, 96, 043503 [Google Scholar]
 Ebrová, I., & Bílek, M. 2020, A&A, 634, A73 [CrossRef] [EDP Sciences] [Google Scholar]
 Eisenstein, D. J., Weinberg, D. H., Agol, E., et al. 2011, AJ, 142, 72 [Google Scholar]
 Erdoǧdu, P., Lahav, O., Huchra, J. P., et al. 2006, MNRAS, 373, 45 [Google Scholar]
 Fairhurst, S. 2009, New J. Phys., 11, 123006 [Google Scholar]
 Fairhurst, S. 2011, Class. Quant. Grav., 28, 105021 [Google Scholar]
 Feeney, S. M., Peiris, H. V., Williamson, A. R., et al. 2019, Phys. Rev. Lett., 122, 061105 [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
 Freedman, W. L., Madore, B. F., Hatt, D., et al. 2019, ApJ, 882, 34 [Google Scholar]
 Hahn, O., Angulo, R. E., & Abel, T. 2015, MNRAS, 454, 3920 [Google Scholar]
 Hawking, S. W., & Israel, W. 1987, Three Hundred Years of Gravitation (Cambridge: Cambridge University Press) [Google Scholar]
 Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19 [Google Scholar]
 Holz, D. E., & Hughes, S. A. 2005, ApJ, 629, 15 [Google Scholar]
 Hotokezaka, K., Nakar, E., Gottlieb, O., et al. 2019, ArXiv eprints [arXiv:1806.10596] [Google Scholar]
 Howlett, C., & Davis, T. M. 2020, MNRAS, 492, 3803 [Google Scholar]
 Hubble, E. 1929, Proc. Natl. Acad. Sci., 15, 168 [Google Scholar]
 Hudson, M. J. 1994, MNRAS, 266, 475 [Google Scholar]
 Hudson, M. J., Smith, R. J., Lucey, J. R., & Branchini, E. 2004, MNRAS, 352, 61 [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Jasche, J., & Lavaux, G. 2019, A&A, 625, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jasche, J., & Wandelt, B. D. 2013, MNRAS, 432, 894 [Google Scholar]
 Jasche, J., Leclercq, F., & Wandelt, B. D. 2015, JCAP, 2015, 036 [Google Scholar]
 Jee, I., Suyu, S., Komatsu, E., et al. 2019, Science, 365, 1134 [Google Scholar]
 Jennings, E. 2012, MNRAS, 427, L25 [NASA ADS] [Google Scholar]
 Jennings, E., Baugh, C. M., & Pascoli, S. 2011, MNRAS, 410, 2081 [NASA ADS] [Google Scholar]
 Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open Source Scientific Tools for Python [Google Scholar]
 Kaiser, N., Efstathiou, G., Saunders, W., et al. 1991, MNRAS, 252, 1 [Google Scholar]
 Knox, L., & Millea, M. 2020, Phys. Rev. D, 101, 043533 [Google Scholar]
 Kolb, E. W., Matarrese, S., Notari, A., & Riotto, A. 2005, Phys. Rev. D, 71, 023524 [Google Scholar]
 Kreisch, C. D., CyrRacine, F.Y., & Doré, O. 2020, Phys. Rev. D, 101, 123505 [Google Scholar]
 Lavaux, G., & Hudson, M. J. 2011, MNRAS, 416, 2840 [Google Scholar]
 Lavaux, G., & Jasche, J. 2016, MNRAS, 455, 3169 [Google Scholar]
 Lavaux, G., Jasche, J., & Leclercq, F. 2019, MNRAS, submitted [Google Scholar]
 Leclercq, F., Jasche, J., & Wandelt, B. 2015, JCAP, 2015, 15 [Google Scholar]
 Leclercq, F., Jasche, J., Lavaux, G., Wandelt, B., & Percival, W. 2017, JCAP, 2017, 049 [Google Scholar]
 Lemaître, G. 1927, Ann. Soc. Sci. Bruxelles, 47, 49 [Google Scholar]
 Lemaître, G. 1931, MNRAS, 91, 483 [NASA ADS] [Google Scholar]
 Lin, M.X., Benevento, G., Hu, W., & Raveri, M. 2019, Phys. Rev. D, 100, 063542 [Google Scholar]
 Lombriser, L., & Lima, N. A. 2017, Phys. Lett. B, 765, 382 [Google Scholar]
 Lombriser, L., & Taylor, A. 2016, JCAP, 1603, 031 [Google Scholar]
 LSST Science Collaboration (Abell, P. A., et al.) 2009, ArXiv eprints [arXiv:0912.0201] [Google Scholar]
 Ma, Y.Z., Branchini, E., & Scott, D. 2012, MNRAS, 425, 2880 [Google Scholar]
 Macaulay, E., Nichol, R. C., Bacon, D., et al. 2019, MNRAS, 486, 2184 [Google Scholar]
 Maggiore, M. 2008, Gravitational Waves: Volume 1: Theory and Experiments, Gravitational Waves (Oxford: Oxford University Press) [Google Scholar]
 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]
 Mooley, K. P., Deller, A. T., Gottlieb, O., et al. 2018a, Nature, 561, 355 [Google Scholar]
 Mooley, K. P., Frail, D. A., Dobie, D., et al. 2018b, ApJ, 868, L11 [Google Scholar]
 Mortlock, D. J., Feeney, S. M., Peiris, H. V., Williamson, A. R., & Nissanke, S. M. 2019, Phys. Rev. D, 100, 103523 [Google Scholar]
 Mukherjee, S., & Wandelt, B. D. 2018, ArXiv eprints [arXiv:1808.06615] [Google Scholar]
 Mukherjee, S., Wandelt, B. D., & Silk, J. 2020a, Phys. Rev. D, 101, 103509 [Google Scholar]
 Mukherjee, S., Wandelt, B. D., & Silk, J. 2020b, MNRAS, 494, 1956 [Google Scholar]
 Neyrinck, M. C., AragónCalvo, M. A., Jeong, D., & Wang, X. 2014, MNRAS, 441, 646 [Google Scholar]
 Nicolaou, C., Lahav, O., Lemos, P., Hartley, W., & Braden, J. 2020, MNRAS, 495, 90 [Google Scholar]
 Nishizawa, A. 2018, Phys. Rev. D, 97, 104037 [Google Scholar]
 Nissanke, S., Holz, D. E., Hughes, S. A., Dalal, N., & Sievers, J. L. 2010, ApJ, 725, 496 [Google Scholar]
 Nissanke, S., Sievers, J., Dalal, N., & Holz, D. 2011, ApJ, 739, 99 [Google Scholar]
 Nissanke, S., Holz, D. E., Dalal, N., et al. 2013a, ArXiv eprints [arXiv:1307.2638] [Google Scholar]
 Nissanke, S., Kasliwal, M., & Georgieva, A. 2013b, ApJ, 767, 124 [Google Scholar]
 Nusser, A., Da Costa, L. N., Branchini, E., et al. 2001, MNRAS, 320, L21 [Google Scholar]
 Oguri, M. 2016, Phys. Rev. D, 93, 083511 [Google Scholar]
 Pan, Y. C., Kilpatrick, C. D., Simon, J. D., et al. 2017, ApJ, 848, L30 [Google Scholar]
 Pardo, K., Fishbach, M., Holz, D. E., & Spergel, D. N. 2018, JCAP, 1807, 048 [Google Scholar]
 Pérez, F., & Granger, B. E. 2007, Comput. Sci. Eng., 9, 21 [Google Scholar]
 Pike, R. W., & Hudson, M. J. 2005, ApJ, 635, 11 [Google Scholar]
 Planck Collaboration XVI. 2014, A&A, 571, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Poisson, E., & Will, C. M. 1995, Phys. Rev. D, 52, 848 [Google Scholar]
 Poulin, V., Smith, T. L., Karwal, T., & Kamionkowski, M. 2019, Phys. Rev. Lett., 122, 221301 [Google Scholar]
 RadburnSmith, D. J., Lucey, J. R., & Hudson, M. J. 2004, MNRAS, 355, 1378 [Google Scholar]
 Reid, M. J., Braatz, J. A., Condon, J. J., et al. 2009, ApJ, 695, 287 [Google Scholar]
 Riess, A. G., Casertano, S., Yuan, W., Macri, L. M., & Scolnic, D. 2019, ApJ, 876, 85 [Google Scholar]
 Sakstein, J., & Jain, B. 2017, Phys. Rev. Lett., 119, 251303 [Google Scholar]
 Saltas, I. D., Sawicki, I., Amendola, L., & Kunz, M. 2014, Phys. Rev. Lett., 113, 191101 [Google Scholar]
 Sasaki, M. 1987, MNRAS, 228, 653 [Google Scholar]
 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]
 Schutz, B. F. 1986, Nature, 323, 310 [Google Scholar]
 Schutz, B. F. 2011, Class. Quant. Grav., 28, 125023 [Google Scholar]
 Scoccimarro, R. 2004, Phys. Rev. D, 70, 083007 [Google Scholar]
 Seto, N., & Kyutoku, K. 2018, MNRAS, 475, 4133 [Google Scholar]
 Shaya, E. J., Tully, R. B., & Pierce, M. J. 1992, ApJ, 391, 16 [Google Scholar]
 Sheth, R. K., & Diaferio, A. 2001, MNRAS, 322, 901 [Google Scholar]
 Springel, V. 2005, MNRAS, 364, 1105 [Google Scholar]
 Springob, C. M., Masters, K. L., Haynes, M. P., Giovanelli, R., & Marinoni, C. 2007, ApJS, 172, 599 [Google Scholar]
 Springob, C. M., Magoulas, C., Colless, M., et al. 2014, MNRAS, 445, 2677 [Google Scholar]
 Turnbull, S. J., Hudson, M. J., Feldman, H. A., et al. 2012, MNRAS, 420, 447 [Google Scholar]
 van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22 [Google Scholar]
 Veitch, J., Mandel, I., Aylott, B., et al. 2012, Phys. Rev. D, 85, 104045 [Google Scholar]
 Verde, L., Protopapas, P., & Jimenez, R. 2013, Phys. Dark Univ., 2, 166 [Google Scholar]
 Verde, L., Treu, T., & Riess, A. G. 2019, Nat. Astron., 3, 891 [Google Scholar]
 Vitale, S., & Chen, H.Y. 2018, Phys. Rev. Lett., 121, 021303 [Google Scholar]
 Wong, K. C., Suyu, S. H., Chen, G. C.F., et al. 2019, MNRAS, 498, 1420 [Google Scholar]
 Yuan, W., Riess, A. G., Macri, L. M., Casertano, S., & Scolnic, D. 2019, ApJ, 886, 61 [Google Scholar]
All Figures
Fig. 1. Display of the correlation coefficient between the nonlinear density and the negative divergence of the nonlinear velocity field, θ ≡ −∇ ⋅ v. The red line is the ensemble average over the posterior distribution of the crosscorrelation and the shaded dark red area is the one standard deviation according to the mean of that crosscorrelation 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 crosscorrelation goes to zero at small scales (high k value). 

In the text 
Fig. 2. Ensemble average powerspectrum of the nonlinear 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 Nbody solver, but starting from a random sample of a Gaussian distribution with a powerspectrum 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 nonlinearities. Intermediate scales (k ≳ 0.1 h Mpc^{−1}, corresponding to ∼60 h Mpc^{−1}) get nonnegligible contributions from nonlinear 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 
Fig. 3. Tests of the BORG algorithm on a Nbody simulation. Topleft panel: center slice of the density contrast computed at a resolution of 32^{3} from the Nbody cosmological simulation used in our test. The simulation itself is obtained with 256^{3} dark matter particles with ΛCDM cosmology as used in Jasche & Lavaux (2019). Topright panel: central slice of grid obtained by assigning objects of the mock catalog to a 32^{3} grid with the nearest grid point filter. Bottomleft 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. Bottomright 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 
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 
Fig. 5. Tests of the BORGPM algorithm on a Nbody simulation and result for the velocity field. We compute the velocity field using the SimplexInCell algorithm using the entire set of particles of the Nbody 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 
Fig. 6. Components (xtopleft, ytopright, 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 Nbody simulation. The points are colorcoded 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. Lowerright 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 
Fig. 7. Global properties of the velocity field that is available from the BORG inference framework applied on the 2M++ compilation and the SDSSIII/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 SDSSIII/BOSS survey. The part at the center is the region from 2M++. 

In the text 
Fig. 8. Nonlinearities (NL) modeled by the BORG dynamical model. Left panel: relative difference between the velocity field derived from the gravity field assuming linearity and the nonlinear 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 nonlinear velocity, (v_{r, BORG} − v_{r, linear}), for each density bin (blue curve). We see here that the nonvirial, nonlinear 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 
Fig. 9. Posterior of the peculiar velocity of NGC 4993. The blue curve displays the large scale flow, v_{p}, 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 
Fig. 10. Posterior of H_{0} 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 H_{0} 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 (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.