Euclid preparation.

As we process information projected into bins of redshift, we consider the galaxy two-point angular correlation function $w(\theta)$, defined as (e.g., Crocce et al. 2011)

with $P_{\ell}$ being the Legendre polynomial and $C(\ell)$ the angular power spectrum defined as

In Eq. (2), $k$ is the wavenumber and $\mathcal{P}_{\Phi}$ stands for the dimensionless power spectrum of the primordial curvature perturbations $\Phi({k})$ that we model with HMCode (Mead et al. 2021) as implemented in the Boltzmann code CAMB (Challinor & Lewis 2011). The term $\Delta_{\ell}$ is the sum of the transfer function of number counts for the galaxy density field $\Delta^{\rm D}_{\ell}$ and the linear RSD contribution $\Delta^{\rm RSD}_{\ell}$ (Chisari et al. 2019)

where the two terms are respectively defined as

and

where $p(z)$ is the normalized galaxy redshift distribution, $b(z)$ is the linear galaxy bias, $z_{\rm max}$ is the maximum redshift of the survey that we fix to 3 for practical purposes, $j_{\ell}$ is the spherical Bessel function of order $\ell$, $j_{\ell}^{\prime\prime}$ its second derivative, and $r(z)$ is the comoving radial distance. In Eqs. (4) and (5), the transfer function $T_{\mathrm{X}}(k,z)$ of a quantity $X$ is defined as the ratio $X(k,z)/\Phi(k)$ so that $T_{\delta}$ is the transfer function of the matter overdensity $\delta(k,z)$ and $T_{\theta}$ is the transfer function of the divergence of the comoving velocity field. Discussed in Lepori et al. (2022) regarding full-shape analyses of photometric galaxy clustering, we checked that the effect of magnification bias over $\alpha$ is smaller than $0.2\,\sigma$ in most redshift bins, which is why no $\Delta^{\rm M}_{\ell}(k)$ term as defined in Eq. (26) of Chisari et al. (2019) is included in the model. Other relativistic effects are completely sub-dominant and not included either (Alonso et al. 2015).

To accurately model large scales at $\ell<220$ needed for BAO analysis, non-Limber integrals are computed using the Fang–Krause–Eifler–MacCrann (FKEM) method described in Fang et al. (2020a), while the faster Limber approximation (Kaiser 1992) is used at $\ell\in[220,10000]$. We checked that the effect of the Limber approximation over $C(\ell)$ is smaller than 1% in this range of $\ell$. Throughout this work, the theoretical $w(\theta)$ is computed using the Core Cosmology Library (Chisari et al. 2019). We refer to $w_{\mathrm{fid}}(\theta)$ when the theoretical $w(\theta)$ is computed with the fiducial cosmological parameters defined in Sect. 4.1.

The two-point angular correlation function is computed from the Flagship Mock Galaxy Catalogue (see Sect. 4) with the Landy–Szalay estimator (Landy & Szalay 1993)

where $N_{\scriptscriptstyle\rm DD}$, $N_{\scriptscriptstyle\rm DR}$, and $N_{\scriptscriptstyle\rm RR}$  are the pair counts where D stands for data and R for a random point. The random catalogues are created by sampling the footprint of the Flagship simulation defined by a HEALPix mask of $N_{\mathrm{side}}=4096$, which is equivalent to an angular resolution of $0.\mathrm{°}014$. We use 30 times as many random points as galaxies, yielding $1.04\times 10^{9}$ points per redshift bin. The angular binning has a resolution of $0.\mathrm{°}1$ and spans between $\theta_{\rm min}=$$$ and $\theta_{\rm max}=\theta_{\scriptscriptstyle\rm BAO}+$$$ where $\theta_{\scriptscriptstyle\rm BAO}$ is

evaluated at the mean redshift of the bin $z_{\mathrm{eff}}$ for the cosmology of the simulation given in Sect. 4.1. The measurement of $w(\theta)$ is performed with the code TreeCorr (Jarvis et al. 2004) and errors are estimated by a jackknife resampling of $N_{\mathrm{patch}}=500$ patches of about $10\,\mathrm{deg}^{2}$ each.

Two approaches are considered to compute the covariance of the two-point angular correlation function. For individual bin analyses, in which one value of the $\alpha$ parameter is fitted for each redshift bin, we use the jackknife covariance matrix built from the data vector measured with TreeCorr while joint analyses use an analytical covariance computed with CosmoCov (Krause & Eifler 2017). The use of an analytical covariance is made necessary by the noise in the covariance of the 13 redshift bins, even with 500 jackknife patches. Increasing to larger number of patches yields very little improvement in that regard.

The jackknife covariance matrix is computed with

where $n$ is the index of the jackknife realization, $w_{ij}\,(\theta)$ is the correlation of redshift bins $i$ and $j$ at angular separation $\theta$, and $\overline{w}_{ij}$ stands for the mean over all jackknife realizations for each angular separation $\overline{w}_{ij}\,(\theta)=\frac{1}{N_{\mathrm{patch}}}\sum_{n=1}^{N_{\mathrm{patch}}}\,w_{ij}^{(n)}\,(\theta)$.

We checked that the jackknife covariance matrix of each individual redshift bin is numerically stable with conditioning numbers ranging from $2.3\times 10^{4}$ for the covariance of the last redshift bin to $2.2\times 10^{5}$ for the first bin. Multiplying each of the jackknife covariances by its inverse yields the identity matrix, as expected from a numerically stable matrix, with a ratio between the diagonal and off-diagonal terms larger than $10^{13}$. In the MCMC analyses of Sect. 5, the inverse of the covariance is corrected by the Hartlap multiplicative factor defined in Hartlap et al. (2007)

where $N_{\mathrm{b}}$ is the number of angular bins in the data vector.

The analytical covariance used in Sect. 5.1 is the sum of the Gaussian and super-sample contributions, leaving aside the connected non-Gaussian term from modes within the footprint. The Gaussian term is computed using CosmoCov as in Krause & Eifler (2017) and Fang et al. (2020b), with a correction for the survey footprint like in Troxel et al. (2018) and without considering the Limber approximation (Fang et al. 2020a)

with $C_{ij}$ the angular power spectrum of redshift bins $i$ and $j$, $\Omega_{\mathrm{\,s}}$ the survey area, $\bar{n}$ the effective number density of galaxies and $i$, $j$, $k$, and $l$ are the indices of the redshift bins. The two-point angular correlation function is then computed from angular power spectra with (Abbott et al. 2024)

with $\overline{\mathrm{P_{\ell}}}$ the Legendre polynomial averaged over angular bins

in which $x=\cos(\theta)$ with $\theta_{\rm min},\theta_{\rm max}$ the lower and upper limits of each angular bin.

As for the super-sample covariance (SSC) contribution to the covariance, it was computed following the fast approximation from Lacasa & Grain (2019) extended to partial-sky in Gouyou Beauchamps et al. (2022) to go beyond the full-sky approximation by taking into account the footprint of the survey

where $\tilde{w}_{ij}\,(\theta)$ is computed as

with $R_{\ell}$ the response of the galaxy power spectrum to variations of the background density $\frac{\partial P_{\mathrm{gal}}(k,z)}{\partial\delta_{\mathrm{b}}}$. The $S_{ijkl}$ matrix element is computed using the implementation from PySSC ¹¹1https://github.com/fabienlacasa/PySSC as

where $W_{i}$ is the kernel of the observable in redshift bin $i$, the kernel being the normalized redshift distribution $p_{i}(z)$ in the case of photometric galaxy clustering. The comoving volume element ${\rm d}V=r^{2}(z)\,({\rm d}r/{\rm d}z)\,{\rm d}z$ is integrated between $z=0$ and $z=z_{\mathrm{max}}=3$. The variance of the background density $\sigma^{2}$ is, for a survey with a window function $\mathcal{W}$ of angular power spectrum $C^{\mathcal{W}}$

where the angular matter power spectrum between redshifts $z_{1}$ and $z_{2}$ is obtained from the 3D matter cross-spectrum $P_{\mathrm{m}}(k\,|\,z_{12})=D\left(z_{1}\right)\,D\left(z_{2}\right)\,P_{\mathrm{m}}(k,z=0)$ with

where $D(z)$ is the linear growth factor, $r_{i}=r(z_{i})$ is the comoving radial distance, $k_{\mathrm{min}}=0.1/r(z_{\mathrm{max}})$, and $k_{\mathrm{max}}=10/r(z_{\mathrm{min}})$.

We use the anafast routine from the HEALPix²²2http://healpix.sf.net library (Zonca et al. 2019; Górski et al. 2005) to compute $C^{\mathcal{W}}$. We make the approximation of a constant $R_{\ell}=5$, which reduces the convolution product to a multiplication. A detailed discussion on the response of the SSC can be found in Euclid Collaboration : Sciotti et al. (2024). The approximation on $R_{\ell}$ used in this work has an impact which does not exceed 0.5% on $\alpha$ and 3% on its uncertainty in all redshift bins, which is expected given the angular scales and redshifts considered.

A fit of the linear galaxy bias $b(z)$ is performed in each redshift bin using the jackknife covariance and the residuals $w(\theta)-b^{2}w_{\mathrm{fid,}b=1}(\theta)$, where $w_{\mathrm{fid,}b=1}$ denotes the theoretical angular correlation function defined in Eq. (1) computed with the Flagship cosmology and a galaxy bias $b=1$. We use scales between $0.\mathrm{°}5$ and $4\mathrm{°}$ in this full shape fit. A third order polynomial is then fitted to the result

which is then used in Eq. (4) to compute the theoretical model of the two-point angular correlation function.

To extract the BAO feature from the photometric sample, we perform the fitting between the measured $w(\theta)$ and a template derived from the theoretical two-point angular correlation function $w_{\mathrm{fid}}(\theta)$ computed for the fiducial cosmology described in Sect. 4.1. The template is defined as

where $\alpha$ is the transverse Alcock–Paczynski parameter quantifying an eventual shift of the BAO peak between the measured $w(\theta)$ and the fiducial $w_{\mathrm{fid}}(\theta)$. The nuisance parameters $B$, $A_{0}$, $A_{1}$, and $A_{2}$ are needed to absorb residual effects like non-linear galaxy bias. Different template parametrisations can be used and we will verify in Sect. 5.5 that the choice of polynomial has minimal impact on the $\alpha$ parameter. The parameter of interest in this work is $\alpha$ and, since the fiducial cosmology used to compute $w_{\mathrm{fid}}(\theta)$ is the same as the simulation, we expect to recover $\alpha=1$. On the contrary, using a different fiducial cosmology to compute $w_{\mathrm{fid}}(\theta)$ should result in $\alpha\neq 1$.

We will consider one set of nuisance parameters per redshift bin, since these parameters can in principle vary with redshift. This represents a total of 53 parameters for the joint analysis and 5 parameters for the analysis of each of the 13 bins. The Markov Chain Monte-Carlo (MCMC) technique is used to quantify the uncertainty on $\alpha$ marginalised over the nuisance parameters. The emcee sampler introduced in Foreman-Mackey et al. (2013) is used with the Gelman–Rubin convergence stopping criterion described in Gelman & Rubin (1992) with a threshold $R_{\mathrm{GR}}-1=0.005$ for analyses on individual redshift bins and $R_{\mathrm{GR}}-1=0.02$ for joint analyses. These thresholds were chosen to stop chains when parameter values and uncertainties reached a plateau. Uniform priors applied to the template parameters are presented in Table 1.

As a comparison to MCMC, we also consider the frequentist approach of profile likelihood to provide constraints in the joint analysis. We obtain the profile likelihood $\Delta\chi^{2}(\alpha)=\chi^{2}(\alpha)-\left.\chi^{2}(\alpha)\right|_{\mathrm{min}}$ by computing the best fit $\chi^{2}$ for each value of $\alpha$ between 0.8 and 1.2 in steps of 0.001. We then fit an 8th-order polynomial to this profile to compute $\Delta\chi^{2}$. The $1\,\sigma$ uncertainty on $\alpha$ is then given by $\Delta\chi^{2}=1$. In the joint analysis, the resolution of the grid of $\alpha$ on which this polynomial is evaluated is increased to use steps of 0.0001 to match the increased constraining power. We also use this frequentist approach to quantify the significance of the BAO detection as described in Sect. 5.1. The code iminuit based on the MINUIT algorithm is used at this effect (Dembinski et al. 2020; James & Roos 1975).