¹

1 Introduction↩︎

Since the discovery of the accelerating expansion of the Universe [1]–[3], type Ia supernovae (SNe Ia), as standardisable candles, have served as a leading probe of the cosmic expansion history. Fits to SN Ia photometric light curves provide their apparent magnitude, stretch and colour, which are used to standardise their peak brightnesses. Together with the redshifts of the SN host galaxy from spectroscopic observations, we obtain the magnitude-redshift relation, from which we can infer the expansion rate of the Universe [4]–[7].

The majority of SN Ia cosmological studies have focused on improving photometric calibration and other systematics to increase accuracy and precision of peak brightness measurements, enabling more accurate inference of the Hubble constant. Typically, the scatter about the underlying magnitude-redshift relation is treated as an uncertainty and propagated into the error budget [8]. Whilst this scatter does include a component originating from the intrinsic brightness variation of SNe Ia, it also contains coherent fluctuations due to the peculiar velocities (PVs) of the SN Ia host galaxies. The dominant contribution to these fluctuations is from peculiar velocities at low redshift, so we focus on redshift range \(z < 0.1\) in this study. There is also a contribution due to gravitational lensing, however this effect is only perceptible in higher-redshift SNe and was not considered in this study.

The peculiar velocity of a galaxy is its motion relative to the Hubble rest frame. Gravitational instability leads to overdensities whose gravitational attraction influences galaxy motion. Peculiar velocities are therefore correlated with the galaxy density field and provide a test of the cosmological model [8], [9]. The velocity field is more sensitive than the density field to matter fluctuations on the largest scales [10], making it a powerful probe of the underlying matter distribution, gravity and dark energy. Velocity fluctuations can be converted to magnitude fluctuations [11], [12] and compared with correlations observed in SN datasets [13]–[18]. In a complementary approach, [19] constrained the cosmological model by studying the moments of the supernova magnitude distribution. Similar information is imprinted in the cross-correlation between SN magnitude fluctuations and large-scale structure, which we consider in this study.

The \(\Lambda\)CDM model of cosmology describes a Universe dominated by dark energy \(\Lambda\) and cold dark matter (CDM), in which gravity is characterised by general relativity at all scales. However, the lack of a theoretical justification for dark energy continues to motivate alternative cosmological models. The growth rate of structure \(f\) is determined by the evolution of density perturbations under the influence of gravity. Peculiar velocities constrain a degenerate combination of \(f\), the growth rate of structure, and \(\sigma_8\), the amplitude of density fluctuations on scales of \(8 \, h^{-1}\) Mpc, serving as a probe of the underlying gravitational physics and offering constraints on modified gravity theories [20]. Because of this degeneracy it is common for peculiar velocity measurements to constrain a combined parameter referred to as the normalised growth rate of structure, \(f \sigma_8\).

Peculiar velocity correlations have been previously measured by a range of “standard candle” probes and methodologies. SNe Ia, as standardisable candles, enable more accurate inference of peculiar velocities because of their smaller intrinsic scatter in standardised peak brightness in comparison to other independent distance indicators such as the Tully-Fisher [21] and Fundamental Plane [22], [23] methods. However, due to the historically limited sample size and sky coverage of existing SN Ia surveys, and the inhomogeneous sky coverage of SN Ia compilations, inference of \(f\sigma_8\) using SNe Ia [15], [17], [24] has been less extensive when compared to the other distance indicators.

Methods for constraining \(f \sigma_8\) using peculiar velocities include: the two-point correlation function, [25]–[31]; the maximum-likelihood fields method, where the covariance between the peculiar velocity and galaxy overdensity fields is modelled analytically [15], [17], [32]–[35]; measuring the momentum power spectrum [36]–[40]; and reconstructions of the peculiar velocity field [24], [41]–[47]. As a whole, these methods have produced growth rate measurements that are largely in agreement with the predictions of the \(\Lambda\)CDM cosmology, though with significant \(\sim 20\%\) errors in most cases. Given the current limited samples of SNe Ia, it is not expected they will currently have constraining power competitive with the methods outlined above, but these constraints will continue to improve as new SNe Ia are detected.

The goal of this paper is to perform new \(f \sigma_8\) measurements using the correlations between SN Ia magnitudes and large-scale structure. We use SNe Ia from the Pantheon+ compilation [48], one of the largest samples of SNe Ia currently available, and for our large-scale structure maps, we use the Bright Galaxy Survey (BGS) component of Data Release 1 (DR1) of the Dark Energy Spectroscopic Instrument [49], [50]. DESI is a galaxy-redshift survey being conducted on the 4m Mayall Telescope, at the Kitt Peak National Observatory. The DESI BGS offers an order-of-magnitude increase in the number of large-scale structure tracers over the previous generation of surveys, resulting in an increase in the significance of the potential correlation signal. This is the first study using the DESI density field in a cross-correlation analysis with SNe Ia. We measure the correlation statistics between these datasets using estimators of the magnitude auto-correlation functions \(\psi_1\) and \(\psi_2\) [51], and the galaxy-magnitude cross-correlation function \(\psi_3\) [52], and test \(\Lambda\)CDM by fitting \(f\sigma_8\) to these measurements. We validate our models using simulations, and forecast how future large-scale structure and SN Ia surveys could improve our constraints on \(f\sigma_8\).

With new generation surveys such as the Zwicky Transient Facility [53], [54] and the Rubin Observatory Legacy Survey of Space and Time [55] currently on-sky or in development, we will observe future SN samples with high cadence and large sky coverage. From this, a more comprehensive and uniform sample of SNe Ia will soon be available for peculiar velocity studies [18], [56]. Local large-scale structure datasets offer direct peculiar velocity measurements and the DESI PV survey plans to measure the peculiar velocities of 186,000 galaxies over its 5 years of operation [57]. This will be further expanded by surveys such as the 4-metre Multi-Object Spectroscopic Telescope [58] Hemisphere Survey [59] and the Australian Square Kilometre Array Pathfinder WALLABY survey [60]. The influx of data from these surveys will improve our distance estimates and combining SN Ia, Tully-Fisher and Fundamental Plane peculiar velocity measurements will provide tighter constraints on \(f\sigma_8\) and, consequently, modified gravity models [61], [62].

Our paper is structured as follows. The theoretical models we use to describe the magnitude fluctuations induced by peculiar velocities are outlined in Sec. 2, and the datasets we utilise from DESI and Pantheon+ are described in Sec. 3. The N-body simulations we created to test the theoretical models and the impact of selection effects on our analysis are discussed in Sec. 4, and the magnitude auto-correlation functions \(\psi_1\) and \(\psi_2\), galaxy-magnitude cross-correlation function \(\psi_3\), and galaxy auto-correlation function \(\xi_{gg}\) are measured in Sec. 5. The growth rate analysis validation using simulations is shown in Sec. 6, and the \(f\sigma_8\) constraint is presented in Sec. 7. We summarise our conclusions and forecast the impact of new large-scale structure and SN Ia surveys in Sec. 8.

2 Magnitude Fluctuation Correlation Theory↩︎

2.1 Growth rate of structure↩︎

The growth rate of structure \(f\) describes the evolution of density perturbations under the influence of gravity. In the linear theory of structure formation, density fluctuations \(\delta\) grow as \(\delta \propto D(a)\), where \(D(a)\) is the linear growth factor, dependent on cosmic scale factor \(a\). The growth rate is then defined as \[f(a) = \frac{d\ln{D(a)}}{d\ln{a}}.\]

Measuring \(f\) allows us to test gravitational physics and constrain modified gravity and dark energy models. The growth rate at redshift \(z\) can be parameterised as \[f(z) = \Omega_m(z)^{\gamma},\] where \(\Omega_m\) is the matter density parameter and \(\gamma\) is the growth index, as shown by [63]. In the \(\Lambda\)CDM model, \(\gamma = {0.55}\) [64]. The value of \(\gamma\) varies in different gravity and dark energy models. For the Dvali, Gabadadze \(\&\) Porrati braneworld model [65], \(\gamma = 0.6875\) [64], and in \(f(R)\) gravity models, which include a non-linear function of the Ricci scalar \(R\) in the Einstein-Hilbert action, the growth rate is dependent on scale [66]. The growth rate modulates the relation between the local peculiar velocity \(\mathbf{v}(\mathbf{x})\) and matter overdensity \(\delta_m(\mathbf{x})\), at some point \(\mathbf{x}\), via the linear theory continuity equation,

\[\mathbf{\nabla} \cdot \mathbf{v}(\mathbf{x}) = -a H f \, \delta_m(\mathbf{x}),\] where \(H\) is the Hubble parameter.

2.2 Relation between magnitude and velocity fluctuations↩︎

In this section, we derive the relation for the magnitude fluctuation imprinted by the velocity of a source [11], [12]. We start from the relation between apparent magnitude \(m\) and luminosity distance \(D_L\) in Mpc, given by the distance modulus equation, \[m = M + 5 \log_{10}{D_L} + 25, \label{eq:dm}\tag{1}\] where \(M\) is the absolute magnitude. From Eq. 1 , we can show that the relation between fluctuations in magnitude \(\delta m\) and luminosity distance \(\delta D_L\) is, \[\delta m = \frac{5}{\ln{10}} \frac{\delta D_L}{D_L} . \label{eq:dmdl}\tag{2}\] The luminosity distance fluctuations \(\delta D_L\) can be determined by relating observations in a perturbed universe to those in a homogeneous Universe (for which quantities will be denoted with an overline). The perturbation in redshift is given by, \[1 + z = (1 + \overline{z}) \left( 1 + \frac{v_r}{c} \right),\] where \(v_r\) is the radial velocity of the source. The relation between the luminosity distance \(D_L\) and angular diameter distance \(D_A\) applies in both perturbed and unperturbed Universes [11], [67], \[\begin{align} D_L(z) &= D_A(z) (1 + z)^2 , \\ \overline{D}_L(\overline{z}) &= \overline{D}_A(\overline{z}) (1 + \overline{z})^2. \end{align}\] Neglecting observer motion, we can write \(D_A(z) = \overline{D}_A(\overline{z})\), hence, \[\frac{D_L(z)}{\overline{D}_L(\overline{z})} = \frac{(1+z)^2}{(1+\overline{z})^2} = \left( 1 + \frac{v_r}{c} \right)^2 \approx 1 + \frac{2v_r}{c}. \label{eq:dl2}\tag{3}\] We now perform a Taylor series expansion about \(z = \overline{z}\), \[\overline{D}_L(z) = \overline{D}_L(\overline{z}) + \frac{\partial \overline{D}_L}{\partial z} \big\vert_{z = \overline{z}} \left( z - \overline{z} \right) .\] Using \(\overline{D}_L(z) = \chi(z) \, (1+z)\) and \(d\chi/dz = c/H(z)\) in terms of radial co-moving co-ordinate \(\chi(z)\), we find, \[\overline{D}_L(z) = \overline{D}_L(\overline{z}) \left\{ 1 + \left[ 1 + \frac{c \, (1 + \overline{z})^2}{H(\overline{z}) \, D_L(\overline{z})} \right] \frac{v_r}{c} \right\} . \label{eq:dl3}\tag{4}\] Substituting Eq. 3 and Eq. 4 into Eq. 2 produces the magnitude fluctuation \(\delta m\) in terms of the radial velocity, which can be conveniently written, \[\delta m = \alpha(z) \times v_r , \label{eq:alpha}\tag{5}\] where \(\alpha\) is the coefficient of proportionality between velocity and magnitude, \[\alpha(z) = \frac{5}{c \, \ln{10}} \left[ 1 - \frac{c \, (1+z)^2}{H(z) \, D_L(z)} \right] . \label{eq:veltomag}\tag{6}\] Using Eq. 5 , we can relate the measurement of magnitude fluctuation correlations to velocity correlations.

2.3 Galaxy and velocity correlation functions↩︎

As shown in Sec. 2.2, magnitude fluctuations can be generated by the peculiar motions of galaxies. By modelling the velocity correlations induced by the growth of large-scale structure, we can predict the observed correlation spectrum of magnitudes. For a Gaussian peculiar velocity field, the two-point correlation tensor between two positions \(A\) and \(B\) is given by [51], [68] as, \[\label{eq:corrtensor} \Psi_{ij}(\mathbf{r}_A, \mathbf{r}_B) = \langle v_i(\mathbf{r}_A) \, v_j(\mathbf{r}_B) \rangle,\tag{7}\] where \(\mathbf{r}\) is a spatial position and \(v_i\) is the peculiar velocity component in directions \(i = \{ x, y, z \}\). For an irrotational, homogeneous and isotropic velocity field with linear velocity perturbations, the velocity correlation tensor can be written as [51], [68], [69], \[\Psi_{ij}(r) = \left[\Psi_{\parallel}(r) - \Psi_{\perp}(r)\right] \, \hat{r}_{Ai} \, \hat{r}_{Bj} + \Psi_{\perp}(r) \, \delta^K_{ij},\] where \(r\) is the separation between positions \(A\) and \(B\), \(\Psi_{\parallel}(r)\) and \(\Psi_{\perp}(r)\) are the correlation functions between components of velocity parallel and perpendicular to the separation vector \(\mathbf{r}\), and \(\delta^K_{ij}\) is the Kronecker delta.

The correlation functions \(\Psi_{\parallel}\) and \(\Psi_{\perp}\) can be written in terms of the matter power spectrum \(P(k)\) as, \[\label{eq:corrperppar} \begin{align} \Psi_{\parallel}(r) &= \frac{H^2 a^2 (f\sigma_8)^2}{2\pi^2} \int \frac{P(k)}{\sigma_{8,\mathrm{fid}}^2} \left[j_0(kr) - 2\frac{j_1(kr)}{kr}\right] \, dk, \\ \Psi_{\perp}(r) &= \frac{H^2 a^2 (f\sigma_8)^2}{2\pi^2} \int \frac{P(k)}{\sigma_{8,\mathrm{fid}}^2} \frac{j_1(kr)}{kr} \, dk, \end{align}\tag{8}\] where \(P(k)\) is the model matter power spectrum as a function of wavenumber \(k\), \(\sigma_{8,\mathrm{fid}}\) is the chosen normalisation of the matter power spectrum, and \(j_\ell\) are spherical Bessel functions of the first kind, \[\begin{align} j_0(x) &= \frac{\sin(x)}{x}, \\ j_1(x) &= \frac{\sin(x)}{x^2} - \frac{\cos(x)}{x}. \end{align}\] The formulation of Eq. 8 illustrates that the measured velocity correlation function may be used to constrain \(f\sigma_8\).

Peculiar velocities are also correlated with galaxy positions, where the cross-correlation function is given by [32], \[\xi_{gv}(r) = -\frac{H a (f \sigma_8) (b \sigma_8)}{2\pi^2} \int dk \, k \, \frac{P(k)}{\sigma_{8,\mathrm{fid}}^2} \, j_1(kr) , \label{eq:galvelcorr}\tag{9}\] where \(b\) is the linear galaxy bias describing how galaxy distribution traces the underlying matter density field, and \(b\sigma_8\) is the the normalised linear galaxy bias. Finally, the galaxy auto-correlation function monopole model (neglecting redshift-space distortions for the moment) is, \[\xi_{gg}(r) = \frac{(b \sigma_8)^2}{2\pi^2} \int dk \, k^2 \, \frac{P(k)}{\sigma_{8,\mathrm{fid}}^2} \, j_0(kr) . \label{eq:galcorr}\tag{10}\] We generate the linear model power spectrum \(P(k)\) using CAMB [70], and use the halofit non-linear matter power spectrum [71] to model the velocity and galaxy correlation functions.

Since we wish to focus on the growth information contained in SN magnitude correlations in this study, we do not consider the galaxy auto-correlation quadrupole, which adds significant additional information through the effect of redshift-space distortions (RSD). We include the effect of RSD on the remaining statistics by scaling the galaxy auto-correlation monopole by a linear RSD factor \((1 + \frac{2}{3}\beta + \frac{1}{5} \beta^2)\), where \(\beta = f/b\), and the galaxy-velocity cross-correlation function by \((1 + \frac{1}{3}\beta)\) [29]. As we are restricting our model fits to large scales, we find that linear bias and RSD models are acceptable, which we test using mock catalogues in Sec. 6.

3 Density field and magnitude fluctuation datasets↩︎

3.1 The Dark Energy Spectroscopic Instrument↩︎

DESI is using the 4m Mayall Telescope, at the Kitt Peak National Observatory, to perform a spectroscopic survey targeting galaxies and quasars to measure the expansion history of the Universe and the growth of large-scale structure [49], [72]–[74]. The DESI Survey is planned to span five years, during which it will observe more than 40 million galaxies and quasars over a footprint of 14,000 square degrees. The survey is divided into four main target classes: Bright Galaxy Survey (BGS), Luminous Red Galaxy Survey (LRG), Emission Line Galaxy Survey (ELG), and Quasar Survey (QSO), which are selected from the DESI Legacy optical imaging surveys [75]. Targets within each observational field are assigned to 5000 optical fibres [76] in the telescope focal plane using a robotic positioner [77] and optical corrector [78], and the observed data are processed by the DESI spectroscopic pipeline [79]. The overall observing strategy is summarised by [80]. DESI has issued an Early Data Release [81], which has been used for scientific validation of the program [82]. In this study, we use a galaxy sample from the publicly available DESI Data Release 1 [83], comprising of spectra obtained during the first year of full-survey observations.

3.2 The Bright Galaxy Survey dataset↩︎

We draw our large-scale structure data from the DESI Bright Galaxy Survey (BGS), which targets a magnitude-limited sample of galaxies with \(r\)-band magnitudes \(14 < r < 19.5\) across redshifts \(z < 0.5\). BGS targets are assigned high priority during DESI bright-time observations. The full selection criteria and validation of the BGS sample are described by [50], resulting in a sample with density 854 deg\(^{-2}\) containing reliable redshift measurements for over 5.5 million galaxies in DR1. The BGS sample consists of a magnitude-limited Bright sample with \(r < 19.5\), and a colour-selected Faint component with \(19.5 < r < 20.175\). For this study, we will only consider the BGS Bright sample, as the Faint sample suffers from complications regarding incompleteness and systematics [50].

Large-scale structure catalogues suitable for cosmological analysis have been built from the redshift and target catalogues for each DESI tracer as described by [84]. The selection function of the sample is defined, along with correction weights designed to compensate for systematic density variations with observational characteristics across the survey. These selection functions yield accompanying unclustered random catalogues which are used in correlation function estimators. Measurements of Baryon Acoustic Oscillations in the clustering pattern of the DR1 BGS sample are presented by [85].

To create a sub-sample overlapping with low-redshift supernovae, we restrict our analysis to the redshift range \(z < 0.1\), resulting in a sample of 578,576 galaxies. The sky positions of the BGS galaxies are shown in Fig. 1 and the redshift distribution of the sample is shown in Fig. 2.

3.3 Pantheon+ dataset↩︎

The Pantheon+ dataset contains 1701 light curves of 1550 distinct spectroscopically-confirmed type Ia supernovae (SNe Ia) in redshift range \(0.001 < z < 2.261\) [48]. A compilation of 18 supernova surveys, the Pantheon+ catalogue standardises SN Ia light curves using a consistent magnitude cross-calibration [86]. The distance modulus \(\mu\) for a SN Ia is determined by fitting the light curves with the SALT2 model originally developed by [87] and updated in [86]. This results in measurements of the peak magnitude or light curve amplitude \(m_B\), the stretch parameter \(x_1\) and light-curve color \(c\). These light-curve fit parameters are related to the distance modulus using the modified Tripp relation [88], [89]: \[\mu = m_B + \alpha x_1 - \beta c - M - \delta_{\mu \mathrm{-bias}} , \label{eq:distmod}\tag{11}\] where the correlation coefficients \(\alpha = 0.14\) and \(\beta=3.1\), \(M\) is the fiducial SN Ia absolute magnitude and \(\delta_{\mu \mathrm{-bias}}\) is the bias correction derived from simulations which account for selection effects and systematics involved in distance recovery [90]. Cosmological constraints from the Pantheon+ dataset are presented by [5].

As magnitude fluctuations are primarily due to peculiar velocities at low redshift, we restrict the SN sample to \(z < 0.1\), which provides 510 unique SNe Ia. The Pantheon+ SN Ia sky positions and redshift distributions are shown in Fig. 1 and Fig. 2, respectively. From the Pantheon+ SN Ia dataset² we used the CMB-frame redshift, right ascension, declination, magnitude, and magnitude error of each SN Ia. We determined the magnitude fluctuation \(\delta m\) in each case by subtracting the fiducial magnitude-redshift relation assuming a flat \(\Lambda\)CDM Universe with \(H_0 = 100 \, h\) km s⁻¹ Mpc⁻¹, \(\Omega_m = 0.3153\), \(\Omega_{\Lambda} = 0.6847\), where the density parameters represent the DESI fiducial cosmology, but do not have a significant effect for our low-redshift sample. We subtract the ensemble mean from the \(\delta m\) values (effectively fitting for the absolute magnitude \(M - 5 \log_{10} h\)).

4 Abacus N-body Simulations↩︎

In this section we summarise the N-body simulation datasets we used to test and validate our methodology, theoretical models, and the impact of supernovae selection effects on our analysis. We do this by utilising mocks matched to the DESI BGS sample, including a realistic population of SN host galaxies. These mocks include realistic tracer selection functions, measurement errors and galaxy bias. We constructed these mock catalogues from the \(z=0.2\) snapshot of the AbacusSummit suite of N-body simulations [91]–[93], the lowest-redshift output, which includes 25 independent boxes generated in the Planck 2018 best-fit flat \(\Lambda\)CDM cosmology using fiducial parameters, \(\Omega_m = 0.3153\), \(\Omega_{\Lambda} = 0.6847\), \(\Omega_b = 0.0493\), \(h = 0.6736\), \(\sigma_8 = 0.8114\) and \(n_s = 0.9649\).

In order to simulate the DESI BGS distribution, dark matter halos were populated with a luminosity-dependent halo occupation distribution model [94], [95], reproducing the comoving number density of galaxies as well as the projected correlation function of DESI BGS data [50]. In addition to positions and velocities, mock galaxies have observed magnitudes and colours that follow a realistic luminosity function including redshift evolution. The full details of the mock production can be found in [96] and references therein.

As we are conducting an analysis of the low-redshift Universe (\(z<0.1\)), we divided the initial 2 \(h^{-1}\) Gpc simulation boxes into \(3^3\) sub-volumes, placing the observer at the centre of each. Given that we have 25 independent Abacus boxes, this yields a total of 675 mock catalogues. While peculiar velocites are correlated on large scales, we assume all sub-volumes to be statistically independent, which is a good approximation for the purposes of this work. We apply observational completeness masks corresponding to the DESI DR1 dataset. Additional properties, such as stellar masses, star formation rates and galaxy sizes are assigned from real data to mock galaxies based on proximity in redshift-magnitude-colour space.

We populated our mock BGS galaxies with SNe Ia that have explosion rates \(R_{Ia}\) dependent on stellar mass \(M_*\) and star-formation rate (SFR), known as the “A+B model” [97]: \[\label{eq:snIa32rate} R_{Ia} = A \times M_*^{n_M} + B \times \mathrm{SFR}^{n_S} ,\tag{12}\] where we assume a linear relation (\(n_M = n_S = 1\)), \(A = 5.3 \times 10^{-14}\) yr^-1.M_^-1 and \(B = 3.9 \times 10^{-4}\) yr⁻¹ ( M yr⁻¹ )⁻¹ [98]. For a given galaxy, the number of SN Ia is a Poisson realisation of the rate multiplied by the number of years of observations. The model for the supernova rate consists of a “delayed” component, with a long delay time driven by the stellar mass of the galaxy, and a “prompt” component, with short delay times caused by the formation of new stars.

4.1 Modelling SN properties for selection effects↩︎

Supernova samples contain observational selection effects which may imprint magnitude- or redshift-dependent trends, presenting potential systematics in the interpretation of magnitude fluctuations as tracing peculiar velocities. In this section we describe how we apply these effects to our simulated Abacus samples to assess their influence on growth rate constraints. We test the influence of SN selection effects by performing growth rate fits with and without applying a range of magnitude selections to our simulated samples. Using Eq. 11 , we model the rest-frame absolute magnitude of each SN Ia (labelled \(i\)) as, \[M^*_i = M_B - \alpha x_{1,i} + \beta c_i + \sigma_{\mathrm{int},i}, \label{eq:absmagstar}\tag{13}\] where the stretch \(x_1\), colour \(c\), and intrinsic scatter \(\sigma_{\mathrm{int}}\) are randomly drawn from statistical distributions, and we use best-fit values of \(\alpha = 0.14\), \(\beta=3.1\) and \(B\)-band absolute magnitude \(M_B = -19.05\) [99]. We rescaled \(M\) to match the Abacus fiducial cosmology using, \[M = -19.05 + 5 \log_{10} \left( \frac{h}{0.7} \right).\]

We modelled the \(x_1\) distribution using a bimodal Gaussian mixture considering the evolution of younger and older SN Ia progenitors as a function of redshift [100], [101]. The \(x_1\) distribution is given by, \[\begin{align} & x_1(z) = \delta(z) \times \mathcal{N}(\mu_1, \sigma^2_1) \, + \\ & (1 - \delta(z)) \times \left[ a \times \mathcal{N}(\mu_1, \sigma^2_1) + (1 - a) \times \mathcal{N}(\mu_2, \sigma^2_2) \right] , \end{align}\] where \(a\) accounts for the relative effect of the younger and older SN Ia progenitors, and the evolving fraction of young SNe Ia is given by \(\delta(z) = (K^{-1} \times (1+z)^{-2.8} + 1)^{-1}\), with coefficient \(K = 0.87\) as outlined by [100]. We used parameters \(\mu_1 = 0.37\), \(\sigma_1 = 0.61\), \(\mu_2 = -1.22\), \(\sigma_2 = 0.56\) and \(a = 0.51\) [101] to determine the stretch parameter for each SN Ia.

The colour parameter \(c\) is modelled as an asymmetric Gaussian distribution from [89] following the low redshift (G10) model, where \(\bar{c} = -0.054\) is the maximum probability colour parameter, \(\sigma_{-} = 0.043\) and \(\sigma_{+} = 0.101\) are the asymmetric Gaussian widths of the distribution. The colour distribution follows, \[\begin{align} P(c) \propto \exp \left[ - \frac{(c - \bar{c})^2}{2 \sigma_{-}^2} \right], \, c \leq \bar{c} , \\ P(c) \propto \exp \left[ - \frac{(c - \bar{c})^2}{2 \sigma_{+}^2} \right], \, c > \bar{c} . \end{align}\] Finally, the intrinsic scatter \(\sigma_{\mathrm{int}}\) is drawn from a Gaussian distribution with dispersion \(\sigma_{M} = 0.12\). After drawing the values of \(x_1\), \(c\) and \(\sigma_{\mathrm{int}}\) for each SN Ia from the distributions described above, the absolute magnitude in Eq. 13 is used to calculate the apparent magnitude in Eq. 1 .

The Pantheon+ compilation consists of 18 individual sub-surveys, resulting in a complex angular and depth selection. For this work, we approximated the angular selection of Pantheon+ using a Healpix map of the distribution of the Pantheon+ dataset and using this map, sub-sampled the mock and random catalogues. Modelling the selection of each individual sub-survey would improve the robustness of this type of analysis and, although beyond the scope of this project, would be worth considering for potential future analysis. However, to test the robustness of our results to the choice of magnitude threshold, we considered applying different magnitude cuts at \(m < 16.0\), \(m < 16.5\), and \(m < 17\), which are representative of the depth of the different samples. This test is hence a first step to consider whether selection effects impact our type of correlation measurement.

5 Correlation function measurements↩︎

In this section we describe our estimators for the auto- and cross-correlations of magnitude fluctuations and galaxy positions, which may be linked to the underlying velocity correlation theory using the models outlined in Sec. 2.

5.1 Magnitude auto-correlation functions↩︎

As discussed in Sec. 2.2, magnitude fluctuations are linked to the underlying velocity fluctuations by the redshift-dependent scaling \(\alpha(z)\) defined in Eq. 6 . Predicting the magnitude fluctuation correlation function requires us to normalise the model velocity correlations by the average product \(\langle \alpha_A \, \alpha_B \rangle\) between a pair of galaxies \(A\) and \(B\) within the separation range, which we also measure during our correlation estimation. In Sec. 6.1 we will test the appropriate separation range to include when fitting the model to the data.

The magnitude correlation between two galaxies depends on the correlation between their radial velocity components. Re-writing Eq. 7 as the correlation for line-of-sight velocities of two galaxies, \(u_A\) and \(u_B\), leads to, \[\begin{align} & \langle u_A(\mathbf{r}_A) \, u_B(\mathbf{r}_B) \rangle = \\ & \Psi_{\perp}(r) \cos\theta_{AB} + \left[ \Psi_{\parallel}(r) - \Psi_{\perp}(r) \right] \cos\theta_A \cos\theta_B , \end{align} \label{eq:corrradial}\tag{14}\] where \(\theta_A\) and \(\theta_B\) are the angles between the lines-of-sight to the two objects and the separation vector, and \(\theta_{AB}\) is the angle subtended by the two lines-of-sight at the origin. Converting Eq. 14 to the corresponding magnitude fluctuations \(\delta m_A\) and \(\delta m_B\) we find, \[\label{eq:corrmag} \langle \delta m_A(\mathbf{r}_A) \, \delta m_B(\mathbf{r}_B) \rangle = \alpha_A \, \alpha_B \langle u_A(\mathbf{r}_A) \, u_B(\mathbf{r}_B) \rangle .\tag{15}\]

The dependence of Eq. 14 on the orientation of the galaxy pair with respect to the line-of-sight shows that the auto-correlation spectrum must be characterised by two separate functions [69]. [68] demonstrated that the information in Eq. 14 can be written in terms of the two correlation statistics \(\psi_1\) and \(\psi_2\). Using magnitude fluctuations as the variable, these may be estimated as, \[\begin{align} \begin{aligned} \tag{16} \hat{\psi}_1(r) ={}& \frac{\sum w_A \, w_B \, \delta m_A \, \delta m_B \, \cos\theta_{AB}}{\sum w_A \, w_B \, \cos^2\theta_{AB}} , \end{aligned} \\ \begin{align} \tag{17} \hat{\psi}_2(r) ={}& \frac{\sum w_A \, w_B \, \delta m_A \, \delta m_B \, \cos\theta_A \, \cos\theta_B}{\sum w_A \, w_B \, \cos\theta_{AB} \, \cos\theta_A \, \cos\theta_B} , \end{align} \end{align}\] where sums are over all pairs of galaxies separated by some fixed distance \(r\), and galaxies are weighted by \(w_A\) and \(w_B\), discussed further below. These estimators for \(\psi_1\) and \(\psi_2\) can be related to the velocity correlation functions \(\Psi_{\parallel}\) and \(\Psi_{\perp}\) using Eqs. 14 , 15 , 16 and 17 to give, \[\langle \psi_1(r) \rangle = \frac{N^{\mathrm{auto}}_{\alpha\alpha}}{N^{\mathrm{auto}}_{gg}} \, \left\{ \mathcal{A}(r) \psi_{\parallel}(r) + \left[ 1 - \mathcal{A}(r) \right] \psi_{\perp}(r) \right\},\] and, \[\langle \psi_2(r) \rangle = \frac{N^{\mathrm{auto}}_{\alpha\alpha}}{N^{\mathrm{auto}}_{gg}} \, \left\{ \mathcal{B}(r) \psi_{\parallel}(r) + \left[ 1 - \mathcal{B}(r) \right] \psi_{\perp}(r) \right\},\] where \(\mathcal{A}(r)\) and \(\mathcal{B}(r)\) are geometry factors dictating the relative contributions of the parallel and perpendicular components of the velocity field to the correlation functions \(\psi_1\) and \(\psi_2\), which are given by [51], \[\mathcal{A}(r) = \frac{\sum w_A \, w_B \, \cos \theta_A \, \cos \theta_B \, \cos \theta_{AB}}{\sum w_A \, w_B \cos^2 \theta_{AB}},\] \[\mathcal{B}(r) = \frac{\sum w_A \, w_B \, \cos \theta_A \, \cos \theta_B \, \cos \theta_{AB}}{\sum w_A \, w_B \, \cos^2\theta_A}.\] The normalisation factor \(N^{\mathrm{auto}}_{gg}\) is evaluated by taking the velocity-velocity pair count in each separation bin weighted by the velocity sample weight \(w_v\) for each galaxy in the pair. The factor \(N^{\mathrm{auto}}_{\alpha\alpha}\) is evaluated by taking the velocity-velocity pair count weighted by \(w_v\) and the values of \(\alpha(z)\) for each galaxy.

The optimal weights for each galaxy are chosen to minimize the statistical error in the 2-point correlation function, balancing sample variance and measurement noise [102], given the varying number density of galaxies across the survey volume. The optimal weight for the galaxy sample is, \[w_{g} = \frac{1}{1 + n_g \, P_g},\] where \(n_g\) is the number density in units of ̉³ Mpc⁻³ of the galaxy sample at the position in question, and the characteristic galaxy power spectrum amplitude is taken as \(P_g = 10^4\) ̉⁻³ Mpc³. The optimal weight for the velocity sample is [29], \[w_{v} = \frac{1}{\sigma^2_m / \alpha(z) + \alpha(z) \, n_v \, P_v},\] where \(n_v\) is the number density of the velocity sample in ̉³ Mpc⁻³, the characteristic velocity power spectrum amplitude is taken as \(P_v = 10^{10}\) ̉⁻³ Mpc³ ( km s⁻¹ )², \(\sigma_m\) is the supernova-specific magnitude error, and \(\alpha(z)\) is the redshift-dependent scaling factor defined in Eq. 6 , evaluated for each galaxy.

5.2 Galaxy-magnitude cross-correlation function↩︎

Converting Eq. 9 to a line-of-sight velocity, the cross-correlation function between a galaxy overdensity at position \(A\) and magnitude fluctuation at position \(B\) is, \[\langle \delta_A (\mathbf{r}_A) \, \delta m_B (\mathbf{r}_B) \rangle = \alpha_B \, \xi_{gv}(r) \cos\theta_B.\] The estimator for the cross-correlation function between the galaxy overdensity and magnitude fluctuations can then be defined following [52] as, \[\hat{\psi_3}(r) = \frac{\sum w_A \, w_B \, \delta m_B \, \cos\theta_B}{\sum w_A \, w_B \, \cos^2\theta_B} . \label{eq:psi3}\tag{18}\] We form the complete cross-correlation estimator by subtracting from Eq. 18 the analogous cross-pair count between the velocity sample and a distribution of unclustered random galaxies with the same distribution as the galaxy data. In this way, the variance of the galaxy-magnitude cross-correlation estimator is significantly reduced [52].

To relate the theoretical galaxy-velocity cross-correlation model from Eq. 9 to the galaxy-magnitude cross-correlation function measurements, we normalise the model as, \[\langle \hat{\psi_3}(r)\rangle = \frac{N^{\mathrm{cross}}_{g\alpha}}{N^{\mathrm{cross}}_{gg}} \, \xi_{gv}(r) .\] The normalisation factor \(N^{\mathrm{cross}}_{gg}\) is evaluated by taking the cross-pair count of the galaxy and velocity samples in each separation bin weighted by the galaxy and velocity optimal weights, \(w_g\) and \(w_v\), respectively. The factor \(N^{\mathrm{cross}}_{g \alpha}\) is evaluated by taking the galaxy-velocity pair count weighted by the optimal weights of each pair and an additional factor of \(\alpha(z)\).

5.3 Galaxy auto-correlation function↩︎

The galaxy auto-correlation function estimator is given by, \[\hat{\xi}_{gg}(r) = \frac{N^2_R}{N^2_D} \frac{D_g D_g(r)}{R_g R_g(r)} - 2 \frac{N_R}{N_D} \frac{D_g R_g(r)}{R_g R_g(r)} + 1, \label{eq:xigg}\tag{19}\] [103], where \(D_g D_g(r)\) and \(R_g R_g(r)\) are the galaxy-galaxy pair count at separation \(r\), for the data and random samples, respectively, and \(D_g R_g(r)\) is the cross-pair count between the samples. \(N_D\) and \(N_R\) are the (weighted) number of galaxies in the data and random sample, respectively and are used here to normalise the pair count. The inclusion of the second term in Eq. 19 decreases the statistical error associated with the distribution of the data with respect to sample boundaries [103].

6 Validation using simulations↩︎

In this section we present the suite of tests performed using N-body simulations to validate the recovery of the fiducial growth rate using our correlation measurements and models.

6.1 Validation by fitting to Abacus simulations (without SN selection effects)↩︎

As described in Sec. 4, the Abacus mock catalogues were constructed by applying a selection function matching the DESI BGS and a maximum redshift of \(z = 0.1\). The SN mock catalogues were matched to the Pantheon+ dataset by sub-sampling the parent catalogues by redshift and sky position, and an intrinsic dispersion \(\sigma_{\mathrm{int}} = 0.126\) was added to the SN magnitudes, which (after the inclusion of peculiar velocities) matches the standard deviation of the observed Pantheon+ magnitudes around the fiducial magnitude-redshift relation.

For our first test, we excluded supernovae selection effects. For each of the 675 mock catalogues we measured the correlation functions \(\left( \psi_1, \psi_2, \psi_3, \xi_{gg} \right)\) in a separation range of \(25 - 100\) ̉⁻¹ Mpc in 5 ̉⁻¹ Mpc bins, where the minimum separation is an estimate of the range of applicability of linear theory (which we test further below), and the maximum separation marks a reduction in the number of available SN pairs. In Fig. 3 we display these statistics averaged over the catalogues. We also used this ensemble of measurements to construct a numerical covariance matrix of the statistics.

To find the best-fit normalised growth rate \(f \sigma_8\) and galaxy bias \(b \sigma_8\), we fit the correlation function measurements using the likelihood function \(\mathcal{L} \propto \exp(-\chi^2/2)\) in terms of the \(\chi^2\) function, \[\chi^2 = \sum_i \sum_j (d_i - m_i)^T \, \left( C^{-1} \right)_{ij} \, (d_j - m_j) , \label{eq:chi2}\tag{20}\] where \(d_i\) and \(m_i\) are concatenations of the four correlation function datasets and models, respectively, and \(C\) is the covariance matrix, and minimising the \(\chi^2\) statistic for each realisation.

The histogram of the best-fit \(f \sigma_8\) and \(b \sigma_8\) values across the realisations is displayed in the top panel of Fig. 4, and the mean best-fit \(f \sigma_8\) values with the \(\chi^2\) and degrees of freedom are shown in Table 1 for different combinations of correlation functions. The distribution of reduced \(\chi^2\) values for the fits to each of the correlation functions is shown in the bottom panel of Fig. 4. All correlations result in acceptable \(\chi^2\) statistics, validating the sufficiency of our covariance and modelling approach. The growth rate fit for the combined correlation function measurements gives \(f \sigma_8 = 0.425^{+0.089}_{-0.150}\) where we quote the mean and standard deviation across 675 realisations. The best-fit galaxy bias parameter was \(b \sigma_8 = 1.042^{+0.065}_{-0.084}\). Agreement between the measurements and model across the four correlation functions demonstrates that our methodology is capable of recovering the fiducial value \(f \sigma_8 = 0.428\) within the statistical error margin, in the presence of realistic survey selection functions, galaxy bias and measurement noise.

In our fiducial analysis we restricted the fitting range of the correlation functions to \(25-100\) ̉⁻¹ Mpc to mitigate the impact of non-linear growth violating our linear modelling assumptions. To determine the impact of this scale cut on the growth rate measurements, we also fit to the combined four correlation functions of the mocks varying the minimum separation in the range \(10-45\) ̉⁻¹ Mpc, finding that the best-fit \(f \sigma_8\) values were consistent across different separation ranges. The variations between each separation range are smaller than the individual statistical uncertainty for the analysis, such that the choice of minimum separation bin leads to a statistically indistinguishable result, as shown in Fig. 5. This demonstrates that our growth rate fits are robust to the choice of fitting range. At the maximum fitting range of 100 ̉⁻¹ Mpc, the SN number density has diminished (Fig. 2) and the low quantity of SN pairs limits the potential signal-to-noise at larger scales.

6.2 Validation including SN selection effect↩︎

To test the impact of SN selection effects on magnitude correlations, we used SN catalogues with and without selection effects applied, as described in Sec. 4.1. We apply magnitude limits of \(m = (16.0, 16.5, 17.0)\) and fit our correlation function models to the 675 Abacus realisations. In each case, we measured the mean difference in the best-fitting \(f\sigma_8\) with and without selection effects applied. For the three magnitude thresholds \(m = (16.0, 16.5, 17.0)\) we found a mean difference in \(f \sigma_8\) of \((0.027, 0.016, 0.0064)\). These offsets are small compared to the statistical errors, demonstrating that SN selection effects do not significantly bias our growth rate measurements at the level of precision of these datasets.

7 Growth rate fits to DESI and Pantheon+↩︎

Having validated our analysis pipeline using mock catalogues, we now apply it to the DESI DR1 and Pantheon+ datasets. In the redshift range \(z < 0.1\), our sample consists of 578,576 DESI BGS galaxies and 510 Pantheon+ SNe Ia. We measured the auto- and cross-correlation statistics between the two samples in the separation range \(25 - 100\) ̉⁻¹ Mpc in 5 ̉⁻¹ Mpc bins. These correlation measurements are shown in Fig. 6.

We fit our correlation models varying \(f \sigma_8\) and \(b \sigma_8\) as before, evaluating \(\chi^2\) using Eq. 20 . We show the marginalised posterior probability distributions for \(f \sigma_8\) and \(b \sigma_8\) from the combined four-statistic fit to the data catalogues in Fig. 7. The best-fit \(f \sigma_8 = 0.384^{+0.094}_{-0.157}\) agrees with the prediction of the Planck cosmology \(f \sigma_8 = 0.429\), and the best-fit bias parameter is \(b \sigma_8 = 0.998^{+0.035}_{-0.036}\). For the best-fit \(f \sigma_8\) and \(b \sigma_8\), the minimum \(\chi^2 = 62.6\), consistent with the expected value from a fit with 58 degrees of freedom. We note that full-shape galaxy clustering using the full DESI BGS sample measured \(f \sigma_8 = 0.38 \pm 0.09\) [104] using redshift-space distortion information, and our constraints, which do not use redshift-space distortion information in the galaxy correlation function, are consistent with this measurement. In Table 2 we list the best-fitting growth rates and minimum \(\chi^2\) values for different combinations of correlation functions, demonstrating that the models provide a good fit to all correlations. The effect of the fitting range on the \(f\sigma_8\) measurements is shown for the DESI+Pantheon correlations in Fig. 5, which demonstrates that the minimum separation does not significantly affect the \(f\sigma_8\) values.

As a measure of the extent to which we “detect” the magnitude correlations in our samples, we considered the \(\Delta \chi^2\) between the best-fitting \(f \sigma_8\) model and a zero-correlation model separately for the \(\psi_1\), \(\psi_2\) and \(\psi_3\) functions, for the ensemble of Abacus mocks and the DESI BGS and Pantheon+ datasets. The results are shown in Fig. 8, with the distribution of mock results displayed as the histograms, and the DESI-Pantheon+ measurement shown as the vertical lines. The differences \(\Delta \chi^2 \sim 4\) for the DESI-Pantheon+ measurements represent a \(\sim 2\)-\(\sigma\) detection for each individual correlation, which accumulates across the different statistics. The results for the ensemble of mocks show that a wide range of detection significances may be obtained. Our results for the real data are consistent with these distributions, resulting in a somewhat more significant detection than the mock average for the auto-correlations \(\psi_1\) and \(\psi_2\), and a somewhat less significant detection for the cross-correlation \(\psi_3\).

This is the first detection of the correlation between the fluctuations in SNe magnitudes and large-scale structure. Our results are consistent with Planck predictions, but due to the large uncertainties we cannot distinguish between \(\Lambda\)CDM and the dynamical dark energy model \(w_0 w_a\)CDM favoured by [105].

7.1 Forecasting improvements in \(f\sigma_8\) constraints with future SN samples↩︎

The statistical significance of these correlations, and the accuracy with which they can constrain the growth rate, will improve as galaxy and SN samples grow in size. To forecast these improvements and demonstrate the potential of this method, we performed Fisher matrix forecasts of the growth rate error for different SN sample sizes, otherwise maintaining the properties of the DESI BGS and Pantheon+ configurations. We consider the combined \((\psi_1, \psi_2, \psi_3)\) statistics varying the growth rate \(f\) and evaluate the Fisher matrix value as, \[\label{eq:fismat} F = \frac{\partial \mathbf{m}^T}{\partial f} \mathbf{C}^{-1} \frac{\partial \mathbf{m}}{\partial f},\tag{21}\] where \(\mathbf{m}\) is the correlation function model and \(\mathbf{C}\) is the covariance matrix. For the purpose of this test, we used Gaussian analytical covariances using the expressions in [69] which include the effects of sample variance, selection functions of the density and velocity tracers, shot noise and velocity measurement errors and curved-sky effects. [69] demonstrate that these expressions are accurate representations of the covariance across statistics and scales, if the selection function of the datasets does not significantly vary on the separation scale in question. We considered the \((\psi_1, \psi_2, \psi_3)\) statistics varying the growth rate \(f\), where the error forecast for \(f\) is given by \(F^{-1/2}\). We considered SN samples of size \(( 100, 300, 1000, 3000 )\), finding fractional errors in the best fit \(f \sigma_8\) to be \(( 49.8\%, 22.5\%, 7.7\%, 3.0\% )\) for these cases. We note that the forecast error is somewhat more optimistic than the result of our analysis of the real data. Two important reasons for this discrepancy are, first, that the forecast is based on an analytical covariance, and the data measurement uses a mock covariance, where the latter should be more reliable. Secondly, a Fisher matrix forecast provides the minimum possible error, and this could highlight that the weighting used in the analysis could be further optimised. Whilst a simplified forecast, the results lie within a factor of 2 of the growth rate measurement using the real sample, and the relative scaling with SNe number should be robust. This analysis demonstrates that future SN surveys have the potential to improve growth rate errors by a further order of magnitude (for more detailed forecasts, see [56] and [18]).

8 Conclusions↩︎

In this study we report a modest detection of the magnitude correlations within current SN samples, and between SNe and local large-scale structure, using the DESI Bright Galaxy Survey and Pantheon+ datasets. These correlations are associated with peculiar velocities induced by local large-scale structure, and may be used to measure the growth rate of structure. The best-fit normalised growth rate and galaxy bias from the correlation fit to the data catalogues are \(f \sigma_8 = 0.384^{+0.094}_{-0.157}\) and \(b \sigma_8 =0.998^{+0.035}_{-0.036}\). This measured growth rate value is consistent with expectations of the \(\Lambda\)CDM cosmological model. At the 1-\(\sigma\) level, our measurement also agrees with the DESI full-shape redshift-space distortion analysis of the BGS sample, which finds \(f \sigma_8 = 0.38 \pm 0.09\) [104].

We used realistic mock catalogues, drawn from the Abacus simulations, to validate our methodology and demonstrate that we are able to recover the fiducial growth rate within the statistical errors in the presence of realistic survey selection functions, galaxy bias and measurement noise. We also demonstrated that the SN selection effects have a negligible impact on the growth rate measurements at the level of precision of this analysis.

The forthcoming combination of the DESI BGS dataset and DESI Peculiar Velocity survey, which contains both Fundamental Plane and Tully-Fisher samples [57], [106], [107], will enable accurate tests of the growth of structure. The next generation of transient surveys such as the Zwicky Transient Facility [53], [54] and Rubin Observatory’s Legacy Survey of Space and Time [55], [108] will provide a 1-2 order-of-magnitude increase in the number of SNe detected, allowing us to achieve unprecedented precision in testing cosmological models using supernovae [18], [56], [109]. This expansion in sample size will necessitate extensions of multiple aspects of our analysis including non-linear modelling of the power spectrum, covariance modelling and the impact of SN selection effects.

Acknowledgements↩︎

We thank Anais Möller for useful discussions of supernovae selection effects, and Greg Aldering for helpful comments on a draft of this paper.

This research was conducted by the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (project number CE230100016) and funded by the Australian Government. AN would like to acknowledge the financial support received through the award of a Research Training Program Stipend scholarship by Swinburne University.

The project leading to this publication has received funding from Excellence Initiative of Aix-Marseille University - A*MIDEX, a French “Investissements d’Avenir” program (AMX-20-CE-02 - DARKUNI).

This material is based upon work supported by the U.S. Department of Energy (DOE), Office of Science, Office of High-Energy Physics, under Contract No. DE–AC02–05CH11231, and by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract. Additional support for DESI was provided by the U.S. National Science Foundation (NSF), Division of Astronomical Sciences under Contract No. AST-0950945 to the NSF’s National Optical-Infrared Astronomy Research Laboratory; the Science and Technology Facilities Council of the United Kingdom; the Gordon and Betty Moore Foundation; the Heising-Simons Foundation; the French Alternative Energies and Atomic Energy Commission (CEA); the National Council of Humanities, Science and Technology of Mexico (CONAHCYT); the Ministry of Science, Innovation and Universities of Spain (MICIU/AEI/10.13039/501100011033), and by the DESI Member Institutions: https://www.desi.lbl.gov/collaborating-institutions. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the U. S. National Science Foundation, the U. S. Department of Energy, or any of the listed funding agencies.

The authors are honored to be permitted to conduct scientific research on I’oligam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation.

Data availability↩︎

Data points for all the figures are available at https://doi.org/10.5281/zenodo.17111172.

Affiliations↩︎

\(^{1}\) Centre for Astrophysics & Supercomputing, Swinburne University of Technology, P.O. Box 218, Hawthorn, VIC 3122, Australia
\(^{2}\) OzGrav: The ARC Centre of Excellence for Gravitational Wave Discovery
\(^{3}\) Aix Marseille Univ, CNRS/IN2P3, CPPM, Marseille, France
\(^{4}\) Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA
\(^{5}\) Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, MA 02215 USA
\(^{6}\) Department of Physics & Astronomy, University of Rochester, 206 Bausch and Lomb Hall, P.O. Box 270171, Rochester, NY 14627-0171, USA
\(^{7}\) Dipartimento di Fisica “Aldo Pontremoli”, Università degli Studi di Milano, Via Celoria 16, I-20133 Milano, Italy
\(^{8}\) INAF-Osservatorio Astronomico di Brera, Via Brera 28, 20122 Milano, Italy
\(^{9}\) Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK
\(^{10}\) Korea Astronomy and Space Science Institute, 776, Daedeokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea
\(^{11}\) Instituto de Fı́sica, Universidad Nacional Autónoma de México, Circuito de la Investigación Cientı́fica, Ciudad Universitaria, Cd. de México C. P. 04510, México
\(^{12}\) Department of Astronomy & Astrophysics, University of Toronto, Toronto, ON M5S 3H4, Canada
\(^{13}\) Department of Physics & Astronomy and Pittsburgh Particle Physics, Astrophysics, and Cosmology Center (PITT PACC), University of Pittsburgh, 3941 O’Hara Street, Pittsburgh, PA 15260, USA
\(^{14}\) University of California, Berkeley, 110 Sproul Hall #5800 Berkeley, CA 94720, USA
\(^{15}\) Departamento de Física, Universidad de los Andes, Cra. 1 No. 18A-10, Edificio Ip, CP 111711, Bogotá, Colombia
\(^{16}\) Observatorio Astronómico, Universidad de los Andes, Cra. 1 No. 18A-10, Edificio H, CP 111711 Bogotá, Colombia
\(^{17}\) Institut d’Estudis Espacials de Catalunya (IEEC), c/ Esteve Terradas 1, Edifici RDIT, Campus PMT-UPC, 08860 Castelldefels, Spain
\(^{18}\) Institute of Cosmology and Gravitation, University of Portsmouth, Dennis Sciama Building, Portsmouth, PO1 3FX, UK
\(^{19}\) Institute of Space Sciences, ICE-CSIC, Campus UAB, Carrer de Can Magrans s/n, 08913 Bellaterra, Barcelona, Spain
\(^{20}\) University of Virginia, Department of Astronomy, Charlottesville, VA 22904, USA
\(^{21}\) Fermi National Accelerator Laboratory, PO Box 500, Batavia, IL 60510, USA
\(^{22}\) Center for Cosmology and AstroParticle Physics, The Ohio State University, 191 West Woodruff Avenue, Columbus, OH 43210, USA
\(^{23}\) Department of Physics, The Ohio State University, 191 West Woodruff Avenue, Columbus, OH 43210, USA
\(^{24}\) The Ohio State University, Columbus, 43210 OH, USA
\(^{25}\) School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia
\(^{26}\) Department of Physics, University of Michigan, 450 Church Street, Ann Arbor, MI 48109, USA
\(^{27}\) University of Michigan, 500 S. State Street, Ann Arbor, MI 48109, USA
\(^{28}\) Department of Physics, The University of Texas at Dallas, 800 W. Campbell Rd., Richardson, TX 75080, USA
\(^{29}\) NSF NOIRLab, 950 N. Cherry Ave., Tucson, AZ 85719, USA
\(^{30}\) Department of Physics, Southern Methodist University, 3215 Daniel Avenue, Dallas, TX 75275, USA
\(^{31}\) Sorbonne Université, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE), FR-75005 Paris, France
\(^{32}\) Department of Astronomy and Astrophysics, UCO/Lick Observatory, University of California, 1156 High Street, Santa Cruz, CA 95064, USA
\(^{33}\) Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95065, USA
\(^{34}\) Departament de Física, Serra Húnter, Universitat Autònoma de Barcelona, 08193 Bellaterra (Barcelona), Spain
\(^{35}\) Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Edifici Cn, Campus UAB, 08193, Bellaterra (Barcelona), Spain
\(^{36}\) Department of Astronomy, The Ohio State University, 4055 McPherson Laboratory, 140 W 18th Avenue, Columbus, OH 43210, USA
\(^{37}\) Institució Catalana de Recerca i Estudis Avançats, Passeig de Lluı́s Companys, 23, 08010 Barcelona, Spain
\(^{38}\) Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, U.K
\(^{39}\) IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France
\(^{40}\) Department of Physics and Astronomy, University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada
\(^{41}\) Perimeter Institute for Theoretical Physics, 31 Caroline St. North, Waterloo, ON N2L 2Y5, Canada
\(^{42}\) Waterloo Centre for Astrophysics, University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada
\(^{43}\) Space Sciences Laboratory, University of California, Berkeley, 7 Gauss Way, Berkeley, CA 94720, USA
\(^{44}\) Instituto de Astrofísica de Andalucía (CSIC), Glorieta de la Astronomía, s/n, E-18008 Granada, Spain
\(^{45}\) Department of Physics and Astronomy, Sejong University, 209 Neungdong-ro, Gwangjin-gu, Seoul 05006, Republic of Korea
\(^{46}\) CIEMAT, Avenida Complutense 40, E-28040 Madrid, Spain
\(^{47}\) National Astronomical Observatories, Chinese Academy of Sciences, A20 Datun Road, Chaoyang District, Beijing, 100101, P. R. China

References↩︎

---

1. \(^*\)E-mail: andnguyen@swin.edu.au↩︎

2. Accessed from https://github.com/PantheonPlusSH0ES/DataRelease/tree/main/Pantheon%2B_Data/4_DISTANCES_AND_COVAR.↩︎