# Measuring the sound horizon and absolute magnitude of SNIa by maximizing the
consistency between low-redshift data sets

Adrià Gómez-Valent1,2<EMAIL_ADDRESS>1 Dipartimento di Fisica,
Università di Roma Tor Vergata, via della Ricerca Scientifica, 1, 00133, Roma,
Italy  2 INFN, Sezione di Roma 2, Università di Roma Tor Vergata, via della
Ricerca Scientifica, 1, 00133 Roma, Italy

###### Abstract

The comoving sound horizon at the baryon drag epoch, $r_{d}$, encapsulates
very important physical information about the pre-recombination era and serves
as a cosmic standard ruler. On the other hand, the absolute magnitude of
supernovae of Type Ia (SNIa), $M$, is pivotal to infer the distances to these
standard candles. Having access to (at least) one of these two quantities is
crucial to measure the Hubble parameter $H_{0}$ from BAO/SNIa data. In this
work we present a new method to measure how long is the cosmic ruler and how
bright are the standard candles independently from the main drivers of the
$H_{0}$ tension, namely by avoiding (i) the use of CMB data; (ii) the
calibration of SNIa in the first steps of the cosmic distance ladder; and
(iii) the assumption of any concrete cosmological model. We only assume that
SNIa can be safely employed as standard candles and $r_{d}$ as a standard
ruler, together with the validity of the Cosmological Principle and the metric
description of gravity, with photons propagating in null geodesics and the
conservation of the photon number. Our method is based on the minimization of
a loss function that represents the level of inconsistency between the low-
redshift data sets employed in this study, to wit: SNIa, BAO and cosmic
chronometers. In our main analysis we obtain: $r_{d}=(146.0^{+4.2}_{-5.1})$
Mpc, $M=-19.362^{+0.078}_{-0.067}$. The former is fully compatible with
Planck’s $\Lambda$CDM best-fit cosmology, but still leaves plenty of room for
new physics before the decoupling, whereas our constraint on $M$ lies closer
to the value preferred by the concordance model, although it is only $\sim
1.4\sigma$ below the SH0ES measurement.

Cosmology: observations – Cosmology: theory – cosmological parameters – Large-
Scale Structure of Universe

##  I Introduction

How old is the universe? How distant are the cosmological objects that we
detect with our telescopes? Our answers to these questions depend crucially on
how accurately we can measure the current universe’s expansion rate, i.e. the
Hubble parameter $H_{0}$, since it sets both, the time and distance scales in
cosmology. Now, after more than nine decades from its first measurement
(Hubble, 1929) several observational sources tell us that its value should
fall in the range $H_{0}\in(65,75)$ km/s/Mpc. However, there is a $5.0\sigma$
tension between the value inferred for this parameter with the classical
cosmic distance ladder by the SH0ES team, $H_{0}=(73.04\pm 1.01)$ km/s/Mpc
(Riess  _et al._ , 2021a), which is almost fully cosmology-independent, and
the one inferred with the cosmic microwave background (CMB) data by the Planck
Collaboration, under the assumption of the flat $\Lambda$CDM model,
$H_{0}=(67.36\pm 0.54)$ km/s/Mpc (Aghanim  _et al._ , 2020). SH0ES measures in
the first steps of the ladder a somewhat higher value of the absolute
magnitude of supernovae of Type Ia (SNIa), $M$, than the one preferred by the
CMB best-fit standard model, and this is what ultimately triggers the $H_{0}$
tension, see e.g. (Camarena and Marra, 2020a; Efstathiou, 2021). Other
observational teams, though, do not find any substantial discrepancy (Freedman
_et al._ , 2019; Freedman, 2021), and it is still possible that systematic
errors play an important role in this story (Efstathiou, 2020; Mortsell  _et
al._ , 2021). Nevertheless, the tension has persisted and grown in a
consistent way in the last years thanks to the gain in precision of the modern
observational facilities employed to explore the local universe (Riess  _et
al._ , 2009, 2016, 2021b) and measure the CMB anisotropies (Hinshaw  _et al._
, 2013; Ade  _et al._ , 2016). Although there exist some internal tensions at
the $\sim 2\sigma$ level in Planck’s data (e.g. between the $\Lambda$CDM best-
fit parameters obtained from low and high multipoles or the amount of lensing
in the TT spectrum), their results seem to be consistent with other CMB
experiments, as WMAP (Hinshaw  _et al._ , 2013), SPT (Aylor  _et al._ , 2017)
and ACT (Aiola  _et al._ , 2020). Moreover, constraints from baryon acoustic
oscillations (BAO) and big bang nucleosynthesis, independent from CMB, lead to
values of the Hubble parameter lying again in the lower range in the context
of the standard model, more in accordance with Planck (Addison  _et al._ ,
2018; Cuceu  _et al._ , 2019). This is also found when use is made of the
inverse cosmic distance ladder (Aubourg  _et al._ , 2015; Cuesta  _et al._ ,
2015; Feeney  _et al._ , 2019; Camarena and Marra, 2020b), which assumes again
standard pre-recombination (and most of the cases also standard late-time)
physics. Interestingly, standard sirens (Abbott  _et al._ , 2017; Palmese  _et
al._ , 2021) and strongly lensed quasars (Denzel  _et al._ , 2021) also allow
us to measure $H_{0}$, but they are still not able to arbitrate the tension.

Cosmological models with a preferred larger critical energy density around the
matter-radiation equality time or an earlier recombination of protons and
electrons before the photon decoupling can lead to an alleviation of the
$H_{0}$ tension, since they allow for lower values of the comoving sound
horizon at the baryon drag epoch, $r_{d}$, what in turn gives room for larger
values of the Hubble parameter, which is needed to keep the good description
of the BAO data and the location of the first peak of the CMB temperature
anisotropies. This is the case of some models that have been proposed to
alleviate the $H_{0}$ tension, based on early dark energy (Poulin  _et al._ ,
2019; Niedermann and Sloth, 2021; Agrawal  _et al._ , 2019; Gómez-Valent  _et
al._ , 2021), modified gravity (Solà Peracaula  _et al._ , 2019, 2020; Braglia
_et al._ , 2020, 2021), primordial magnetic fields (Jedamzik and Pogosian,
2020), varying atomic constants Liu  _et al._ (2020); Sekiguchi and Takahashi
(2021), and running vacuum models (Solà Peracaula  _et al._ , 2021).

The measurement of $r_{d}$ has been carried out using a plethora of techniques
in the past:

  * •

From CMB data alone or in combination with other data sets, assuming concrete
cosmological models. See e.g. (Verde  _et al._ , 2017a; Aghanim  _et al._ ,
2020; Benisty and Staicova, 2021) and the references of the previous
paragraph.

  * •

Using the local (distance ladder) $H_{0}$ value as an anchor at $z=0$,
together with BAO and SNIa data, assuming the $\Lambda$CDM (Cuesta  _et al._ ,
2015) or using cubic splines to reconstruct the Hubble function (Bernal  _et
al._ , 2016; Aylor  _et al._ , 2019). Alternatively, it is also possible to
use the local calibration of the absolute magnitude of SNIa from SH0ES, again
with BAO and SNIa data Camarena and Marra (2020b).

  * •

Fixing the scale with cosmic chronometers (CCH), employed in combination with
BAO and SNIa, either reconstructing the shape of $H(z)$ with linearly-
interpolated values that are left free in the fit (Heavens  _et al._ , 2014;
Verde  _et al._ , 2017b), with a Multi-Task Gaussian Process (Haridasu  _et
al._ , 2018), or using alternative parametrizations Cai  _et al._ (2022).

  * •

Model-independent determinations of the product $H_{0}r_{d}$ are also possible
using BAO and SNIa data alone L’Huillier and Shafieloo (2017); Shafieloo  _et
al._ (2018); Camarena and Marra (2020b).

We propose here a new method that lets us measure $r_{d}$ and $M$ under very
mild theoretical assumptions (which reduce mainly to the validity of the
Cosmological Principle and the metric description of gravity), and without
making use of the two data sets driving the $H_{0}$ tension, namely the CMB
and the calibration of the SNIa performed in the first steps of the distance
ladder by SH0ES. Our method is close in spirit to those employed in Heavens
_et al._ (2014); Verde  _et al._ (2017b); Haridasu  _et al._ (2018), since we
also use data on BAO, SNIa and CCH, and do not assume a particular
cosmological model. Here we take, though, a different approach. While in
Heavens  _et al._ (2014); Verde  _et al._ (2017b) the authors used a spline
for $H(z)$ and reconstructed it, allowing the values of $H(z_{i})$ at several
redshift nodes $z_{i}$ to vary in the fitting analysis (together with $M$,
$r_{d}$ and the curvature density parameter $\Omega^{0}_{k}$), here we only
vary the last three. Our method is also different from the one presented in
(Haridasu  _et al._ , 2018; Cai  _et al._ , 2022). In this work we do not aim
to find the shape of $H(z)$ preferred by the low-redshift data. We measure the
size of the cosmic ruler and the SNIa intrinsic brightness by minimizing a
loss function that tells us what is the level of inconsistency between the
BAO, SNIa and CCH data sets for each triplet of values
$\vec{\theta}=(M,r_{d},H_{0}^{2}\Omega^{0}_{k})$. Every $\vec{\theta}$ can be
used to translate the BAO+SNIa+CCH data into a list of cosmological distances
and values of the Hubble rate at different redshifts, with their corresponding
covariance matrix. The difference between the distances and Hubble rates
inferred from these data sets is clearly a function of $\vec{\theta}$ and, in
particular, is very sensitive to the calibrators $M$ and $r_{d}$. In absence
of unaccounted systematic errors in the data we should expect the calibrators
$M$ and $r_{d}$ to minimize these differences. Thus, we can associate a larger
probability to a vector $\vec{\theta}$ that leads to a lower statistical
tension between data sets. In this work we employ the index of inconsistency
(IOI) proposed in Lin and Ishak (2017a) to quantify the level of disagreement
between the transformed data sets and obtain the posterior distribution for
the three parameters contained in $\vec{\theta}$. The IOI has been previously
applied in Lin and Ishak (2017a, b); Garcia-Quintero  _et al._ (2019);
Gonzalez  _et al._ (2021) to estimate the tension between Gaussian posterior
distributions of parameters that are obtained from the fitting analyses of the
same model, but using alternative data sets. Hence, it was used as a measure
of the existing tensions between data sets in the context of concrete models,
e.g. the $\Lambda$CDM. Here we deal only with data, i.e. no cosmological model
nor parametrization is assumed, with the aim to reduce as much as we can the
degree of model dependence of our results. We demand the minimization of the
inconsistency between the SNIa, BAO and CCH data sets through the optimal
choice of the calibrators contained in $\vec{\theta}$. Our method has not been
employed before in the literature. Model-independent approaches as the one
proposed and studied in this paper will allow us to further test the viability
of the candidate models aiming to loosen the $H_{0}$ tension, as well as the
calibration of the SNIa in the first steps of the distance ladder.

The paper is organized as follows. In Sec. II we describe the data sets that
we have employed in our analysis. In Sec. III we explain our statistical
method. Our results are presented in Sec. IV, and the conclusions in Sec. V.

##  II Data sets and observables

These are the data sets that we have employed in our study.

###  II.1 Supernovae of Type Ia

The expression for the apparent magnitude $m(z)$ of standardized SNIa reads,

$m(z)=M+25+5\log_{10}\left(\frac{D_{L}(z)}{1\,{\rm Mpc}}\right)\,,$ (1)

with $M$ the absolute magnitude and $D_{L}(z)$ the luminosity distance. In a
Friedmann-Lemaître-Robertson-Walker universe the latter takes the following
form,

$D_{L}(z)=\frac{c(1+z)}{\sqrt{\Omega_{k}^{0}H_{0}^{2}}}\sinh\left(\sqrt{\Omega_{k}^{0}H_{0}^{2}}\int_{0}^{z}\frac{dz^{\prime}}{H(z^{\prime})}\right)\,,$
(2)

where $\Omega_{k}^{0}=-kc^{2}/(R_{0}H_{0})^{2}$ is the curvature density
parameter, with $k=0,-1,+1$ for a flat, open and closed universes,
respectively. $R_{0}$ is a constant with units of length that can be
interpreted as the current radius of curvature in a closed universe.

In this work we employ the SNIa contained in the Pantheon compilation (Scolnic
_et al._ , 2018). For practical purposes we opt to use the 40 binned data
points provided in 111http://github.com/dscolnic/Pantheon, with their
corresponding covariance matrix, which includes the statistical and systematic
uncertainties. The bins span the redshift range $z\in[0.014,1.6123]$.

###  II.2 Baryon acoustic oscillations

Baryons and photons were tightly coupled through electromagnetic interactions
during the pre-recombination era. The intense fight between radiation pressure
and gravity in the photon-baryon plasma generated sound waves that left an
imprint in the distribution of baryons when the universe was cool enough for
the CMB photons to escape and start their travel towards us. The maximum
distance traveled by this wave, the sound horizon at the baryon drag epoch, is
prompted in the distribution of matter in the universe. It manifests as a peak
in the matter two-point correlation function or as wiggles in the matter power
spectrum (Cole  _et al._ , 2005; Eisenstein  _et al._ , 2005). Galaxy surveys
use this characteristic length, $r_{d}$, as a standard ruler with respect to
which they can measure cosmological distances at various redshifts. The latter
can be employed to constrain cosmological models in a quite robust way (Bernal
_et al._ , 2020). Their constraints are given either in terms of the dilation
scale $D_{V}$,

$\frac{D_{V}(z)}{r_{d}}=\frac{1}{r_{d}}\left[D_{M}^{2}(z)\frac{cz}{H(z)}\right]^{1/3}\,,$
(3)

with $D_{M}=(1+z)D_{A}(z)$ being the comoving angular diameter distance, or by
splitting (when possible) the transverse and line-of-sight BAO information,
providing data on $D_{A}(z)/r_{d}$ and $H(z)r_{d}$ separately, with some
degree of correlation. In any metric theory of gravity with photons traveling
on null geodesics and conservation of the photon number, the Etherington
relation (Etherington, 1933) applies,

$D_{A}(z)=\frac{D_{L}(z)}{(1+z)^{2}}\,$ (4)

with $D_{L}(z)$ given by Eq. (2). This expression is useful to translate
luminosity distances into angular diameter distances, and viceversa.

We employ the following BAO data points:

  * •

$D_{V}/r_{d}$ at $z=0.122$ provided in (Carter  _et al._ , 2018), which
combines the dilation scales previously reported by the 6dF Galaxy Survey
(6dFGS) (Beutler  _et al._ , 2011) at $z=0.106$ and the one obtained from the
Sloan Digital Sky Survey (SDSS) Main Galaxy Sample at $z=0.15$ (Ross  _et al._
, 2015).

  * •

The anisotropic BAO data measured by BOSS using the LOWZ ($z=0.32$) and CMASS
($z=0.57$) galaxy samples (Gil-Marín  _et al._ , 2017).

  * •

The dilation scale measurements by WiggleZ at $z=0.44,0.60,0.73$ (Kazin  _et
al._ , 2014).

  * •

$D_{A}/r_{d}$ at $z=0.81$ measured by the Dark Energy Survey Year 1 (DESY1)
(Abbott  _et al._ , 2019). We will also study the impact of substituting this
point by the more recent measurement from DESY3 at the effective redshift
$z=0.835$ (Abbott  _et al._ , 2021), which is in $2.3\sigma$ tension with the
Planck prediction assuming the $\Lambda$CDM.

  * •

The anisotropic BAO data from the extended BOSS Data Release 16 (DR16) quasar
sample at $z=1.48$ (Neveux  _et al._ , 2020; Hou  _et al._ , 2020).

We avoid the use of the anisotropic BAO data obtained from the Ly$\alpha$
absorption and quasars of the final data release (SDSS DR16) of eBOSS at
$z=2.334$ (du Mas des Bourboux  _et al._ , 2020) because it falls out of the
measurement ranges of SNIa and CCH, see Sec. II.1 and II.3. The full BAO data
vector and associated covariance matrix is provided in Appendix A.

###  II.3 Cosmic chronometers

Spectroscopic dating techniques of passively evolving galaxies, i.e. galaxies
with old stellar populations and low star formation rates, have become a good
tool to obtain observational values of the Hubble function at redshifts
$z\lesssim 2$ (Jiménez and Loeb, 2002). These measurements do not rely on any
particular cosmological model, although are subject to other sources of
systematic uncertainties, as to the ones associated with the modeling of
stellar ages, see e.g. Moresco  _et al._ (2020), which is carried out through
the so-called stellar population synthesis (SPS) techniques, and also with a
possible contamination due to the presence of young stellar components in
quiescent galaxies (López-Corredoira  _et al._ , 2017; López-Corredoira and
Vazdekis, 2018; Moresco  _et al._ , 2018). Given a pair of ensembles of
passively evolving galaxies at two different redshifts it is possible to infer
$dz/dt$ from observations under the assumption of a concrete SPS model and
compute $H(z)=-(1+z)^{-1}dz/dt$. Thus, cosmic chronometers allow us to obtain
the value of the Hubble function at different redshifts, contrary to other
probes which do not directly measure $H(z)$, but integrated quantities as e.g.
luminosity distances.

In this study we use the 24 data points on $H(z)$ from CCH in the redshift
range $z\in[0.07,1.965]$ reported in (Jiménez  _et al._ , 2003; Stern  _et
al._ , 2010; Moresco  _et al._ , 2012; Zhang  _et al._ , 2014; Moresco, 2015;
Moresco  _et al._ , 2016; Ratsimbazafy  _et al._ , 2017; Borghi  _et al._ ,
2021). See also Fig. 4 in Appendix C. More concretely, in our main analyses we
make use of the processed sample provided in Table 2 of (Gómez-Valent, 2019),
but adding the data point of (Borghi  _et al._ , 2021) and removing the ones
of (Simon  _et al._ , 2005)222Serious concerns about the statistical analysis
carried out in (Simon  _et al._ , 2005) have been recently raised in Kjerrgren
and Mortsell (2021). Thus, we prefer to omit at the moment the use of these
CCH data.. Our resulting CCH data set is robust, since it introduces
corrections accounting for the systematic errors mentioned above. In addition,
we also study the effect of the covariance matrix of the CCH data given in
Moresco  _et al._ (2020). We compute it with the code
333https://gitlab.com/mmoresco/CCcovariance, cf. Sec. IV and the comments in
Appendix C.

##  III The method

If the SNIa, BAO and CCH data sets described in the previous section are
unbiased they should be, of course, consistent with each other. The level of
consistency, though, depends on the values of the parameters contained in the
vector $\vec{\theta}=(M,r_{d},H_{0}^{2}\Omega_{k}^{0})$, since they are used
to translate the original BAO and SNIa data into a set of cosmological
distances and Hubble rates. Our goal in this paper is to infer the
distribution of these parameters by maximizing the consistency of the data
sets under consideration. To do so we need first to find or define a function
$L(\vec{\theta})$ that lets us quantify the degree of inconsistency between
the data sets for every $\vec{\theta}$, and then use it to build a
distribution that maps the latter to concrete probability values. The function
$L(\vec{\theta})$ should include the contribution of the inconsistency
estimates of all pairs of data sets. $L(\vec{\theta})$ can be thought of as a
loss function, since for wide enough priors the distribution that we build
from it will be maximized for the value of $\vec{\theta}$ that minimizes $L$.
We construct it making use of the index of inconsistency proposed in Lin and
Ishak (2017a), as follows,

$L(\vec{\theta})=L_{1}(M,r_{d})+L_{2}(M,H_{0}^{2}\Omega_{k}^{0})\,,$ (5)

with

$L_{1}(M,r_{d})={\rm IOI}[{\rm BAO,SNIa+CCH}]$ (6)

and

$L_{2}(M,H_{0}^{2}\Omega_{k}^{0})={\rm IOI}[{\rm SNIa,CCH}]\,.$ (7)

(6) is the index of inconsistency between the BAO data and the string
SNIa+CCH, which accounts for the inconsistencies between the pairs (BAO,SNIa)
and (BAO,CCH). On the other hand, (7) is the IOI between the SNIa and the CCH
data sets. The two-experiment IOI is defined as follows (Lin and Ishak,
2017a),

${\rm IOI}[i,j]=\frac{1}{2}\delta^{T}(C^{(i)}+C^{(j)})^{-1}\delta\,,$ (8)

where $C^{(i)}$ is the covariance matrix of the ith data set and
$\delta=\mu^{(i)}-\mu^{(j)}$ is the difference between the data vectors of the
two data sets under consideration, i.e. data sets $i$ and $j$. The IOI is a
generalization of the Mahalanobis distance (Mahalanobis, 1936), and strictly
speaking it is reliable only for Gaussian-distributed data sets444Fortunately,
this is the case in the current study, in very good approximation. We will
comment on this explicitly when needed.. The specific dependence of $L_{1}$
and $L_{2}$ on the parameters of $\vec{\theta}$ and all the details about how
to build them will be duly explained below, in Secs. III.1 and III.2.

Now the question is how to build the probability distribution out of (5). We
consider the following form,

$P(\vec{\theta})=\pi(\vec{\theta})\times\mathcal{N}e^{-L(\vec{\theta})}\,,$
(9)

where $\pi(\vec{\theta})$ is the prior distribution of the parameters and
$\mathcal{N}$ is a constant normalization factor. The motivation for the
exponential term is that given a pair of calibrated data sets $i=1$ and $j=2$
the value of IOI$[1,2]$ is exactly the same as the one that would be obtained
by considering a fake theoretical point located at $\mu^{(1)}$ and fake data
centered at $\mu^{(2)}$ with covariance matrix $C^{(1)}+C^{(2)}$. Different
values of $\vec{\theta}$ lead to a different pair of calibrated data sets, of
course, but their associated IOI can be mapped to the one obtained from a
different fake theoretical point with exactly the same fake data as considered
before. Let us express this in more mathematical terms. Given a reference fake
data vector $\tilde{\mu}$ with covariance matrix $\tilde{C}$ it is always
possible to perform a change of coordinates such that

$(C^{(i)}+C^{(j)})=B^{\rm T}\tilde{C}B\,,$ (10)

with $B$ the transformation matrix associated to the aforesaid change of
coordinates, which will obviously depend on the parameters entering the
covariance matrices $C^{(i)}$ and $C^{(j)}$. Using this we can do

$\displaystyle(\mu^{(i)}-\mu^{(j)})^{\rm T}$
$\displaystyle(C^{(i)}+C^{(j)})(\mu^{(i)}-\mu^{(j)})=$ (11)
$\displaystyle(\mu^{(i)}-\mu^{(j)})^{\rm T}B^{\rm
T}\tilde{C}B(\mu^{(i)}-\mu^{(j)})\,.$ (12)

Thus, if we define

$\mu^{(i,j)}-\tilde{\mu}\equiv B(\mu^{(i)}-\mu^{(j)})\,,$ (13)

we can rewrite (8) as follows,

${\rm IOI}[i,j]=\frac{1}{2}(\mu^{(i,j)}-\tilde{\mu})^{\rm
T}\tilde{C}(\mu^{(i,j)}-\tilde{\mu})\,.$ (14)

One can map every set of data $(i,j)$ to a particular fake theoretical point
$\mu^{(i,j)}$, considering a data vector $\tilde{\mu}$ with covariance matrix
$\tilde{C}$. The only element that depends on $\vec{\theta}$ in the last
formula is the fake theoretical point $\mu^{(i,j)}$. It is then natural to
consider (9), as in standard $\chi^{2}$ analyses. Moreover, if data are
unbiased and our minimal theoretical assumptions hold (9) leads to the correct
estimation of the calibrators, with their uncertainties decreasing for an
increasing number of data points. This can be easily understood. Imagine that
we add unbiased data to our problem, let us say with a similar constraining
power as the preceding data set. If we compute the ratios between the
probability density at the peak $\hat{\theta}$ and at a point $\tilde{\theta}$
away from it before and after adding the new data we find that these ratios
read: $R_{\rm old}=e^{L(\tilde{\theta})-L(\hat{\theta})}$, $R_{\rm new}\sim
e^{2(L(\tilde{\theta})-L(\hat{\theta}))}$, so $R_{\rm new}/R_{\rm old}\sim
R_{\rm old}>1$. This means that after the addition of the new data the point
in parameter space located far away from the peak is less preferred than
before when compared to the one at which (9) is maximized, what in turn leads
to lower uncertainties for the parameters contained in $\vec{\theta}$.

Once we have the probability distribution (9) we can compute the confidence
regions for the parameters of $\vec{\theta}$ either by evaluating our
distribution on a three-dimensional grid or by sampling it with a Monte Carlo
method. We opt to do the latter, see the details in Sec. III.3.

Figure 1: This figure shows the histograms obtained for three out of the
eleven elements of $\mu^{\rm BAO}$ (first row) and $\mu^{\rm SNIa+CCH}$
(second row), obtained using $10^{4}$ points in the sampling process and using
$M=-19.454$ and $r_{d}=145$ Mpc. The values of the dilation scale
$D_{V}(z=0.60)$ and the angular diameter distance $D_{A}(z=0.57)$ are
expressed in Mpc, whereas the Hubble function $H(z=0.32)$ is in km/s/Mpc. As
mentioned in the main text, the distributions are very well approximated by
Gaussians with the mean and standard deviation inferred from the sample. The
latter are also plotted with solid lines, in blue. See Sec. III.1 for details.

###  III.1 IOI between BAO and SNIa+CCH

Given a value of $r_{d}$ we can rewrite the BAO constraints on
$\\{D_{V}(z)/r_{d},D_{A}(z)/r_{d},H(z)r_{d}\\}$ of Sec. II.2 (cf. (23)-(24) in
Appendix A) as constraints on $\\{D_{V}(z),D_{A}(z),H(z)\\}$, with a data
vector

$\mu^{\rm BAO}(r_{d})=\begin{pmatrix}D_{V}(z=0.122)\\\ D_{V}(z=0.44)\\\
D_{V}(z=0.60)\\\ D_{V}(z=0.73)\\\ D_{A}(z=0.32)\\\ D_{A}(z=0.57)\\\
D_{A}(z=0.81)\\\ D_{A}(z=1.48)\\\ H(z=0.32)\\\ H(z=0.57)\\\
H(z=1.48)\end{pmatrix}\,,$ (15)

and its associated covariance matrix $C^{\rm BAO}(r_{d})$. Both depend on
$r_{d}$. The resulting distribution is Gaussian in very good approximation, as
we have explicitly checked for different values of the sound horizon, cf. Fig.
1. In order to compute the IOI between the BAO and SNIa+CCH data sets we need
to build the analogous quantities from the SNIa+CCH joint data set, i.e.
$\mu^{\rm SNIa+CCH}$ and $C^{\rm SNIa+CCH}$. As we do not have SNIa and CCH
data at the exact BAO redshifts, we need first to compute the extrapolated
values of $D_{A}(z)$ and $H(z)$ at these $z$’s using some technique that
allows us to get the vectors

$\mu^{\rm SNIa}(M)=\begin{pmatrix}D_{A}(z=0.122)\\\ D_{A}(z=0.44)\\\
D_{A}(z=0.60)\\\ D_{A}(z=0.73)\\\ D_{A}(z=0.32)\\\ D_{A}(z=0.57)\\\
D_{A}(z=0.81)\\\ D_{A}(z=1.48)\end{pmatrix}$ (16)

and

$\mu^{\rm CCH}=\begin{pmatrix}H(z=0.122)\\\ H(z=0.44)\\\ H(z=0.60)\\\
H(z=0.73)\\\ H(z=0.32)\\\ H(z=0.57)\\\ H(z=1.48)\end{pmatrix}\,,$ (17)

with their individual covariances, $C^{\rm SNIa}(M)$ and $C^{\rm CCH}$. Notice
that $\mu^{\rm SNIa}$ and $C^{\rm SNIa}$ depend on $M$, since we can only get
angular diameter distances from the SNIa data after fixing the absolute
magnitude of the supernovae. Using all these ingredients we can build the
final mean and covariance matrix $\mu^{\rm SNIa+CCH}$ and $C^{\rm SNIa+CCH}$
to be employed (together with $\mu^{\rm BAO}$ and $C^{\rm BAO}$) in the
computation of $L_{1}$ (6).

For a given value of $M$, and making use of Eqs. (2) and (4), we can translate
the SNIa apparent magnitudes $m(z)$ into data on $D_{A}(z)$ at the redshifts
of the (binned) SNIa data. By sampling this distribution and using e.g. a
cubic interpolation method, it is easy to infer the Gaussian distribution for
the $D_{A}$’s at the redshifts of interest, i.e. those specified in (16). In
order to be more efficient, though, it is better actually to split $D_{A}$ as
follows,

$D_{A}(z)=10^{-M/5}B(z)\,{\rm Mpc}\quad{\rm with}\quad
B(z)=\frac{10^{\frac{m(z)-25}{5}}}{(1+z)^{2}}\,,$ (18)

and sample the distribution of $m$’s to generate the mean ($\mu_{B}$) and
covariance matrix ($C_{B}$) for the $B$’s at the redshifts of (16) before
starting the Monte Carlo, since this part of $D_{A}$ does not depend on $M$
and hence can be employed at each step of the sampling process. This
distribution is Gaussian too, see Appendix B. Given a value of $M$, it is easy
to obtain then the mean and covariance for $D_{A}$ from the distribution of
$B(z)$. We just need to do:

$\mu^{\rm SNIa}=10^{-M/5}\mu_{B}\quad;\quad C^{\rm SNIa}=10^{-2M/5}C_{B}\,.$
(19)

We proceed in a similar way to obtain $\mu^{\rm CCH}$ (17) and the
corresponding covariance matrix $C^{\rm CCH}$. In this case the result does
not depend on any of the parameters contained in $\vec{\theta}$. Therefore, we
can employ the same vector $\mu^{\rm CCH}$ and covariance $C^{\rm CCH}$ in
each step of the Monte Carlo routine. We reconstruct the shape of $H(z)$ using
the data described in Sec. II.3 and Gaussian Processes (GP), with a Gaussian
kernel. In appendix C we describe the GP reconstruction of the Hubble function
and provide the resulting $\mu^{\rm CCH}$ and $C^{\rm CCH}$, and also comment
on several technical aspects for the interested reader.

The obtention of $\mu^{\rm SNIa+CCH}(M)$ and $C^{\rm SNIa+CCH}(M)$ from
($\mu^{\rm SNIa}(M),C^{\rm SNIa}$(M)) and ($\mu^{\rm CCH},C^{\rm CCH}$) is
straightforward. It can be done through a simple sampling of the two
multivariate Gaussians. The result is in very good approximation Gaussian too,
cf. again Fig. 1. Equipped with these tools we can finally evaluate $L_{1}$
(6), which is of course a function of the calibrators $M$ and $r_{d}$.

Loss function | $M$ |  $r_{d}$ [Mpc] |  $H_{0}^{2}\Omega_{k}^{0}$ [(km/s/Mpc)2] | $\Omega_{k}^{0}$
---|---|---|---|---
$L_{1}(M,r_{d})$ | $-19.355\pm 0.095$ | $145.5^{+5.4}_{-6.0}$ | - | -
$L(M,r_{d},0)$, flat universe | $-19.374\pm 0.080$ | $146.5^{+4.8}_{-5.7}$ | 0 | 0
$L(M,r_{d},H_{0}^{2}\Omega_{k}^{0})$ | $-19.362^{+0.078}_{-0.067}$ | $146.0^{+4.2}_{-5.1}$ | $-36\pm 490$ | $-0.01\pm 0.10$
$\mathcal{L}(M,r_{d},H_{0}^{2}\Omega_{k}^{0})$ | $-19.363\pm 0.076$ | $146.7\pm 5.1$ | $-24\pm 490$ | $0.00\pm 0.10$
$L_{c}(M,r_{d},H_{0}^{2}\Omega_{k}^{0})$ | $-19.399\pm 0.098$ | $148.3\pm 6.6$ | $-38\pm 490$ | $-0.01\pm 0.10$
Table 1: Means and uncertainties at the $1\sigma$ c.l. of the fitting
parameters obtained with the loss functions listed in Sec. III.3. The
constraints on $\Omega^{0}_{k}$ are computed by breaking the degeneracy in the
$H_{0}-\Omega_{k}^{0}$ plane with the Gaussian prior $H_{0}=(70.72\pm 6.44)$
km/s/Mpc (or $70.36\pm 5.69$ km/s/Mpc when we consider the covariance matrix
of the CCH data) obtained from our reconstruction with GP (cf. Appendix C).
See the comments in Sec. IV and Fig. 2.

###  III.2 IOI between SNIa and CCH

We compute the index of inconsistency between the SNIa and CCH data sets by
first noting that for a given pair $(M,H_{0}^{2}\Omega_{k}^{0})$ we can
translate a particular value of $m(z_{i})$ into a value of the Hubble function
$H(z_{i})$, if we also know the derivative of the apparent magnitude at that
redshift. This becomes evident if we perform the derivative of Eq. (2) with
respect to the redshift, and isolate $H(z)$,

$H(z)=\frac{\left[\left(\frac{c(1+z)}{D_{L}(z)}\right)^{2}+H_{0}^{2}\Omega_{k}^{0}\right]^{1/2}}{\frac{\ln(10)}{5}\frac{\partial
m}{\partial z}-\frac{1}{1+z}}\,,$ (20)

where the luminosity distance can be written in terms of a function
$\tilde{B}(z)$ that does not depend on $M$,

$D_{L}(z)=10^{-\frac{(M+25)}{5}}\tilde{B}(z)\,{\rm Mpc}\quad{\rm
with}\quad\tilde{B}(z)=10^{\frac{m(z)}{5}}\,.$ (21)

Formula (20) is very useful to translate the SNIa data into data on the Hubble
rate and establish the link with the CCH data. We can sample the distribution
of $m$’s at the SNIa redshifts and obtain from it the distribution of
$\tilde{B}$’s and $\partial m/\partial z$’s at those redshifts at which we
have CCH data555We exclude the CCH data point at $z=1.965$ Moresco (2015)
because we have no SNIa data at this high redshift., e.g. making use again of
a cubic interpolation method. The resulting distribution is Gaussian. The
derivative of the apparent magnitude can be computed numerically, using finite
differences with $\Delta z=0.01$. Finally, given a pair
$(M,H_{0}^{2}\Omega_{k}^{0})$ we can construct from the latter the
distribution of values of the Hubble function at the CCH redshifts using (20).
It is obviously a function of $M$ and $H_{0}^{2}\Omega_{k}^{0}$. The
computation of $L_{2}$ (7) is at this stage very easy because we already have
all the ingredients to apply (8). In order to sample physically motivated
values of the product $H_{0}^{2}\Omega_{k}^{0}$ we use the Gaussian prior

$\pi(H_{0}^{2}\Omega_{k}^{0})=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(H_{0}^{2}\Omega_{k}^{0})^{2}}{2\sigma^{2}}}\,,$
(22)

with $\sigma=500$ (km/s/Mpc)2. The latter is motivated by the constraint on
$H_{0}$ obtained from our GP reconstruction, $H_{0}=(70.72\pm 6.44)$ km/s/Mpc
(cf. Appendix C), and also by the ones on $\Omega_{k}^{0}$ obtained under some
CMB data sets (Aghanim  _et al._ , 2020; Handley, 2021; Di Valentino  _et al._
, 2019), which allow values of $\Omega_{k}^{0}\sim-0.1$ at $\sim 1.5\sigma$
c.l. Much tighter constraints on the curvature parameter are derived when CMB
lensing, SNIa and/or BAO data are also considered (Aghanim  _et al._ , 2020;
Efstathiou and Gratton, 2020), but we want to proceed as model-independently
as possible. This is why we choose a wide prior for this parameter, but still
forbidding values that clearly fall out of the region allowed by the CMB. We
center the Gaussian prior at $\Omega_{k}^{0}=0$. In this way we keep a
reasonable upper bound on the absolute value of this parameter, without
favoring a particular sign. As we will discuss in Sec. IV,
$H_{0}^{2}\Omega_{k}^{0}$ has a very low impact on our constraints on $M$ and
$r_{d}$, and its posterior is basically dominated by the prior.

###  III.3 Monte Carlo analyses

We sample the distributions (9) built with the following functions:

  * •

$L_{1}(M,r_{d})$, Eq. (6).

  * •

$L(M,r_{d},H_{0}^{2}\Omega_{k}^{0}=0)$ assuming a flat universe.

  * •

$L(M,r_{d},H_{0}^{2}\Omega_{k}^{0})$.

  * •

$L(M,r_{d},H_{0}^{2}\Omega_{k}^{0})$, but using the data point from DESY3
(Abbott  _et al._ , 2021) instead of DESY1 (Abbott  _et al._ , 2019). We call
this loss function $\mathcal{L}(M,r_{d},H_{0}^{2}\Omega_{k}^{0})$.

  * •

$L(M,r_{d},H_{0}^{2}\Omega_{k}^{0})$, but considering the covariance matrix of
the CCH data Moresco  _et al._ (2020). We denote this loss function
$L_{c}(M,r_{d},H_{0}^{2}\Omega_{k}^{0})$.

We make use of the Metropolis–Hastings algorithm (Metropolis  _et al._ , 1953;
Hastings, 1970). For the calibrators we use the flat priors
$M\in[-19.7,-19.0]$ and $r_{d}\in[110,170]$ Mpc, which are much wider than
their associated uncertainties and hence have no impact on their posterior
distributions. The results obtained in our Monte Carlo analyses are presented
and discussed in the next section.

Figure 2: Contour plots at $1\sigma$ and $2\sigma$ c.l. in all the planes of
the ($M,r_{d},H_{0}^{2}\Omega_{k}^{0}$)-parameter space, obtained with the
loss function $L(M,r_{d},H_{0}^{2}\Omega_{k}^{0})$, cf. Sec. III.3. $r_{d}$ is
expressed in Mpc and $H_{0}^{2}\Omega_{k}^{0}$ in (km/s/Mpc)2. The vertical
green bands correspond to the (1$\sigma$ and 2$\sigma$) constraint on the
absolute magnitude of SNIa obtained by SH0ES from the first steps of the
distance ladder Riess  _et al._ (2021a), $M=-19.253\pm 0.027$. Our constraint
on $H_{0}^{2}\Omega_{k}^{0}$ is basically dominated by the prior (22). In the
absence of this prior we would obtain much more elongated contours in the
$y$-direction both in the $(M,H_{0}^{2}\Omega_{k}^{0})$ and
$(r_{d},H_{0}^{2}\Omega_{k}^{0})$ planes. We have made use of GetDist Lewis
(2019) to analyze the Markov chains and generate this figure. See Sec. IV for
further comments.

##  IV Results

Our constraints on $M$, $r_{d}$, $H_{0}^{2}\Omega^{0}_{k}$ and
$\Omega^{0}_{k}$ are shown in Table I. In Fig. 2 we plot the confidence
regions at $1$ and $2\sigma$ c.l. in all the relevant planes of the parameter
space, together with the corresponding one-dimensional posteriors. The central
values for $M$ and $r_{d}$ are compatible in all the scenarios analyzed in
this study, even when the loss function is only built with the
IOI(BAO,SNIa+CCH), i.e. when $L=L_{1}$. In this case we find $M=-19.355\pm
0.095$ and $r_{d}=(145.5^{+5.4}_{-6.0})$ Mpc. If we penalize also the degree
of inconsistency between SNIa and CCH by considering the loss function
$L=L_{1}+L_{2}$ we find, though, that the uncertainties on $M$ and $r_{d}$
decrease by $\sim 20\%$, yielding $M=-19.362^{+0.078}_{-0.067}$ and
$r_{d}=(146.0^{+4.2}_{-5.1})$ Mpc. These are our main results. Notice that the
inclusion of $L_{2}$ is important to study also the compatibility between the
CCH and SNIa data sets, which is not automatically ensured with the use of
$L_{1}$ alone. It allows us to put more stringent and more precise limits on
the calibrators, i.e. to better draw the borders of the allowed region in
parameter space. The absolute magnitude of SNIa is lower than the SH0ES
measurement, at $\sim 1.4\sigma$ c.l. (cf. Fig. 2). This result is compatible
with Efstathiou (2020); Renzi and Silvestri (2020). Using the value of
$M-5\log_{10}\left(H_{0}{\rm Mpc}/c\right)\sim-1.185$ obtained by SH0ES from
the analysis of the SNIa in the Hubble flow we obtain the following values of
$H_{0}$ associated to our $M-\sigma_{M}$, $M$ and $M+\sigma_{M}$,
respectively: $H_{0}=(67.3,69.5,72.0)$ km/s/Mpc. The central value of the
Hubble parameter is very close to the distance ladder measurement obtained
with the TRGB calibration Freedman  _et al._ (2019); Freedman (2021). However,
it lies at $<1\sigma$ from the Planck result obtained under the standard model
Aghanim  _et al._ (2020) and at $\sim 1\sigma$ from the SH0ES value (Riess
_et al._ , 2021a). This departure is not significant, but it will be
interesting to revisit this calculation in the future to check whether this
difference grows. If so, and if not caused by systematics in the data, it
could indicate a flaw in our minimal (although still model-dependent)
assumptions Macpherson and Heinesen (2021); Krishnan  _et al._ (2021).
Similarly, if the constraints on the calibrators $M$ and $r_{d}$ obtained with
our method from different subsets of future CCH+BAO+SNIa data exhibit some
level of tension this could be due, again, to systematic biases in the data or
to some incorrect theoretical assumptions.

On the other hand, the comoving sound horizon agrees with the $\Lambda$CDM-
based inference by Planck, although the one and two sigma bands also encompass
much lower values of $r_{d}$, as those found e.g. in early dark energy and
modified gravity scenarios that have been explored in the literature to
alleviate (in greater or lesser extent) the $H_{0}$ tension. The correlation
coefficient in the $M-r_{d}$ plane is negative, as expected, since larger
values of $M$ lead also to larger values of the Hubble function and,
therefore, $r_{d}$ needs to decrease in order to calibrate appropriately the
BAO ruler.

The error bars of both, $M$ and $r_{d}$, are a bit smaller if a flat universe
is taken for granted, of course, just because we reduce by one the
dimensionality of our parameter space. It is also worth to mention that the
change of the DESY1 BAO data point Abbott  _et al._ (2019) by the DESY3 one
Abbott  _et al._ (2021), which is in $2.3\sigma$ tension with the Planck
$\Lambda$CDM best-fit cosmology, does not alter our results significantly. The
impact of considering the covariance matrix of the CCH data from Moresco  _et
al._ (2020) is not important either, and our constraints on the curvature
parameter are weak and essentially dominated by the prior (cf. again Table I).
We find $\Omega_{k}^{0}=-0.01\pm 0.10$.

##  V Conclusions

The absolute magnitude of supernovae of Type Ia, $M$, and the comoving sound
horizon, $r_{d}$, are the anchors of the direct and inverse cosmic distance
ladders, respectively, and therefore play a very important role in the
discussion of the $H_{0}$ tension. Many models aiming to solve the latter
modify the physics of the pre-recombination era and predict a value of $r_{d}$
which is much lower than the one preferred by the $\Lambda$CDM when it is
constrained with CMB data. In this paper we have employed low-$z$ data to
constrain $M$ and $r_{d}$ under minimal model assumptions (which reduce to the
Cosmological Principle and an underlying metric theory of gravity) and also
avoiding the use of any input coming from the main drivers of the $H_{0}$
tension, similarly to what was done by other groups in the past Heavens  _et
al._ (2014); Verde  _et al._ (2017b); Haridasu  _et al._ (2018). We have
applied a novel method, though, based on the minimization of a loss function
that quantifies the degree of inconsistency between the BAO, SNIa and cosmic
chronometer data sets, which is built from the IOI estimator proposed by Lin
and Ishak Lin and Ishak (2017a). It is only a function of $M$, $r_{d}$ and
$H_{0}^{2}\Omega_{k}^{0}$. Our constraints read $M=-19.362^{+0.078}_{-0.066}$
and $r_{d}=(146.0^{+4.2}_{-5.1})$ Mpc at $1\sigma$ c.l. The former is slightly
lower ($1.4\sigma$) than the value inferred in the first steps of the cosmic
distance ladder by SH0ES, $M=-19.253\pm 0.027$ Riess  _et al._ (2021a),
whereas the comoving sound horizon is not very constrained by our principle of
consistency. It is fully compatible with the standard model value from the
TT,TE,EE+lowE+lensing analysis by Planck Aghanim  _et al._ (2020),
$r_{d}=(147.09\pm 0.26)$ Mpc, but still leaves plenty of room for new physics
Poulin  _et al._ (2019); Niedermann and Sloth (2021); Agrawal  _et al._
(2019); Gómez-Valent  _et al._ (2021); Solà Peracaula  _et al._ (2019, 2020);
Braglia  _et al._ (2020, 2021); Jedamzik and Pogosian (2020); Liu  _et al._
(2020); Sekiguchi and Takahashi (2021); Solà Peracaula  _et al._ (2021); Verde
_et al._ (2017a). With the advent of future data e.g. from Euclid and LSST,
our method will provide tighter constraints on these relevant parameters,
allowing us to study the viability of cosmological models with non-standard
pre-recombination physics and to test the agreement between the two ends of
the cosmic ladder in a quite model-independent way.

Acknowledgements

The author is funded by the Istituto Nazionale di Fisica Nucleare (INFN)
through the project “Dark Energy and Modified Gravity Models in the light of
Low-Redshift Observations” (n. 22425/2020). He is grateful to the Institute
for Theoretical Physics (ITP) Heidelberg for letting him use its computational
facilities from remote, and to Dr. Elmar Bittner for his valuable technical
help in the use of the ITP computational resources. He would also like to
thank Prof. Joan Solà Peracaula for feedback on this manuscript, and Dr.
Javier de Cruz Pérez for pointing out the existence of the CCH covariance
matrix. Finally, the author wants to express his gratitude to the anonymous
referee for their insightful suggestions.

##  Appendix A BAO data

These are the full data vector and covariance matrix of the BAO data described
in Sec. II.2:

$\vec{d}_{\rm BAO}=\begin{pmatrix}D_{V}(z=0.122)/r_{d}\\\
D_{V}(z=0.44)/r_{d}\\\ D_{V}(z=0.60)/r_{d}\\\ D_{V}(z=0.73)/r_{d}\\\
D_{A}(z=0.32)/r_{d}\\\ D_{A}(z=0.57)/r_{d}\\\ D_{A}(z=0.81)/r_{d}\\\
D_{A}(z=1.48)/r_{d}\\\ r_{d}H(z=0.32)\\\ r_{d}H(z=0.57)\\\
c/(r_{d}H(z=1.48))\end{pmatrix}=\begin{pmatrix}3.65\\\ 11.55\\\ 14.95\\\
16.93\\\ 6.60\\\ 9.39\\\ 10.75\\\ 12.18\\\ 11549\\\ 14021\\\
13.23\end{pmatrix}\,,$ (23)

$\mathcal{D}_{\rm
BAO}=\setcounter{MaxMatrixCols}{11}\begin{pmatrix}0.0144&0&0&0&0&0&0&0&0&0&0\\\
-&4.812&-2.465&1.037&0&0&0&0&0&0&0\\\ -&-&3.770&-1.587&0&0&0&0&0&0&0\\\
-&-&-&3.650&0&0&0&0&0&0&0\\\ -&-&-&-&0.0179&0&0&0&28.929&0&0\\\
-&-&-&-&-&0.0106&0&0&0&13.827&0\\\ -&-&-&-&-&-&0.1849&0&0&0&0\\\
-&-&-&-&-&-&-&0.101&0&0&0.0057\\\ -&-&-&-&-&-&-&-&148099&0&0\\\
-&-&-&-&-&-&-&-&-&50717&0\\\ -&-&-&-&-&-&-&-&-&-&0.2195\end{pmatrix}\,.$ (24)

The product $r_{d}H(z)$ is given in km/s. The other quantities are
dimensionless. We do not specify the lower half of the covariance matrix (it
is of course symmetric). When the data point from DESY3 Abbott  _et al._
(2021) is used instead of the one from DESY1 Abbott  _et al._ (2019), we
change $D_{A}(z=0.81)/r_{d}=10.75\pm 0.43$ by $D_{A}(z=0.835)/r_{d}=10.31\pm
0.28$.

Figure 3: Histograms of $B(z_{i})$ (18) at the redshifts $z_{i}$ of interest,
see (25). They have been obtained from the sampling of $10^{4}$ points with
the procedure explained in Sec. III.1. We also plot the Gaussian distributions
built making use of (25). One can see that the distributions are almost
perfectly Gaussian.

##  Appendix B Gaussian distribution for $B(z)$

The distribution for the values $B(z_{i})$ (18) at the redshifts $z_{i}$ of
interest, needed to compute $\mu^{\rm SNIa}$ and $C^{\rm SNIa}$ (see Sec.
III.1), is Gaussian (cf. Fig. 3). The corresponding mean and covariance matrix
read,

$\mu_{B}=\begin{pmatrix}B(z=0.122)\\\ B(z=0.44)\\\ B(z=0.60)\\\ B(z=0.73)\\\
B(z=0.32)\\\ B(z=0.57)\\\ B(z=0.81)\\\
B(z=1.48)\end{pmatrix}=\begin{pmatrix}0.06094\\\ 0.15869\\\ 0.18705\\\
0.20048\\\ 0.13150\\\ 0.18374\\\ 0.21156\\\ 0.22756\end{pmatrix}\quad
C_{B}=10^{-8}\begin{pmatrix}39.45&0.26&3.91&-6.15&2.36&1.60&-9.17&-23.26\\\
-&410.30&6.99&37.58&8.21&33.21&0.72&-49.87\\\
-&-&460.39&70.55&-10.73&364.85&72.72&113.22\\\
-&-&-&784.12&-26.46&119.12&291.36&51.58\\\
-&-&-&-&168.67&-15.79&-33.62&-83.11\\\ -&-&-&-&-&521.70&96.06&-112.98\\\
-&-&-&-&-&-&906.54&716.72\\\ -&-&-&-&-&-&-&10869.4\end{pmatrix}\,.$ (25)

##  Appendix C Gaussian Processes to reconstruct the Hubble function

Gaussian processes (GPs) can be thought of as Gaussian distributions over
functions (Rasmussen and Williams, 2006). They can be used to reconstruct the
most probable underlying continuous function describing a set of Gaussian-
distributed data and obtain its associated confidence regions without assuming
any parametrization for the aforesaid function. One has to assume, though, a
concrete kernel, which is in charge of controlling the correlations between
the points of the reconstructed function. A GP is defined by two objects: its
mean function, $\mu(z)$, and its two-point covariance function
$\mathcal{C}(z,z^{\prime})$,

$\xi(z)\sim\mathcal{GP}\left(\mu(z),\mathcal{C}(z,z^{\prime})\right)\,.$ (26)

Any realization $\xi(z)$ of a GP is a continuous curve, and it is of course
possible to compute the probability of finding a realization inside any region
$\xi(z)\pm\Delta\xi(z)$. The covariance function $\mathcal{C}$ is defined as
follows:

Figure 4: Reconstructed shape of the Hubble function using Gaussian Processes
under different kernels and CCH data sets. $H(z)$ is expressed in km/s/Mpc.
The data points employed in our main analysis are shown in red, whereas those
from (Simon  _et al._ , 2005) appearing in the lower left plot are in green.
The dotted lines are the mean curves of the reconstructed Hubble function and
the solid ones are the borders of the corresponding $2\sigma$ confidence
regions. In all plots we show the results obtained with the Gaussian kernel
(28) performing the full Bayesian analysis for the hyperparameters (in black).
We compare it with the reconstructed Hubble function (in blue) obtained by
minimizing the marginal likelihood (31) and using: the Gaussian kernel (28)
(upper left plot); the Cauchy kernel (29) (upper right plot); and the Gaussian
kernel, considering also the data from (Simon  _et al._ , 2005) (lower left
plot). The values of $H_{0}$ obtained for these three cases read $(66.16\pm
6.02,67.05\pm 6.26,66.02\pm 5.31)$ km/s/Mpc, respectively, at $1\sigma$ c.l.,
whereas for the main scenario (in black) we obtain $(70.72\pm 6.44)$ km/s/Mpc.
The orange curve in the lower right plot is the reconstruction obtained as the
black one, but using the covariance matrix of the CCH data (see Sec. II.3). In
this case we obtain $H_{0}=(70.36\pm 5.69)$ km/s/Mpc.

$\mathcal{C}(z_{1},z_{2})=\left\\{\begin{array}[]{ll}\mathcal{K}(z_{1},z_{2})&\mbox{if
}z_{1}{\rm\,or}\,z_{2}\,\mbox{are not among the data points}\\\
\mathcal{D}(z_{1},z_{2})+\mathcal{K}(z_{1},z_{2})&\mbox{otherwise}\end{array}\right.$
(27)

where $\mathcal{D}$ is the covariance matrix of the data points and
$\mathcal{K}$ is the kernel function. Two of the most famous kernels are:

  * •

The Gaussian kernel:

$\mathcal{K}(z_{1},z_{2})=\sigma_{f}^{2}e^{-\frac{1}{2}\left(\frac{z_{1}-z_{2}}{l_{f}}\right)^{2}}\,.$
(28)

  * •

The Cauchy kernel:

$\mathcal{K}(z_{1},z_{2})=\frac{\sigma_{f}^{2}l_{f}}{(z_{1}-z_{2})^{2}+l_{f}^{2}}\,.$
(29)

$\sigma_{f}$ and $l_{f}$ are the so-called hyperparameters of the kernel
function. The first one controls the uncertainties’ size, whereas the second
regulates the scope of the correlations in $z$, i.e. for distances
$|z_{1}-z_{2}|\gg l_{f}$ the values of the function at $z_{1}$ and $z_{2}$ are
very uncorrelated. The reconstructed function depends on the values of the
hyperparameters and, in principle, also on our choice of the kernel. We will
study the impact of the latter below. For now, let us see how to select
$\sigma_{f}$ and $l_{f}$ assuming a concrete kernel. In the GP philosophy, our
data set is conceived just as a realization of the Gaussian process. The
hyperparameters are usually chosen so as to maximize the probability of the GP
to produce our data set. If we marginalize the GP (26) over the points at
$z^{*}$ (with no data) we get the following multivariate normal distribution,

$\vec{\xi}\sim\mathcal{N}\left(\\{\mu_{i}(\tilde{z}_{i})\\},\mathcal{C}\right)\,,$
(30)

where $i=1,...,N$, with $N$ being the dimension of the vector of data points
$\vec{y}\equiv\\{\tilde{z}_{i},y_{i}\\}$ at our disposal.
$\mu_{i}(\tilde{z}_{i})$ can be set e.g. to $\vec{0}$ $\forall{i}$, since the
result is almost insensitive to this. Thus, the hyperparameters will be
obtained upon the minimization of

$-2\ln\mathcal{L}(\sigma_{f},l_{f})=N\ln(2\pi)+\ln|\mathcal{C}(\sigma_{f},l_{f})|+\vec{y}^{T}\mathcal{C}^{-1}(\sigma_{f},l_{f})\vec{y}\,,$
(31)

with $\mathcal{L}$ being here the marginal likelihood and $|\mathcal{C}|$ the
determinant of $\mathcal{C}$. Using (26) one can compute now the conditional
probability of finding a given realization of the Gaussian process in the case
in which $\xi(\tilde{z}_{i})=y_{i}(\tilde{z}_{i})$. The resulting mean and
covariance function extracted from the conditioned GP read, respectively,

$\bar{\xi}(z^{*})=\sum_{i,j=1}^{N}\mathcal{C}^{-1}(\tilde{z}_{i},\tilde{z}_{j})y(\tilde{z}_{j})\mathcal{K}(\tilde{z}_{i},z^{*})\,,$
(32) ${\rm
Cov}(z^{*}_{1},z^{*}_{2})=\mathcal{K}(z_{1}^{*},z_{2}^{*})-\sum_{i,j=1}^{N}\mathcal{C}^{-1}(\tilde{z}_{i},\tilde{z}_{j})\mathcal{K}(\tilde{z}_{i},z_{1}^{*})\mathcal{K}(\tilde{z}_{j},z_{2}^{*})\,.$
(33)

A complete Bayesian analysis would demand to marginalize over the
hyperparameters, rather than optimizing Eq. (31). We prefer to consider in our
main analysis the full distribution for $(l_{f},\sigma_{f})$. In Fig. 4 we
show a plot comparing the reconstructed shape of $H(z)$ obtained with the
Gaussian kernel (28) by both, optimizing and performing the complete Bayesian
analysis. We also compare these results with those obtained with the
optimization procedure using the Cauchy kernel (29), and using the optimized
Gaussian kernel, but adding the data points from (Simon  _et al._ , 2005) that
we have excluded in our main analysis. In addition, we also study the impact
of the covariance matrix of the CCH data Moresco  _et al._ (2020). Figure 4
allows us to check that: (i) the impact of using the full Bayesian approach is
not huge, but it is safer to propagate the uncertainties on the
hyperparameters; (ii) our reconstructed Hubble function is stable under the
change of the kernel; (iii) the effect that causes the removal of the data
points from Simon  _et al._ (2005) is sizable, but the reconstructed shape is
still compatible at $1\sigma$ c.l. with the one employed by us in the
evaluation of $L(\vec{\theta})$ (5); (iv) the impact of the covariance matrix
of the CCH data is small. For previous reconstructions of $H(z)$ with CCH data
see e.g. (Busti  _et al._ , 2014; Verde  _et al._ , 2014; Yu  _et al._ , 2018;
Haridasu  _et al._ , 2018; Gómez-Valent and Amendola, 2018).

We use our reconstructed Hubble function to obtain the data vector $\mu^{\rm
CCH}$ and its corresponding covariance matrix $C^{\rm CCH}$, which are needed
to compute $L_{1}$ (6) (see Sec. III.1 for details). In our main analyses we
employ

$\mu^{\rm CCH}=\begin{pmatrix}H(z=0.122)\\\ H(z=0.44)\\\ H(z=0.60)\\\
H(z=0.73)\\\ H(z=0.32)\\\ H(z=0.57)\\\
H(z=1.48)\end{pmatrix}=\begin{pmatrix}75.55\\\ 88.20\\\ 94.60\\\ 99.81\\\
83.41\\\ 93.40\\\ 129.94\end{pmatrix}\,,\qquad C^{\rm
CCH}=\begin{pmatrix}24.89&7.57&-1.12&-8.15&14.10&0.51&-47.92\\\
-&11.33&13.11&14.50&9.95&12.78&21.69\\\ \ -&-&20.21&25.95&7.75&18.88&58.11\\\
-&-&-&35.28&5.93&23.80&88.34\\\ -&-&-&-&11.55&8.17&-5.03\\\
-&-&-&-&-&17.74&51.22\\\ -&-&-&-&-&-&270.45\\\ \end{pmatrix}\,,$ (34)

whereas in the case in which we use the covariance from Moresco  _et al._
(2020) we consider

$\mu^{\rm CCH}=\begin{pmatrix}H(z=0.122)\\\ H(z=0.44)\\\ H(z=0.60)\\\
H(z=0.73)\\\ H(z=0.32)\\\ H(z=0.57)\\\
H(z=1.48)\end{pmatrix}=\begin{pmatrix}74.91\\\ 86.79\\\ 92.77\\\ 97.63\\\
82.31\\\ 91.65\\\ 125.65\end{pmatrix}\,,\qquad C^{\rm
CCH}=\begin{pmatrix}24.54&16.79&12.91&9.75&19.72&13.64&-8.51\\\
-&22.06&24.68&26.80&20.08&24.19&38.99\\\ \
-&-&30.63&35.41&20.25&29.52&63.10\\\ -&-&-&42.42&20.37&33.80&82.76\\\
-&-&-&-&19.97&20.21&21.01\\\ -&-&-&-&-&28.52&58.57\\\ -&-&-&-&-&-&197.21\\\
\end{pmatrix}\,.$ (35)

The data vectors and covariance matrices are expressed in km/s/Mpc and $({\rm
km/s/Mpc})^{2}$, respectively.

## References

  * Hubble (1929) E. Hubble, Proc. Nat. Acad. Sci. 15, 168 (1929).
  * Riess  _et al._ (2021a) A. G. Riess _et al._ , (2021a), arXiv:2112.04510 [astro-ph.CO] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), arXiv:1807.06209 [astro-ph.CO] .
  * Camarena and Marra (2020a) D. Camarena and V. Marra, Phys. Rev. Res. 2, 013028 (2020a), arXiv:1906.11814 [astro-ph.CO] .
  * Efstathiou (2021) G. Efstathiou, Mon. Not. Roy. Astron. Soc. 505, 3866 (2021), arXiv:2103.08723 [astro-ph.CO] .
  * Freedman  _et al._ (2019) W. L. Freedman _et al._ ,  (2019), 10.3847/1538-4357/ab2f73, arXiv:1907.05922 [astro-ph.CO] .
  * Freedman (2021) W. L. Freedman, Astrophys. J. 919, 16 (2021), arXiv:2106.15656 [astro-ph.CO] .
  * Efstathiou (2020) G. Efstathiou, (2020), arXiv:2007.10716 [astro-ph.CO] .
  * Mortsell  _et al._ (2021) E. Mortsell, A. Goobar, J. Johansson, and S. Dhawan, (2021), arXiv:2105.11461 [astro-ph.CO] .
  * Riess  _et al._ (2009) A. G. Riess _et al._ , Astrophys. J. 699, 539 (2009), arXiv:0905.0695 [astro-ph.CO] .
  * Riess  _et al._ (2016) A. G. Riess _et al._ , Astrophys. J. 826, 56 (2016), arXiv:1604.01424 [astro-ph.CO] .
  * Riess  _et al._ (2021b) A. G. Riess, S. Casertano, W. Yuan, J. B. Bowers, L. Macri, J. C. Zinn, and D. Scolnic, Astrophys. J. Lett. 908, L6 (2021b), arXiv:2012.08534 [astro-ph.CO] .
  * Hinshaw  _et al._ (2013) G. Hinshaw _et al._ (WMAP), Astrophys. J. Suppl. 208, 19 (2013), arXiv:1212.5226 [astro-ph.CO] .
  * Ade  _et al._ (2016) P. A. R. Ade _et al._ (Planck), Astron. Astrophys. 594, A13 (2016), arXiv:1502.01589 [astro-ph.CO] .
  * Aylor  _et al._ (2017) K. Aylor _et al._ (SPT), Astrophys. J. 850, 101 (2017), arXiv:1706.10286 [astro-ph.CO] .
  * Aiola  _et al._ (2020) S. Aiola _et al._ (ACT), JCAP 12, 047 (2020), arXiv:2007.07288 [astro-ph.CO] .
  * Addison  _et al._ (2018) G. E. Addison, D. J. Watts, C. L. Bennett, M. Halpern, G. Hinshaw, and J. L. Weiland, Astrophys. J. 853, 119 (2018), arXiv:1707.06547 [astro-ph.CO] .
  * Cuceu  _et al._ (2019) A. Cuceu, J. Farr, P. Lemos, and A. Font-Ribera, JCAP 10, 044 (2019), arXiv:1906.11628 [astro-ph.CO] .
  * Aubourg  _et al._ (2015) E. Aubourg _et al._ , Phys. Rev. D 92, 123516 (2015), arXiv:1411.1074 [astro-ph.CO] .
  * Cuesta  _et al._ (2015) A. J. Cuesta, L. Verde, A. Riess, and R. Jiménez, Mon. Not. Roy. Astron. Soc. 448, 3463 (2015), arXiv:1411.1094 [astro-ph.CO] .
  * Feeney  _et al._ (2019) S. M. Feeney, H. V. Peiris, A. R. Williamson, S. M. Nissanke, D. J. Mortlock, J. Alsing, and D. Scolnic, Phys. Rev. Lett. 122, 061105 (2019), arXiv:1802.03404 [astro-ph.CO] .
  * Camarena and Marra (2020b) D. Camarena and V. Marra, Mon. Not. Roy. Astron. Soc. 495, 2630 (2020b), arXiv:1910.14125 [astro-ph.CO] .
  * Abbott  _et al._ (2017) B. P. Abbott _et al._ (LIGO Scientific, Virgo, 1M2H, Dark Energy Camera GW-E, DES, DLT40, Las Cumbres Observatory, VINROUGE, MASTER), Nature 551, 85 (2017), arXiv:1710.05835 [astro-ph.CO] .
  * Palmese  _et al._ (2021) A. Palmese, C. R. Bom, S. Mucesh, and W. G. Hartley, (2021), arXiv:2111.06445 [astro-ph.CO] .
  * Denzel  _et al._ (2021) P. Denzel, J. P. Coles, P. Saha, and L. L. R. Williams, Mon. Not. Roy. Astron. Soc. 501, 784 (2021), arXiv:2007.14398 [astro-ph.CO] .
  * Poulin  _et al._ (2019) V. Poulin, T. L. Smith, T. Karwal, and M. Kamionkowski, Phys. Rev. Lett. 122, 221301 (2019), arXiv:1811.04083 [astro-ph.CO] .
  * Niedermann and Sloth (2021) F. Niedermann and M. S. Sloth, Phys. Rev. D 103, L041303 (2021), arXiv:1910.10739 [astro-ph.CO] .
  * Agrawal  _et al._ (2019) P. Agrawal, F.-Y. Cyr-Racine, D. Pinner, and L. Randall, (2019), arXiv:1904.01016 [astro-ph.CO] .
  * Gómez-Valent  _et al._ (2021) A. Gómez-Valent, Z. Zheng, L. Amendola, V. Pettorino, and C. Wetterich, Phys. Rev. D 104, 083536 (2021), arXiv:2107.11065 [astro-ph.CO] .
  * Solà Peracaula  _et al._ (2019) J. Solà Peracaula, A. Gómez-Valent, J. de Cruz Pérez, and C. Moreno-Pulido, Astrophys. J. Lett. 886, L6 (2019), arXiv:1909.02554 [astro-ph.CO] .
  * Solà Peracaula  _et al._ (2020) J. Solà Peracaula, A. Gómez-Valent, J. de Cruz Pérez, and C. Moreno-Pulido, Class. Quant. Grav. 37, 245003 (2020), arXiv:2006.04273 [astro-ph.CO] .
  * Braglia  _et al._ (2020) M. Braglia, M. Ballardini, W. T. Emond, F. Finelli, A. E. Gumrukcuoglu, K. Koyama, and D. Paoletti, Phys. Rev. D 102, 023529 (2020), arXiv:2004.11161 [astro-ph.CO] .
  * Braglia  _et al._ (2021) M. Braglia, M. Ballardini, F. Finelli, and K. Koyama, Phys. Rev. D 103, 043528 (2021), arXiv:2011.12934 [astro-ph.CO] .
  * Jedamzik and Pogosian (2020) K. Jedamzik and L. Pogosian, Phys. Rev. Lett. 125, 181302 (2020), arXiv:2004.09487 [astro-ph.CO] .
  * Liu  _et al._ (2020) M. Liu, Z. Huang, X. Luo, H. Miao, N. K. Singh, and L. Huang, Sci. China Phys. Mech. Astron. 63, 290405 (2020), arXiv:1912.00190 [astro-ph.CO] .
  * Sekiguchi and Takahashi (2021) T. Sekiguchi and T. Takahashi, Phys. Rev. D 103, 083507 (2021), arXiv:2007.03381 [astro-ph.CO] .
  * Solà Peracaula  _et al._ (2021) J. Solà Peracaula, A. Gómez-Valent, J. de Cruz Pérez, and C. Moreno-Pulido, EPL 134, 19001 (2021), arXiv:2102.12758 [astro-ph.CO] .
  * Verde  _et al._ (2017a) L. Verde, E. Bellini, C. Pigozzo, A. F. Heavens, and R. Jiménez, JCAP 04, 023 (2017a), arXiv:1611.00376 [astro-ph.CO] .
  * Benisty and Staicova (2021) D. Benisty and D. Staicova, Astron. Astrophys. 647, A38 (2021), arXiv:2009.10701 [astro-ph.CO] .
  * Bernal  _et al._ (2016) J. L. Bernal, L. Verde, and A. G. Riess, JCAP 10, 019 (2016), arXiv:1607.05617 [astro-ph.CO] .
  * Aylor  _et al._ (2019) K. Aylor, M. Joy, L. Knox, M. Millea, S. Raghunathan, and W. L. K. Wu, Astrophys. J. 874, 4 (2019), arXiv:1811.00537 [astro-ph.CO] .
  * Heavens  _et al._ (2014) A. Heavens, R. Jiménez, and L. Verde, Phys. Rev. Lett. 113, 241302 (2014), arXiv:1409.6217 [astro-ph.CO] .
  * Verde  _et al._ (2017b) L. Verde, J. L. Bernal, A. F. Heavens, and R. Jiménez, Mon. Not. Roy. Astron. Soc. 467, 731 (2017b), arXiv:1607.05297 [astro-ph.CO] .
  * Haridasu  _et al._ (2018) B. S. Haridasu, V. V. Luković, M. Moresco, and N. Vittorio, JCAP 10, 015 (2018), arXiv:1805.03595 [astro-ph.CO] .
  * Cai  _et al._ (2022) R.-G. Cai, Z.-K. Guo, S.-J. Wang, W.-W. Yu, and Y. Zhou, Phys. Rev. D 105, L021301 (2022), arXiv:2107.13286 [astro-ph.CO] .
  * L’Huillier and Shafieloo (2017) B. L’Huillier and A. Shafieloo, JCAP 01, 015 (2017), arXiv:1606.06832 [astro-ph.CO] .
  * Shafieloo  _et al._ (2018) A. Shafieloo, B. L’Huillier, and A. A. Starobinsky, Phys. Rev. D 98, 083526 (2018), arXiv:1804.04320 [astro-ph.CO] .
  * Lin and Ishak (2017a) W. Lin and M. Ishak, Phys. Rev. D 96, 023532 (2017a), arXiv:1705.05303 [astro-ph.CO] .
  * Lin and Ishak (2017b) W. Lin and M. Ishak, Phys. Rev. D 96, 083532 (2017b), arXiv:1708.09813 [astro-ph.CO] .
  * Garcia-Quintero  _et al._ (2019) C. Garcia-Quintero, M. Ishak, L. Fox, and W. Lin, Phys. Rev. D 100, 123538 (2019), arXiv:1910.01608 [astro-ph.CO] .
  * Gonzalez  _et al._ (2021) J. E. Gonzalez, M. Benetti, R. von Marttens, and J. Alcaniz, JCAP 11, 060 (2021), arXiv:2104.13455 [astro-ph.CO] .
  * Scolnic  _et al._ (2018) D. M. Scolnic _et al._ , Astrophys. J. 859, 101 (2018), arXiv:1710.00845 [astro-ph.CO] .
  * Cole  _et al._ (2005) S. Cole _et al._ (2dFGRS), Mon. Not. Roy. Astron. Soc. 362, 505 (2005), arXiv:astro-ph/0501174 [astro-ph] .
  * Eisenstein  _et al._ (2005) D. J. Eisenstein _et al._ (SDSS), Astrophys. J. 633, 560 (2005), arXiv:astro-ph/0501171 [astro-ph] .
  * Bernal  _et al._ (2020) J. L. Bernal, T. L. Smith, K. K. Boddy, and M. Kamionkowski, Phys. Rev. D 102, 123515 (2020), arXiv:2004.07263 [astro-ph.CO] .
  * Etherington (1933) I. Etherington, Philos. Mag. 15, 761 (1933).
  * Carter  _et al._ (2018) P. Carter, F. Beutler, W. J. Percival, C. Blake, J. Koda, and A. J. Ross, Mon. Not. Roy. Astron. Soc. 481, 2371 (2018), arXiv:1803.01746 [astro-ph.CO] .
  * Beutler  _et al._ (2011) F. Beutler, C. Blake, M. Colless, D. H. Jones, L. Staveley-Smith, L. Campbell, Q. Parker, W. Saunders, and F. Watson, Mon. Not. Roy. Astron. Soc. 416, 3017 (2011), arXiv:1106.3366 [astro-ph.CO] .
  * Ross  _et al._ (2015) A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden, and M. Manera, Mon. Not. Roy. Astron. Soc. 449, 835 (2015), arXiv:1409.3242 [astro-ph.CO] .
  * Gil-Marín  _et al._ (2017) H. Gil-Marín, W. J. Percival, L. Verde, J. R. Brownstein, C.-H. Chuang, F.-S. Kitaura, S. A. Rodríguez-Torres, and M. D. Olmstead, Mon. Not. Roy. Astron. Soc. 465, 1757 (2017), arXiv:1606.00439 [astro-ph.CO] .
  * Kazin  _et al._ (2014) E. A. Kazin _et al._ , Mon. Not. Roy. Astron. Soc. 441, 3524 (2014), arXiv:1401.0358 [astro-ph.CO] .
  * Abbott  _et al._ (2019) T. M. C. Abbott _et al._ (DES), Mon. Not. Roy. Astron. Soc. 483, 4866 (2019), arXiv:1712.06209 [astro-ph.CO] .
  * Abbott  _et al._ (2021) T. M. C. Abbott _et al._ (DES), (2021), arXiv:2107.04646 [astro-ph.CO] .
  * Neveux  _et al._ (2020) R. Neveux _et al._ , Mon. Not. Roy. Astron. Soc. 499, 210 (2020), arXiv:2007.08999 [astro-ph.CO] .
  * Hou  _et al._ (2020) J. Hou _et al._ , Mon. Not. Roy. Astron. Soc. 500, 1201 (2020), arXiv:2007.08998 [astro-ph.CO] .
  * du Mas des Bourboux  _et al._ (2020) H. du Mas des Bourboux _et al._ , Astrophys. J. 901, 153 (2020), arXiv:2007.08995 [astro-ph.CO] .
  * Jiménez and Loeb (2002) R. Jiménez and A. Loeb, Astrophys. J. 573, 37 (2002), arXiv:astro-ph/0106145 [astro-ph] .
  * Moresco  _et al._ (2020) M. Moresco, R. Jiménez, L. Verde, A. Cimatti, and L. Pozzetti, Astrophys. J. 898, 82 (2020), arXiv:2003.07362 [astro-ph.GA] .
  * López-Corredoira  _et al._ (2017) M. López-Corredoira, A. Vazdekis, C. M. Gutiérrez, and N. Castro-Rodríguez, Astron. Astrophys. 600, A91 (2017), arXiv:1702.00380 [astro-ph.GA] .
  * López-Corredoira and Vazdekis (2018) M. López-Corredoira and A. Vazdekis, Astron. Astrophys. 614, A127 (2018), arXiv:1802.09473 [astro-ph.CO] .
  * Moresco  _et al._ (2018) M. Moresco, R. Jiménez, L. Verde, L. Pozzetti, A. Cimatti, and A. Citro, Astrophys. J. 868, 84 (2018), arXiv:1804.05864 [astro-ph.CO] .
  * Jiménez  _et al._ (2003) R. Jiménez, L. Verde, T. Treu, and D. Stern, Astrophys. J. 593, 622 (2003), arXiv:astro-ph/0302560 [astro-ph] .
  * Stern  _et al._ (2010) D. Stern, R. Jiménez, L. Verde, M. Kamionkowski, and S. A. Stanford, JCAP 1002, 008 (2010), arXiv:0907.3149 [astro-ph.CO] .
  * Moresco  _et al._ (2012) M. Moresco _et al._ , JCAP 1208, 006 (2012), arXiv:1201.3609 [astro-ph.CO] .
  * Zhang  _et al._ (2014) C. Zhang, H. Zhang, S. Yuan, T.-J. Zhang, and Y.-C. Sun, Res. Astron. Astrophys. 14, 1221 (2014), arXiv:1207.4541 [astro-ph.CO] .
  * Moresco (2015) M. Moresco, Mon. Not. Roy. Astron. Soc. 450, L16 (2015), arXiv:1503.01116 [astro-ph.CO] .
  * Moresco  _et al._ (2016) M. Moresco, L. Pozzetti, A. Cimatti, R. Jiménez, C. Maraston, L. Verde, D. Thomas, A. Citro, R. Tojeiro, and D. Wilkinson, JCAP 1605, 014 (2016), arXiv:1601.01701 [astro-ph.CO] .
  * Ratsimbazafy  _et al._ (2017) A. L. Ratsimbazafy, S. I. Loubser, S. M. Crawford, C. M. Cress, B. A. Bassett, R. C. Nichol, and P. Väisänen, Mon. Not. Roy. Astron. Soc. 467, 3239 (2017), arXiv:1702.00418 [astro-ph.CO] .
  * Borghi  _et al._ (2021) N. Borghi, M. Moresco, and A. Cimatti, (2021), arXiv:2110.04304 [astro-ph.CO] .
  * Gómez-Valent (2019) A. Gómez-Valent, JCAP 05, 026 (2019), arXiv:1810.02278 [astro-ph.CO] .
  * Simon  _et al._ (2005) J. Simon, L. Verde, and R. Jiménez, Phys. Rev. D71, 123001 (2005), arXiv:astro-ph/0412269 [astro-ph] .
  * Kjerrgren and Mortsell (2021) A. A. Kjerrgren and E. Mortsell, (2021), arXiv:2106.11317 [astro-ph.CO] .
  * Mahalanobis (1936) P. Mahalanobis, Proceedings of the National Institute of Science of India 2, 49 (1936).
  * Handley (2021) W. Handley, Phys. Rev. D 103, L041301 (2021), arXiv:1908.09139 [astro-ph.CO] .
  * Di Valentino  _et al._ (2019) E. Di Valentino, A. Melchiorri, and J. Silk, Nature Astron. 4, 196 (2019), arXiv:1911.02087 [astro-ph.CO] .
  * Efstathiou and Gratton (2020) G. Efstathiou and S. Gratton, Mon. Not. Roy. Astron. Soc. 496, L91 (2020), arXiv:2002.06892 [astro-ph.CO] .
  * Metropolis  _et al._ (1953) N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, and E. Teller, Journal of Chemical Physics 21, 1087 (1953).
  * Hastings (1970) W. Hastings, Biometrika 57, 97 (1970).
  * Lewis (2019) A. Lewis,  (2019), arXiv:1910.13970 [astro-ph.IM] .
  * Renzi and Silvestri (2020) F. Renzi and A. Silvestri, (2020), arXiv:2011.10559 [astro-ph.CO] .
  * Macpherson and Heinesen (2021) H. J. Macpherson and A. Heinesen,  (2021), 10.1103/PhysRevD.104.109901, [Erratum: Phys.Rev.D 104, 109901 (2021)], arXiv:2103.11918 [astro-ph.CO] .
  * Krishnan  _et al._ (2021) C. Krishnan, R. Mohayaee, E. O. Colgáin, M. M. Sheikh-Jabbari, and L. Yin, Class. Quant. Grav. 38, 184001 (2021), arXiv:2105.09790 [astro-ph.CO] .
  * Rasmussen and Williams (2006) C. Rasmussen and C. Williams,  _Gaussian Processes for Machine Learning_ (MIT Press, Massachusetts, 2006).
  * Busti  _et al._ (2014) V. C. Busti, C. Clarkson, and M. Seikel, Mon. Not. Roy. Astron. Soc. 441, 11 (2014), arXiv:1402.5429 [astro-ph.CO] .
  * Verde  _et al._ (2014) L. Verde, P. Protopapas, and R. Jiménez, Phys. Dark Univ. 5-6, 307 (2014), arXiv:1403.2181 [astro-ph.CO] .
  * Yu  _et al._ (2018) H. Yu, B. Ratra, and F.-Y. Wang, Astrophys. J. 856, 3 (2018), arXiv:1711.03437 [astro-ph.CO] .
  * Gómez-Valent and Amendola (2018) A. Gómez-Valent and L. Amendola, JCAP 04, 051 (2018), arXiv:1802.01505 [astro-ph.CO] .