The variance of the locally measured Hubble parameter explained with different estimators

TL;DR

This paper compares different estimators for the variance of local Hubble constant measurements, highlighting how their differences impact the interpretation of observational discrepancies.

Contribution

It clarifies how estimator choice affects variance calculations of Hubble constant measurements in perturbation theory and N-body simulations.

Findings

N-body simulation estimators give larger weights to distant sources.

Different estimators lead to different variance estimates.

Both approaches produce variances too small to explain the Hubble tension.

Abstract

We study the expected variance of measurements of the Hubble constant, $H_0$, as calculated in either linear perturbation theory or using non-linear velocity power spectra derived from $N$-body simulations. We compare the variance with that obtained by carrying out mock observations in the N-body simulations, and show that the estimator typically used for the local Hubble constant in studies based on perturbation theory is different from the one used in studies based on N-body simulations. The latter gives larger weight to distant sources, which explains why studies based on N-body simulations tend to obtain a smaller variance than that found from studies based on the power spectrum. Although both approaches result in a variance too small to explain the discrepancy between the value of $H_0$ from CMB measurements and the value measured in the local universe, these considerations are important in light of the percent determination of the Hubble constant in the local universe.

Figures

Tables

Table 1. Table 1 : This table gives the spread in each of the two estimates of the Hubble constant considered in this paper, H loc mean superscript subscript 𝐻 loc mean H_{\textrm{loc}}^{\textrm{mean}} ( top ) and H loc lsq superscript subscript 𝐻 loc lsq H_{\textrm{loc}}^{\textrm{lsq}} ( bottom ). In both cases, the spread is calculated based on three different distributions of sources, and using four different techniques: using mock observations with the observed sources chosen among the halos in the simulation; using mock observations with the sources chosen among the grid points onto which the velocity field has been interpolated; and using the velocity power spectrum (PS) obtained from respectively linear theory and from the N-body simulation.

	Top-hat ($R=256\text{Mpc/h}$)	CfA3+OLD	Union 2.1
———————— ${\sigma }_{\mathrm{mean}}$ ————————
Mock observations of halos	0.49%	1.86%	0.81%
Mock observations of grid points	0.47%	1.77%	0.77%
Linear velocity PS	0.41%	1.72%	0.71%
N-body velocity PS	0.40%	1.62%	0.68%
———————— ${\sigma }_{\mathrm{lsq}}$ ————————
Mock observations of halos	0.39%	0.96%	0.30%
Mock observations of grid points	0.39%	0.93%	0.27%
Linear velocity PS	0.34%	0.86%	0.23%
N-body velocity PS	0.33%	0.82%	0.22%

Equations

$$H_0^{\mathrm{Planck}+\mathrm{BAO}}=67.6\pm 0.6\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$$

$$H_0^{\mathrm{local}}=73.24\pm 1.74\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}.$$

$$H_{\text{loc}}^{\text{mean}}\equiv ⟨v_r/r⟩.$$

$$H_{\text{loc}}^{\text{lsq}}\equiv ⟨v_{r}r⟩/⟨r^2⟩.$$

$$H_{\text{loc}}^{\text{lsq}}\equiv \frac{⟨\frac{v_r}{r}{r}^2⟩}{⟨r^2⟩},$$

$$H_{\text{loc}}^{\text{mean}}=H_0+\frac{1}{N}\sum _i\frac{{\text{v}}_i\cdot ({\text{r}}_i-{\text{r}}_0)}{|{\text{r}}_i-{\text{r}}_0{|}^2}.$$

$$H_{\text{loc}}^{\text{mean}}=H_0+\int d^3\text{r}\frac{\text{v}(\text{r})\cdot (\text{r}-{\text{r}}_0)}{|\text{r}-{\text{r}}_0{|}^2}W(\text{r}-{\text{r}}_0).$$

$$\text{v}(\text{r})=\frac{1}{(2\pi {)}^3}\int d^3\text{k}{\text{v}}_{\text{k}}{e}^{-i\text{k}\cdot \text{r}},$$

$${\text{v}}_{\text{k}}=\frac{f{H}_0}{ik}{\delta }_{\text{k}}\hat{\text{k}},$$

$$\frac{H_{\text{loc}}^{\text{mean}}-H_0}{H_0}=\frac{f}{(2\pi {)}^3}\int d^3\text{k}{\delta }_{\text{k}}{e}^{-i\text{k}\cdot {\text{r}}_0}\int d^3\text{r}\frac{\text{k}}{i{k}^2}\cdot \frac{\text{r}-{\text{r}}_0}{|\text{r}-{\text{r}}_0{|}^2}W(\text{r}-{\text{r}}_0)e^{-i\text{k}\cdot (\text{r}-{\text{r}}_0)}.$$

$$W_{\text{top-hat}}(\text{x})=\{\begin{array}{ll}\frac{1}{4/3\pi R^3},& |\text{x}|\le R,\\ 0,& |\text{x}|>R.\end{array}$$

$${\mathcal{L}}_{\text{top-hat}}^{\text{mean}}(kR)$$

$$=\frac{1}{ik}{\int }_0^{R}dx\frac{3x}{4\pi R^3}{\int }_0^{2\pi }d\varphi {\int }_0^{\pi }d\theta \sin (\theta )\cos (\theta )e^{-ikx\cos (\theta )}$$

$$=\frac{3}{(kR{)}^3}\left(\sin (kR)-{\int }_0^{kR}dx\frac{\sin (x)}{x}\right).$$

$${\sigma }_{\text{mean}}^2$$

$$=\frac{f^2}{2{\pi }^2{R}^2}{\int }_0^{\infty }dkP(k){\left[kR{\mathcal{L}}^{\text{mean}}(kR)\right]}^2$$

$${\sigma }_{\text{mean}}^2$$

$$H_{\text{loc}}^{\text{lsq}}=H_0+\frac{\int d^3\text{r}\text{v}(\text{r})\cdot (\text{r}-{\text{r}}_0)W(\text{r}-{\text{r}}_0)}{\int d^3\mathbf{r}|\mathbf{r}-{\mathbf{r}}_0{|}^{2}W(\mathbf{r}-{\mathbf{r}}_0)}.$$

$${\mathcal{L}}^{\text{lsq}}(kR)\equiv \frac{\int d^3\text{x}\text{k}\cdot \text{x}W(x)e^{-i\text{k}\cdot \text{x}}}{i{k}^2\int d^3\mathbf{x}{x}^{2}W(x)}.$$

$$P_v(k)=⟨\mathbf{v}(\mathbf{k})\cdot {\mathbf{v}}^{\ast }(\mathbf{k})⟩,$$

$$\mathbf{v}(\mathbf{k})\cdot {\mathbf{v}}^{\ast }(\mathbf{k})=$$

$$P_v=k^{-2}\left(P_{\theta }+P_w\right).$$

$$\theta (\mathbf{k})=\int d^3\text{r}\nabla \cdot \mathbf{v}(\mathbf{r})e^{-i\mathbf{k}\cdot \mathbf{r}}=i\mathbf{k}\cdot \mathbf{v}(\mathbf{k})=\delta (\mathbf{k})H_{0}f,$$

Full text

1 Introduction

The Hubble parameter, $H$, measures the expansion rate of the Universe through the time derivative of the scale factor, $H\equiv\frac{\dot{a}}{a}$. Its current value, $H_{0}$, is one of the most fundamental cosmological observables. Historically the measurement of $H_{0}$ has been of great importance in e.g. establishing the age of the Universe. However, in recent years precision observations of the Cosmic Microwave Background (CMB) and large scale structure (LSS) have allowed for an indirect determination of $H_{0}$ with a precision exceeding that of current direct measurements. Intriguingly the value of $H_{0}$ inferred from CMB and LSS observations assuming a flat $\Lambda$CDM model [1]

[TABLE]

is significantly lower than the value obtained from direct measurement in the local universe [2]

[TABLE]

There have been several suggestions for the origin of this discrepancy. It could be caused by systematic issues in either the value inferred from the CMB [3, 4] or in the local measurements [5, 6] – or both. However, as the considerations raised in these papers have been reexamined, and as more data is coming in, the discrepancy does not seem to be getting smaller [7, 8, 2]. Since the value of $H_{0}$ inferred from CMB and LSS observations is inferred using $\Lambda$CDM-like models, this could be an indication that the true cosmological model is different from $\Lambda$CDM.

An alternative possibility which has been discussed extensively is that the local universe is described by different parameters than the global ones. For example the presence of a large local underdensity – a so-called Hubble bubble – could lead to a significantly higher local value of $H_{0}$ relative to the global value [9, 10, 11, 12].

A very interesting question is whether the variance predicted in local measurements of $H_{0}$ by $\Lambda$CDM models could actually account for the discrepancy. This has been studied in the literature several times using different approaches, resulting in somewhat differing conclusions. Calculations based on perturbation theory, such as [13, 14, 15, 16, 17, 18], point to a variance in the local Hubble parameter of 2-3%, for distributions of sources similar to the one used in the most recent measurements of the local Hubble parameter within $z=0.1$. In contrast, studies based on N-body simulations result in a variance of less than 1% in the same redshift range [19, 20, 21]. This discrepancy between the variance deduced from perturbation theory and from N-body simulations is the subject of this paper.

In this work, we estimate the variance in the local Hubble constant using both mock observations in a large N-body simulation and perturbation theory, and compare the results. In order to bridge the gap between the two methods, we also use a hybrid approach in which the variance is determined with perturbation theory but using the velocity power spectrum from the N-body simulation. In section 2, we describe the simulation used for the study, and explain how the mock observations are carried out. We also introduce two different estimators for $H_{0}$ – as we will see, a large part of the difference between the results from studies based on respectively N-body simulations and perturbation theory stems from the fact that they use two different estimators. In section 3, we derive the equations which describe the variance in $H_{0}$ in linear perturbation theory for each of the two estimators. In section 4, we present the velocity power spectrum from the N-body simulation. Our results are presented in section 5, and finally we conclude in section 6.

2 Mock observations in N-body simulations

The most precise measurements of the local value of $H_{0}$ are based on measurements of the distances to galaxies, inferred from the observed luminosities of type Ia supernovae, together with the spectroscopic redshifts of their host galaxies which provide the radial velocities [22, 2]. An intuitive way to mimic such measurements is by carrying out mock observations in synthetic data sets, which can for example be generated from the halo catalogs from N-body simulations. In this section, we will describe such an analysis.

2.1 Simulation

The mock observations are based on an N-body simulation with cosmological parameters given by $(\Omega_{b},\Omega_{m})=(0.048,0.26)$, $(h,\sigma_{8})=(0.68,0.85)$ and a spectral index of $n_{s}=1$. The simulation is performed using a modified version of the GADGET-2 code [23], with initial conditions generated using a code described in [24] based on transfer functions computed using CAMB111http://camb.info/ [25]. The simulation is run in a box of side length $2048\,\,\text{Mpc/h}$, containing $1024^{3}$ dark matter particles. It is initiated at a redshift of $z=49$ and run until the present time. In this paper, we only consider the snapshot at $z=0$. We use the halo finder ROCKSTAR [26] to generate the halo catalog from which sources for the local determination of $H_{0}$ will be selected.

2.2 Observers and observations

We carry out mock observations for 1,000 observers distributed throughout the simulation volume. Three different distance-distributions of sources are used for estimating the Hubble constant. The first two are chosen to be similar to the distributions used in actual observations; these are the CfA3+OLD collection of supernovae, which is close to the one used in [22], and the Union 2.1 collection, which covers the same redshift range as the one used in [2]. In order to obtain source distributions correponding to each of these for the observers in the simulation, the volume around each observer is divided into bins of equal width – as illustrated in figure 1 – and within each bin the number of sources dictated by the distance distribution is chosen among the halos in the bin. We also consider a third distribution, in which the sources are chosen randomly among all the halos in the simulation within a radius of $256\,\,\text{Mpc/h}$ of the individual observers; this corresponds to a top-hat distribution, and it peaks at high redshift, since the volumes of concentric shells grow with distance. Note that no measurements of the local value of the Hubble constant exists for such a distribution.

Observers are placed at random positions throughout the simulation volume. In e.g. [19, 21], it was shown that the observer position causes a bias in the measured Hubble constant, since observers placed in massive halos will tend to measure a low value of $H_{0}$ due to local inflow, whereas observers placed at random positions – which will tend to be located in voids as these take up a larger part of the simulation volume – tend to measure a value of $H_{0}$ slightly higher than the true value. However, even though the choice of the observer positions causes a bias in the measured value of the Hubble constant, we find that it does not have a significant effect on its spread for any of the redshift distributions of sources considered in this study. That is, placing the observers in massive halos or at random positions results in approximately the same variance in the locally measured Hubble constant.

When the sources for the mock observations are chosen among the halos in the simulation, overdense regions are sampled more than voids. In order to determine how this affects the measured variance, we also carry out a set of mock observations in which the sources are chosen randomly from a regular grid onto which the peculiar velocities in the simulation has been interpolated.

2.3 Estimating the Hubble constant

For each observer, we determine the Hubble constant using two different approaches which result in slightly different values for $H_{0}$. In one, the local Hubble constant is obtained by calculating $v_{r}/r$ for each of the observed sources, and then taking the mean of all these values, i.e.

[TABLE]

A different calculation of the local Hubble constant is obtained by using the least-squares estimator for the slope of the relationship $v_{r}=H_{0}r$. This is given by

[TABLE]

We can compare this to $H_{\textrm{loc}}^{\textrm{mean}}$ by rewriting $H_{\textrm{loc}}^{\textrm{lsq}}$ as

[TABLE]

which shows that the least squares estimate is equivalent to the mean of the individual slopes with each of these weighted by the square of the distance to the source.

The expectation values for $H_{\textrm{loc}}^{\textrm{mean}}$ and $H_{\textrm{loc}}^{\textrm{lsq}}$, i.e. the values obtained by taking the mean for a large number of observers in the simulation, are not exactly the same. Since $H_{\textrm{loc}}^{\textrm{mean}}$ gives relatively more weight to nearby sources, local flows affect this estimator more than $H_{\textrm{loc}}^{\textrm{lsq}}$. The previously mentioned tendency for the observers to measure a value of the local Hubble constant which is higher than the true value, due to the fact that they tend to be located in underdense regions, is therefore somewhat more pronounced for $H_{\textrm{loc}}^{\textrm{mean}}$ than for $H_{\textrm{loc}}^{\textrm{lsq}}$. The effect is largest for the CfA3+OLD distribution of sources, in the case where the sources are chosen among the halos in the simulation; in this case, the average values of the Hubble constant among the observers in the simulation for each of the two estimators are $\left\langle H_{\textrm{loc}}^{\textrm{mean}}\right\rangle=100.41\,\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}\,{\rm h}$ and $\left\langle H_{\textrm{loc}}^{\textrm{lsq}}\right\rangle=100.16\,\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}\,{\rm h}$.

In figure 1, the mock observation and subsequent determination of the Hubble constant is illustrated for an example observer. From the Hubble diagram for the specific observer it is difficult to make out the difference between the two estimates. To amplify the difference, we therefore also plot $v_{r}-H_{0}r$ against $r$.

3 Estimating the variance in $H_{0}$ from the velocity power spectrum

The expected variance in the Hubble constant caused by over- and underdensitites in the Universe, and the resulting peculiar velocities, can be estimated from the power spectrum of density fluctuations. Below, we derive the expressions for the local variance of both $H_{\textrm{loc}}^{\text{mean}}$ and $H_{\textrm{loc}}^{\textrm{lsq}}$. We first derive the relationship between the velocity power spectrum and $H_{\textrm{loc}}^{\mathrm{mean}}$, as this is the simpler of the two cases.

3.1 The local variation in $H_{\textrm{loc}}^{\textrm{mean}}$

Consider an observer at a comoving position $\textbf{r}_{0}$, making a local estimate of the Hubble constant based on observations of $N$ galaxies at positions $\textbf{r}_{i}$ and with peculiar velocities $\textbf{v}_{i}$. In this section it is assumed that the observer estimates the Hubble constant from the observed sources by taking the average of $v_{r}^{i}/r^{i}$, where $i$ indicates the individual sources - this is the approach used in [27, 13, 14], and it is similar to that used in [16, 17] (which, however, also incorporates non-linear and relativistic effects to a certain degree). This is equivalent to

[TABLE]

We can replace the individual velocities $\textbf{v}_{i}$ by the velocity field, $\textbf{v}(\textbf{r})$, and the sum by an integral over a window function $W(\textbf{r}-\textbf{r}_{0})$, which describes the distribution of sources:

[TABLE]

The velocity field is generated by inhomogeneities in the density field, $\rho(\textbf{r})$, typically expressed through the over-density parameter $\delta(\textbf{r})=(\rho(\textbf{r})-\rho_{0})/\rho_{0}$, where $\rho_{0}$ is the mean density. The velocity and the over-density fields can be expressed in terms of their Fourier transforms, i.e.

[TABLE]

and likewise for $\delta(\textbf{r})$. In linear perturbation theory, the relationship between these fields is [28, 29]

[TABLE]

where $f$ is the growth function, which we approximate as $f=\Omega_{m}^{0.6}$, and $\hat{\textbf{k}}$ is a unit vector pointing in the direction of k. Inserting this in the integral of equation 3.2, it is found that the fractional deviation of the local Hubble constant from $H_{0}$ is given by

[TABLE]

The last integral represents the window function determining how density perturbations of different sizes affect the local Hubble constant. As a simple example, consider the case where $W$ is a top-hat window with radius $R$, i.e.

[TABLE]

In this case, the window function reduces to a relatively simple expression in Fourier-space [27, 13]. By making the substitution $\textbf{x}=\textbf{r}-\textbf{r}_{0}$, and choosing the coordinate system with the $z$-axis parallel to k, such that $\textbf{k}\cdot\textbf{x}=kx\cos(\theta)$, the integral can be evaluated as:

[TABLE]

For more realistic redshift distributions, such as those describing the CfA3+OLD and Union 2.1 collections of type Ia supernovae, an analytic expression for the window function cannot be obtained, and the integral has to be evaluated numerically.

In linear perturbation theory, the density field $\delta_{k}$, and consequently the deviations in the Hubble constant it produces, is a Gaussian random variable with mean [math] and a variance, $\sigma^{2}_{\textrm{mean}}\equiv\left\langle\left(\frac{H^{\textrm{mean}}_{\textrm{loc}}-H_{0}}{H_{0}}\right)^{2}\right\rangle$, that can be determined from equation 3.5 as

[TABLE]

where we have exploited that $\int d^{3}\textbf{x}e^{-i\textbf{k}\cdot\textbf{x}}=(2\pi)^{3}\delta^{3}_{\textrm{D}}(\textbf{k})$, with $\delta^{3}_{\textrm{D}}$ being the 3-dimensional Dirac delta function, and $P(k)=\langle|\delta_{\textbf{k}}|^{2}\rangle$ the matter power spectrum. It follows from equation 3.4 that the velocity power spectrum is related to the matter power spectrum as $P_{v}(k)=f^{2}H_{0}^{2}P(k)/k^{2}$. Hence, equation 3.11 for the variations in the Hubble constant can also be written in terms of the velocity power spectrum:

[TABLE]

In figure 3, we show the window function, $\left[k^{2}R\mathcal{L}^{\textrm{mean}}(kR)\right]^{2}$ corresponding to the top-hat, the CfA3+OLD and the Union 2.1 redshift distributions.

Based on the velocity power spectrum from CAMB, we calculate and plot the variance in the Hubble constant calculated as in equation 3.12. The resulting 68%, 95% and 99% confidence regions as a function of the radius, $R$, of the top-hat window function are shown in figure 2. We compare it to the equivalent result from the mock observations in the N-body simulation.

3.2 The local variation in $H_{\textrm{loc}}^{\textrm{lsq}}$

The calculation of the window function determining the local variations in the Hubble constant can be repeated for the least squares estimate of the Hubble constant, as given in equation 2.3, by replacing equation 3.2 by

[TABLE]

This leads to a slightly different set of window functions, $\mathcal{L}^{\textrm{lsq}}$, given by

[TABLE]

The variance in the least squares estimates of $H_{0}$, $\sigma^{2}_{\textrm{lsq}}\equiv\left\langle\left(\frac{H^{\textrm{lsq}}_{\textrm{loc}}-H_{0}}{H_{0}}\right)^{2}\right\rangle$, can then be found by simply replacing $\mathcal{L}^{\textrm{mean}}$ by $\mathcal{L}^{\textrm{lsq}}$ in equation 3.12. We show the corresponding window functions, $\left[k^{2}R\mathcal{L}^{\textrm{lsq}}(kR)\right]^{2}$, in figure 3. Using the least squares estimates moves the window functions towards lower $k$-values, corresponding to larger scales, due to the larger weight given to distant sources. The top-hat window, in which the more distant sources already contribute much more than nearby sources, is almost unchanged.

We note that for each of the considered distributions, only modes with $0.005\,\,\text{h/Mpc}\lesssim k\lesssim 0.1\,\,\text{h/Mpc}$ – which are well resolved in the simulation used for this study – contribute significantly to the variance in the Hubble constant.

4 A hybrid approach: Velocity power spectrum from N-body simulations

We can replace the linear prediction for the velocity power spectrum by that measured in the N-body simulation in which the mock observations are carried out. By using the velocity power spectrum measured from the N-body simulation in the calculation of the variance in $H_{\textrm{loc}}^{\textrm{mean}}$ and $H_{\textrm{loc}}^{\textrm{lsq}}$, and comparing with the variance predicted by linear theory, we can determine how much the variance is affected by the non-linear evolution of the velocity field. In contrast to the matter power spectrum, the velocity power spectrum decreases when non-linear effects of clustering are included. Therefore, the variance measured in N-body simulations can be expected to be lower than that determined in linear perturbation theory. However, this effect is only significant at relatively small scales ($k\gtrsim 0.1\,\,\text{h/Mpc}$), which only contribute very little to the variance in the measured Hubble constant for each of the distributions of sources studied here.

4.1 The non-linear velocity power spectrum

From the velocity field we construct its divergence, $\theta=\nabla\cdot\mathbf{v}$, and curl, $\mathbf{w}=\nabla\times\mathbf{v}$. The velocity field can then be written as $\mathbf{v}=\nabla^{-2}\nabla\theta+(\nabla\times)^{-2}\nabla\times\mathbf{w}$. Fourier transforming we get $\mathbf{v}=\frac{i}{k^{2}}(-\mathbf{k}\theta+\mathbf{k}\times\mathbf{w})$. The velocity power spectrum is then given by

[TABLE]

with

[TABLE]

The second and third terms vanish through cyclic permutation, and the last term can be written as $k^{2}\mathbf{w}\cdot\mathbf{w}^{*}-(\mathbf{k}\cdot\mathbf{w})(\mathbf{k}\cdot\mathbf{w}^{*})$, where the latter term vanishes since the curl field is divergence free. Defining $P_{\theta}=\theta\theta^{*}$ and $P_{w}=\mathbf{w}\cdot\mathbf{w}^{*}$ we get

[TABLE]

The normalization of the divergence and velocity power spectra is given by the requirement that $\theta(\mathbf{k})=\delta(\mathbf{k})H_{0}f$ in the linear limit. This follows from the definition of the divergence power spectrum, since

[TABLE]

where the last equality only holds in linear theory.

All power spectra are normalized to the matter power spectrum from CAMB, which is shown in figure 4 along with the matter power spectrum from the N-body simulation. Since for the velocity power spectrum, non-linear effects are non-negligible even at the largest scales in the N-body simulation at $z=0$, the normalization of $\theta$ is instead determined from a snapshot at $z=49$. In figure 4, the divergence, curl, and total velocity power spectra from the N-body simulation are shown, as well as the linear velocity power spectrum obtained from CAMB. The non-linear velocity power spectrum has been found by interpolating the N-body particle velocities to a $1024^{3}$ regular grid using the interpolation method of [30].

5 Results

In the previous sections, we presented two different estimators for $H_{0}$: $H_{\textrm{loc}}^{\textrm{mean}}$ and $H_{\textrm{loc}}^{\textrm{lsq}}$, and explained how their variances could be calculated from either mock observations or from the velocity power spectrum. In table 1, we show the results corresponding to each of the window functions shown in figure 3. The calculation is both carried out for the velocity power spectrum predicted from linear perturbation theory, and for the velocity power spectrum measured in the N-body simulation. The mock observations are also carried out in two different ways, by selecting the observed sources among the halos in the simulation and among the grid points from which the velocity power spectrum is calculated. This results in a slightly smaller variance, as expected since peculiar velocities will tend to be smaller in voids than in overdense regions.

Even when the sources are chosen among the grid points, the spread in the Hubble constant found from the mock observations is 10-20% higher than that deduced from the velocity power spectrum. The finite thickness of the bins and the randomness in the selection of sources described in section 2 implies that there will be some variations between the distribution of sources used by the individual observers. Therefore, we expect the variance obtained from the mock observations to be slightly larger than that found from the calculations based on the power spectrum. By increasing the number of observed sources in the mock observations we have confirmed that the variance does indeed decrease by several percent when the velocity field is sampled more densely.

In figure 5 we show how the variance depends on the highest value of $k$, $k_{\textrm{max}}$, included in the analysis. As expected from the window functions, the value of $k_{\textrm{max}}$ has the largest impact on the variance associated with the CfA3+OLD supernova distribution, which peaks at the largest $k$-value, i.e. the smallest scales. It has nearly no effect on the top-hat window, for which only the largest scales in the simulation contribute significantly to the variance.

6 Discussion and conclusions

We have calculated the variance in the local value of $H_{0}$ for two different estimators and three different distributions of sources. We have considered the differences caused by the non-linear evolution of the velocity power spectrum, as well as the significance of how the velocity field is sampled in the N-body simulation. In the calculations based on perturbation theory, it was found that replacing the linear prediction for the velocity power spectrum by the fully non-linear velocity power spectrum measured in the N-body simulation only decreases the variance by 0.1% or less. Similarly, changing the way the velocity field is sampled in the mock observations changes the variance by less than 0.1%.

In contrast, the estimator one uses to obtain the Hubble constant has a strong influence on the variance obtained, with the estimator typically used in studies based on the power spectrum, $H_{\textrm{loc}}^{\textrm{mean}}$, resulting in a variance about twice as large as the variance obtained using the estimator typically used in studies based on N-body simulations, $H_{\textrm{loc}}^{\textrm{lsq}}$. We have demonstrated that there is reasonable agreement between the variances obtained from perturbation theory and mock observations in N-body simulations if the same estimator is used in both cases.

Neither of the two estimators considered in this study exactly correspond to how the Hubble constant is determined in actual observations. In real observations, greater weights are given to distant sources due to the fitting procedure, but this is most likely counteracted by the larger uncertainties in the redshifts and velocities of these sources compared to nearby ones. Moreover, actual surveys do not sample the whole sky, and hence have anisotropic window functions, which causes an additional variance, especially at large redshifts [19, 31]. An interesting and relevant extension of this study could take this approach even further, and calculate the exact window function corresponding to the estimator of $H_{0}$ used in modern determinations of the local Hubble constant, such as the one used in [2]. As we have shown, this could be based either on mock observations or on the velocity power spectrum. Further, by basing the measurement of the velocity power spectrum on the halos in the simulation, each weighted by their type Ia supernova rate (which could e.g. be calculated from the star formation histories as in [32]), the velocity power spectrum could be probed in the same way as it is in supernova surveys. This would allow for a thorough understanding of what exactly is measured in surveys of the local value of the Hubble constant, and how much this quantity can be expected to vary within any given cosmological model.

Acknowledgments

We thank Dominik Schwarz for discussions that inspired this work and for useful comments.

Parts of the numerical work for this project have been done using computer resources from the Center for Scientific Computing Aarhus.

This work was supported by the Villum Foundation.