Field-level inference of H₀ from simulated type Ia supernovae in a local Universe analogue

ABSTRACT

Two particular challenges face type Ia supernovae (SNeIa) as probes of the expansion rate of the Universe. One is that they may not be fair tracers of the matter velocity field and the second is that their peculiar velocities distort the Hubble expansion. Although the latter has been estimated at |${\lesssim} 1.5~{{\ \rm per\ cent}}$| for |$z>0.023$|⁠, this is based either on constrained linear or unconstrained (random) non-linear velocity simulations. In this paper, we address both challenges by incorporating a physical model for the locations of supernovae, and develop a Bayesian Hierarchical Model that accounts for non-linear peculiar velocities in our local Universe, inferred from a Bayesian analysis of the 2M++ spectroscopic galaxy catalogue. With simulated data, the model recovers the ground truth value of the Hubble constant |$H_0$| in the presence of peculiar velocities including their correlated uncertainties arising from the Bayesian inference, opening up the potential of including lower redshift SNeIa to measure |$H_0$|⁠. Ignoring peculiar velocities, the inferred |$H_0$| increases minimally by |${\sim} 0.4 \pm 0.5$| km s|$^{-1}$| Mpc|$^{-1}$| in the range |$0.023\ {<}\ z<0.046$|⁠. We conclude it is unlikely that the |$H_0$| tension originates in unaccounted-for non-linear velocity dynamics.

1 INTRODUCTION

The emergence of data-constrained cosmological simulations equipped with galaxy formation models enables the study of supernova (SN) physics in realistic reconstructions of the local large-scale structure for the first time (e.g. Forero-Romero et al. 2011; Nuza et al. 2014; Sorce et al. 2021). Unlike random N-body simulations, constrained ones enable the reconstruction of the large-scale structure from initial conditions that reproduce the observed clustering of real galaxies in the Universe. Here, we simulate SNeIa from individual galaxy properties in such an analogue of the local Universe. We demonstrate a Bayesian hierarchical framework that enables the study of the impact of non-linear peculiar velocities on |$H_0$| and the inference of |$H_0$| from SNeIa at |$z<0.023$|⁠, which are typically discarded. The latter is enabled by the fact that realistic peculiar velocity reconstructions in the local Universe are available.

We build upon the SIBELIUS-DARK (McAlpine et al. 2022; Sawala et al. 2022) cosmological simulation, whose initial conditions were constrained with 2M++ galaxies (Lavaux & Hudson 2011) through the BORG algorithm (Jasche & Wandelt 2013; Lavaux & Jasche 2016; Jasche & Lavaux 2019), out to a distance of 200 Mpc (⁠|$z=0.046$|⁠). Whilst SIBELIUS-DARK is a dark-matter only cosmological simulation, it is equipped with GALFORM (Lacey et al. 2016), a semi-analytic galaxy formation and evolution model that emulates baryonic physics. The SIBELIUS-DARK simulation provides a catalogue of simulated galaxies, each with complete star-formation histories, which we use to derive SN rates.

SNe are clustered in the cosmic large-scale structure (Carlberg et al. 2008; Mannucci et al. 2008; Tsaprazi et al. 2022) as they typically reside in star-forming galaxies. The latter are less clustered than typical galaxies, therefore the clustering pattern of SNe may differ from such galaxies (Molero et al. 2021). At the same time, more complex intergalactic physics affects the emergence of SNe. Earlier studies assumed the distribution of SNe on the sky to be uniform (Goobar et al. 2002; Feindt et al. 2019), whereas more recent ones simulate type Ia supernovae in individual haloes (Carreres et al. 2023), or from galaxy properties (Vincenzi et al. 2021; Lokken et al. 2023). The catalogues presented here contain SNe simulated from individual galaxy properties and distributed in a realistic large-scale structure.

We first employ the mock SN catalogues we generate to study the effect of peculiar velocities induced by our local large-scale structure on the derived Hubble constant, |$H_0$|⁠, inferred from the measured redshifts and distance moduli of simulated SNeIa. The argument that local inhomogeneities affect the measured expansion dates back to Riess et al. (1998). Ever since, a large number of studies have investigated the effect providing evidence in favour or against a special configuration of the cosmic large-scale structure that could result in the emergence of a discrepancy (Riess et al. 2016, 2022) between measurements of |$H_0$| in the local and distant Universe (e.g. Wojtak et al. 2014; Carrick et al. 2015; Odderskov, Koksbang & Hannestad 2016; Kenworthy, Scolnic & Riess 2019; Böhringer, Chon & Collins 2020; Kazantzidis & Perivolaropoulos 2020; Sedgwick et al. 2021; Camarena et al. 2022; Carr et al. 2022; Castello, Högås & Mörtsell 2022; Peterson et al. 2022; Mc Conville & Ó Colgáin 2023; Poulin, Smith & Karwal 2023; Carreres et al. 2025; Giani et al. 2024; Hollinger & Hudson 2025), even in the absence of a local void (Jasche & Lavaux 2019).

Our study adds to the above investigations by (a) considering our specific large-scale structure instead of random N-body simulations which provides us with knowledge of peculiar velocities, (b) using non-linear peculiar velocities at the locations of simulated SNe along with their non-linear correlations through a constrained gravity model, (c) exploiting correlations with the full large-scale structure, instead of individual structures that gravitationally influence the local Universe dynamics, (d) looking at the locations of SNeIa generated according to galaxy evolution. Further, since SIBELIUS-DARK emulates an analogue of our Local Group, our results naturally account for effects of our specific environment on the observer and therefore, potential deviations from the Copernican Principle (Camarena et al. 2023). The demonstration presented in this study paves the way for a future analysis of observed SNIa data at |$z<0.023$|⁠.

This paper is structured as follows: In Section 2, we describe the SNIa rate model and the Bayesian Hierarchical Model for the inference of |$H_0$|⁠. In Section 3, we present our findings on the effect of peculiar velocities on |$H_0$| measurements and further validate the model on simulated SNIa data with assumed observational uncertainties at |$z<0.023$|⁠. Finally, in Section 4 we summarize our conclusions and provide an outlook on upcoming developments.

2 METHOD

Our analysis is based on the SIBELIUS-DARK catalogue of simulated galaxies (McAlpine et al. 2022). The catalogue was generated within the dark matter density and velocity field of the SIBELIUS dark-matter only simulation (Sawala et al. 2022) with the application of the GALFORM galaxy formation and evolution model (Lacey et al. 2016). The SIBELIUS-DARK catalogue contains 4527 471 hosts with stellar masses |$M>10^7 \mathrm{ M}_\odot$| at |$z=0$|⁠. This mass range includes all SNIa host galaxies (Toy et al. 2023, fig. 1). The GALFORM star formation histories are informed about galaxy mergers, cooling, feedback, photoionization, as well as the evolution of cosmic metallicity and mass–metallicity correlations. The gravity model through which the constrained initial conditions of the SIBELIUS simulation were inferred, further gives us access to non-linear peculiar velocities at the field level (3D reconstruction). Since SIBELIUS (Sawala et al. 2022) is a simulation constrained with the initial conditions inferred from the 2M++ catalogue (Jasche & Lavaux 2019), our analysis carries over all the conclusions reported therein. Specifically, the 2M++ large-scale structure reconstruction we use shows no evidence of a large-scale underdensity which could account for the Hubble tension (Jasche & Lavaux 2019). Therefore our conclusions below are not driven by a local void which distorts the observed cosmic expansion. Below, we describe our modelling of SN occurrence in each simulated galaxy.

2.1 SNIa rate modelling

SNIa rates, |$R_\mathrm{IA}$|⁠, are typically modelled as (e.g. Sullivan et al. 2006; Graur, Bianco & Modjaz 2015; Wiseman et al. 2021)

where |$\tau$| is the time following a star-formation burst, |$t_0$| corresponds to the lookback time to a galaxy’s redshift, |$t_\mathrm{f}$| indicates the lookback time to the epoch at which the first star formation burst occurred, |$\Psi$| is the star formation history and |$\Phi$| the delay time distribution. The latter indicates the time between the formation of a stellar population and the explosion of the first SNeIa (Freundlich & Maoz 2021). It can be modelled as a power law (Wiseman et al. 2021)

where |$t_\mathrm{p}$| is the time required for a population of white dwarves to form after a star formation burst. We take |$t_\mathrm{f}$| to be the lookback time to the redshift at which |$\Psi$| becomes non-zero for the first time in each galaxy. |$\Psi$| is provided by SIBELIUS-DARK for every galaxy in the simulation volume out to |$z=20$|⁠. We adopt |$A = 2.11 \times 10^{-13}$| SNe |$\mathrm{ M}_\odot ^{-1}$| yr|$^{-1}$|⁠, |$\beta = -1.13$|⁠, |$t_\mathrm{p} = 40$| Myr (Wiseman et al. 2021). We rescale all rates by |${\sim} 3$| to match the mean global volumetric rate reported by Perley et al. (2020): (2.35 |$\pm$| 0.24) |$\times 10^4$| SNeIa yr|$^{-1}$| Gpc|$^{-3}$|⁠. The need for rescaling likely arises due to a combination of factors. GALFORM was calibrated using data sets other than SNIa ones, therefore a deviation from observed rates is in principle expected. Further, a mild underestimation of mass was reported in SIBELIUS-DARK in the range |${\sim} 5\times 10^9{\!-\!}10^{11}\,\mathrm{ M}_\odot$|⁠, where most SNIa hosts reside, which explains the need to apply a rescaling factor |${>}1$| to match the observed SNIa production efficiency.

The underestimation of the global volumetric rate affects almost the entire population systematically, so the relative object-by-object weights we recover, which are of interest in our analysis, are credible. Such rescalings have been performed in hydrodynamical simulations when comparing the derived SNIa rates to observations (fig. 3, Schaye et al. 2015). Our analysis is by construction insensitive to global rate rescalings, since it is only affected by the SNIa host velocity and distance distribution. As we will demonstrate in Section 3, it is also on average insensitive to the per-host rate modelling given the local Universe realization. It is the subject of future work to test various SNIa rate and galaxy formation models against observations, since it has been shown that the assumptions on the star-formation history may significantly affect astrophysical predictions (Briel et al. 2022).

Finally, we model the occurrence of SNe in a galaxy in a given period of time |$\Delta t$|⁠, as a Poisson process (Wiseman et al. 2021)

where |$N_\mathrm{s}$| is the number of SN explosions per galaxy, and |$\lambda _\mathrm{SN} = R_\mathrm{IA} \Delta t$| the intensity of the Poisson process. In order to obtain the SN mock catalogues, we sample from the above distribution. A simulated Poisson SN sample around the Coma cluster along with 2M++ and simulated SIBELIUS galaxies is shown in Fig. 1.

2.2 Hubble constant inference

We will use realizations drawn from equation (3) to validate a framework for the incorporation of peculiar velocities in inference of |$H_0$| and quantify the effect of our local velocity environment on |$H_0$|⁠. The latter impact is given by the fractional Hubble uncertainty

where |$H_{0,\mathrm{fid}}$| represents a fiducial value and |$H_0$| the Hubble constant inferred if peculiar velocities are properly accounted for. The fractional uncertainty can also be written as the intercept of the Hubble diagram, |$\alpha _\mathrm{B}$| (Riess et al. 2022; Scolnic et al. 2023; Carreres et al. 2025)

where |$M_\mathrm{B}$| is the absolute magnitude of SNeIa in the B band. We write equation (5) once for the fiducial and once for the inferred Hubble constant and subtract the former from the latter, assuming that |$M_\mathrm{B}$| is fixed. Recasting the resulting terms into the fractional Hubble variation with the help of equation (4), we arrive at a proxy for the fractional variation in the Hubble constant

which we will use with |$\Delta H/H_0$| interchangeably according to the above equation, since it represents the difference of the inferred from the fiducial Hubble constant. While the simulation extends in |$0\,{<}\,z<0.046$|⁠, we will truncate it to |$z>0.023$| in which the magnitude–redshift relation is computed for the purposes of investigating the impact of peculiar velocities on |$H_0$| (Efstathiou 2021). We will also consider SNeIa at |$z<0.023$| to demonstrate that our framework allows us to account for peculiar velocities for SNeIa which are typically discarded from |$H_0$| analyses. This is the only redshift cut applied for part of our analysis, for comparison with typical cuts made in the literature to avoid excessive peculiar velocity corrections. The SIBELIUS-DARK catalogue is otherwise redshift complete. Note that the Bayesian hierarchical model that we present would allow lower redshift SNeIa to be used in an analysis of |$H_0$|⁠.

2.3 Non-linear velocity covariance matrix estimation

Here, we use our knowledge of the non-linear peculiar velocity field from the 2M++ Bayesian reconstruction (Jasche & Lavaux 2019) to reduce the impact of peculiar velocities. The Bayesian reconstruction provided a posterior of the initial conditions for structure formation which were evolved to present-day density and velocity fields that reproduce the observed clustering of the 2M++ spectroscopic galaxy catalogue. The velocity fields were estimated using a Cloud-in-Cell kernel (e.g. Hockney & Eastwood 2021; Mukherjee et al. 2021). The minimum scale at which we trust the physical modelling of peculiar velocities in the 2M++ reconstruction is 2.65 Mpc h⁻¹. For a real data application, the velocity dispersion on smaller scales must be modelled externally.

Carreres et al. (2025) argued that linear peculiar velocity correlations can affect the derived |$\alpha _\mathrm{B}$|⁠. We take into account both the non-linear velocity correlations and the varying quality of the 2M++ reconstruction as a function of location in the reconstruction box by accounting for the full non-linear peculiar velocity covariance matrix. We consider |$N_\mathrm{sim} \sim 240$| non-linear velocity fields of the 2M++ reconstruction at the locations of simulated SNeIa for each realization of SN positions (Jasche & Lavaux 2019). This provides us with |$N_\mathrm{sim} \times N_\mathrm{SN}$| velocity values from which we build an |$N_\mathrm{SN} \times N_\mathrm{SN}$| covariance matrix per SNIa data set, which we denote by |$C_\mathrm{2M++}$| in what follows. The covariance matrix is then estimated as

where |$v^i$| is the velocity field at the locations of SNe in realization i and |$\bar{v}$| is the mean velocity field across all |$N_\mathrm{sim}$| realizations at the SN locations. The number of velocity fields, |$N_\mathrm{sim}$|⁠, is determined by the data availability and is statistically adequate for the purposes of our analysis, because within that number, all uncertainties in the reconstruction have been captured by the velocity covariance matrices. Our covariance matrices are robust to this number. The velocity fields are created by forward modelling from initial conditions that are sampled from the BORG posterior. The variability in the samples arises from various sources, such as Poisson shot noise of galaxies and prior ranges for parameters. At the same time, each velocity field, having been constructed with a gravity model, accounts for physical correlations between any points. We choose to simulate |$N_\mathrm{SN}$| SNe, such that |$N_\mathrm{SN} < N_\mathrm{sim}-2$|⁠, otherwise the estimated covariance matrix is not invertible. We generate multiple SNIa Poisson samples which satisfy the above condition. We keep only those that contain at most 1 SN per voxel, such that information is not repeated in the covariance matrix. These two conditions guarantee that we end up with invertible covariance matrices for each SN realization. In our setting, these conditions are satisfied for |$N_\mathrm{SN} \sim 200$| SNe (within Poisson noise) per data set.

The impact of peculiar velocities on |$H_0$| in observations is typically inferred from distance modulus |$\hat{\mu }$| and redshift estimates |$\hat{z}$|⁠, which are known with some noise |$\sigma _\mathrm{\mu }$| and |$\sigma _\mathrm{z}$|⁠, respectively. The error on distance moduli can be significant and therefore cannot simply be added in quadrature to redshift uncertainties as is typically done within the errors assumed here. More importantly, at |$z<0.023$| where the contribution of peculiar velocities to the observed redshift is significant, the dependence of distance moduli on the peculiar velocities must be self-consistently modelled along each line of sight. In this section, we advance the formalism in Roberts et al. (2017) to infer the intercept of the Hubble diagram from redshift and distance moduli data, by considering their uncertainties along both their corresponding axes, instead of adding them in quadrature. This requires a Bayesian Hierarchical Model that accounts for the known non-linear 2M++ peculiar velocities in |$H_0$| inference.

2.4 Hubble constant inference

We will consider a set of observed redshifts, |$\hat{z}$|⁠, observed distance moduli, |$\hat{\mu }$|⁠, and the three-dimensional reconstruction of velocities from which we obtain samples of peculiar redshifts at any cosmological redshift, |$z_\mathrm{c}$|⁠. These data are given in the presence of uncertainties which render the exact true value of peculiar redshifts, |${z}_\mathrm{p}$|⁠, and cosmological redshifts unknown. We characterize the posterior peculiar redshifts with means |$\bar{z}_\mathrm{p}$| and covariance |$C_{2\mathrm{ M}++}$| from the 2M++ reconstruction. As all quantities except for |$\alpha _\mathrm{B}$| are vectors, we will omit writing them in bold. Marginalizing over |$z_\mathrm{p}$| and |$z_\mathrm{c}$|⁠, the posterior may be written as (note that this is a multidimensional integral over all the SN redshifts)

where |$\pi (\alpha _\mathrm{B})$| is the prior on |$\alpha _\mathrm{B}$| (which we will take to be uniform), and we have assumed that the estimated redshifts and distance moduli conditioned on the true values are independent. We have removed variables where there is no dependence with |$\cancel {\phantom{.}}$|⁠. The distance modulus is given by |$\bar{\mu }=\mu _\mathrm{c}+10\log (1+z_\mathrm{p})-5\alpha _\mathrm{B}$|⁠, |$\mu _\mathrm{c}$| being the fiducial distance modulus at any given cosmological redshift |$z_\mathrm{c}$|⁠. The Doppler boosting term is a small correction, and its dependence on the true |$z_\mathrm{p}$| complicates the integral enormously and prevents a necessary simplification, so we approximate this dependence by replacing |$z_\mathrm{p}$| by |$\bar{z}_\mathrm{p}$|⁠, i.e. we take

This assumption has very little impact on our analysis, as our final posterior on |$H_0$| is insensitive to it. We simplify the posterior in equation (10)

We also neglect a small correction due to selection on the estimated redshifts |$\hat{z}$|⁠. Retaining only the relevant dependencies, the integral simplifies to

Reordering the integrals, we arrive at

In what follows, we will formulate the integral over |$z_\mathrm{p}$| analytically to ensure tractability in high-dimensional spaces. Following the derivation in Sellentin & Heavens (2016), |$p(z_\mathrm{p}|\bar{z}_\mathrm{p})$| the likelihood marginalized over the uncertain true covariance matrix is a multivariate Student-t likelihood (see also Percival et al. 2022). However, this does not allow analytic marginalization over |$z_\mathrm{p}$|⁠, so we instead use the Hartlap-corrected Gaussian approximation with covariance C, which provides a debiased estimate of a covariance matrix that has been estimated with uncertainty, as the one in 2M++. The new covariance reads (Hartlap, Simon & Schneider 2007), |$C = C_\mathrm{2M++}{N_\mathrm{sim}}/({N_\mathrm{sim}-N_\mathrm{SN}-1})$|⁠. The mean is |$\bar{z}_\mathrm{p}$|⁠. For the redshift likelihood |$p(\hat{z}|z_\mathrm{p},z_\mathrm{c})$|⁠, we will assume independent Gaussian errors with covariance |$\Sigma =$| diag(⁠|$\sigma _\mathrm{z}^2$|⁠), and mean |$\bar{z}=(1+z_\mathrm{c})(1+z_\mathrm{p})-1$| (Davis et al. 2011). For shortness of notation, we define

which becomes

where

By rearranging the terms in the first term in the integrand above, we arrive at

where

Substituting |$u=z_\mathrm{p}-\bar{z}_\mathrm{p}$| and |$v=\tilde{z}-\bar{z}_\mathrm{p}$| above and evaluating the n-dimensional Gaussian integral gives

Using the Woodbury (1950) formula and the matrix determinant theorem we find

Then, our final posterior reads

where |$\bar{q} = \hat{z}-[(1+z_\mathrm{c})[1+\bar{z}_\mathrm{p}(z_\mathrm{c})]-1]$|⁠, |$\tilde{C} = C\times {\rm diag}(1+z_\mathrm{c})^2$|⁠, |$\Sigma _\mathrm{\mu }=$| diag(⁠|$\sigma _\mathrm{\mu }^2$|⁠), and we have assumed a |$z_\mathrm{c}^2$| prior on cosmological redshifts, motivated by the volume effect in a low-redshift survey. We find that deviations from |$z_\mathrm{c}^2$| are minimal and have insignificant impact on our conclusions in our redshift range. This is an integral of the same dimensionality as the |$\hat{\mu }$| vector. The flowchart of the method is shown in Fig. 2.

In the above formulation, we have further assumed that the change in |$H_0$| leaves the peculiar velocity reconstruction invariant. It is unlikely that changes in the reconstructed velocity field with |$H_0$| affect our inference. In the regions of highest sensitivity of the density field to |$H_0$| (fig. 2, Kostić et al. 2022), changes of |$5$| km s|$^{-1}$| Mpc|$^{-1}$|⁠, would result in a change in overdensity of no more than |$\Delta \delta \approx 0.25$| at the resolution of the 2M++ inference, in filaments where the density is typically |$\delta \sim 6$| at this resolution (fig. 9, Jasche & Lavaux 2019), so the effect is likely to be small. Moreover the peculiar velocity field has a longer correlation length than the dark matter density field, and outside the filaments, there is very little sensitivity to the |$H_0$| assumed in the reconstruction (Kostić et al. 2022).

The dimensionality of the integral in equation (23) is the length of the SNIa data vector, and such high-dimensional integrals can be extremely challenging to compute, even with nested sampling (Skilling 2004; Skilling 2006). Attempting to compute this integral with the dynesty nested sampler (Feroz, Hobson & Bridges 2009; Handley, Hobson & Lasenby 2015a, b; Speagle 2020; Koposov et al. 2024), led to reasonable posteriors on |$H_0$|⁠, which could not, none the less, become robust to the sampler’s tunable parameters. Further, the reported uncertainties were underestimated. Therefore, due to the challenging nature of computing this high-dimensional noisy integral, we will proceed with an analytical approach.

We will expand |$\mu$| in a linear Taylor expansion, which will yield a Gaussian posterior (Laplace approximation). In what follows, we will derive an analytically tractable form for the |$H_0$| posterior, where we will approximate the integrand as a multivariate Gaussian, beginning with

where

and

We will expand around the |$\tilde{z}$| that minimizes |$\bar{q}^2$|⁠, such that

where |$\mu ^{\prime }_i = \frac{\partial \mu _i}{\partial z}$| at |$z=\tilde{z}_i$|⁠. We minimize |$\bar{q}^2$| over a redshift grid at the lowest resolution that guarantees no new voxels are added to the line of sight to each SNIa when increasing the resolution further, in order to avoid artefacts resulting from the discrete nature of the velocity field. Therefore, we arrive at

where

The above expansion is likely to be inadequate for photometric redshift errors, but for typical spectroscopic uncertainties a linear Taylor expansion about |$\tilde{z}$| is sufficient. Setting |$\Lambda = \mathrm{diag}(1/\mu ^{\prime })$|⁠, we arrive at

This is a Gaussian integral in |$z_\mathrm{c}$|⁠, if we approximate the |$z_\mathrm{c}^2$| prior by |$\tilde{z}^2$|⁠, and extend the lower integration limits to |$-\infty$|⁠. Thus, we arrive at

where

Applying the Woodbury formula, we arrive at the simple expression

where

Finally, the integral can be simplified to

where |$\hat{\beta }$| is the maximum a posteriori (MAP) value of |$5\alpha _\mathrm{ B}$| as we assume a uniform prior

and

and |$\mathbb {I}$| is a vector of 1s. The above may be written in an easily interpretable form

i.e. the MAP value is a weighted sum of estimates |$\hat{\beta }_i$| from each SNIa (with weights, |$w_i$|⁠, which reduce to minimum variance weights if the velocity field is uncorrelated).

2.5 Mock data generation

In this work, we will validate the above framework on self-consistent mock data, use it to explore the impact of peculiar velocities on |$H_0$| with a statistically principled treatment of observational uncertainties and demonstrate that SNeIa at |$z<0.023$| need not be discarded in an |$H_0$| inference, since their peculiar velocities can be accounted for. For this purpose, we first create simulated redshifts and distance moduli. The mock data generation is as follows

• Draw a Poisson realization of SNe from equation (1).

• Grid their hosts in three dimensions using Nearest Grid Point assignment and the known cosmological redshifts, |$z_\mathrm{c}^\mathrm{t}$|⁠, right ascension and declination from the SIBELIUS-DARK catalogue.

• The 2M++ non-linear peculiar velocity at those locations will be |$v_\mathrm{r}$| and the corresponding peculiar redshift will be |$\bar{z}_\mathrm{p}^\mathrm{t}=v_\mathrm{r}/c$|⁠.

• The observed peculiar redshifts are assumed to have been drawn from a multivariate Gaussian with mean |$\bar{z}_\mathrm{p}^\mathrm{t}$| and covariance C. This represents our assumption that the true peculiar velocities are known within the uncertainty with which the velocity reconstruction could constrain them. The simulation velocities are |$\bar{z}_\mathrm{p}^\mathrm{t}$|⁠.

• The observed distance moduli are drawn from a multivariate Gaussian with mean |$\mu _\mathrm{c}(H_{0,\mathrm{true}},z_\mathrm{c}^\mathrm{t})+10\log {(1+\bar{z}_\mathrm{p}^\mathrm{t})}$| and covariance |$\Sigma _\mathrm{\mu }$| = diag(⁠|$\sigma _\mathrm{\mu }^2$|⁠).

• The observed redshifts are drawn from a multivariate Gaussian with mean |$(1+z_\mathrm{c}^\mathrm{t})(1+\bar{z}_\mathrm{p}^\mathrm{t})-1$| and covariance |$\Sigma$|  = diag(⁠|$\sigma _\mathrm{z}^2$|⁠), where |$\bar{z}_\mathrm{p}^\mathrm{t}$| is the 2M++ peculiar redshift at |$z_\mathrm{c}^\mathrm{t}$|⁠.

We assume |$\sigma _\mathrm{\mu }=0.13$| mag (Gris et al. 2023, and references therein). For the redshift covariance, we assume a diagonal matrix with uncertainties |$\sigma _\mathrm{z} = 0.001$| for all SNeIa. While this uncertainty can be reduced by considering host redshifts (Aubert et al. 2025), here we assume this larger scatter to demonstrate the robustness of the model to it. For the generation of mock data, we use the one 2M++ peculiar velocity realization that has the same initial conditions as the ones with which SIBELIUS was constrained. For our analysis, we use the mean 2M++ peculiar velocity across all realizations. This is because only one among the 2M++ realizations will correspond to the velocity field in the Universe, but we can only infer the velocity field within observational uncertainties. Further, for the analysis of the data, we assume the 2M++ covariance matrix estimated at |$\tilde{z}$|⁠, whereas for the simulated SNeIa velocities we have assumed the velocity covariance at |$z_\mathrm{c}^\mathrm{t}$|⁠.

3 RESULTS

3.1 SNIa rate modelling

In Fig. 3, we show our results on the derived SNIa rates. In (a), we show the SNIa host star-formation rates at |$z=0$|⁠, in comparison to those of the entire galaxy sample. The latter presents a strong bimodality corresponding to passive and star-forming galaxies. The bimodality is seen in observations of galaxy colours (Cano-Díaz et al. 2019). As shown in low-redshift SNIa observations in Childress et al. (2013, fig. 13), the distribution peaks in the range |$-1< \log _{10}$|(SFR/[|$\mathrm{ M}_\odot$| yr|$^{-1}])<1$|⁠, in accordance with our results. The error bars represent the uncertainty from 300 Poisson realizations. In (b), we plot the SNIa rate in star-formation rate bins. We compare with the corresponding plot from the low-redshift SDSS-II SN Survey (Smith et al. 2012, fig. 4) (‘S + 12’). Our rates are consistent with the linear fit. In (c), we plot the SNIa hosts on the SFR–|$M_*$| plane. We consider a host to be passive if |$\log _{10}$|(sSFR/yr|$^{-1}){<}-12$|⁠, weakly star-forming if |$-12<\log _{10}$|(sSFR/yr|$^{-1}){<}-9.5$|⁠, and strongly star-forming if |$\log _{10}$|(sSFR/yr|$^{-1})>-9.5$|⁠. The SNIa hosts follow the same trends as Smith et al. (2012, fig. 2). In Figs 3c and d, we add the ‘split 1’ line to visually envelope the distribution of strongly star-forming hosts in the above study. This is the dashed line in Smith et al. (2012, fig. 2). ‘split 2’ represents the dashed-dotted line in the above figure. Note that Smith et al. (2012) assigned random SFRs to passive hosts, whereas we maintain the values by GALFORM. Finally, in (d) we show the SNIa rates as a function of stellar mass, compared to Smith et al. (2012, fig. 3). The overall trends agree within error bars, with the exception of the highest mass bin of strongly star-forming galaxies in (c), which contains very few members for a statistically meaningful report on the moments of the SFR distribution. We therefore consider the rates to be sufficiently realistic for the purposes of our analysis here, as we will show in what follows.

3.2 Bayesian hierarchical framework

For the first part of our analysis, we investigate the effect of the specific configuration of the non-linear peculiar velocities in our local Universe on |$H_0$|⁠. For this purpose, we generate mock data with non-linear peculiar velocities added and wrongly ignore them in our analysis. We run the analysis across |${\sim} 100$| mock data sets of |${\sim} 200$| SNeIa each at |$0.023<z<0.046$| and |${\sim} 100$| noise realizations for each data set, beyond which our results do not change significantly. This is because we want to probe the effect of the local Universe on the average SNIa sample. We present our results in Fig. 4. We find that, on average, the configuration of the local large-scale structure gives rise to |$\langle \Delta H/H_0\rangle = 0.006 \pm 0.008$| (⁠|$0.4 \pm 0.5$| km s|$^{-1}$| Mpc|$^{-1}$|⁠) in the presence of the observational uncertainties assumed here.

In the second part of our analysis, we probe what drives this mild positive change on |$H_0$|⁠. Equation (39) (along with its associated weight) provides a proxy for the contribution of each source to the |$H_0$| posterior for each data set. We therefore derive |$\beta _i$| and |$w_i$| for each source in the mock data sets. We isolate all sources at |$\pm 10^{\circ }$| in right ascension and declination from known structures in the local Universe and report the average effect |$\langle \Delta H/H_0 \rangle _\mathrm{patch} = \sum _i w_i \beta _i/(5\sum _j w_j)$|⁠, where i runs over the sources in each structure and j runs over all sources in the sample. We choose the querying radius of |$\pm 10^{\circ }$| to minimize overlap of the patches. This configuration yields |$(4.2, 1.43, 1.75, 1.45, 3.34, 1.94, 2.65, 0.93, 1.42, 0.71)$| sources per patch in the Shapley, Coma, Leo, Pavo-Indus, Hercules, Perseus-Pisces, Hydra-Centaurus, Ophiuchus, Boötes, and Horologium-Reticum superclusters, respectively. The number of SNe in each data set is limited by the number of available velocity fields, preventing tighter error bars. Upcoming reconstructions will provide more velocity fields with more advanced velocity modelling, allowing a more precise and accurate determination of the effect of local superstructures individually. We find that the largest positive contribution to |$H_0$| at |$0.023{<}z<0.046$| is from SNeIa in the direction of the Hydra-Centaurus and Shapley superclusters, although the evidence is weak. SNeIa which on average exhibit |$\langle \Delta H/H_0 \rangle _\mathrm{patch}>0$| have predominantly positive peculiar velocities, whereas those associated with |$\langle \Delta H/H_0 \rangle _\mathrm{patch}<0$| have mainly negative peculiar velocities. Given the diminishing impact of peculiar velocities on |$H_0$|⁠, the above results are consistent with the value reported by Odderskov et al. (2016) who placed the impact of non-linear velocities in random N-body simulations at |${\sim} 0.3~{{\ \rm per\ cent}}$| on |$H_0$| in the wider redshift range |$0.01<z<0.1$|⁠, ignoring velocity correlations.

Odderskov et al. (2016) explored the impact of peculiar velocities on |$H_0$| in random simulations assuming the SNIa redshift distribution of an observed sample and one resulting from rates proportional to the halo mass, finding that the rate choice affected the derived |$H_0$|⁠. To test this statement in a constrained simulation and assess the sensitivity of the first part of our analysis to the SNIa rate modelling, we repeat our analysis assuming uniform rates across the simulated galaxies instead of the rates in equation (1). Assuming uniform rates, we find |$\langle \Delta H/H_0 \rangle = 0.006 \pm 0.008$|⁠, i.e. no change on average with respect to the results where we derived SNIa rates from the galaxy star-formation histories. We therefore conclude that the per-host SNIa rate modelling is unlikely to significantly change our conclusions for a typical SNIa sample, as it does not significantly impact the host velocity and distance distributions at the level of the uncertainties assumed here. The latter is also corroborated by the fact that the 1-point velocity statistics show no evidence of relative SNIa host velocity bias with respect to galaxies. Our findings suggest that it is the location of SNIa hosts in the large-scale structure which predominantly drives the derived |$\langle \Delta H/H_0 \rangle >0$|⁠, rather than the preferential occurrence of SNeIa in star-forming hosts which occupy specific velocity environments. Therefore, once a constrained simulation of a cosmological volume has been obtained, accurate SNIa rate modelling, which is often challenging to perform, may be of secondary importance to the study of velocity dynamics. Further work is required at the level of 2-point statistics and beyond to probe the SNIa velocity bias.

The statistical framework for |$H_0$| inference presented here can be used to extend the SNIa redshift range to |$z<0.023$|⁠, which are typically discarded, as the effects of peculiar velocities are larger. To demonstrate this, we generate data with |$H_0=H_{0,\mathrm{true}}=67.77$| km s|$^{-1}$| Mpc|$^{-1}$|⁠, wrongly assume in our analysis that |$H_0^\mathrm{fid}=63$| km s|$^{-1}$| Mpc|$^{-1}$| and infer |$H_0$|⁠. The choice of |$H_0^\mathrm{fid}$| is arbitrary and our results do not depend on it, because we infer the offset of |$H_0$| from the truth. In Fig. 5, we show the |$H_0$| posterior for 9 random data sets with |${\sim} 165$| SNeIa each, according to the current number of sources in the Zwicky Transient Facility (ZTF) Bright Transient Survey (BTS) Sample Explorer (Fremling et al. 2020; Perley et al. 2020). This was the number of SNe in the BTS sample explorer with the appropriate cuts at |$z<0.023$|⁠. A similar number of SNeIa suitable for our analysis is available from the Pantheon + SH0ES data set (Riess et al. 2022; Scolnic et al. 2022). We show only 9 mocks for the purposes of presentation, but 100 were checked and found consistent with the ground truth. The posteriors are consistent with |$H_{0,\mathrm{true}}$|⁠. The posteriors are also consistent with the ground truth when we set |$H_0^\mathrm{fid}>H_{0,\mathrm{true}}$|⁠. The consistency persists even when we generate data with the SIBELIUS peculiar velocities, which have been generated by a full gravity solver which potentially produces more non-linear velocities than the BORG quasi-linear velocity model. The SIBELIUS velocities are at the scale of individual hosts, so they, in principle, account for velocity dispersion to some degree. However, we present our results using data generated with the 2M++ velocities, as only the velocity covariance matrix and mean velocity field from the 2M++ inference are available.

4 DISCUSSION & CONCLUSIONS

In this work, we presented a more complete treatment of peculiar velocities for |$H_0$| inference at the level of the three-dimensional field, marginalizing over the SNeIa’s unknown cosmological redshifts and peculiar velocities. A number of studies have corrected the observed redshifts for peculiar velocities (e.g. Carr et al. 2022), and the correction is uncertain and should be marginalized over (equation 10, Peterson et al. 2022). In our approach, all assumed contributions to the observed redshifts (peculiar velocities and errors, redshift errors) are forward-modelled and accounted for in the distance modulus likelihood. Further, instead of assuming a typical dispersion of 250 km s⁻¹ for all sources, we utilize our knowledge of the constrained non-linear velocity field at 2.65 Mpc h⁻¹ resolution to provide improved velocity uncertainties and their correlations, accounting for the varying quality of the reconstruction at each location (Carr et al. 2022). The overall difference can be understood by comparing our equation (14) to Peterson et al. (equation 10, 2022). Instead of each SNIa being associated with a peculiar velocity estimate, it is associated with a full velocity posterior along its line of sight, which accounts for the quality of the reconstruction, redshift errors, and non-linear velocity correlations. All these uncertainties are part of the Bayesian modelling of the distance modulus likelihood.

Comparing the simulated SNIa host galaxy properties to observations in a similar redshift range recovers good agreement up to a mild rescaling factor, which could originate in the semi-analytic nature of GALFORM and to which our analysis is robust. We found an insignificant increase of |${\sim} 0.4 \pm 0.5$| km s|$^{-1}$| Mpc|$^{-1}$| on |$H_0$| on average in the range |$0.023<z<0.046$|⁠. The upper redshift cut is imposed by the redshift extent of the SIBELIUS-DARK simulation. This increase is weak, predominantly in the direction of the Shapley and Hydra-Centaurus superclusters. Given the diminishing impact of peculiar velocities on |$H_0$| at higher redshifts, we conclude that it is unlikely that the origin of the Hubble tension is in the assumptions on the velocity dynamics or the specific configuration of our local Universe. This conclusion is reinforced by Wu & Huterer (2017), who found a sample variance error on |$H_0$| at 0.3 km s|$^{-1}$| Mpc|$^{-1}$| in the range |$0.023<z<0.15$| using unconstrained simulations and observers hosted by haloes of varying properties.

Repeating the analysis assuming uniform SNIa rates instead of rates that depend on the star-formation history makes no significant difference to our main conclusions. This suggests that the per-host SNIa rate modelling is of secondary importance to the study of SNIa velocity dynamics, indicating that the positive |$\Delta H/H_0$| is driven by the locations of galaxies in the large-scale structure, rather than the particular star-forming locations of SNeIa. These results further indicate that once a constrained simulation of a cosmological volume is available, precise modelling of the SNIa rate – often challenging – could be less critical for studying velocity dynamics at SNIa locations. The |$H_0$| errors reported here are likely expected to increase further in the case of a real-data application, where one might want to include the effects of velocity dispersion and peculiar velocities sourced outside the 2M++ volume, but already broadly agree in magnitude with the inference of the impact of linear velocities on |$H_0$| from ZTF DR2 data (Carreres et al. 2025). This implies that the inclusion of non-linear velocity correlations leaves the conclusions on the impact of velocities on |$H_0$| qualitatively invariant for a typical SNIa sample with respect to the Carreres et al. (2025) analysis, which assumed a linear velocity covariance in the range |$0.023<z<0.06$|⁠. However, for the purposes of an accurate quantitative estimation or inferring |$H_0$| from SNeIa within our local flow (⁠|$z<0.023$|⁠), accounting for non-linear velocity correlations is necessary. This is because at |$z<0.023$| where contributions to the observed redshifts from peculiar velocities are significant, uncertainties which are not accounted for in the |$H_0$| posterior self-consistently may, in principle, introduce biases. Our method accounts both for non-linear velocity correlations and the 2M++ survey uncertainties self-consistently, combined with state-of-the-art velocity modelling (Stiskalek et al. 2025). It remains to be checked against real data whether the more accurate velocity modelling in the BORG 2M++ reconstruction will yield results on |$H_0$| at |$z<0.023$| which will be consistent with existing studies. Therefore, the effectiveness of the Bayesian treatment to deal with peculiar velocity effects allows the use, in conjunction with posterior samples of the local peculiar velocity field, of low-redshift SNe that are normally discarded from Hubble constant analyses. The uncertainties in the presented |$H_0$| posteriors are dominated by the number of SNeIa used, so extending the minimum redshift of |$H_0$| analyses to |$z\sim 0$| will reduce the variance of the |$H_0$| posteriors presented here.

ACKNOWLEDGEMENTS

We thank the anonymous referee for the useful comments, which improved the quality of the manuscript. We thank Stuart McAlpine for the SIBELIUS-DARK catalogue, Jens Jasche and Guilhem Lavaux for the 2M++ BORG inference data and their comments on the draft. We thank Josh Speagle for his insights on the performance of dynesty in high dimensions. ET would like to thank Ariel Goobar for his comments on the draft, Bruce Bassett, Boris Leistedt, and Michelle Lochner for useful discussions. This work was supported by STFC through Imperial College Astrophysics Consolidated Grant ST/W000989/1. This research utilized the HPC facility supported by the Research Computing Service at Imperial College London. ET further acknowledges support from the Centre for Cosmological Studies Balzan Fellowship, that contributed to the successful completion of this work. This work has been done within the Aquila Consortium (https://www.aquila-consortium.org).

DATA AVAILABILITY

Data products underlying this article can be made available upon reasonable request to the corresponding author.

REFERENCES