# Baryon acoustic oscillations in thin redshift shells from BOSS DR12 and
eBOSS DR16 galaxies

Ranier Menote1 and Valerio Marra,2,3,4
1PPGCosmo, Universidade Federal do Espírito Santo, 29075-910, Vitória, ES,
Brazil
2Núcleo de Astrofísica e Cosmologia & Departamento de Física, Universidade
Federal do Espírito Santo, 29075-910, Vitória, ES, Brazil
3INAF – Osservatorio Astronomico di Trieste, via Tiepolo 11, 34131, Trieste,
Italy
4IFPU – Institute for Fundamental Physics of the Universe, via Beirut 2,
34151, Trieste, Italy

###### Abstract

In an age of large astronomical datesets and severe cosmological tensions, the
case for model independent analyses is compelling. We present a set of 14
baryon acoustic oscillations measurements in thin redshift shells with $3\%$
precision that were obtained by analyzing BOSS DR12 and eBOSS DR16 galaxies in
the redshift range $0.32<z<0.66$. Thanks to the use of thin shells, the
analysis is carried out using just redshifts and angles so that the fiducial
model is only introduced when considering the mock catalogs, necessary for the
covariance matrix estimation and the pipeline validation. We compare our
measurements, with and without supernova data, to the corresponding
constraints from Planck 2018, finding good compatibility. A Monte Python
module for this likelihood is available at github.com/ranier137/angularBAO.

###### keywords:

large-scale structure of Universe – cosmology: observations – cosmological
parameters – cosmology: theory

††pagerange: Baryon acoustic oscillations in thin redshift shells from BOSS
DR12 and eBOSS DR16 galaxies–A

##  1 Introduction

Observational cosmology is being characterized by ever larger datasets, but
also ever worse discrepancies in the predictions of its standard model
(Perivolaropoulos & Skara, 2021), as remarked by the $5\sigma$ tension between
late and early determinations of the Hubble constant (Riess et al., 2021a).
This calls for model independent methods for analyzing observational data.

Model independent observational constraints are generally weaker than the ones
that adopt assumptions to reduce the degrees of freedom of a problem.
Sufficiently large dataset are thus necessary. Common hypotheses regard the
metric and content of the universe. Regarding the former, one usually assumes
large-scale homogeneity and isotropy, and so the FLRW metric. Regarding the
latter, raw data are usually analyzed withing the framework of the
$\Lambda$CDM model, often fixing its parameters to the fiducial standard model
values. This means that model independent methods allow one to test deviations
from the FLRW metric (Camarena et al., 2022) and the $\Lambda$CDM paradigm (Di
Valentino et al., 2021) in a more consistent way. In particular, model
independent analyses could be useful to understand the causes of the standard
model tensions, exactly because it is easier, with less assumptions, to narrow
down the problem.

Here, we present a set of 14 baryon acoustic oscillations (BAO) measurements
in thin redshift shells with $3\%$ precision that were obtained by analyzing
BOSS DR12 and eBOSS DR16 galaxies in the redshift range $0.32<z<0.66$. Thanks
to the use of thin shells (Sanchez et al., 2011), the analysis is carried out
using just redshifts and angles so that the fiducial model is only introduced
when considering the mock catalogs, necessary for the covariance matrix
estimation and the pipeline validation.111See Marra & Chirinos Isidro (2019)
for the estimation of the radial BAO signal using only angles and redshifts.

We compare these angular BAO measurements, with and without supernova data, to
the corresponding constraints from Planck 2018, finding good compatibility.
BAO measurements in thin redshift shells were also obtained by the
Observatório Nacional Group (ON Group, Carvalho et al., 2020, 2016; Alcaniz et
al., 2017), which, however, found a tension with respect to CMB results. We
also combine the angular BAO and supernova data with the SH0ES prior on $M$.
We find that it causes a strong $4.5\sigma$ tension in the $r_{d}$-$H_{0}$
plane with respect to Planck 2018: it is the $H_{0}$ crisis. A Monte Python
module for this likelihood is available at github.com/ranier137/angularBAO.

This paper is organized as follows. In Section 2 we describe the data, and in
Section 3 our methodology. Our set of 14 BAO measurements in thin redshift
shells is obtained in Section 4 and cosmological inference is carried out in
Section 5. We conclude in Section 6.

##  2 Data

The Sloan Digital Sky Survey (SDSS) is an international scientific
collaboration that has created the most detailed three-dimensional maps of the
Universe. This project was divided into 4 phases, SDSS-I (2000-2005), SDSS-II
(2005-2008), SDSS-III (2008-2014), and SDSS-IV (2014-2020). This work deals
with the last Data Release (DR) 12 of the Baryon Oscillation Spectroscopic
Survey (BOSS, Dawson et al., 2013) and the DR16 of the extended Baryon
Oscillation Spectroscopic Survey (eBOSS, Dawson et al., 2016), subsets of
SDSS-III and SDSS-IV, respectively. Fig. 1 shows their redshift distributions.
We considered the catalogs of luminous red galaxies (LRGs) that cover the
northen sky. We do not consider the southern sky because the smaller number of
galaxies does not allow for a robust determination of the BAO signal in thin
redshift slices. DR12-North offers a catalog with about 953 thousand LRGs with
redshift up to $z=0.8$, while DR16-North contains about 256 thousand LRGs
within the redshift range $0.6<z<1.0$.

In order to estimate the angular correlation function we use random catalogs
from BOSS and eBOSS about 50 and 20 times larger than the corresponding real
ones, respectively. In order to estimate the covariance matrix and the
robustness of the method we use 1000 mock catalogs from BOSS (Patchy-Mocks-
DR12NGC-COMPSAM V6C, Kitaura et al., 2016; Rodríguez-Torres et al., 2016) and
eBOSS (EZmock_eBOSS_LRGpCMASS_NGC, Zhao et al., 2021).

Figure 1: Redshift distributions of the luminous red galaxies from BOSS DR12
and eBOSS DR16. We considered the redshift range $0.32<z<0.66$.

##  3 Method

As said earlier, we wish to estimate the projected correlation function
$w(\theta)$ in thin redshift shells, which is related to the 3D correlation
function $\xi(s)$ via:

$\displaystyle w(\theta,z)$ $\displaystyle\xrightarrow{\text{thin-shell
limit}}\xi\big{(}s(\theta,z),z\big{)}\,,$ (1) $\displaystyle s(\theta,z)$
$\displaystyle=2\,(1+z)d_{A}(z)\sin\frac{\theta}{2}\,,$ (2)

so that the BAO peak is at the scale:

$\displaystyle\theta_{\rm BAO}(z)=\frac{r_{d}}{(1+z)d_{A}(z)}\,,$ (3)

where $r_{d}$ is the sound horizon at the drag epoch and $d_{A}$ is the
angular diameter distance. Our goal is to estimate $\theta_{\rm BAO}$.

Here, we focus on the BAO peak of the angular correlation function, that is,
we consider only the background expansion. Alternatively, one could use the
correlation function multipoles to constrain also the growth of structures,
see Wang et al. (2018). Furthermore, in order to avoid converting redshifts
and angles via a fiducial model, one could constrain the growth of the large-
scale structure via the angular power spectrum as proposed by Tanidis & Camera
(2021); Camera et al. (2018).

For later use, considering the late-time flat $\Lambda$CDM model, we list the
relevant background quantities:

$\displaystyle r(z)$ $\displaystyle=\int_{0}^{z}\frac{c\,{\rm
d}\bar{z}}{H(\bar{z})}\,,$ (4) $\displaystyle{H^{2}(z)\over H_{0}^{2}}$
$\displaystyle=\Omega_{m}(1+z)^{3}+1-\Omega_{m}\,,$ (5) $\displaystyle
d_{A}(z)$ $\displaystyle=(1+z)^{-1}r(z)\,,$ (6) $\displaystyle d_{L}(z)$
$\displaystyle=(1+z)\,r(z)\,,$ (7) $\displaystyle m(z)$
$\displaystyle=5\log_{10}\left[\frac{d_{L}(z)}{10\text{pc}}\right]+M\,,$ (8)
$\displaystyle\mu(z)$ $\displaystyle=m(z)-M\,,$ (9)

where we defined the comoving distance, the Hubble rate, the angular diameter
distance, the luminosity distance, the apparent magnitude and the distance
modulus, respectively. The parameter $M$ is the absolute magnitude of the
source.

###  3.1 Angular Correlation Function

The $i$-th galaxy of the catalog has redshift $z_{i}$ and angular position
$\beta_{i}=(\alpha_{i},\delta_{i})$, where $\alpha_{i}$ and $\delta_{i}$ are
its right ascension and declination, respectively. Considering two galaxies
within the redshift bin $[z-\frac{\delta z}{2},z+\frac{\delta z}{2}]$, their
angular separation $\gamma$ is given by:

$\displaystyle\gamma(\beta_{i},\beta_{j})=\arccos\Big{(}$
$\displaystyle\sin(\delta_{i})\sin(\delta_{j})+$
$\displaystyle\cos(\delta_{i})\cos(\delta_{j})\cos(\alpha_{i}-\alpha_{j})\Big{)}\,,$
(10)

with which one can estimate the two-point angular correlation function (ACF)
by counting pairs in various angular bins.

The catalogs give statistical weights associated to each galaxy in order to
correct for possible observational biases. Galaxies are then counted according
to the total weight (Anderson et al., 2013):

$\displaystyle w=w_{{}_{\rm FKP}}w_{\rm sys}\left(w_{\rm rf}+w_{\rm
cp}-1\right)\,,$ (11)

where $w_{\rm sys}=w_{\rm star}w_{\rm seeing}$ is the total systematic weight,
$w_{\rm rf}$ is the redshift failure weight, and $w_{\rm cp}$ is the close
pair weight. The FKP weights are meant to optimally weight the galaxies
(Feldman et al., 1994). In the present case of thin redshift shells they have
at most a 10% variation. Then, the total count of pairs separated by an
angular scale within the interval
$I_{k}\equiv[\theta_{k},\theta_{k}+\Delta\theta[$ is given by
$\sum_{i,j>i}\phi_{\theta_{k}}(\beta_{i},\beta_{j})$, where the counting
function is defined according to:

$\displaystyle\phi_{\theta}(\beta_{i},\beta_{j})=\begin{cases}w_{i}w_{j}&\theta\leq\gamma(\beta_{i},\beta_{j})<\theta+\Delta\theta\\\
0,&\text{otherwise}\end{cases}\,.$ (12)

Considering the real catalog, the normalized count of pairs of galaxies
relative to an angular separation within the interval $I_{k}$ is:

$\displaystyle DD(\theta_{k})\equiv
2\frac{\sum_{i,j>i}w_{i}w_{j}}{\left(\sum_{l}w_{l}\right)^{2}-\sum_{l}w_{l}^{2}}\text{\hskip
2.84526pt for \hskip 2.84526pt}\gamma_{ij}\in I_{k}\,,$ (13)

so that $\sum_{k}DD(\theta_{k})=1$.

In order to estimate the correlation function we adopt the optimal estimator
by Landy & Szalay (1993). Therefore, we need the normalized counts relative to
the random catalog:

$\displaystyle RR(\theta_{k})\equiv
2\frac{\sum_{i,j>i}w_{i}w_{j}}{\left(\sum_{l}w_{l}\right)^{2}-\sum_{l}w_{l}^{2}}\text{\hskip
2.84526pt for \hskip 2.84526pt}\gamma_{ij}\in I_{k}\,,$ (14)

and also the cross count between the real and random catalog:

$\displaystyle
DR(\theta_{k})=\frac{\sum_{i,j}w_{i}w_{j}}{\sum_{l,m}w_{l}w_{m}}\text{\hskip
2.84526pt for \hskip 2.84526pt}\gamma_{ij}\in I_{k}\,,$ (15)

where the $i$-galaxy belongs to the real catalog while the $j$-galaxy to the
random one. Finally, the angular correlation function is estimated according
to:

$\displaystyle
w(\theta_{k})=\frac{DD(\theta_{k})-2DR(\theta_{k})+RR(\theta_{k})}{RR(\theta_{k})}\,.$
(16)

###  3.2 Covariance matrix

The data analysis requires an accurate estimation of the covariance between
$w(\theta_{k})$ and $w(\theta_{k^{\prime}})$. We then use Eq. (16) to compute
the angular correlation function for the 1000 mock catalogs described in
Section 2. The covariance matrix is then estimated according to:

$\displaystyle
S_{kk^{\prime}}=\frac{1}{N\\!\\!-\\!\\!1}\\!\sum_{m=1}^{N}[w_{m}(\theta_{k})\\!-\\!\overline{w}(\theta_{k})][w_{m}(\theta_{k^{\prime}})\\!-\\!\overline{w}(\theta_{k^{\prime}})]\,,$
(17)

where $w_{m}(\theta_{k})$ is the correlation function for the $m$-th mock
catalog, $\overline{w}(\theta_{k})$ is the mean, and $N=1000$ is the number of
mock catalogs. As the estimation of the covariance matrix is based on mock
catalogs, it depends to some extent on the fiducial cosmology that was
adopted. This introduces a model-dependent component in the analysis.

To avoid using model-dependent mock catalogs, the covariance matrix could be
estimated by splitting the complete dataset into many sub-samples. While the
galaxy density from BOSS and eBOSS is not high enough to carry out such an
analysis, future datasets by DESI (Aghamousa et al., 2016)222www.desi.lbl.gov
and Euclid (Amendola et al., 2018)333www.euclid-ec.org could be suitable.
Alternatively, one could adopt the method developed by Angulo & White (2010);
Contreras et al. (2020) in order to change the cosmology of a halo sample and
study the impact of the choice of a fiducial model on cosmological inference.

###  3.3 Phenomenological model

Because of gravitational interaction, it is expected that the correlation of
galaxy positions is larger at smaller scales, then decreasing, like a power
law, towards larger scales, reaching eventually homogeneity ($w\approx 0$).
However, we should see, at intermediate scales, a sudden increase in the
correlation, marking the BAO feature. Following Sanchez et al. (2011), we use
the following phenomenological model in order to constrain the BAO peak
position directly from the data:

$w(\theta)=A+B\theta^{\nu}+Ce^{-\frac{(\theta-\theta_{\rm
fit})^{2}}{2\sigma^{2}}}\,.$ (18)

We consider angular bins with $\Delta\theta=0.3^{\circ}$, which are fine
enough so as to resolve the BAO feature, but wide enough in order to reduce
shot noise. We compute the correlation functions in the range
$1.5^{\circ}<\theta<7^{\circ}$.

In the above parametrization, the only physically interesting parameters are
the BAO position $\theta_{\rm fit}$ and the BAO signal strength $C$. The other
parameters are nuisance parameters, to be marginalized over, that are
necessary to reach a satisfactory fit to the data and whose correlations with
$\theta_{\rm fit}$ will affect the final constraint on the BAO signal. We
tested more complicated templates for the overall shape of the correlation
function. We examined the double exponential template of Sanchez et al.
(2013), one with a double Gaussian and also one with a coupling between the
exponential and Gaussian terms. We found that data does not justify them, as
these templates do not provide better results relative to their increased
number of parameters.

###  3.4 Detection of the BAO signal and priors

The chi-squared function is defined as

$\displaystyle\chi^{2}=\sum_{ij}\left[w_{i}-w(\theta_{i})\right]S^{-1}_{ij}\left[w_{j}-w(\theta_{j})\right]\,,$
(19)

where $w_{i}$ are the measurements, $w(\theta_{i})$ is the model of Eq. (18)
and $S_{ij}$ is the covariance matrix of Eq. (17).

We adopt wide-enough flat priors on the parameters $A$, $B$ and $\nu$ that
describe the correlation function without BAO feature. Regarding the
parameters describing the BAO peak, we adopt a flat prior on $C$ and
$\theta_{\rm fit}$, but an informative flat prior on $\sigma$:
$0.15^{\circ}<\sigma<1^{\circ}$. The lower bound is justified by the binning
width that we adopted as the narrowest Gaussian peak we can detect satisfies
$4\sigma\sim 2\Delta\theta$, that is, the 4$\sigma$ extension of the Gaussian
would cover two angular bins (three points $\theta_{k}$). The upper limit
comes from theoretical expectations, fitting the template of Eq. (18) to Eq.
(1). Ideally, one should use a left-bounded prior for $\sigma$ but the number
of galaxies is not high enough to permit such an analysis. The angular scale
that maximizes the marginalized posterior on $\theta_{\rm fit}$ is denoted by
$\theta_{\rm mc}$.

As said earlier, $C$ parametrizes the BAO signal strength and it is expected
to be positive. Shot noise, however, could produce false signals of positive
or negative intensity. Therefore, we use the posterior on $C$ in order to
assess the detection of the BAO peak. More precisely, we define the detection
strength according to:

$\displaystyle S_{\rm BAO}=\int_{0}^{\infty}f(C)\,{\rm d}C\,,$ (20)

where $f(C)$ is the posterior on $C$. We consider only measurements that
satisfy $S_{\rm BAO}>0.95$. In other words, we consider only measurements for
which the null hypothesis that the BAO feature is absent is excluded at more
than 95% confidence level.

###  3.5 Biases in the estimation of the BAO peak

We now discuss the possible sources of bias that affect the determination of
the angular scale of the BAO peak. The major sources of bias are due to
projection effects, modeling and parametrization assumptions, and redshift
space distortions. We will discuss their magnitude and how we took them into
account in the next Sections.

####  3.5.1 Projection bias

An important bias that needs to be considered is the one caused by the finite
width of the redshift bins. The number of galaxies in the BOSS and eBOSS
catalogs is not sufficient in order to consider a redshift shell of negligible
thickness. Projection effects will then displace the BAO peak to smaller
angles. Sanchez et al. (2011) showed that this shift depends on the bin width,
its redshift and, to a smaller degree, on the underlying cosmological model.

One could then proceed in two ways. The first approach is to deliver the
$\theta_{\rm BAO}(z)$ value relative to the bin $\Delta z$ that was used. In
other words, when comparing data to model, one should correct the theoretical
BAO angle of Eq. (3) according to the projection effect relative to the bin
$\Delta z$, redshift and cosmology that are considered. This would require
estimating the correlation function from the power spectrum after applying the
relevant selection function for every point of the parameter space. The second
approach is to deliver the $\theta_{\rm BAO}(z)$ value relative to $\Delta
z=0$ so that cosmological inference can be carried out by simply considering
Eq. (3). We adopt the second approach as it offers a much simpler likelihood.

Here, we try to estimate the projection bias directly from the data, without
adopting a cosmological model. More precisely, for a given redshift $z$ we
estimate the correlation function for the following range of bin widths
$\delta z_{i}$:

$\displaystyle\\{\delta
z_{i}\\}\\!\\!=\\!\\!\\{0.005,0.0075,0.01,0.0125,0.015,0.0175,0.02\\},$ (21)

and study how $\theta_{\rm mc}$ depends on $\delta z$. Note that adjacent
values of $z$ are displaced by 0.02 in order to avoid overlap.

As the highest value of $\delta z$ is much smaller than unity, we can expand
the dependency of $\theta_{\rm mc}$ on $\delta z$ in a Taylor series up to
second order. Based on geometrical arguments,444Because of the finite
thickness of the bin, the projected BAO scale is proportional to
$\cos\epsilon\sim 1-\epsilon^{2}/2$, where $\epsilon\propto\delta z$ is the
angle with respect to the tangent to the redshift shell. also confirmed by the
results of Sanchez et al. (2011), the first derivative of this correction at
$\delta z=0$ is null so that the expansion has only two free parameters:

$\displaystyle\theta(z,\delta z)=\theta_{0}(z)+E(z)\delta z^{2}\,.$ (22)

In order to correct for projection, we need to determine the function $E(z)$.
We adopt a linear model:

$\displaystyle E(z)=E_{0}+E_{1}z\,,$ (23)

but we found that the data does not justify $E_{1}$, which we then set to
zero.

The function $\theta_{0}(z)$ is meant to absorb the cosmology dependence of
$\theta_{\rm BAO}(z)$. To this end we model it as a piece-wise function:

$\displaystyle\theta_{0}(z)=\begin{cases}\dots&\dots\\\ \theta_{0,j}&\text{ if
}z=z_{j}\\\ \dots&\dots\end{cases}\,.$ (24)

Each redshift is given its own parameter in order not to force an ansatz on
the redshift dependence of $\theta_{0}(z)$.

We then adopt the following chi-squared function:

$\displaystyle\chi^{2}(\\{\theta_{0,j}\\},E_{0})$
$\displaystyle=\sum_{j}\chi_{j}^{2}(\theta_{0,j},E_{0})\,,$ (25)
$\displaystyle\chi_{j}^{2}(\theta_{0,j},E_{0})$
$\displaystyle=\sum_{ik}\big{[}\theta_{{\rm
mc},ij}-\theta_{ji}\big{]}\Sigma_{j,ik}^{-1}\big{[}\theta_{{\rm
mc},kj}-\theta_{jk}\big{]}\,,$ (26)

where $\theta_{ji}=\theta(z_{j},\delta z_{i})$, $i$ labels the thickness
$\delta z_{i}$ and $j$ the redshift $z_{j}$. We include only measurements for
which $S_{\rm BAO}>0.95$ so that the covariance matrix $\Sigma$ ranges from
$7\times 7$ to a single entry. The covariance matrix is obtained via:

$\displaystyle\Sigma_{j}=\text{diag}\big{(}\\{\sigma_{{\rm
mc},ij}\\}\big{)}\,R_{{\rm mc},j}\,\text{diag}\big{(}\\{\sigma_{{\rm
mc},ij}\\}\big{)}\,,$ (27)

where $\sigma_{{\rm mc},ij}$ are the uncertainties from the posterior on
$\theta_{{\rm mc},ij}$ and $R_{{\rm mc},j}$ is the correlation matrix relative
to the $\\{\delta z_{i}\\}$. We estimate the correlation matrix $R_{{\rm
mc},j}$ by minimizing the $\chi^{2}$ of Eq. (19) for the 1000 mock catalogs,
keeping only the fits with $C>0$. Note that, for different redshift $z_{j}$,
the $\theta_{{\rm mc},ij}$ are independent because we adopted non-overlapping
redshift slices.

Summarizing, the constraints on $\theta_{\rm BAO}(z)$ relative to $\Delta z=0$
are given by the posteriors on $\theta_{0,j}$ that are obtained from the chi-
squared function of Eq. (25).

####  3.5.2 Inference bias

To test the validity of the template of Eq. (18) we fit it to the mean of the
correlation functions from the 1000 mocks. In this case we use the chi-squared
function of Eq. (19) with the covariance matrix relative to the mean, that is,
$S/1000$. In order to neglect projection effects we consider the thinnest
shell with $\delta z=0.005$. We found that across the various redshift bins
this inference bias $B$ is approximately constant and of magnitude555The same
result was found when considering the projection effect.:

$\displaystyle B=\frac{\theta_{\rm mc}}{\theta_{\rm fid}}-1=(3.2\pm 0.1)\%\,.$
(28)

We checked that this bias can be reduced with a more complicated model, which,
however, would not be justified by the actual data: the extra parameters would
lead to poor constraints on most of the parameters of the model.

Considering this bias, the final measurement of the BAO peak that we report in
this work is given by the expression:

$\displaystyle\theta_{\rm BAO}=\frac{\theta_{0}(z)}{1+B}$ (29)

The error on $B$ is included in the error budget.

####  3.5.3 Further bias and systematic errors

As pointed out by Sanchez et al. (2011), there are other sources of bias and
systematic errors that are relative to the approach adopted in this analysis.
The first one is associated to the choice of the angular interval that is used
to measure the angular correlation function. Moving, for example, the starting
point to smaller angles, the fitting procedure will be exposed to higher
correlation and this may shift the recovered position of the BAO peak. This
parametrization bias $\delta\theta_{\rm par}$ is around $1\%$.

Then, the nonlinear gravitational growth of structure, scale dependent and
non-local bias, and redshift-space distortions also introduce a systematic
bias $\delta\theta_{\rm rsd}$ in the determination of $\theta_{BAO}$, which is
also around 1% (Sanchez et al., 2011). See Prada et al. (2016) for a thorough
discussion of possible systematic shifts and damping in baryon acoustic
oscillations due to the effects above.

###  3.6 Total error budget

The total uncertainty in the determination of the BAO scale includes all the
sources of error that were previously discussed:

$\displaystyle\sigma_{\rm sys}$ $\displaystyle=\sqrt{\delta\theta_{\rm
par}^{2}+\delta\theta_{\rm rsd}^{2}+\delta\theta_{B}^{2}}\,,$ (30)
$\displaystyle\sigma_{\rm BAO}$ $\displaystyle=\sqrt{\sigma_{\rm
stat}^{2}+\sigma_{\rm sys}^{2}}\,,$ (31)

where $\sigma_{\rm stat}$ is the uncertainty from the posterior on
$\theta_{0}$, $\delta\theta_{\rm par}=\delta\theta_{\rm
rsd}=0.01\,\theta_{0}$, and $\delta\theta_{B}\simeq 0.001\,\theta_{0}$. The
full covariance matrix is given by:

$\displaystyle\Sigma_{\rm BAO}=\text{diag}\big{(}\\{\sigma_{{\rm
BAO},j}\\}\big{)}\,R_{\rm BAO}\,\text{diag}\big{(}\\{\sigma_{{\rm
BAO},j}\\}\big{)}\,,$ (32)

where $R_{\rm BAO}$ is the correlation matrix relative to the posterior on
$\\{\theta_{0,j}\\}$, see Eq. (25).

###  3.7 MCMC exploration

Summarizing, our analysis consists of two steps.666In Appendix A we discuss a
method with only one step. First, we estimate the constraints on
$\\{\theta_{{\rm mc},ij}\\}$ via Eq. (19) for all the redshift slices and bin
widths. As we will see, this amounts to a total of $17\times 7=119$ MCMC
evaluations. Second, we explore the posterior relative to Eq. (25) in order to
obtain the constraints on $\\{\theta_{0,j}\\}$ and $E_{0}$.

For the posterior exploration we use the numerical code emcee (Foreman-Mackey
et al., 2013) and for the visualization of the chains the code getdist (Lewis,
2019). The convergence of the chains is tested following the method by Goodman
& Weare (2010), which consists in estimating the effective number of
independent samples using the integrated autocorrelation time (which is the
mean number of steps needed before the chain forgets where it started). In the
Supplementary Materials we provide plots and convergence tests for all the
redshift slices. In the MCMC exploration, we start the chain from the best-fit
model in order to avoid the burn-in phase.

##  4 BAO measurements

$z$ | $\theta_{\rm BAO}$ | $\sigma_{\rm stat}$ | $\sigma_{\rm sys}$ | $\sigma_{\rm BAO}$
---|---|---|---|---
0.35 | 5.80 | 0.063 | 0.085 | 0.106
0.37 | 6.07 | 0.103 | 0.089 | 0.136
0.39 | 5.89 | 0.067 | 0.086 | 0.109
0.41 | 5.30 | 0.137 | 0.078 | 0.157
0.43 | 4.87 | 0.060 | 0.071 | 0.093
0.45 | 4.52 | 0.134 | 0.066 | 0.150
0.47 | 4.69 | 0.111 | 0.069 | 0.131
0.49 | 4.69 | 0.041 | 0.069 | 0.080
0.51 | 4.65 | 0.088 | 0.068 | 0.112
0.53 | 4.03 | 0.067 | 0.059 | 0.089
0.55 | 3.56 | 0.058 | 0.052 | 0.078
0.57 | 4.36 | 0.081 | 0.064 | 0.103
0.61 | 3.78 | 0.056 | 0.055 | 0.079
0.63 | 3.90 | 0.057 | 0.057 | 0.080
Table 1: This table summarizes our final BAO measurements from the luminous red galaxies of BOSS DR12 and eBOSS DR16. The analysis that led to these constraints did not assume a fiducial cosmological model, except for the part relative to the galaxy mocks. The errors are approximately 3%. The corresponding correlation matrix $R_{\rm BAO}$ is given in Table 2. | $\theta_{35}$ | $\theta_{37}$ | $\theta_{39}$ | $\theta_{41}$ | $\theta_{43}$ | $\theta_{45}$ | $\theta_{47}$ | $\theta_{49}$ | $\theta_{51}$ | $\theta_{53}$ | $\theta_{55}$ | $\theta_{57}$ | $\theta_{61}$ | $\theta_{63}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$\theta_{35}$ | 1 | 0.19 | 0.35 | 0.12 | 0.09 | 0.08 | 0.05 | 0.21 | 0.05 | 0.06 | 0.06 | 0.06 | 0.33 | 0.07
$\theta_{37}$ |  | 1 | 0.26 | 0.1 | 0.07 | 0.07 | 0.04 | 0.16 | 0.03 | 0.05 | 0.05 | 0.05 | 0.26 | 0.05
$\theta_{39}$ |  |  | 1 | 0.17 | 0.13 | 0.11 | 0.07 | 0.28 | 0.07 | 0.08 | 0.08 | 0.08 | 0.46 | 0.09
$\theta_{41}$ |  |  |  | 1 | 0.04 | 0.04 | 0.02 | 0.1 | 0.02 | 0.03 | 0.03 | 0.03 | 0.17 | 0.03
$\theta_{43}$ |  |  |  |  | 1 | 0.03 | 0.03 | 0.08 | 0.01 | 0.02 | 0.02 | 0.02 | 0.12 | 0.03
$\theta_{45}$ |  |  |  |  |  | 1 | 0.02 | 0.06 | 0.01 | 0.02 | 0.02 | 0.01 | 0.11 | 0.02
$\theta_{47}$ |  |  |  |  |  |  | 1 | 0.04 | 0.01 | 0.02 | 0.01 | 0.02 | 0.06 | 0.02
$\theta_{49}$ |  |  |  |  |  |  |  | 1 | 0.04 | 0.05 | 0.05 | 0.05 | 0.27 | 0.06
$\theta_{51}$ |  |  |  |  |  |  |  |  | 1 | 0.02 | 0.01 | 0.01 | 0.06 | 0.01
$\theta_{53}$ |  |  |  |  |  |  |  |  |  | 1 | 0.02 | 0.01 | 0.08 | 0.01
$\theta_{55}$ |  |  |  |  |  |  |  |  |  |  | 1 | 0.02 | 0.07 | 0.01
$\theta_{57}$ |  |  |  |  |  |  |  |  |  |  |  | 1 | 0.08 | 0.02
$\theta_{61}$ |  |  |  |  |  |  |  |  |  |  |  |  | 1 | 0.09
$\theta_{63}$ |  |  |  |  |  |  |  |  |  |  |  |  |  | 1
Table 2: Correlation matrix $R_{\rm BAO}$ relative to the 14 BAO measurements
of Table 1. The entry $\theta_{35}$ refers to the $\theta_{\rm BAO}$
estimation at $z=0.35$ and similarly for the other redshifts. Figure 2:  Fit
of the model of Eq. (23) via Eq. (25). Only the measurements with $S_{\rm
BAO}>0.95$ are considered, see Eq. (20).

We consider the redshift range $0.32<z<0.66$ that has the highest density of
galaxies, see Fig. 1. As we use bins with widths up to $\delta z=0.02$, we
consider the non-overlapping bin centers $0.33,0.37,\dots,0.65$. As said
earlier, we only consider measurements with $S_{\rm BAO}>0.95$. None of the
measurements relative to the shells at $z=0.33,0.59,0.65$ pass this quality
cut and are thus absent from the following results, leaving 14 BAO
measurements.

Figure 2 shows the best-fit de-projection model for the redshift slices and
bin widths that satisfy $S_{\rm BAO}>0.95$. We find $E_{0}=-281\pm 138$. Our
final BAO results, corrected for projection and modeling biases via Eq. (29),
with total uncertainties from Eq. (31) are given in Table 1. The correlation
matrix $R_{\rm BAO}$ is given in Table 2. The covariance matrix is built via
Eq. (32).

The evolution of our measurements with the redshift is shown in Fig. 3. Also
shown are the results that were obtained by the Observatório Nacional (ON)
Group (Carvalho et al., 2016; Alcaniz et al., 2017; Carvalho et al., 2020),
which used a similar approach to measure the angular BAO scale. The ON Group
found a tension with the results from the CMB; from Fig. 3 it seems that our
determinations are closer to the CMB expectation. We will discuss this in more
details in the next Section.

Figure 3:  Top panel: Redshift evolution of the BAO measurements of Table 1 as
compared with the evolution as predicted by the $\Lambda$CDM model with the
parameters estimated by Planck 2018 (Aghanim et al., 2018, Table 2, last
column). The gray band corresponds to the $1\sigma$ uncertainty on $r_{d}\,h$.
Bottom panel: as above, but normalized with respect to the $\Lambda$CDM/Planck
2018 expectation.

##  5 Cosmological analysis

We now use our new BAO measurements of Table 1 to constrain the flat
$\Lambda$CDM model, whose prediction for $\theta_{\rm BAO}$ is given by Eq.
(3). The chi-squared function is:

$\displaystyle\chi^{2}_{\rm BAO}$ $\displaystyle(H_{0},\Omega_{m},r_{d})$ (33)
$\displaystyle=\sum_{ij}[\theta_{\rm BAO}(z_{i})-\theta_{{\rm
BAO},i}]\Sigma^{-1}_{{\rm BAO},ij}[\theta_{\rm BAO}(z_{j})-\theta_{{\rm
BAO},j}].$

In order to obtain competitive constraints it is useful to include Type Ia
Supernovae (SNe) in the analysis. We adopt the Pantheon dataset, consisting of
1048 supernovae spanning the redshift range $0.01<z<2.3$ (Scolnic et al.,
2018). The corresponding chi-squared function is:

$\displaystyle\chi^{2}_{\rm
sne}(H_{0},\Omega_{m},M)\\!=\\!\\!\sum_{ij}[m_{i}\\!-\\!m(z_{i})]S^{-1}_{\text{sne},ij}[m_{j}\\!-\\!m(z_{j})],$
(34)

where the apparent magnitudes $m_{i}$, redshifts $z_{i}$ and covariance matrix
$S_{\rm sne}$ are from the Pantheon catalog (considering both statistical and
systematic errors).

Finally, we will also consider the local prior on the SNe absolute magnitude
$M$ relative to the $H_{0}$ determination by Riess et al. (2021b), as derived
by Camarena & Marra (2020, 2021):

$M_{\rm R21}=-19.244\pm 0.037\text{ mag}\,.$ (35)

The chi-squared function is:

$\displaystyle\chi^{2}_{\rm R21}(M)=\frac{(M-M_{\rm
R21})^{2}}{\sigma^{2}_{M}}\,.$ (36)

###  5.1 Angular BAO alone

First, we only consider the angular BAO estimates. Alone, they can only
constrain the combination $r_{d}h$, where $H_{0}=h\,H_{100}$ and
$H_{100}=100\frac{\rm Km/s}{\rm Mpc}$. Therefore, we adopt the following chi-
squared function:

$\displaystyle\tilde{\chi}^{2}_{\rm BAO}(r_{d}h,\Omega_{m})=\chi^{2}_{\rm
BAO}(H_{100},\Omega_{m},r_{d}h)\,.$ (37)

Figure 4 shows the constraints from our dataset, the ON Group, and Planck 2018
(Aghanim et al., 2018, Table 2, last column). Our constraints are well
compatible with the ones from Planck 2018, but in tension with the ones from
the ON Group.

Figure 4: Marginalized constraints (68% and 95% credible regions) from the BAO
measurements of Table 1 for the flat $\Lambda$CDM model using the chi-squared
function of Eq. (37). Also shown are the constraints from the ON Group
(Carvalho et al., 2020, 2016; Alcaniz et al., 2017) and Planck 2018 (Aghanim
et al., 2018).

###  5.2 Angular BAO + SNe

Next, we include SN data. Without the use of a prior on $M$ (see next
Section), the parameters $H_{0}$ and $M$ are degenerate. Therefore, it is
convenient to adopt the variable $\hat{M}\equiv M-5\log_{10}h$ so that we
adopt the following chi-squared function:

$\displaystyle\tilde{\chi}^{2}_{\rm sne}(\Omega_{m},\hat{M})=\chi^{2}_{\rm
sne}(H_{100},\Omega_{m},\hat{M})\,.$ (38)

The total chi-squared function will be then:

$\displaystyle\chi^{2}(r_{d}h,\Omega_{m},\hat{M})=\tilde{\chi}^{2}_{\rm
BAO}(r_{d}h,\Omega_{m})+\tilde{\chi}^{2}_{\rm sne}(\Omega_{m},\hat{M}).$ (39)

Figure 5 shows the corresponding constraints together with the ones from
Planck 2018. Also in this case we find good agreement.

Figure 5: Marginalized constraints (68% and 95% credible regions) from the BAO
measurements of Table 1 and the Pantheon supernovae (Scolnic et al., 2018) for
the flat $\Lambda$CDM model using the chi-squared function of Eq. (39). Also
shown are the constraints from Planck 2018 (Aghanim et al., 2018).

###  5.3 Angular BAO + SNe + SH0ES

Finally, we include the local prior on $M$ of Eq. (35). This breaks the
degeneracy between $M$ and $H_{0}$, and this breaks the degeneracy between
$H_{0}$ and $r_{d}$. The total chi-squared function is then:

$\displaystyle\chi^{2}(h,r_{d},\Omega_{m},M)=$ $\displaystyle\chi^{2}_{\rm
BAO}(H_{0},\Omega_{m},r_{d})$ (40) $\displaystyle+\chi^{2}_{\rm
sne}(H_{0},\Omega_{m},M)+\chi^{2}_{\rm R21}(M)\,.$

Figure 6 shows the corresponding constraints together with the ones from
Planck 2018. We see that the SH0ES prior on $M$ causes a strong tension in the
$r_{d}$-$H_{0}$ plane with respect to Planck 2018: it is the $H_{0}$ crisis.

Figure 6:  Marginalized constraints (68% and 95% credible regions) from the
BAO measurements of Table 1, the Pantheon supernovae and the local calibration
of the supernova absolute magnitude for the flat $\Lambda$CDM model using the
chi-squared function of Eq. (40). Also shown are the constraints from Planck
2018.

###  5.4 Tension between angular BAO and CMB

We show in Table 3 the tension between the constraints from our angular BAO
measurements and the ones from Planck 2018 (Aghanim et al., 2018, Table 2,
last column). In order to quantity the tension in the $r_{d}$-$H_{0}$ plane,
we adopt the index of inconsistency (IOI) (Lin & Ishak, 2017):

$\displaystyle\sqrt{2\text{IOI}}$ $\displaystyle\equiv\sqrt{\delta^{T}{(C_{\rm
BAO}+C_{\rm P18})}^{-1}\delta\,}\,,$ (41) $\displaystyle\delta$
$\displaystyle=\\{r_{d}^{\rm BAO}-r_{d}^{\rm P18},\;H_{0}^{\rm BAO}-H_{0}^{\rm
P18}\\}\,,$

where $C$ are the covariance matrices on $r_{d}$ and $H_{0}$ from the analysis
of Figure 6 and $\delta$ is the difference vector. Note that this estimator
assumes Gaussianity and that the posteriors on $r_{d}$ and $H_{0}$ are very
close to Gaussian.

We find a good agreement between the $r_{d}h$ determination from Planck and
the ones using our BAO measurements, alone or with supernova data. When we add
the local prior on the supernova magnitude we find instead a strong
disagreement at the $4.5\sigma$ level.

Analysis |  $r_{d}h$ [Mpc] |  | Tension
---|---|---|---
Planck 2018 | $99.23\pm 0.94$ |  | -
Angular BAO | $97.5\pm 3.6$ |  | $0.5\sigma$
Angular BAO+SNe | $101.2\pm 1.1$ |  | $1.3\sigma$
|  $r_{d}$ [Mpc] | $H_{0}\left[\frac{\rm Km/s}{\rm Mpc}\right]$ | Tension
Planck 2018 | $147.21\pm 0.23$ | $67.66\pm 0.42$ | -
Angular BAO+SNe+M | $137.6\pm 2.6$ | $73.5\pm 1.3$ | $4.5\sigma$
Table 3: Constraints on $r_{d}h$, $r_{d}$ and $H_{0}$ for the analyses carried
out in this paper. The tension in the $r_{d}$-$H_{0}$ plane is estimated via
the index of inconsistency of Eq. (41).

##  6 Conclusions

Motivated by the standard model tensions (for a recent review, see
Perivolaropoulos & Skara, 2021), we analyzed BOSS DR12 and eBOSS DR16 galaxies
in thin redshift shells, obtaining a set of 14 BAO measurements with $3\%$
precision. These results are weakly model dependent as the fiducial model is
only introduced when considering the mock catalogs, necessary for the
covariance matrix estimation and the pipeline validation.

We find good compatibility with Planck 2018, also when combining these angular
BAO measurements with supernova data. Once we add the SH0ES prior on $M$ we
find a strong $4.5\sigma$ tension in the $r_{d}$-$H_{0}$ plane with respect to
Planck 2018. This model-independent analysis further highlights the $H_{0}$
crisis. A Monte Python module for this likelihood is available at
github.com/ranier137/angularBAO.

BAO measurements in thin redshift shells were also obtained by the
Observatório Nacional Group, which, however, found a tension with respect to
CMB results. This raised questions regarding the consistency of CMB and BAO
measurements. According to our results, which are relative to the latest SDSS
catalogs, CMB and BAO constraints are compatible.

## Acknowledgements

It is a pleasure to thank Armando Bernui and Gabriela C. Carvalho for useful
discussions and comments. RM thanks FAPES for financial support. VM thanks
CNPq and FAPES for partial financial support. This project has received
funding from the European Union’s Horizon 2020 research and innovation
programme under the Marie Skłodowska-Curie grant agreement No 888258. This
work also made use of the Virgo Cluster at Cosmo-ufes/UFES, which is funded by
FAPES and administrated by Renan Alves de Oliveira.

Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the
Participating Institutions, the National Science Foundation, and the U.S.
Department of Energy Office of Science. The SDSS-III web site is sdss3.org.
SDSS-III is managed by the Astrophysical Research Consortium for the
Participating Institutions of the SDSS-III Collaboration including the
University of Arizona, the Brazilian Participation Group, Brookhaven National
Laboratory, Carnegie Mellon University, University of Florida, the French
Participation Group, the German Participation Group, Harvard University, the
Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA
Participation Group, Johns Hopkins University, Lawrence Berkeley National
Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for
Extraterrestrial Physics, New Mexico State University, New York University,
Ohio State University, Pennsylvania State University, University of
Portsmouth, Princeton University, the Spanish Participation Group, University
of Tokyo, University of Utah, Vanderbilt University, University of Virginia,
University of Washington, and Yale University.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

The massive production of all MultiDark-Patchy mocks for the BOSS Final Data
Release has been performed at the BSC Marenostrum supercomputer, the Hydra
cluster at the Instituto de Fısica Teorica UAM/CSIC, and NERSC at the Lawrence
Berkeley National Laboratory. We acknowledge support from the Spanish MICINNs
Consolider-Ingenio 2010 Programme under grant MultiDark CSD2009-00064, MINECO
Centro de Excelencia Severo Ochoa Programme under grant SEV- 2012-0249, and
grant AYA2014-60641-C2-1-P. The MultiDark-Patchy mocks was an effort led from
the IFT UAM-CSIC by F. Prada’s group (C.-H. Chuang, S. Rodriguez-Torres and C.
Scoccola) in collaboration with C. Zhao (Tsinghua U.), F.-S. Kitaura (AIP), A.
Klypin (NMSU), G. Yepes (UAM), and the BOSS galaxy clustering working group.

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Aghamousa et al. (2016) Aghamousa A., et al., 2016, [1611.00036].
  * Aghanim et al. (2018) Aghanim N., et al., 2018, [1807.06209].
  * Alcaniz et al. (2017) Alcaniz J. S., Carvalho G. C., Bernui A., Carvalho J. C., Benetti M., 2017, Fundam. Theor. Phys., 187, 11, [1611.08458].
  * Amendola et al. (2018) Amendola L., et al., 2018, Living Rev. Rel., 21, 2, [1606.00180].
  * Anderson et al. (2013) Anderson L., et al., 2013, Mon. Not. Roy. Astron. Soc., 427, 3435, [1203.6594].
  * Angulo & White (2010) Angulo R. E., White S. D. M., 2010, Mon. Not. Roy. Astron. Soc., 405, 143, [0912.4277].
  * Camarena & Marra (2020) Camarena D., Marra V., 2020, Phys. Rev. Res., 2, 013028, [1906.11814].
  * Camarena & Marra (2021) Camarena D., Marra V., 2021, Mon. Not. Roy. Astron. Soc., 504, 5164, [2101.08641].
  * Camarena et al. (2022) Camarena D., Marra V., Sakr Z., Clarkson C., 2022, Mon. Not. Roy. Astron. Soc., 509, 1291, [2107.02296].
  * Camera et al. (2018) Camera S., Fonseca J., Maartens R., Santos M. G., 2018, Mon. Not. Roy. Astron. Soc., 481, 1251, [1803.10773].
  * Carvalho et al. (2016) Carvalho G. C., Bernui A., Benetti M., Carvalho J. C., Alcaniz J. S., 2016, Phys. Rev. D, 93, 023530, [1507.08972].
  * Carvalho et al. (2020) Carvalho G. C., Bernui A., Benetti M., Carvalho J. C., de Carvalho E., Alcaniz J. S., 2020, Astropart. Phys., 119, 102432, [1709.00271].
  * Contreras et al. (2020) Contreras S., Angulo R. E., Zennaro M., Aricò G., Pellejero-Ibañez M., 2020, Mon. Not. Roy. Astron. Soc., 499, 4905, [2001.03176].
  * Dawson et al. (2013) Dawson K. S., et al., 2013, AJ, 145, 10, [1208.0022].
  * Dawson et al. (2016) Dawson K. S., et al., 2016, AJ, 151, 44, [1508.04473].
  * Di Valentino et al. (2021) Di Valentino E., et al., 2021, Class. Quant. Grav., 38, 153001, [2103.01183].
  * Feldman et al. (1994) Feldman H. A., Kaiser N., Peacock J. A., 1994, ApJ, 426, 23, [astro-ph/9304022].
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, Publ. Astron. Soc. Pac., 125, 306, [1202.3665].
  * Goodman & Weare (2010) Goodman J., Weare J., 2010, Communications in Applied Mathematics and Computational Science, 5, 65.
  * Kitaura et al. (2016) Kitaura F.-S., et al., 2016, MNRAS, 456, 4156, [1509.06400].
  * Landy & Szalay (1993) Landy S. D., Szalay A. S., 1993, ApJ, 412, 64\.
  * Lewis (2019) Lewis A., 2019, [1910.13970].
  * Lin & Ishak (2017) Lin W., Ishak M., 2017, Phys. Rev., D96, 023532, [1705.05303].
  * Marra & Chirinos Isidro (2019) Marra V., Chirinos Isidro E. G., 2019, Mon. Not. Roy. Astron. Soc., 487, 3419, [1808.10695].
  * Perivolaropoulos & Skara (2021) Perivolaropoulos L., Skara F., 2021, [2105.05208].
  * Prada et al. (2016) Prada F., Scóccola C. G., Chuang C.-H., Yepes G., Klypin A. A., Kitaura F.-S., Gottlöber S., Zhao C., 2016, Mon. Not. Roy. Astron. Soc., 458, 613, [1410.4684].
  * Riess et al. (2021a) Riess A. G., et al., 2021a, [2112.04510].
  * Riess et al. (2021b) Riess A. G., Casertano S., Yuan W., Bowers J. B., Macri L., Zinn J. C., Scolnic D., 2021b, Astrophys. J. Lett., 908, L6, [2012.08534].
  * Rodríguez-Torres et al. (2016) Rodríguez-Torres S. A., et al., 2016, MNRAS, 460, 1173, [1509.06404].
  * Sanchez et al. (2011) Sanchez E., et al., 2011, Mon. Not. Roy. Astron. Soc., 411, 277, [1006.3226].
  * Sanchez et al. (2013) Sanchez E., Alonso D., Sanchez F. J., Garcia-Bellido J., Sevilla I., 2013, Mon. Not. Roy. Astron. Soc., 434, 2008, [1210.6446].
  * Scolnic et al. (2018) Scolnic D., et al., 2018, Astrophys. J., 859, 101, [1710.00845].
  * Tanidis & Camera (2021) Tanidis K., Camera S., 2021, [2107.00026].
  * Wang et al. (2018) Wang Y., Zhao G.-B., Chuang C.-H., Pellejero-Ibanez M., Zhao C., Kitaura F.-S., Rodriguez-Torres S., 2018, Mon. Not. Roy. Astron. Soc., 481, 3160, [1709.05173].
  * Zhao et al. (2021) Zhao C., et al., 2021, Mon. Not. Roy. Astron. Soc., 503, 1149, [2007.08997].

##  Appendix A Alternative analysis

For completeness, we present an alternative analysis which is, however,
computationally challenging. Let us consider, for each redshift $z_{j}$, the
vectors $\tilde{\theta}_{i}$ and $\tilde{w}_{j,i}$ of $7\times 19=133$ data
points, the 7 correlation functions at the various $\delta z_{i}$, with the
corresponding 19 values of $\theta$. Let us compute the $133\times 133$
covariance matrix $Z_{j}$ from the mocks. The model be:

$\displaystyle w_{j}(\tilde{\theta})$
$\displaystyle=\begin{cases}\dots&\dots\\\ w_{ji}(\theta)&\text{ if }\delta
z=\delta z_{i}\\\ \dots&\dots\end{cases}\,,$ (42) $\displaystyle
w_{ji}(\theta)$
$\displaystyle=A_{ji}+B_{ji}\theta^{\nu_{ji}}+C_{ji}e^{-\frac{(\theta-\Theta_{ji})^{2}}{2\sigma_{ji}^{2}}}\,,$
(43) $\displaystyle\Theta_{ji}$ $\displaystyle=\theta_{0,j}+E(z_{j})\delta
z_{i}^{2}\,,$ (44) $\displaystyle E(z_{j})$
$\displaystyle=E_{0}+E_{1}z_{j}\,.$ (45)

The chi-squared function is:

$\displaystyle\chi^{2}$ $\displaystyle=\sum_{j=1}^{15}\chi_{j}^{2}\,,$ (46)
$\displaystyle\chi_{j}^{2}$
$\displaystyle=\sum_{i=1,k=1}^{133,133}\left[\tilde{w}_{j,i}-w_{j}(\tilde{\theta}_{i})\right]Z^{-1}_{j,ik}\left[\tilde{w}_{j,k}-w_{j}(\tilde{\theta}_{k})\right]\,.$
(47)

In this way, one constrains the 17 BAO values $\theta_{0,j}$ via one single
analysis. One can then apply quality cuts via $S_{{\rm BAO},ji}$. The drawback
is that it involves $(5\times 7\times 17+2)=597$ nuisance parameters and,
therefore, a 614-dimensional parameter space.