The impact of anisotropic sky-sampling on the Hubble constant in numerical relativity

The simulations we use were performed using the Einstein Toolkit¹¹1https://einsteintoolkit.org (ET; Löffler et al., 2012; Zilhão & Löffler, 2013); an open source numerical relativity (NR) code which has been adapted and used for cosmological simulations (e.g. Bentivegna, 2016; Bentivegna & Bruni, 2016; Macpherson et al., 2017, 2018, 2019). We used FLRWSolver (Macpherson et al., 2017) to generate the initial data as an Einstein-de Sitter (EdS) space-time with linear perturbations following the matter power spectrum output from CLASS²²2http://class-code.net (Lesgourgues, 2011) at redshift $z=1000$ (with default Planck parameters). The initial data is evolved using the McLachlan thorn ML_BSSN which implements the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) method for NR, alongside the GRHydro thorn for hydrodynamics evolution. GRHydro does not allow for a pure dust (i.e., pressure $P=0$) evolution, so instead we choose a negligible pressure with respect to the mass density, i.e. we set $P\ll\rho$ (which was shown in Macpherson et al., 2017, to be sufficient for matching a dust FLRW evolution).

We evolve the simulations from the initial redshift $z_{\rm ini}=1000$ to $z\approx 0$, with the final slice defined by the change in average volume of the entire domain. We use the effective scale factor $a_{\mathcal{D}}\equiv V_{\mathcal{D}}^{1/3}$ where the volume is $V_{\mathcal{D}}=\int_{\mathcal{D}}\sqrt{\gamma}\,d^{3}X$, $\gamma$ is the determinant of the spatial metric of the hypersurfaces, and $\mathcal{D}$ (in this case) is the entire simulation domain. The final slice is chosen to be the one with $a_{\mathcal{D}}\approx 1$ where the initial value is $a_{\mathcal{D},{\rm ini}}=1/(1+z_{\rm ini})$. This corresponds to a change in effective scale factor of $10^{3}$.

We choose a harmonic-type gauge for the evolution of the lapse function (and zero shift vector) such that while our simulations are close to the initial EdS background our simulation coordinate time coincides with the conformal time, $\eta$. We wish to emphasise that there is no enforced background cosmology and the evolution of the initial data is completely general using NR. Also of importance to mention here is the fact that the simulations contain no dark energy due to the nature of the ET thorns used for the evolution. Including $\Lambda$ is an important development of using the ET for cosmology, and is the focus of ongoing work. For further details on our numerical method we refer the reader to Macpherson et al. (2017) and Macpherson et al. (2019).

We use a simulation with box side length $L=1024\,h^{-1}$ Mpc and numerical resolution $N=256$, where the total number of grid cells is $N^{3}$. To ensure that small-scale structure is sampled with sufficiently many grid cells, we apply a sharp cut-off to the initial power spectrum below scales with wavelength $\sim 10$ grid cells for both simulations. While structure can (and does) grow beneath this scale when the structure reaches the nonlinear regime, removing this from the initial data greatly reduces our numerical error at late times (see also Giblin et al., 2016). Importantly, to ensure numerical convergence of our results we perform a second, lower-resolution simulation which samples the same physical scales. In Appendix A we use these two simulations to show that our main results are robust to changes in resolution.

Throughout this work, we quote length units in $h^{-1}$ Mpc. This $h$ enters our simulations only via the setting of initial data; arising in the units of the CLASS power spectrum. We choose $h=0.45$ in generating the power spectrum; corresponding to the EdS model which is the background for our initial data. However, we wish to clarify that this value of $h$ is not necessarily a good fit to the Hubble expansion at late times in the simulation since this is not enforced via the use of a background FLRW expansion. An important consequence of this is that constraints on $H_{0}$ from our simulations can thus vary from the typical value of $H_{0}=100\,h$ km/s/Mpc due to local inhomogeneity and/or anisotropy.

The Hubble expansion that we find when smoothing over the entire $z\approx 0$ spatial slice is $H_{0}=100.001\,h$ km/s/Mpc (with $h=0.45$, see Macpherson et al., 2018). This tells us that, on average, the simulation matches the predicted Hubble expansion of the EdS model used as a background for the initial data to $0.001\%$. In this work, we will constrain $H_{0}/h$ for each observer, such that a constraint which is consistent with the “true” simulation $H_{0}$ would be 100.001 km/s/Mpc.

We perform our ray tracing analysis as post-processing of the simulation data using the mescaline code (Macpherson, 2023). We randomly place observers on the $z\approx 0$ spatial slice of the simulation and propagate light rays back in time through subsequent spatial slices by solving the geodesic equation using the simulated metric tensor. This allows us to calculate the photon energy, and thus redshift, along the path of the light ray. Additionally, we solve the Jacobi matrix equation to obtain the angular diameter distance along each light ray in the simulation. Details of the specific equations solved and tests of the software we use can be found in Macpherson (2023). Mescaline has been shown to produce distances and redshifts to better than $0.01\%$ accuracy.

We propagate rays which cover a range directions on the observer’s skies. Thus, the maximum redshift we can reach without introducing spurious correlations coincides with the co-moving length of the simulation domain. We are interested in reproducing the low-redshift Pantheon SNe sample used in the determination of $H_{0}$, the maximum redshift of which is $z=0.15$, which lies within the domain of our main and lower-resolution simulations.

We randomly position 25 observers in the simulation domain. We choose lines of sight for each observer coinciding with the directions of SNe in the Pantheon catalog, which we describe in the next section.

Figure 1 shows all-sky maps for the luminosity distance from ray-tracing for one example observer. We show the relative difference between the ray-traced distance in our NR simulation, $D_{L}$, and the EdS distance, $D_{L,{\rm FLRW}}$, for the same redshift (calculated along each individual line of sight, since $z$ varies across the sky as well). The top panel shows the lower-bound of the redshift interval we choose to generate mock data, namely $\bar{z}\approx 0.023$, where $\bar{z}$ is the mean redshift across the sky. The bottom panel shows the upper bound of the redshift interval of $\bar{z}\approx 0.15$. At $\bar{z}\approx 0.023$, our ray-traced distance is different from the FLRW distance we use to fit the data by $\sim$20–40%, and at $\bar{z}\approx 0.15$ our distance is $\sim$5–6% different. The differences we see here can be attributed to inhomogeneites in the space-time of the simulation which affect the path of incoming photons and hence the luminosity distance.

Our observers and ‘sources’ are chosen to be co-moving with the simulation fluid in this work. The redshift everywhere along the geodesic is calculated from the observed energy at that point, $E\equiv k^{\mu}u_{\mu}$. Here, $k^{\mu}$ is the photon 4–momentum and $u^{\mu}$ is the observer’s 4–velocity. Choosing our observers and ‘sources’ to be co-moving with the fluid means identifying $u^{\mu}$ with the 4–velocity of the fluid flow at that position in the simulation.

In SNe constraints of the Hubble constant, it is common to apply peculiar velocity (PV) corrections to measured redshifts (Peterson et al., 2021) as well as a boost due to our own velocity as inferred from the CMB dipole (Sullivan & Scott, 2021). PV corrections account for motion at a variety of scales: peculiar motion of galaxies within their group, inter-group dynamics, as well as larger-scale coherent flow corrections.

We do not make any PV corrections to the redshifts calculated from our simulations. This is because our simulations sample structures down to a minimum scale of $8\,h^{-1}$ Mpc. Additionally, structures below $\sim 40\,h^{-1}$ Mpc have been artificially removed from the initial data for purposes of minimising numerical error (see Section 2.1). Thus, structure on scales beneath $\sim 40\,h^{-1}$ Mpc is significantly reduced with respect to expectations based on the EdS model. Due to this—as well as the fact our simulations rely on a fluid approximation instead of particles—we do not form bound structures. We thus restrict our simulations to scales above that of the largest bound objects. Most of the PV corrections to standard analyses occur beneath the resolution scale of our simulations, and are thus not necessary.

An interesting avenue would be to study whether coherent flow corrections could reduce the scatter we find in Section 3.4. However, we leave this analysis to future work and continue without any velocity corrections to our data.

We make a mock catalog of objects based on the 237 Pantheon SNe³³3https://github.com/dscolnic/Pantheon in the redshift range $z\in[0.023,0.15]$. We initialise the ray-tracer to trace geodesics in the directions of this sample of objects. To do this, we use the measured right ascension (RA) and declination (dec) of each SN to generate Cartesian $(x,y,z)$ coordinate directions. These direction vectors are subsequently used to generate initial data for the photon 4–momentum $k^{\mu}$ (see Section 6.3 of Macpherson, 2023, for full details on setting initial data).

Coordinate systems are oriented with respect to some reference direction, in the case for our observations this is typically the Earth’s poles. However, there is no a priori preferred orientation of an observers coordinate system. For this reason, we choose five random orientations of the coordinate axes for each observer to ultimately generate five different samples of 237 objects for each observer. To perform the rotations we randomly choose three Euler angles and use the Rotation class in scipy to apply the transform to the $(x,y,z)$ coordinates of each SNe.

Since the distribution of Pantheon SNe is anisotropic, as we vary the coordinate orientation we might expect our observers to infer different values of $H_{0}$. This comes from the underlying anisotropy of the space-time expansion nearby the observer due to local structures. In particular, if an observers anisotropic distribution of objects samples the underlying expansion in a particular way we might expect them to measure a higher or lower $H_{0}$ than the true “background” expansion in their environment.

We propagate rays along all lines of sight until the redshift reaches $z=0.15$. The resulting data consists of discretised geodesic paths for each line of sight, namely, we have $(z,D_{L})$ data for many points along each observed direction. We then take the measured redshifts of the Pantheon SNe and find the closest match from the simulation data along the corresponding geodesic.

We test the impact of an anisotropic sky sampling by also constraining $H_{0}$ using an isotropic sample of the same size. For the isotropic sample, we use ray-tracing data with a distribution of objects across each observers sky coinciding with HEALPix indices with $N_{\rm side}=32$. From this catalog, we randomly select 237 lines of sight and match the Pantheon redshift distribution for these objects. We thus have an equivalent number of objects and distribution in redshift as described above, the only change is the distribution of objects on our observer’s skies. We use the same observers as in the mock Pantheon samples for this analysis.

The top panel of Figure 2 shows the five random rotations of the Pantheon SNe directions we use for all observers. In the bottom panel of Figure 2 we show an example of randomly-drawn directions (with the same number of objects and redshift distribution as the top panel) which we use to compare to the constraints using the anisotropic distributions of objects in the top panel.

Figure 3 shows the redshift distribution of our synthetic catalog (red step histogram) and the Pantheon SNe (grey shaded histogram) after making our chosen redshift cuts of $z\in[0.023,0.15]$. We show the distribution of 237 objects for only one observer and one rotation of coordinates (corresponding to the blue circles in Figure 2).

Local inferences of the Hubble constant using low-redshift data Riess et al. (such as 2022); Freedman et al. (such as 2019) make use of a Taylor series expansion of the luminosity distance in redshift within the FLRW models. This is known as FLRW cosmography (Visser, 2004) and we begin by writing the luminosity distance, $d_{L}$, of an object as

where the coefficients of the expansion can be written in terms of the familiar cosmological parameters

Here, $q_{0}$ is the deceleration parameter, $j_{0}$ is the jerk parameter, and $\Omega_{k,0}$ is the spatial curvature parameter. All parameters are evaluated at redshift zero, as indicated by the subscript.

In our constraints, we translate the luminosity distance $D_{L}$ calculated from the ray tracing into the distance modulus $\mu$, defined as

We do this so that we can use the Pantheon covariance matrix in our analysis, which is defined for the distance moduli of the SNe.

We place constraints on $H_{0}$ by modelling a $\chi^{2}$ logarithm of the likelihood, ${\rm ln}\mathcal{L}=-0.5\chi^{2}$ where

$\Delta=\mu_{\rm sim}-\mu_{\rm FLRW}$ is the difference vector between simulation distance moduli, $\mu_{\rm sim}$, and FLRW distance moduli, $\mu_{\rm FLRW}$, for the same redshift calculated using (1) and (4). In (5), $C_{SN}$ is the covariance matrix of the Pantheon SNe on which our mock catalog is built. In this work, we use the statistical covariance only and neglect the systematic covariance for our mock data. The statistical covariances for the Pantheon SNe are of order $\sim 10^{-3}-10^{-2}$, whereas the error in our simulated luminosity distances is maximum $10^{-4}$ (see Appendix E.5.2 of Macpherson, 2023), so we neglect the latter in the covariance matrix.

We use the Python package emcee⁴⁴4https://emcee.readthedocs.io to find the posterior distributions. We adopt a uniform prior in $H_{0}$ with range $H_{0}/h\in[50,150]$.

The impact of inhomogeneity in the matter distribution on the observed distance-redshift relation has been studied for decades (see, e.g. Bonvin et al., 2006; Clarkson et al., 2012; Kaiser & Peacock, 2016; Fleury et al., 2017, and references therein). Inhomogeneity induces fluctuations about the best-fit FLRW model in the Hubble diagram, which has been studied in the context of perturbation theory (e.g. Camarena & Marra, 2018; Ben-Dayan et al., 2014), Newtonian simulations (e.g. Odderskov et al., 2014; Wu & Huterer, 2017), as well as in NR simulations (Macpherson et al., 2018). Usually, SNe with redshift $z<0.023$ are excluded from the analysis in attempt to minimise the impact of inhomogeneity on the resulting $H_{0}$ constraints (Riess et al., 2016). However, even at this scale the value of $H_{0}$ can depend on the observer’s local environment (e.g. Marra et al., 2013; Macpherson et al., 2018).

In this work, we are interested in making some assessment on how the inferred value of $H_{0}$ depends on the observer’s local density contrast. The minimum redshift of the observational sample might be thought of as an effective “smoothing scale” of the data; washing out any large variances in distances due to small-scale inhomogeneities. Thus, we want a measure of the observer’s local environment at this chosen scale. To this end, we smooth the density field, $\rho$, at each observer’s position within a sphere of radius $z\approx 0.023$. In the EdS model, this corresponds to a co-moving radius of $\sim 70\,h^{-1}$ Mpc. We use a Gaussian kernel with $3\sigma=70\,h^{-1}$ Mpc centered on each observer’s position to smooth the density field and estimate the local over-density on these scales, namely $\langle\delta_{o}\rangle_{70}$. We can then make a comparison of this quantity to the inferred value of $H_{0}$ to estimate the impact of local inhomogeneity on the Hubble expansion.

Since we are primarily interested in how inhomogeneities and anisotropies can impact the inference of $H_{0}$, in order to get the tightest constraints we set $q_{0}=0.5$ and $j_{0}=\Omega_{k}=0$. These values coincide with the EdS model—which we find to match our simulation averages on large scales to within 1% (see also Macpherson, 2019).

Figure 4 shows the distance-redshift relation for one example mock catalog (one observer and one rotation of their coordinate system). Cyan crosses are the mock $(\mu,z)$ data generated from the simulation, the dashed black curve is the cosmographic relation using the best-fit $H_{0}$ value to the data, and the solid grey curve is the cosmographic relation using the global $H_{0}$ value from the simulation.

Figure 5 shows posterior distributions on $H_{0}$ for an example set of six observers in the $N=256$ resolution simulation. Each panel shows constraints from six different mock surveys for a single observer. Thick black contours in each panel are the constraints from the isotropic distribution of objects (similar to the bottom panel in Figure 2) and the thin coloured contours are the constraints from the five rotations of Pantheon SNe shown in the top panel of Figure 2. These examples demonstrate the variance among observers that we see here, in particular the sometimes many-$\sigma$ change in $H_{0}$ when sampling the sky in a different way.

Figure 6 shows the best-fit values of $H_{0}$ as a function of local over-density $\langle\delta_{o}\rangle_{70}$ (see Section 3.3). Light blue squares indicate constraints from all five rotations for each of the 25 observers, while dark blue triangles are the constraints from an isotropic distribution of objects. The horizontal dotted grey line is the mean over all light blue squares, with a value of $H_{0}=100.231\pm 0.129\,h$ km/s/Mpc. The black dashed horizontal line is the globally-averaged simulation value of $H_{0}=100.001\,h$ km/s/Mpc. We see a weak correlation between the constrained value of $H_{0}$ and the local density in a region of $\sim 70\,h^{-1}$ Mpc surrounding the observer. As expected, observers in under-dense regions tend to infer a higher $H_{0}$ by about $1\,h$ km/s/Mpc and observers in over-dense regions infer a lower $H_{0}$ by less than $1\,h$ km/s/Mpc.

In particular, the dark blue triangles show a stronger correlation than the light blue squares. In these constraints, the anisotropic distribution of SNe is removed so we are mostly left with the impact of the inhomogeneities on the isotropic component of $H_{0}$. Unless we have a very isotropic distribution of objects, the local density within a sphere of $z\approx 0.023$ surrounding the observer does not necessarily indicate whether they will infer a higher or lower $H_{0}$. The anisotropy of distances across their sky often plays a more important role.