Probing the Distance Duality Relation with Machine Learning and Recent Data

Abstract

The distance duality relation (DDR) relates two independent ways of measuring cosmological distances, namely the angular diameter distance and the luminosity distance. These can be measured with baryon acoustic oscillations (BAO) and Type Ia supernovae (SNe Ia), respectively. Here, we use recent DESI DR1, Pantheon+, SH0ES and DES-SN5YR data to test this fundamental relation. We employ a parametrised approach and also use model-independent Generic Algorithms (GA), which are a machine learning method where functions evolve loosely based on biological evolution. When we use DESI and Pantheon+ data without Cepheid calibration or big bang nucleosynthesis (BBN), there is a $2\sigma$ violation of the DDR in the parametrised approach. Then, we add high-redshift BBN data and the low-redshift SH0ES Cepheid calibration. This reflects the Hubble tension since both data sets are in tension in the standard cosmological model $\Lambda$CDM. In this case, we find a significant violation of the DDR in the parametrised case at $6\sigma$. Replacing the Pantheon+ SNe Ia data by DES-SN5YR, we find similar results. For the model-independent approach, we find no deviation in the uncalibrated case and a small deviation with BBN and Cepheids which remains at 1$\sigma$. This shows the importance of considering model-independent approaches for the DDR.

1 Introduction

The Etherington [1] or distance duality relation (DDR) relates the angular diameter distance with the luminosity distance. It holds in metric theories of gravity that have conservation of the photon number and where photons travel along unique null geodesics. Thus, testing this relation can rule out extensions to the standard cosmological model $\Lambda$CDM or it can be a strong indicator for new physics. Deviations from the DDR can be expressed with a function that quantifies this discrepancy. The DDR has been extensively tested and has gained renewed interest recently [2, 3, 4].

The DDR could be violated if there are deviations from a metric theory of gravity, if photons do not travel along null geodesics or if fundamental constants vary in time. The violation photon number conservation would break the DDR and could be caused, i.e. by attenuation through interstellar gas, dust or plasma [5]. Measuring the DDR can also provide limits on axion-like particles which change photon number as discussed in [6, 7]. These are scalars beyond the standard model of particle physics that couple to photons Thus, in presence of external magnetic fields, photons could convert into axion-like particles[8]. The necessary conditions can possibly be provided by intergalactic magnetic fields.

Type Ia Supernovae (SNe Ia) are a measurement for the luminosity distance and can then be combined with other measurements related to the angular diameter distance. This has been done with measurements of the Hubble expansion rate from luminous red galaxies to constrain cosmic opacity and new physics [6]. Galaxy clusters can also be used for a measurement on the angular diameter distance to compare it with SNe Ia data [9]. A common measurement of the angular diameter distance are baryon acoustic oscillations (BAO). In [10], Pantheon+ SNe Ia [11] and a compilation of BAO data by [12] were used to test different parametrisations of DDR violations.

Other than BAO, the DDR can also be tested using SNe Ia combined with strong gravitational lensing. No statistically significant violation has been found with the Pantheon+ data set [13, 14]. This has been confirmed using compact radio quasars for the angular diameter distance, also in combination with Pantheon+ [15].

It is important to test the DDR in a model-independent way that does not assume a certain parametrisation that is not necessarily physically motivated. The authors of [16] have conducted a model-independent test of the DDR utilising Pantheon SNe Ia [17] and eBOSS DR16 quasar data. BAO have been used in [18] with Monte Carlo methods to constrain the DDR violation agnostically at different redshifts. Strong gravitational lensing with SNe Ia has also been tested model-independently, yielding no deviation [4]. However, in [2], they use data from the Sunyaev–Zeldovich effect for galaxy clusters finding a $2$–$3$ $\sigma$ discrepancy depending on the cluster data.

Here we use genetic algorithms (GA), which are machine learning methods. GA have already been used in other cosmological contexts, e.g. as a null test for the cosmological constant [19], for the dark energy equation of state [20, 21, 22] or other consistency tests including the Om statistic and the $r_{0}$ test for dark sector interactions [23]. They also have been employed to reconstruct the Hubble expansion history [24] and to test curvature [25] or the coupling of dark energy to the electromagnetic sector [26].

The DDR has already been tested with GA using older data sets, e.g. in [21]. The authors in [27] used GA with Pantheon SNe Ia and various BAO surveys (6dFGS [28], SDSS [29], BOSS CMASS [30], WiggleZ [31], MGS [32], BOSS DR12 [33], DES [34], Ly-$\alpha$ [35], SDSS DR14 LRG [36] and quasar observations [37]). So far, no studies have found significant violation.

Other models outside the DDR have been tested with DESI DR1 [38] and Pantheon+ [11] and alternatively DES-SN5YR [39]. In [40], the authors reconstruct the expansion history with the Om statistic and the evolution of the total equation of state. They find a $2\sigma$ discrepancy for Pantheon+ and an over $3\sigma$ discrepancy with DES-SN5YR data. An extension of the Friedman–Lemaître–Robertson–Walker (FLRW) metric has been tested on this data, finding a discrepancy to $\Lambda$CDM greater than $1\sigma$ [41].

This work provides a model-independent approach for the DDR in addition to a frequentist approach using new data sets, namely Pantheon+, SH0ES [11], DES-SN5YR [39] and DESI DR1 [38]. The frequentist method minimises the $\chi^{2}$ statistic numerically in every point over all other parameters whose contours are not considered. In a next step, we use GA as our model-independent approach. This is a machine learning approach that tries to find the best-fitting function, here based also on the $\chi^{2}$ statistic.

In addition to the uncalibrated SNe Ia and BAO analysis, we calibrate our data sets in such a way that they reflect the Hubble tension which exists in $\Lambda$CDM. For a summary of the Hubble tension, see e.g. [42]. To put the DDR into the context of this tension, we calibrate the SNe Ia with SH0ES yielding $H_{0}=73.5\pm 1.1\,$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{Mpc}}^{-1}$$ [43] and the DESI BAO with the big bang nucleosynthesis (BBN) value for $\omega_{\mathrm{b}}$ resulting in $H_{0}=68.53\pm 0.80\,$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{Mpc}}^{-1}$$ for the parametrised approach.

In Sect. 2, we introduce and derive the DDR. Sect. 3 describes the SNe Ia and BAO data sets we employ and details the Cepheid calibration. The parametrisation with the $\epsilon$ parameter is presented in Sect. 4 and the model-independent GA approach in Sect. 5. We show our results for both approaches in Sect. 6 before we summarise and conclude in Sect. 7.

2 Distance Duality Relation

The luminosity distance $d_{\rm{L}}$, which can be measured by SNe Ia, compares flux $F$ and absolute luminosity L of an object.

$$d_{\rm{L}}=\sqrt{\frac{\mathrm{L}}{4\,\pi\,F}}\>.$$ (2.1)

If we assume a flat universe, we can express the flux with the comoving distance $d_{\mathrm{M}}$

$$F=\frac{L}{4\,\pi\,a_{0}^{2}\,d_{\mathrm{M}}^{2}\,(1+z)^{2}}\>.$$ (2.2)

Here, $a_{0}$ is the scale factor of the universe today. The factor $(1+z)^{2}$ comes from the fact that the photon energy is proportional to $\frac{1}{a(t)}\propto(1+z)$ coming from the dilation of the wavelength. Additionally, the rate at which photons arrive, which is measured by the luminosity, is also diluted by $\frac{1}{a(t)}$. Using this, we obtain the luminosity distance as a function of the comoving distance

$$d_{\rm{L}}=a_{0}\,(1+z)\,d_{\mathrm{M}}\>.$$ (2.3)

On the other hand, the angular diameter distance can be obtained using trigonometry by comparing the proper length $l$ and the angular size $\theta$ of an object. To do this, we have to know the proper length of the object which is the case for BAO, making it a standard ruler.

$$d_{\rm{A}}=\frac{l}{\theta}\>.$$ (2.4)

In the FLRW metric, a proper length is equal to the product of the scale factor, the comoving distance and the angular size, yielding $l=a(t)\,d_{\mathrm{M}}\,\theta$ [44]. Due to the expansion of the universe, distances are redshifted by the factor $1+z=\lambda(t_{0})/\lambda(t)=a_{0}/a(t)$ where $\lambda$ is the physical size. This leads us to:

$$\frac{l}{\theta}=a(t)\,d_{\mathrm{M}}=\frac{a_{0}}{(1+z)}\,d_{\mathrm{M}}\>.$$ (2.5)

The DDR now follows from comparing both distance measures

$$d_{\rm{A}}=\frac{d_{\rm{L}}}{(1+z)^{2}}\>.$$ (2.6)

To test this empirically, we parametrise deviations from the DDR by defining the following phenomenological function

$$\eta(z)=\frac{d_{\rm{L}}(z)}{d_{\rm{A}}(z)\,(1+z)^{2}}\>.$$ (2.7)

Thus, if $\eta(z)=1$, there is no discrepancy, which serves as the null hypothesis.

3 Data

3.1 Supernovae Type Ia

The Pantheon+ SNe Ia data set [11] contains 1701 Type Ia supernovae. In the data set, they provide the Hubble Diagram Redshift with which we compute the luminosity distance. This is then rescaled using the heliocentric redshift for the likelihood. For the magnitude, we use the Tripp 1998 [45] corrected magnitude. The measured apparent magnitude is then compared with the theoretical one calculated from the luminosity distance:

$$m_{\rm{th}}(z,\Omega_{\mathrm{m}},H_{0})=M_{0}+5\,\log_{10}\left[\frac{d_{\rm{L}}(z,\,\Omega_{\rm m},\,H_{0})}{\mathrm{Mpc}}\right]+25\;.$$ (3.1)

Thus, our $\chi^{2}$ statistic depends on the parameters $\Omega_{\mathrm{m}}$, $H_{0}$ and $M_{0}$. We calculate it with the Fisher matrix $F_{ij}$ which is the inverse of the covariance matrix [46].

$$\chi^{2}_{\rm{Panth}}=(m_{\rm{th}}-m_{\rm{data}})_{i}F_{ij}(m_{\rm{th}}-m_{\rm{data}})_{j}\,.$$ (3.2)

When using the Pantheon+ data set without SH0ES, a redshift cut ($z>0.01$) has to be made due to systematic errors in the data [43]. Furthermore, in this uncalibrated case, we marginalise over the absolute magnitude $M_{0}$ with the expression taken from Appendix C of [47].

3.1.1 Cepheid Calibration

When using the Pantheon+ SNe Ia catalogue in conjunction with SH0ES, it is possible to calibrate all SNe Ia with Cepheids. In the cases where a Cepheid has been found in the same galaxy as the SNe Ia, the apparent magnitude will be calculated with the absolute distance from the Cepheid:

$$m(z)=M_{0}+d_{\mathrm{abs},\mathrm{SH0ES}}.$$ (3.3)

Every SN with the calibration flag has this associated apparent magnitude instead of the one using the luminosity distance. This calibration effectively fixes the absolute magnitude $M_{0}$ and thus sets the Hubble parameter $H_{0}$ to the combined Pantheon+ and SH0ES value of $73.5\pm 1.1\,$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{Mpc}}^{-1}$$ [43].

3.2 Baryon Acoustic Oscillations

In the early universe, before the decoupling of baryons from photons, acoustic waves were propagating in the tightly coupled fluid with the same phase. After decoupling, at the drag epoch, these waves were frozen-in and are the basis for matter fluctuations which grow with the expansion of the universe. This depends on the sound horizon $r_{\rm s}$ at the drag redshift $z_{\rm d}$ which length is theoretically well-predicted using early-time physics. It is obtained by integrating over the sounds speed $c_{\mathrm{s}}$ from the drag redshift shortly after decoupling up to infinity. We write it in this form to show that it is proportional to $1/H_{0}$ which will be important for our analysis.

$$r_{\rm d}=r_{\rm s}(z_{\rm d})=\frac{1}{H_{0}}\int_{z_{\rm d}}^{\infty}\frac{c_{\mathrm{s}}(z)}{H(z)/H_{0}}\mathrm{d}z\;.$$ (3.4)

We use the DESI DR1 data release [38] where they provide seven redshift bins from $0.1<z<4.2$ using the bright galaxy survey, luminous red galaxies, emission line galaxies, quasars and Lyman-$\alpha$ data. Since they measure angular separation $\Delta\theta$, they provide a ratio of a distance divided by the sound horizon at drag epoch $r_{\rm d}$. They use the approximation by [48] which we adopt as well. It assumes $\Lambda$CDM and standard early-time physics to estimate the $r_{\rm d}$ given the physical matter and baryon densities $\omega_{\mathrm{m}}$ and $\omega_{\mathrm{b}}$ and the effective number of neutrinos $N_{\rm{eff}}$:

$$r_{\rm d}=147.05\,\left(\frac{\omega_{\mathrm{m}}}{0.1432}\right)^{-0.23}\,\left(\frac{N_{\rm{eff}}}{3.04}\right)^{-0.1}\,\left(\frac{\omega_{\mathrm{b}}}{0.02236}\right)^{-0.13}\,\mathrm{Mpc}\;.$$ (3.5)

The distance in each bin is either the comoving distance $d_{\rm M}(z)$ and the Hubble distance $d_{\rm H}(z)$ or the angle-averaged distance $d_{\rm V}(z)$ in bins with low signal-to-noise ratio:

$$d_{\rm V}(z)=\left[z\,d_{\rm M}(z)^{2}\,d_{\rm H}(z)\right]^{1/3}\;.$$ (3.6)

In the $\Lambda$CDM case, the comoving distance is the angular diameter distance multiplied by $1+z$:

$$d_{\rm M}=d_{\rm A}\,(1+z)\;.$$ (3.7)

From Eq. 3.4, we know that $r_{\rm d}\propto H_{0}^{-1}$. The same is true for the above distances $d_{X}$. This means that the ratio $d_{X}/r_{\rm d}$ is independent of the value of $H_{0}$. In practice, we need to compute the function $d_{\rm A}(z)$ from which we derive $d_{\rm M}(z)$ and $d_{\rm H}(z)$ and $d_{\rm V}(z)$. In the next step, the ratio is computed with $r_{\rm d}$. The $\chi^{2}$ is computed analogously to the SNIa case, with its respective Fisher matrix $\tilde{F}$.

$$\chi^{2}_{\rm{DESI}}=(v_{\rm{th}}-v_{\rm{data}})_{i}\,\tilde{F}_{ij}\,(v_{\rm{th}}-v_{\rm{data}})_{j}\;,$$ (3.8)

where $v_{i}=d_{X,i}/r_{\rm d}$. In the standard case without BBN, we minimise over the value of $r_{\rm d}$ to compute the minimum $\chi^{2}$.

4 Parametrisation with the $\epsilon$-Parameter

Since the deviations grow larger at higher redshifts, we parameterise $\eta(z)$ with a constant exponential parameter $\epsilon$ which is common in the literature [6, 49, 10]:

$$\eta(z)=(1+z)^{\epsilon}.$$ (4.1)

We test this parametrisation of $\eta(z)$ using the frequentist approach. In every calculated point, the $\chi^{2}$ statistic is calculated using the data sets, minimising over all other parameters. For the BAO data, we need to know the evolution of $H(z)$. For the GA, this is computed from the angular diameter distance $d_{\rm A}=d_{\rm L}/[(1+z)^{2}\eta(z)]$ where $d_{\rm L}$ and $\eta(z)$ are the resulting GA functions.

$$d_{\rm L}(z,\Omega_{\mathrm{m}})=a_{0}\,(1+z)\,\eta(z)\,\int_{0}^{z}\>\frac{c}{a_{0}\,H(z^{\prime},\Omega_{\mathrm{m}})}\,\mathrm{d}z^{\prime}\>.$$ (4.2)

Now, we can derive both sides to find $H(z)/H_{0}$.

$$\frac{H(z,\Omega_{\mathrm{m}})}{H_{0}}=\left[\frac{\mathrm{d}}{\mathrm{d}z}\left(\frac{d_{\rm L}(z,\Omega_{\mathrm{m}})}{c/H_{0}(1+z)\eta(z,\epsilon)}\,\right)\right]^{-1}\;.$$ (4.3)

From these, we can derive the other distances needed for the BAO. The $\chi^{2}$ is independent of $H_{0}$ since it cancels out as discussed in 3.2.

Due to the $\epsilon$-parametrisation, only one data set is rescaled with $\epsilon$. We choose BAO in this case, so the SNe Ia likelihood does not depend on $\epsilon$ here. To calculate the $\chi^{2}$ for the BAO, we fit $d_{\rm A}(\Omega_{\mathrm{m}},\epsilon)$ to the data using the modified DDR:

$$d_{\rm A}(z,\Omega_{\mathrm{m}},\epsilon)=d_{\rm L}(z,\Omega_{\mathrm{m}})\cdot(1+z)^{-2-\epsilon}.$$ (4.4)

With this, we can compute a profile likelihood for $\Omega_{\mathrm{m}}$ and $\epsilon$.

4.1 BBN and Planck values

The DESI data by itself cannot constrain $H_{0}$, as detailed in Sect. 3.2. However, we can use Eq. (3.5) to compute the sound horizon at the drag redshift. This depends on the value of the baryon density $\omega_{\mathrm{b}}=\Omega_{\mathrm{b}}\,h^{2}$, which can be computed from the BBN. Using data about primordial Deuterium and Helium abundances yields $\omega_{\mathrm{b}}=0.02218\pm 0.00055$ [50]. When this is added as a fixed parameter to the BAO data from DESI, it fixes $r_{\rm d}$, meaning we can compute $H_{0}$. The uncertainty on $\omega_{\mathrm{b}}$ is very small, therefore we assume that there is no significant correlation with $\Omega_{\mathrm{m}}$ or $\epsilon$. As can be seen in the constraints from the DESI DR1 data release, the information from BBN in practice does not change the value or the error of $\Omega_{\mathrm{m}}$, as can be seen in [38] in Table 3.

4.2 SH0ES magnitude for DES-SN5YR

The DES-SN5YR data are not directly coupled to a Cepheid data set. So, to incorporate the SH0ES value, we calculate the absolute magnitude $M_{0}$ for the minimum $\chi^{2}$ fixing $\Omega_{\mathrm{m}}=\Omega_{\mathrm{m},\mathrm{DES}}$ [39] and $H_{0}=H_{0,\mathrm{SH0ES}}$ [43]. We repeat this for adding and subtracting the $H_{0}$ error from Pantheon+SH0ES to construct an additional $\chi^{2}$ for the value of $M_{0}$.

$$\chi^{2}_{\mathrm{SH0ES}}(M_{0})=\left(\frac{M_{0}-M_{0,\mathrm{SH0ES}}}{\sigma_{M_{0},\mathrm{SH0ES}}}\right)^{2}\;.$$ (4.5)

This essentially calibrates the SNe Ia with the SH0ES absolute magnitude which in turn determines $H_{0}$.

5 Genetic Algorithms

GA are a Machine Learning method that finds functions to fit a data set without previous parametrisation. Here, we will briefly outline how they function. For a more detailed introduction, we refer the reader to Ref. [20].

To find the best-fitting function, GA have a random initial population of functions according to a specified grammar. We use polynomials here because they are differentiable and sufficient to describe our data. From the initial population, in each generation the functions are evaluated for how well they match the data using a $\chi^{2}$-function which is the sum of the SNe Ia and BAO $\chi^{2}$ statistic:

$$\chi^{2}_{\mathrm{tot}}=\chi^{2}_{\mathrm{DESI}}+\chi^{2}_{\mathrm{Pantheon+SH0ES}}.$$ (5.1)

(a) Crossover operation

(b) Mutation operation

Figure 1: Illustration of GA operations, figure taken from [51].

The functions with the lowest $\chi^{2}$ can generate offspring functions, and the rest do not propagate to the next generation. The remaining functions "reproduce", with a preselected probability called the crossover rate, see Fig. 1a. Here, both parent functions combine parts of them, which represents the DNA in the analogy. These parts are then summed. Analogously, the mutation operation (see Fig. 1b) is selected with a probability called the mutation rate. If that is the case, one parameter of the function will randomly change at the time of "reproduction" to the new generation. For each of the 10 samples, the function with the lowest $\chi^{2}$ is the final candidate of the sample. From these 10, we select the one with the lowest $\chi^{2}$ as the resulting best-fit function.

We compare the resulting functions of the GA with $\Lambda$CDM for each case. When we calculate $H(z)$ in the $\Lambda$CDM case from the Friedmann equation, we also consider radiation and neutrinos.

$$\frac{H^{2}(z)}{H_{0}^{2}}=a^{-3}\Omega_{\mathrm{m}}\left(1+a_{\rm{eq}}/a\right)+\left[1-\Omega_{\mathrm{m}}\left(1+a_{\rm{eq}}\right)\right]f(a,w_{0},w_{\rm a})$$ (5.2)

Here, $f$ is a function that depends on $w_{0}$ and $w_{a}$ from the Chevallier–Polarski–Linder parametrisation of dark energy, i.e. $w_{\rm{DE}}(a)=w_{0}+w_{a}(1-a)$ [52, 53]. This factor depends on the ordinary hypergeometric function ${}_{3}F_{2}(a_{i};b_{j};x)$ with argument $x$

$$f(a,w_{0},w_{a})=a^{-3\,(1+w_{0}+w_{a})}\exp[-3\,w_{a}+3\,a\,w_{a3}\,F_{2}(1,1,0;2,2;a)].$$ (5.3)

The scale factor at matter–radiation equality $a_{\rm{eq}}$ depends on the density of radiation $\Omega_{\gamma}$ and ultra-relativistic neutrinos $\Omega_{\nu}$:

$$a_{\rm{eq}}=\frac{\Omega_{\gamma}+\Omega_{\nu}}{\Omega_{\mathrm{m}}}.$$ (5.4)

For the effective number of neutrinos, we use $N_{\rm{eff}}=3.044$ [54].

The GA tries to fit two functions, namely the luminosity distance $d_{\rm L}$ and the DDR deviation $\eta(z)$, from which the angular distance is calculated. In the first three cases, we either marginalise over the absolute SNe Ia magnitude $M_{0}$ or minimise the $\chi^{2}$ by varying $r_{\rm d}$. This implies that the GA can rescale the distance functions by any value without changing the $\chi^{2}$. Any global factor is reabsorbed in the marginalisation or the minimisation, respectively. Thus, we choose the following dimensionless parametrisation for the two GA functions $f_{\mathrm{GA}}$ and $g_{\mathrm{GA}}$:

$$\eta(z)=(1+z)^{\epsilon(z)}=(1+z)^{1/10\,f_{\mathrm{GA}}(z)},$$ (5.5)

$$\frac{H_{0}}{c}\,d_{\rm L}(z)=z[1+z\,g^{2}_{\mathrm{GA}}(z)].$$ (5.6)

Using these results, we can calculate $d_{\rm A}$ using Eq. 2.7 and $H(z)$ using Eq. 4.3.

Figure 2: $\chi^{2}$ of the GA functions per generation for the standard DESI+Pantheon+ case.

In Fig. 2, we show the evolution of the GA over 300 generations. All 10 populations approach the $\chi^{2}$ of $\Lambda$CDM which is at $\chi^{2}=1785$ where we minimised over $H_{0}$ and $\Omega_{\mathrm{m}}$. This is represented by the dashed line. In this figure, we used the standard Pantheon+ case for which 5 of the 10 samples have a lower $\chi^{2}$ than $\Lambda$CDM at the end.

For the model-independent GA, it is necessary to use the unparametrised sound horizon at the drag redshift $r_{\rm d}$. We cannot add only the BBN value since this would imply other free parameters. Thus, we calculate this quantity not only with BBN but with additional Planck 2018 data. We fix the sound horizon to a value using the Planck values [55] $\omega_{\rm{m}}$ and $H_{0}$ in addition to the BBN.

For the case using Cepheid calibration and the sound horizon with BBN andPlanck values, we do not predict the distances in a dimensionless way as before, since in this case it is possible for the GA to predict both the amplitude of $d_{\rm L}$ and $d_{\rm A}$. This emphasises the violation of the DDR due to the Hubble tension, since we predict $d_{\rm L}$ instead of $H_{0}\cdot d_{\rm L}$. However, it also increases the computing time for the GA since machine learning algorithms generally work better with values closer to unity [56]. So, even after rescaling our functions, the GA needs 1000 generations in this case to converge compared to 300 in the standard case.

5.1 Error Calculation

To estimate the GA error, we generate multiple samples using different random seeds. This affects the initial population of the GA. From these samples of different resulting functions $\eta(z)$, we select the 68.27% of functions closest to our best-fit solution. This would correspond to the 1$\sigma$ error in the Gaussian case, but we find that our results are not symmetric in every case. This is taken from 100 GA samples, calculated in redshift intervals of $0.1$ and then interpolated to smooth it out.

6 Results

6.1 Parametrised Sampling

6.1.1 Pantheon+ with DESI DR1

(a) Contour plot for Pantheon+ (yellow), DESI (blue) and their combination (red).

(b) Contour plot for Pantheon+ (yellow), DESI+BBN (blue) and their combination (red).

(c) Contour plot for Pantheon+SH0ES (yellow), DESI (blue) and their combination (red).

(d) Contour plot for Pantheon+SH0ES (yellow), DESI+BBN (blue) and their combination (red).

Figure 3: Evaluation of the profile likelihood of $\Omega_{\mathrm{m}}$ and $\epsilon$ with Pantheon+ and DESI DR1 data for different cases.

We sample the parameters $\epsilon$ (defined in Eq. 4.1) and $\Omega_{\mathrm{m}}$ using a profile likelihood approach and sampling a grid since this does not require much computational time. The profile likelihood contour plot for the standard case is shown in Fig. 3a. The best-fit point is $\epsilon=-0.019$ with a distance from $\epsilon=0$ of $1.54\sigma$.

Pantheon+ finds a value of $\Omega_{\mathrm{m}}=0.334\pm 0.018$ [43] and DESI DR1 computes a lower value of $\Omega_{\mathrm{m}}=0.295\pm 0.0015$. Since this constitutes a discrepancy of $2.2\sigma$, it is not surprising that we also find an $\epsilon$ different form $0$ using these data sets.

In a next step, we add data sets to our likelihoods that are in tension with respect to the Hubble parameter $H_{0}$ to reflect this tension in a violation of the DDR. As detailed in Sect. 4.1, for the second case, we add BBN information to the BAO data, which makes the likelihood sensitive to $H_{0}$. The result is plotted in Fig. 3b. This is extremely similar to Fig. 3a, since the marginalisation over the absolute SNe Ia magnitude $M_{0}$ can compensate for a value of $H_{0}$ that would not match the SNe Ia data in the calibrated case.

For the third case, we calibrate the Pantheon+ SNe Ia using Cepheids with the SH0ES data set, as was also done in the Pantheon+ analysis [43] but we do not use BBN information. This is shown in Fig. 3c. Here as well, there is not much visible change, compared to Fig. 3a, because the BAO likelihood is minimised over $r_{\rm d}\,h$, essentially leaving $H_{0}$ adjustable to the SH0ES value.

Lastly, we combine the last two cases and use both the Cepheid calibration and the BBN value of $\omega_{\mathrm{b}}$. Now, both likelihoods are sensitive to $H_{0}$ but prefer a different value. This, in turn, creates a violation of the DDR which can be seen in Fig. 3d. For the combined contour, the best-fit is at $\epsilon=-0.040$ and $6.2\sigma$ away from 0. DESI+BBN finds $H_{0}=68.52\pm 0.80\;\rm{km}\;\rm{s}^{-1}\;\rm{Mpc}^{-1}$, whereas Pantheon+SH0ES lies at $H_{0}=73.5\pm 1.1\;\rm{km}\;\rm{s}^{-1}\;\rm{Mpc}^{-1}$. This constitutes a $4.5\sigma$ difference. Together with the discrepancy in $\Omega_{\mathrm{m}}$, it makes sense that we find a value of $\epsilon$ that is far from $\epsilon=0$. Here, we find a $\chi^{2}=1569$ in $\Lambda$CDM and a much improved value of $\chi^{2}=1541$ when allowing for $\epsilon\neq 0$.

Exclusion of the $d_{\rm H}$ BAO data

In the DESI BAO data release, per redshift bin, they report either both the comoving distance $d_{\rm M}$ and the Hubble distance $d_{\rm H}$ or only the volume-averaged distance $d_{\rm V}$. To test for example the conservation of photons, it is useful to incorporate all distance measurements into the analysis. This is common in the literature, see e.g. [27, 2, 10]. To be more conservative, it is also possible to only use the data points of the comoving distance $d_{\rm M}=(1+z)d_{\rm A}$ since this is the transverse distance. We recalculate the earlier case with BBN and Cepheid calibration (see Fig. 3d) with only these values to show this difference in Fig. 4. The BAO contours are much larger in this case because only 5 instead of 12 data points from DESI are used. This also increases the combined contour. However, the preferred value is $\epsilon=-0.065$, which is lower than in Fig. 3d and still $4.4\sigma$ apart from the DDR. We note that to focus on this aspect, it would help to add other BAO data sets to DESI. In the rest of this work, we will use all 12 DESI data points for the GA and refer this more conservative test to future work.

Figure 4: Contour plot of $\Omega_{\mathrm{m}}$ and $\epsilon$ for Pantheon+SH0ES (yellow), DESI+BBN (blue) and their combination (red) using only the comoving distance measurements $d_{\rm M}$ from DESI.

6.1.2 DES-SN5YR with DESI DR1

We also calculate the profile likelihood showing contour plots for $\Omega_{\mathrm{m}}$ and $\epsilon$, replacing the SNe Ia with the DES-SN5YR data set [39]. The results are plotted in Fig. 5. All four contours are similar to those for Pantheon+ in Fig. 3. For the standard DES-SN5YR in Fig. 5a, the best-fit value for $\epsilon$ is $2.5\sigma$ apart from $0$ compared to $2.2\sigma$ for Pantheon+. This is due to DES-SN5YR preferring a higher value of $\Omega_{\mathrm{m}}$, making it less compatible with DESI. When adding either only BBN (Fig. 5b) or only the $\chi^{2}$ from Eq. 4.5 incorporating the $H_{0}$ value from Pantheon+SH0ES (Fig. 5c), the contours barely change with respect to the first case). In the last case, combining BBN and the SH0ES $\chi^{2}$, we find a strong violation of the DDR at $6.5\sigma$ compared to $6.2\sigma$ from Pantheon+. When comparing the $\chi^{2}$ values, we find a $\chi^{2}=1693$ in $\Lambda$CDM and a much smaller value of $\chi^{2}=1657$ when allowing for $\epsilon\neq 0$.

(a) Contour plot for DES-SN5YR (green), DESI (blue) and their combination (brown).

(b) Contour plot for DES-SN5YR (green), DESI+BBN (blue) and their combination (brown).

(c) Contour plot for DES-SN5YR+SH0ES (green), DESI (blue) and their combination (brown).

(d) Contour plot for DES-SN5YR+SH0ES (green), DESI+BBN (blue) and their combination (brown).

Figure 5: Evaluation of the profile likelihood of $\Omega_{\mathrm{m}}$ and $\epsilon$ with DES-SN5YR and DESI DR1 data for different cases.

6.2 GA Best-Fit Functions

We test five different cases with the GA, four using the Pantheon+ data set combined with SH0ES corresponding to the parametrised approach from Fig. 3 and one using the DES-SN5YR data set for the SN corresponding to Fig. 5a. In all cases, we use DESI data for the BAO.

6.2.1 Pantheon+ with DESI DR1

(a) GA Deviation from $\Lambda$CDM in percent for the angular diameter distance (red), the luminosity distance (green) and the Hubble parameter (purple).

(b) Deviation from the standard DDR $\eta(z)$ for the GA (in solid red) with the associated errors (shaded in orange) compared to $\Lambda$CDM (in dashed black), i.e. no deviation.

Figure 6: Results of the GA using Pantheon+ and DESI DR1 data.

When applying the GA to the uncalibrated SNe Ia and BAO data, we obtain the results of the GA for the angular diameter distance, the luminosity distances, and the Hubble parameter. We show the difference of each of these functions with respect to the best-fit $\Lambda$CDM in percent in Fig. 6a. For the angular diameter distance, no discrepancy to $\Lambda$CDM can be seen. The luminosity distance deviates slightly below $\Lambda$CDM as is visible in the DDR deviation as well. The Hubble parameter also shows no deviation from the standard model. The GA result for the DDR violation $\eta(z)=d_{\rm L}(z)/[d_{\rm A}(z)(1+z)^{2}]$ is shown in Fig. 6b where $\eta(z)=1$ corresponds to no violation. This is consistent with the deviation of $\epsilon=-0.0192$ we find in the parametrised approach in Sect. 6.1. However, here in the unparametrised approach this change constitutes less than $0.5\sigma$ using our sampled error.

BBN and Planck Values

As described in Sect. 4.1, we show the deviation from the DDR using uncalibrated SNe Ia and setting $\omega_{\mathrm{b}}$ to BBN and $H_{0}$ and $\omega_{\mathrm{m}}$ to the Planck 2018 values in In Fig. 7. The residuals with respect to $\Lambda$CDM are slightly larger than in the standard case but still below 5%, see Fig. 7a. The evolution of $\eta(z)$ looks very similar to 6b from the standard analysis. This was expected since DESI DR1 is mostly in agreement with Planck 2018 and the absolute SN magnitude is marginalised over, which essentially leaves the amplitude of the luminosity distance as a free parameter. This can then be adjusted to fit together with the amplitude of the BAO data.

(a) GA Deviation from $\Lambda$CDM in percent for the angular diameter distance (red), the luminosity distance (green) and the Hubble parameter (purple).

(b) Deviation from the standard DDR $\eta(z)$ for the GA (in solid red) with the associated errors (shaded in orange) compared to $\Lambda$CDM (in dashed black), i.e. no deviation.

Figure 7: Results of the GA using Pantheon+ and DESI DR1 with Planck values.

Cepheid Calibration

(a) GA Deviation from $\Lambda$CDM in percent for the angular diameter distance (red), the luminosity distance (green) and the Hubble parameter (purple).

(b) Deviation from the standard DDR $\eta(z)$ for the GA (in solid red) with the associated errors (shaded in orange) compared to $\Lambda$CDM (in dashed black), i.e. no deviation.

Figure 8: Results of the GA using using Pantheon+SH0ES and DESI DR1.

When enabling the Cepheid calibration (see Sect. 3.1.1), the GA finds the results shown in Fig. 8. Here, the deviations from $\Lambda$CDM are larger than in Fig. 6 and Fig. 7 but they remain below 5% for the residuals in Fig. 8a and the error on $\eta(z)$ is smaller than 1$\sigma$.

Cepheid Calibration and Planck Values

Here, we test the GA for the case in which the SNe Ia data reflect the Pantheon+ with SH0ES $H_{0}$ value of $H_{0}=73.5\pm 1.1\,$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{Mpc}}^{-1}$$ and at the same time the BAO reflect the Planck 2018 [55] values of $h,\,\omega_{\mathrm{m}}$ and $\omega_{\mathrm{b}}$. Together with the DESI BAO this yields $H_{0}=68.53\pm 0.80\,$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{Mpc}}^{-1}$$.

(a) GA Deviation from $\Lambda$CDM in percent for the angular diameter distance (red), the luminosity distance (green) and the Hubble parameter (purple).

(b) Deviation from the standard DDR $\eta(z)$ for the GA (in solid red) with the associated errors (shaded in orange) compared to $\Lambda$CDM (in dashed black), i.e. no deviation.

Figure 9: Results of the GA using using Pantheon+SH0ES and DESI DR1 with Planck values.

Fig. 9a again shows the results of the GA for the angular diameter distance, the luminosity distance, and the Hubble parameter. The angular diameter distance found by the GA is much higher than in $\Lambda$CDM. In contrast to the standard case, the $H(z)$ predicted by the GA is much lower than in the standard model at higher redshifts. This is expected since at lower $z$, the SNe Ia data with the higher $H_{0}$ dominate while for higher redshifts BAO+BBN data dominate which prefer a lower $H_{0}$. The deviation from the DDR can be seen in Fig. 9b. This shows a much stronger violation than in the standard case, see Fig. 6b. This is consistent with our results in the parametrised approach that found $\epsilon=-0.040$. These results reflect the Hubble tension present in the combined data sets very well. However, with this model-independent approach the error remains at $1\sigma$.

6.2.2 DES-SN5YR with DESI DR1

We also implement the DES-SNYR5 [39] likelihood and use the GA on this data set in combination with DESI DR1. In Fig. 10, the GA result for both distances is shown. Since DES-SN5YR covers a redshift range up to 1.13, which is roughly half of the Pantheon+ range, we shrink the $\eta$ axis range to $0.95$–$1.05$. There is no discrepancy to be seen with respect to $\Lambda$CDM. This is mirrored in Fig. 10b which also shows the sampled error for $\eta(z)$. The error for the GA lies slightly below $2\%$ in this case which is much lower than for the Pantheon+ case. For this SNe Ia data set, we only consider this standard case since the results in the parametrised approach are very similar, and the GA including the error calculation is computation intensive.

(a) GA Deviation from $\Lambda$CDM in percent for the angular diameter distance (red), the luminosity distance (green) and the Hubble parameter (purple).

(b) Deviation from the standard DDR $\eta(z)$ for the GA (in solid red) with the associated errors (shaded in orange) compared to $\Lambda$CDM (in dashed black), i.e. no deviation.

Figure 10: Results of the GA using using DES-SN5YR and DESI DR1.

7 Conclusion

This work tests the cosmic DDR using new data sets and provides a new way of looking at the Hubble tension. We use a parametrised approach and a model-agnostic GA to constrain the violation of the cosmic DDR for the Pantheon+ or the DES-SN5YR SNe Ia in combination with the DESI DR1 BAO data set. For the parametrised, model-dependent approach, we calculate the profile likelihood to find contours in the $\Omega_{\mathrm{m}}$–$\epsilon$ plane, where $\epsilon$ quantifies deviations from the DDR given by Eq. 4.1. We find that the best-fit point is $\epsilon=0.019$ and $1.5\sigma$ apart from 0. The reported $\Omega_{\mathrm{m}}$ values from both surveys are over $2\sigma$ apart, which explains this value.

When calibrating the SNe Ia with Cepheids using SH0ES and fixing the value of $\omega_{\mathrm{b}}$ to the value from BBN, we find the best-fit for $\epsilon$ at $-0.040$ which is $6.2\sigma$ from $0$. This is due to the fact that both likelihoods are sensitive to the Hubble constant in this combination and the discrepancy of the two data sets lies at $4.5\sigma$. With the DES-SN5YR data set, we find very similar results when combining with BBN and SH0ES, yielding a $6.5\sigma$ discrepancy.

Then, we employ GA to test this in a model-independent way. For the standard, uncalibrated case, we find a very slight difference from $\eta(z)=1$ which indicates no DDR deviations. To calculate the error, we run the GA 100 times for each case and select 68.27% of functions that are closest to our best-fit for an estimation of the $1\sigma$ error, analogously to the Gaussian error. In the standard analysis, this error is much larger than the DDR discrepancy.

In a next step, the BBN value for $\omega_{\mathrm{b}}$ and Planck values for $H_{0}$ and $\omega_{\mathrm{m}}$ are assumed in the sound horizon at drag redshift when computing the DESI BAO likelihood. This does not significantly change the result for $\eta(z)$. The third case considers SH0ES Cepheid-calibrated SNe Ia with the standard DESI likelihood. There, we find a slightly stronger deviation that starts at much earlier redshifts.

When combining the last two cases, i.e. BBN and Planck values for BAO and Pantheon+SH0ES, we run the analysis using the dimension-full luminosity and angular diameter distances. In this case, we see a much stronger violation of the DDR which is also reflected in the $H(z)$ calculated from the GA result. This is confirmed when looking at the sampled error. At redshifts below 1.4, the DDR violation is close to our sampled error corresponding to $1\sigma$.

We also consider DES-SN5YR using GA. There, the maximum redshift lies at $1.13$, which is less than half of Pantheon+. In that case, there is no preference for a DDR violation. Since the parametrised results for both data sets are extremely similar, we do not test all cases with the GA since they are computationally expensive and leave this to future work.

In conclusion, we show that the DDR is slightly violated in the uncalibrated parametrised approach and more strongly violated when adding BBN and the SH0ES calibration, reflecting the Hubble tension. This holds for both considered SNe Ia data sets. When using a model-independent approach, the DDR is confirmed within 1$\sigma$. This emphasises the importance of model-independent approaches for possible DDR violations.