© The Authors 2026

1. Introduction

The Euclid satellite will generate two main surveys. A survey of about 10⁷ galaxies with very precise spectroscopic redshifts and a second survey of about 1.5 × 10⁹ galaxies with less precise photometric redshifts, but which also contains galaxy shapes for a weak lensing shear analysis (Euclid Collaboration: Mellier et al. 2025). One of the main summary statistics of this second survey will be the 3 × 2 angular power spectra of the lensing shear, the galaxy number counts, and their cross-correlation (Euclid Collaboration: Blanchard 2020; Tutusaus et al. 2020).

In order to be able to exploit the data optimally, the calculation of the theoretical power spectra to be compared with the data must be both fast and accurate. For example, when performing a Markov chain Monte Carlo (MCMC) analysis for cosmological parameter inference, many theoretical spectra will need to be calculated, as quickly as possible, without biasing results, something that is currently not feasible to do via the full calculation. To date, mainly the Limber approximation (which can provide an increase in speed of the order of 100 times with respect to the full calculation in CLASS, for example) has been used to determine power spectra from photometric surveys (for a recent topical paper, see Leonard et al. 2023). While the Limber approximation is excellent when determining the power spectrum of the lensing potential via shear measurements (Kitching et al. 2017; Kilbinger et al. 2017), it is known to perform poorly for galaxies distributed narrowly in the radial direction (e.g. Simon 2007). For this paper our aim was to assess the particular impact of this failure on the results set to come from Euclid, both in the first and final data releases. This was done by determining the errors incurred, relative to the full calculation, via the Limber approximation and another approximation, the flat-sky approximation, which improves accuracy at large scales, while also speeding up the computation with respect to the full calculation if using an optimised code. While they do not use the flat sky approximation, BLAST (Beyond Limber Angular power Spectra Toolkit) from Chiarenza et al. (2024) and SwiftC_ℓ (Reymond et al. 2026) are examples of highly optimised codes for calculating the power spectrum.

The paper is structured as follows. In the next section we present a brief introduction to the flat-sky approximation and the Limber approximation and make contact with the relevant literature. In Sect. 3 we discuss the numerical implementation of the approximations used in this work, and in Sect. 4 we present results for the photometric survey of Euclid. In Sect. 5 we discuss our results.

1.1. Notation.

We calculated angular power spectra of different variables in different redshift bins. Power spectra of variables A and B at redshifts z and z′ calculated with the full expression (by using numerical codes such as CAMB Challinor & Lewis 2011 or CLASS Di Dio et al. 2013) are denoted C_ℓ^AB(z, z′) or, when integrated over the redshift bins i and j, they are denoted C_ℓ^AB(i, j). When neither A nor B is present in the superscript, C_ℓ(i, j) refers to the total observed power spectrum, including correlations of all necessary variables. The corresponding notations for the flat-sky approximation and Limber approximation will be ^FC_ℓ^AB(z, z′) and ^LC_ℓ^AB(z, z′) or, when integrated over redshift bins i and j, we denote them ^FC_ℓ^AB(i, j) and ^LC_ℓ^AB(i, j), respectively.

We work in a spatially flat background Friedmann-Lemaître universe with conformal time denoted η. As spatial curvature is certainly very small (see Planck Collaboration VI 2020) this is sufficient for our purpose. While it is fairly straightforward to extend the expressions for the angular power spectrum of galaxy number counts and its approximations to curved cosmologies, the expressions are more complicated and beyond the scope of our study. The background metric is then given by

$$\begin{array}{c}\hfill \mathrm{d}{s}^2=a^2(\eta )[-c^2\mathrm{d}{\eta }^2+\mathrm{d}\mathit{x}·\mathrm{d}\mathit{x}],\end{array}$$(1)

with c the vacuum speed of light. The scale factor is currently normalised to one: a(η₀)≡a₀ = 1. The Hubble factor is defined as $H=\frac{1}{a}\frac{\mathrm{d}a}{\mathrm{d}t}$, where t (dt = a dη) is the cosmic time.

2. The flat-sky and Limber approximations for the angular power spectrum of galaxy number counts

The flat-sky approximation has been investigated in the past (see Datta et al. 2007; White & Padmanabhan 2017; Castorina & White 2018; Jalilvand et al. 2020; Matthewson & Durrer 2021; Gao et al. 2024). While it has been found to be in excellent agreement with the full calculation of angular power spectra for very slim redshift bins (see Matthewson & Durrer 2021), an accurate extension to photometric bins has been developed in Gao et al. (2024). Here we present this extension (the flat-sky approximation) along with a second approach (the Limber approximation) as means to approximate the angular power spectrum for an arbitrary observable in the sky (see also Euclid Collaboration: Tessore et al. 2025). We applied both these approximations to the angular power spectrum of galaxy number counts.

We first considered an arbitrary cosmological observable A(x, η) at redshift z which is observed in the sky in direction $\widehat{\mathit{n}}$. If r(z) is the comoving distance out to redshift z, in terms of the variables $\widehat{\mathit{n}}$ and z, we have

$$\begin{array}{c}\hfill A(\widehat{\mathit{n}},z)=A(r(z)\widehat{\mathit{n}},{\eta }_0-r(z)/c),\end{array}$$(2)

where η₀ is the present time. This function is typically expanded on the sphere using spherical harmonics (if A is a complex spin-s field, like the shear, spin-weighted spherical harmonics are used),

$$\begin{array}{c}\hfill A(\widehat{\mathit{n}},z)={\displaystyle \sum _{\ell ,m}}{a}_{\ell m}(z)Y_{\ell m}(\widehat{\mathit{n}}).\end{array}$$(3)

Assuming statistical homogeneity and isotropy, only expectation values of equal ℓ and m do not vanish and their value depends only on ℓ. These define the angular power spectrum of two variables $A(\widehat{\mathit{n}},z)$ and $B(\widehat{\mathit{n}}\prime ,z\prime )$, i.e.

$$\begin{array}{c}\hfill ⟨a_{\ell m}^A(z)a_{{\ell }^{\prime }{m}^{\prime }}^{B^{\ast }}(z^{\prime })⟩={\delta }_{\ell {\ell }^{\prime }}^{\mathrm{K}}{\delta }_{m{m}^{\prime }}^{\mathrm{K}}{C}_{\ell }^{\mathit{AB}}(z,z^{\prime }).\end{array}$$(4)

For A = B, this is called the auto-correlation power spectrum, which is always positive definite, while for A ≠ B it is referred to as a cross-correlation spectrum. If $A(\widehat{\mathit{n}},z)$ is of the simple form of Eq. (2), it can be expressed its angular power spectrum in terms of the 3D power spectrum in Fourier space as (see e.g. Durrer 2020)

$$\begin{array}{c}\hfill C_{\ell }^{\mathit{AB}}(z,z^{\prime })=\frac{2}{\pi }{\displaystyle {\int }_0^{\infty }}\mathrm{dln}k{k}^3{P}_{\mathit{AB}}(k,z,z^{\prime })j_{\ell }(kr)j_{\ell }(k{r}^{\prime }).\end{array}$$(5)

Here, r and r′ are the comoving distances out to redshifts z and z′, respectively, and j_ℓ is the spherical Bessel function of order ℓ, whilst P_AB(k, z, z′) is the 3D unequal-time power spectrum of A and B, namely

$$\begin{array}{c}\hfill ⟨A(\mathit{k},z)B^{\ast }({\mathit{k}}^{\prime },z^{\prime })⟩={(2\pi )}^3{\delta }^{\mathrm{D}}(\mathit{k}-{\mathit{k}}^{\prime })P_{\mathit{AB}}(k,z,z^{\prime }).\end{array}$$(6)

If the variables A or B also have an intrinsic dependence on the direction $\widehat{\mathit{n}}$, the spherical Bessel functions in Eq. (5) may be replaced by more complicated expressions. For example, for redshift space distortions we find d²j_ℓ/dx² instead of j_ℓ(x) (see Bonvin & Durrer 2011, for details on all terms appearing in the number counts). Henceforth, we refer to various terms with A, B being replaced with the shorthand D for density, RSD for redshift space distortions, and L for lensing.

More realistically, especially for photometric surveys, we do not measure the C_ℓ at a precise redshift but integrated over redshift windows, w_i(z). Eq. (5) then becomes

$$\begin{array}{c}\hfill C_{\ell }^{\mathit{AB}}(i,j)={\displaystyle {\int }_0^{\infty }}\mathrm{d}z\mathrm{d}{z}^{\prime }{w}_i(z)w_j(z^{\prime })C_{\ell }^{\mathit{AB}}(z,z^{\prime }),\end{array}$$(7)

where, depending on the observable, j_ℓ in C_ℓ^AB(z, z′) may have to be replaced by its first or second derivative, as noted before in the case of redshift space distortions.

The Limber approximation now consists of replacing the integral over k in Eq. (5) by a Dirac-delta function such that (Limber 1954; Kaiser 1992)

$$\begin{array}{cc}\hfill C_{\ell }^{\mathit{AB}}(z,z^{\prime })& \simeq \frac{{\delta }^{\mathrm{D}}(r-r^{\prime })}{r^2}{P}_{\mathit{AB}}(\frac{\ell +1/2}{r},z,z)\hfill \\ \hfill & ={\delta }^{\mathrm{D}}(z-z^{\prime })\frac{H(z)}{c{r}^2}{P}_{\mathit{AB}}(\frac{\ell +1/2}{r},z,z),\hfill \end{array}$$(8)

hence

$$\begin{array}{c}\hfill {}^{\mathrm{L}}{C}_{\ell }^{\mathit{AB}}(i,j)\simeq {\displaystyle {\int }_0^{\infty }}\mathrm{d}z{w}_i(z)w_j(z)\frac{H(z)}{c{r}^2(z)}{P}_{\mathit{AB}}(\frac{\ell +1/2}{r},z,z).\end{array}$$(9)

In this way, the heavily oscillating integral over k and the double integral over z and z′ can be replaced by a single integral over z. If A = B, the integrand is even a manifestly non-negative function. This provides an enormous speed-up leading to excellent results for weak lensing which has a very broad window function at multipoles ℓ ≳ 20 (Kilbinger et al. 2017). However, it is well known that for galaxy number counts, especially in slim redshift bins, the Limber approximation differs by a factor 2 and more (depending on the bin width) from the true result for ℓ ≲ 100 (Di Dio et al. 2014, 2019; Fang et al. 2020; Martinelli et al. 2022). Here we compared it with the flat-sky approximation for the photometric survey of Euclid.

We now introduce the flat-sky approximation: it approximates the sky by a plane and the pair (ℓ,m) is replaced by a 2D vector ℓ. We assume that $\widehat{\mathit{n}}$ is close to some reference direction $\widehat{\mathit{e}}$, and we write $\widehat{\mathit{n}}=\widehat{\mathit{e}}+\mathit{\alpha }$. The amplitudes a_ℓm are then replaced by 2D Fourier transforms:

$$\begin{array}{cc}\hfill a^A(\mathit{\ell },z)& =\frac{1}{2\pi }{\displaystyle {\int }_{{\mathbb{R}}^2}}{\mathrm{d}}^2\alpha A(\mathit{\alpha },z){\mathrm{e}}^{-\mathrm{i}\mathit{\alpha }·\mathit{\ell }},\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill A(\mathit{\alpha },z)& =\frac{1}{2\pi }{\displaystyle {\int }_{{\mathbb{R}}^2}}{\mathrm{d}}^2\ell a^A(\mathit{\ell },z){\mathrm{e}}^{\mathrm{i}\mathit{\alpha }·\mathit{\ell }}.\hfill \end{array}$$(11)

Inserting the 3D Fourier representation for A, we find

$$\begin{array}{c}\hfill A(\mathit{\alpha },z)=\frac{1}{{(2\pi )}^3}{\displaystyle {\int }_{{\mathbb{R}}^3}}{\mathrm{d}}^{3}kA(\mathit{k},z){\mathrm{e}}^{\mathrm{i}\mathit{k}·(\widehat{\mathit{e}}+\mathit{\alpha })r(z)}.\end{array}$$(12)

Setting

$$\begin{array}{cc}\hfill \mathit{k}& =k_{\Vert }\widehat{\mathit{e}}+{\mathit{k}}_{\perp }\hfill \end{array}$$(13)

and comparing Eqs. (12) and (10), we obtain

$$\begin{array}{cc}\hfill a^A(\mathit{\ell },z)& =\frac{1}{{(2\pi )}^2}{\displaystyle {\int }_{{\mathbb{R}}^3}}{\mathrm{d}}^{3}kA(\mathit{k},z){\mathrm{e}}^{\mathrm{i}\mathit{k}·(\widehat{\mathit{e}}+\mathit{\alpha })r(z)}{\delta }^{\mathrm{D}}[\mathit{\ell }-r(z){\mathit{k}}_{\perp }]\hfill \\ \hfill & =\frac{1}{{(2\pi )}^2{r}^2(z)}{\displaystyle {\int }_{-\infty }^{+\infty }}\mathrm{d}{k}_{\Vert }A(\mathit{k}(k_{\Vert },\ell ,z),z){\mathrm{e}}^{\mathrm{i}{k}_{\Vert }r(z)},\hfill \end{array}$$(14)

where k(k_∥, ℓ, z) is given by Eq. (13) with k_⊥ = ℓ/r(z). This expression for k_⊥ is reminiscent of the Limber approximation, but overall the flat-sky approximation is distinct, since it retains the parallel component of k. If the redshifts are similar, we may replace z and z′ by $\overline{z}$, defined by $r(\overline{z})=\sqrt{r(z)r(z\prime )}$, in ℓ/r(z), such that

$$\begin{array}{c}\hfill ⟨a^A(\mathit{\ell },z)a^{B\ast }({\mathit{\ell }}^{\prime },z^{\prime })⟩\simeq \frac{1}{2\pi }\\ \hfill \times {\displaystyle {\int }_{{\mathbb{R}}^3}}{\mathrm{d}}^{3}k{P}_{\mathit{AB}}(k,\overline{z}){\mathrm{e}}^{\mathrm{i}{k}_{\Vert }(r-r^{\prime })}{\delta }^{\mathrm{D}}[\mathit{\ell }-r(\overline{z}){\mathit{k}}_{\perp }]{\delta }^{\mathrm{D}}[{\mathit{\ell }}^{\prime }-r(\overline{z}){\mathit{k}}_{\perp }]\\ \hfill ={\delta }^{\mathrm{D}}(\mathit{\ell }-{\mathit{\ell }}^{\prime })\frac{1}{2\pi r^2(\overline{z})}{\displaystyle {\int }_{-\infty }^{\infty }}\mathrm{d}{k}_{\Vert }{P}_{\mathit{AB}}(k,\overline{z},\overline{z}){\mathrm{e}}^{\mathrm{i}{k}_{\Vert }(r-r^{\prime })},\end{array}$$(15)

in which, the imaginary part vanishes and the real part, gives

$$\begin{array}{c}\hfill {}^{\mathrm{F}}{C}_{\ell }^{\mathit{AB}}(z,z^{\prime })=\frac{1}{\pi r^2(\overline{z})}{\displaystyle {\int }_0^{\infty }}\mathrm{d}{k}_{\Vert }{P}_{\mathit{AB}}(k,\overline{z},\overline{z})cos[k_{\Vert }(r-r^{\prime })].\end{array}$$(16)

It should be noted that we would not recover the physically relevant factor δ^D(ℓ − ℓ′), were it not for the fact that we had set $r=r(\overline{z})=r\prime$ in the expression for k_⊥. This is due to the fact that A(α, z) and A(α, z′) live on spheres with different radii, which modifies the relation between k_⊥ and ℓ. On the other hand, in the full calculation, see Eqs. (4) and (7), statistical isotropy ensures δ_ℓ ℓ′^K and δ_m m′^K. Fortunately, the oscillations of the cosine in Eq. (16) significantly suppress the signal for different redshifts so that this inaccuracy is irrelevant because of the shape of the power spectrum, which decreases in amplitude at small k, see Gao et al. (2023, 2024).

The power spectra from redshift bins with window functions w_i and w_j are thus given by

$$\begin{array}{cc}\hfill {}^{\mathrm{F}}{C}_{\ell }^{\mathit{AB}}(i,j)& \simeq {\displaystyle {\int }_0^{\infty }}\frac{\mathrm{d}z\mathrm{d}{z}^{\prime }}{\pi r^2(\overline{z})}{w}_i(z)w_j(z^{\prime })\hfill \\ \hfill & \times {\displaystyle {\int }_0^{\infty }}\mathrm{d}{k}_{\Vert }{P}_{\mathit{AB}}(k,\overline{z},\overline{z})cos[k_{\Vert }(r-r^{\prime })],\hfill \\ \hfill & ={\displaystyle {\int }_0^{\infty }}\frac{\mathrm{d}z\mathrm{d}{z}^{\prime }}{\pi r^2(\overline{z})}{w}_i(z)w_j(z^{\prime })\hfill \\ \hfill & \times {\displaystyle {\int }_{\ell /r(\overline{z})}^{\infty }}\mathrm{d}k{P}_{\mathit{AB}}(k,\overline{z},\overline{z})cos[k_{\Vert }(r-r^{\prime })],\hfill \end{array}$$(17)

where the change in integral bounds follows from the change of variable, given the relation between k and k_∥,

$$\begin{array}{c}\hfill k=\sqrt{k_{\Vert }^2+\frac{{\ell }^2}{r^2(\overline{z})}}·\end{array}$$(18)

In the literature, including in Gao et al. (2024), the significantly better ‘recalibrated’ flat-sky approximation is generally employed, where ℓ is replaced by ℓ + 1/2 in the above expression for k.

Even though the number of integrals over the power spectrum P_AB(k,z,z′), for the flat-sky expression in Eq. (18), is the same as in Eq. (7), for r = r′ (the value with the dominant contribution to $C_{\mathit{ij}}^{\mathit{AB}}$) the ℓ-dependence has simply been absorbed into the boundary of the integral, reducing the complexity of the integration. In addition, the integral no longer contains the heavily oscillating product of spherical Bessel functions. More importantly, as described in Gao et al. (2024), this form of the integral allows the k-integral to be precomputed and stored, which significantly saves computation time in applications where the angular power spectra from many combinations of different cosmological parameter values are required, such as parameter estimation and MCMC analysis. From the flat-sky approximation we may further obtain the Limber expression simply by setting k_∥ = 0 and replacing the k-integral by δ^D(r − r′) = δ^D(z − z′)H(z)/c.

We used this flat-sky expression to calculate the C_ℓ for galaxy number counts. In addition to density fluctuations, we also included redshift space distortions and relativistic effects, for example, the lensing magnification, which arise in observations made on the past light cone – see Yoo et al. (2009), Bonvin & Durrer (2011), and Challinor & Lewis (2011) for a full derivation of all terms in the exact solution. We refer to these as the contributions from observational effects. For the density we have

$$\begin{array}{c}\hfill A^{\mathrm{D}}(\mathit{k},z)=b(z)D(\mathit{k},z),\end{array}$$(19)

where b(z) is the linear galaxy bias, that is determined via simulations (see Sect. 3.2.1), and D(k, z) is the matter density contrast in Fourier space. The redshift space distortions are given by

$$\begin{array}{c}\hfill A^{\mathrm{RSD}}(\mathit{k},z)=\frac{k_{||}^2}{k\mathcal{H}(z)}V(\mathit{k},z),\end{array}$$(20)

where k⁻¹ V(k, z) is the Fourier transform of the velocity potential, and $\mathcal{H}=\frac{1}{a}\frac{\mathrm{d}a}{\mathrm{d}\eta }=Ha$ is the conformal Hubble factor. The next important contribution to the number counts comes from the lensing magnification and is given by

$$\begin{array}{c}\hfill A^{\mathrm{L}}(\mathit{k},z,L)=(2-5s(z,L))\kappa ,\end{array}$$(21)

where κ is the lensing convergence and contains an integral over redshift, which is well approximated by the Limber approximation. The function s(z, L) is the logarithmic slope of the galaxy number density, n(z, L), as a function of luminosity, L,

$$\begin{array}{c}\hfill s(z,L)=\frac{2}{5}\frac{\partial lnn(z,L)}{\partial lnL}·\end{array}$$(22)

We used the result obtained by the Flagship simulation (see Sect. 3.2.1; for a more detailed representation of all Limber and flat-sky expressions and their cross-correlations, see Matthewson & Durrer 2021). In this analysis we included the appropriate terms for these three contributions in each of the approximations of the full angular power spectra. The remaining relativistic contributions in Bonvin & Durrer (2011) are negligible and, therefore, set to zero in both approximations.

3. Implementation, performance, and simulations

In this section we first describe the implementation of our power spectra computations under different approximations. We then clarify the methodology that has been used to compare them, and finally present the different Euclid scenarios that were studied in this work.

3.1. Computation of the power spectra

As laid out in the previous section, we considered three alternative methods to compute the power spectra of galaxy number counts; that is, the full computation (Eq. (7)), the flat-sky approximation (Eq. (18)), and the Limber approximation (Eq. (9)). Starting with the full computation case, we calculated the galaxy number count angular power spectra with the CLASS Boltzmann solver (Blas et al. 2011) to linear order, including all contributions and relativistic terms. For the flat-sky approximation, we made use of the GZCphysics code described in Gao et al. (2024), further modified by us to use CLASS instead of the Boltzmann solver CAMB. We checked that the flat-sky results from the above paper are in agreement with our analogous calculations. We then extended the code to include lensing magnification, and to accept arbitrary (not just Gaussian) windows in redshift, needed for the application to Euclid. Following this, we implemented the Euclid photometric survey’s simulated tomographic bins. Finally, we computed the power spectra with the Limber approximation using CAMB via the CosmoSIS framework (Zuntz et al. 2015). This choice was made to facilitate the control of the scale at which the Limber approximation is used in the computation which is not straightforward when simply using the Limber implementation available in CLASS or CAMB.

Our main goal in this analysis is to compare the performance of the flat-sky approximation and Limber approximation with respect to the full computation. Given that the differences between the approximations and the full result are only relevant at low multipoles, we restricted our analysis to the largest scales (ℓ ≤ 300) and considered a linear matter power spectrum. We note that even if there is some leakage from higher k modes due to integrated terms like lensing magnification, when applying a multipole cut, there is no need to consider accurate recipes to account for the non-linearities in the matter power spectrum at the scales relevant for this analysis.

3.2. Euclid specifications

In this work we extracted the specifications that correspond to the Euclid photometric survey from the Flagship simulation for two different scenarios: the first data release (DR1) and the third and final data release (DR3). Below, both the simulation and the data releases are briefly described.

3.2.1. The Flagship simulation

The Flagship galaxy mock (Euclid Collaboration: Castander et al. 2025) is a simulated catalogue of galaxies built from one of the largest N-body dark matter simulations ever run, with a box size of 3600 h⁻¹ Mpc and a particle mass resolution of 10⁹ h⁻¹ M_⊙. The dark matter simulation was run with PKDGRAV3 (Potter et al. 2017), and haloes were identified with the rockstar algorithm (Behroozi et al. 2013). Haloes were then populated with galaxies using a combination of the halo occupation distribution and sub-halo abundance matching techniques (Carretero et al. 2015). The galaxy mock was calibrated to reproduce several observations, including the luminosity function (Blanton et al. 2003, 2005a), the clustering of galaxies (Zehavi et al. 2011), and the colour-distribution as a function of the absolute magnitude in the r band (Blanton et al. 2005b) of Sloan Digital Sky Survey data.

The cosmological parameters for this simulation are given in Table 1 along with the ones used in this paper, which are quite similar, but not exactly the same. Since the only simulation output used here was the galaxy distribution as a function of redshift in each tomographic bin, these small differences in cosmological parameters are not relevant. What is relevant is that we used the same galaxy redshift distribution, linear galaxy bias, and magnification bias for all three calculations. The catalogue can be accessed via the CosmoHub platform¹ (Tallada et al. 2020; Carretero et al. 2017).

3.2.2. DR1 settings

In order to study the performance of the flat-sky approximation and Limber approximation for Euclid, we first considered the initial data release. We followed the approach presented in Euclid Collaboration: Mellier et al. (2025) and selected a galaxy sample from the Flagship galaxy mock, selecting a conservative magnitude cut at I_E < 23.5. This choice of a fairly bright sample decreases the presence of systematic effects that introduce spurious clustering signals. Once the sample had been defined, we classified galaxies into six equi-populated bins, depending on the photometric redshift assigned to each galaxy. The galaxy distributions were then built by making a histogram using the bins in observed redshift associated to each object. We note that these were derived assuming Legacy Survey of Space and Time (LSST)-like photometry, which will not be available for DR1 and, as such, these DR1 settings correspond to an optimistic scenario. The main specifications are summarised in Table 2. In this case, we measured the corresponding galaxy bias defined in Eq. (20), and the magnification bias defined in Eq. (23) for this sample and their values are provided in Table 3. The galaxy distributions are also shown in Fig. 1.

Fig. 1. Number density bins simulated for the photometric surveys of Euclid, in 6-bin and 13-bin equi-populated configurations. The numbering convention we used for the ith bin is shown in the coloured blocks above the maxima of each bin.

3.2.3. DR3 settings

In addition to the initial Euclid data release, we also considered the performance of the flat-sky approximation and Limber approximation for the full data sample of Euclid at the end of the mission. We considered 13 equi-populated bins with a higher (dimmer) magnitude cut at I_E < 24.5. The final sample contains a total number density of 24.3 galaxies per arcmin², and we considered the galaxy and magnification biases measured for this sample in Euclid Collaboration: Mellier et al. (2025). The main specifications are summarised in Table 2 and the galaxy distributions are shown in Fig. 1.

3.3. Assessment of performance

For both survey configurations considered, we wished to assess the accuracy of each of the approximations in recovering the full calculation result, taking into account the uncertainty on the measurement of the power spectra – which, crucially, is largest at small multipoles, and depends on the survey configuration under consideration.

To compare the accuracy of each approximation in the context of real measurements, we simulated noisy power spectra, based on the full CLASS spectra and the theoretical covariance matrix for the planned tomographic bins. In this way, we obtained a set of simulated mock spectra that encompasses the expected spread of measurements, subject to the relevant observational uncertainties (see Tables 2 and 3). Specifically, we took the values of the full spectrum for some subset of multipoles ℓ as the mean, around which we drew 1000 random noise realisations which are normally distributed with a Gaussian covariance (i.e. we neglect non-Gaussian contributions from the connected four-point function). We used the SpaceBorne code² (Euclid Collaboration: Sciotti et al. 2024) to estimate an analytic Gaussian covariance for the different spectra considered computed using the full spectrum calculation from CLASS. This Gaussian covariance takes the form

$$\begin{array}{cc}\hfill \mathrm{Cov}\left[C_{\ell }^{\mathit{GG}}(i,j)C_{{\ell }^{\prime }}^{G^{\prime }{G}^{\prime }}(k,l)\right]=\frac{1}{2\ell (\ell +1)f_{\mathrm{sky}}\mathrm{\Delta }\ell }{\delta }_{\ell {\ell }^{\prime }}^K& \\ \hfill \times [(C_{\ell }^{G{G}^{\prime }}(i,k)+N_{\ell }^{G{G}^{\prime }}(i,k))(C_{{\ell }^{\prime }}^{G{G}^{\prime }}(j,l)+N_{{\ell }^{\prime }}^{G{G}^{\prime }}(j,l))& \\ \hfill +(C_{\ell }^{G{G}^{\prime }}(i,l)+N_{\ell }^{G{G}^{\prime }}(i,l))(C_{{\ell }^{\prime }}^{G{G}^{\prime }}(j,k)+N_{{\ell }^{\prime }}^{G{G}^{\prime }}(j,k))],\end{array}$$(23)

see Equation (138) of Blanchard (2020), where f_sky is the fraction of the sky observed, Δℓ is the width of the multipole bins, δ_ℓℓ′^K is the Kronecker delta, G indicates the full observed galaxy clustering power spectrum, and indices i, j, k, l run over the tomographic bins. The noise term is given by the shot noise

$$\begin{array}{c}\hfill N_{\ell }^{\mathit{GG}}(i,j)=\frac{1}{{\overline{n}}_i}{\delta }_{\mathit{ij}}^K,\end{array}$$(24)

with ${\overline{n}}_i$ the galaxy surface density per bin in units of inverse steradians.

By comparing to each other the χ² distributions, formed from the residuals of each method’s result with respect to every one of the spectra in the set of mock observations, we could then obtain an estimate of how far each of the calculations might effectively lie from an observed spectrum, and the frequency with which this occurs.

The χ² for a particular approximation of the spectra ^XC_ℓ, X ∈ {L, F}, with respect to one of the noisy realisations, ${\widehat{C}}_{\ell }$, is computed as

$$\begin{array}{c}\hfill {\chi }_X^2={\displaystyle \sum _{n=1}^{12}}{\left({\mathit{r}}_n^X\right)}^{\mathsf{T}}{({\mathrm{\Sigma }}^{-1})}_n{\mathit{r}}_n^X,\end{array}$$(25)

where the sum is over values corresponding to 12 bins distributed evenly in the logarithmic space of ℓ, between ℓ_min = 2 and ℓ_max = 300. In Eq. (26), ${\mathit{r}}_n^X$ is the vector of the residuals,

$$\begin{array}{c}\hfill {\mathit{r}}_n^X{\mathrm{vech}}_{\mathit{ij}}\left[{\widehat{C}}_n(i,j){-}^X{C}_n(i,j)\right],\end{array}$$(26)

with vech_ij the half-vectorisation operator³, acting on redshift indices of the symmetric tomographic matrix, and ${\mathrm{\Sigma }}_n^{-1}$ is the inverse of the covariance matrix in a given multipole bin, ${\mathrm{\Sigma }}_n=⟨{\mathit{r}}_n^X{({\mathit{r}}_n^X)}^T⟩$.

The covariance (and thus also the mock realisations of the spectra) includes all possible redshift correlations and is calculated for both survey configurations described above. It is often useful to consider the χ² in the context of the number of degrees of freedom (DOF). In this case, there are N = 12 chosen covariance bins, which effectively act as our ‘data points’, and n redshift bins, leading to n(n + 1)/2 independent spectra of equal and unequal redshift correlations. Thus, we can expect DOF_DR1 = N × n × (n + 1)/2 = 252 for DR1 and DOF_DR3 = 1092 for DR3.

4. Results

We now present the main results of our analysis. The aim of this section is to quantify the performance of the flat-sky approximation and Limber approximation (for the latter, see also Euclid Collaboration: Cardone et al. 2025; Euclid Collaboration: Joudaki et al. 2025) in the context of Euclid (both for DR1 and DR3 set-ups), assessing how both perform compared to the exact solution.

4.1. Power spectrum recovery

In Figs. 2 and 3 we plot results for the tomographic angular power spectra for the first bin of the DR1 and DR3 configurations, respectively. The true result taken from a full-sky calculation of CLASS is displayed, along with the results using the Limber and flat-sky approximations. The relative error with respect to the full calculation at each multipole is shown in the lower panels.

Fig. 2. Upper panel: Equal redshift angular power spectrum for the lowest redshift bin of Euclid, using six equi-populated bins. We show the full calculation result from CLASS as a black solid line, and the Limber approximation result as a purple dashed line. The green dashed line corresponds to the recalibrated version of the flat-sky approximation. Lower panel: Relative error (in %) of each approximation from the full calculation result.

Fig. 3. Upper panel: Equal redshift angular power spectrum for the lowest redshift bin of Euclid, using 13 equi-populated bins. We show the full calculation result from CLASS as a black solid line, and the Limber approximation result as a purple dashed line. The green dashed line corresponds to the recalibrated version of the flat-sky approximation. Lower panel: Relative error (in %) of each approximation from the full calculation result.

We calculated only until ℓ = 300 since, on smaller scales than this, the Limber approximation performs sufficiently well to recover the full calculation result without bias and it is always well within the grey shaded 1σ error band for the survey specifications considered here. We also include an estimate of the error on the observed angular power spectrum, shown as 1σ contours in light grey. In addition, we show an example of the observational effect contributions to the angular power spectrum of DR1 in Fig. 4, for the full calculation from CLASS.

Fig. 4. Equal redshift angular power spectrum for the fifth redshift bin of Euclid (DR1). We show the full calculation result from CLASS as a black solid line, and the Limber approximation result as a purple dashed line. The green dashed line corresponds to the recalibrated version of the flat-sky approximation. The various grey lines represent the contributions, from observational effects, to the exact solution (calculated using CLASS). The L-L contribution in this correlation is too insignificant to be visible on these axes.

Figures 5 and 6, instead, show the relative errors of the approximations for various auto- and cross-bin correlations of the DR3 configuration. In general, the relative errors grow, for unequal redshift correlations between increasingly disparate bins, due to the sensitivity required to detect the smaller signals. In addition, the relative errors do not vary much for equal-z with increasing redshift, though the Limber approximations for the density and redshift terms become worse. This can be understood as a fixed angular scale corresponding to larger co-moving separations as redshift increases, and thus a poorer recovery of the spectrum by the Limber approximation, which completely neglects k_∥ – an approximation that becomes worse when k_⊥ = ℓ/r(z) gets smaller⁴. We use this to inform our choice of correlations shown, and in Fig. 5, for the case of DR3, we include correlations in redshift where Limber first becomes significantly different to the full result (Bin 4 × Bin 4), where it reaches one of the most significant departures (Bin 8 × Bin 8), as well as the first bin and last bin equal-z correlations.

Fig. 5. Relative errors in the equal redshift angular power spectrum for redshift bin 1, 4, 8, and 13 of Euclid, using 13 equi-populated bins (DR3). We show the full calculation result from CLASS as a black solid line, and the Limber approximation result as a purple dashed line. The green dashed line corresponds to the recalibrated version of the flat-sky approximation. The relative error associated with the Gaussian covariance calculated at 1σ for the given survey configuration is shown as the grey contour. We note that the y-axis ranges are [ − 50, 50], in contrast to Figs. 2 and 3.

Fig. 6. Same as in Fig. 5, but for un-equal redshift bin correlations 1 × 2, 4 × 5, 8 × 9, and 11 × 12 from DR3. The near vertical excursions in the 11 × 12 correlation is due to a zero crossing of the true correlation spectrum.

Overall, the results from the DR1 configuration show that the improved accuracy of the flat-sky approximation does not result in any useful improvement compared to the Limber approximation in the context of the expected observational error, i.e. it will not be possible to measure the spectra accurately enough to distinguish between the two approximation results. This is true across the entire redshift range, and for correlations between any two redshift bins, though we do not show all these spectra here, for the sake of brevity. This is expected, since the DR1 footprint will be about 10–15% of the full area covered by DR3, meaning that the largest scales, where the discrepancy between Eqs. (7), (9), and (18) is more apparent, will be undersampled with respect to the final data release.

However, in the DR3 configuration, the situation is different. Starting as low as the fourth redshift bin (see second panel from the top in Fig. 5), the Limber approximation fails to recover the spectrum well enough to fall within the 1σ contour. This is also true for the equal-redshift spectra in all of the higher redshift bins, and occurs at relatively low ℓ, where we would expect Limber to encounter difficulty in approximating the spectrum. As in the DR1 case, the Limber approximation of neighbouring bin cross-correlations are not, for the most part, significantly different to the flat-sky approximation results, falling within the 1σ contour for the survey in all but one case (Bin 11 × Bin 12, which is itself merely due to a zero crossing in the full calculation correlation spectrum).

If we consider the individual contributions to the spectra from observational effects, it may provide a hint as to why, in particular, the Limber approximation is insufficient for the full Euclid data release. For the DR1 case (see e.g. Fig. 4), the largest contributions⁵ are D-D followed by D-L (more important at higher redshift, and correlations between large redshift separated bins). L-L also increases at higher redshift, but is still only of negligible secondary significance with respect to the former two. As expected for a photometric survey, RSDs are not very relevant. RSD-RSD is always smaller than D-D and for each correlation the largest contributions from D-RSD and D-L are of comparable magnitude, with D-RSD being more important at lower ℓ and closer z correlations.

The DR3 case (not shown) is very similar. The major difference is that due to the slimmer redshift bins, the RSD-RSD contribution can be comparable to D-D and D-RSD in the regimes where they are important, i.e. at large scales in correlations of nearby bins (Tanidis & Camera 2019; Euclid Collaboration: Tanidis et al. 2024). The increased overall inaccuracy of the Limber approximation is likely caused by the failure to accurately recover the RSD terms, as a result.

4.2. Simulation results

We then assessed which, if either, of the two approximations is best to model Euclid data, using a χ² test. This was done by computing the distance between the correct computation of the tomographic harmonic-space power spectra, as per Eq. (7), and the results obtained with either approximation Eqs. (9) and (18). The advantage of this approach is that it takes into account all the possible correlations between different bins, across a range of multipoles ℓ, allowing for the evaluation of the performance as a whole, which is not as easy through simple inspection of the individual spectra alone.

As outlined in Sect. 3.3, we calculated the χ² distribution of each approximation, given the expected error budget of each survey configuration, arising from 1000 Gaussian realisations of measurement noise. Specifically, we adopted each simulation as a synthetic, noisy data set, and computed the χ_X² for both approximations, as per Eq. (26). The further apart these distributions are from the one produced for the CLASS full-sky result, χ_true², the more sensitive the observations are to differences between the calculated spectra and the full calculation.

We examined the distribution of the true model to verify that the framework of our analysis is returning a reasonable fit to the data, and to give a benchmark for the expected width of the distribution in the case of realistic data. In order to judge whether the model is a good enough description of the data, we can determine the number of degrees of freedom and check whether the mean of the χ² distribution lies close to this, i.e. ${\overline{\chi }}^2/\mathrm{DOF}\sim 1$. The overlap in distributions, or Δχ² > 0, can be associated with a p-value, which acts as a frequentist measure of the probability that each distribution can achieve the χ² of the full calculation, subject to data uncertainties. Since we ran only 1000 noisy realisations, our smallest possible p-value is 0.001. Our null hypothesis is that the Δχ² distribution for a particular approximation is consistent with zero. That is, the mean of the χ² distribution of the approximation is not significantly different to (in this particular case, greater than) that of the full calculation. At the 95% confidence limit, the null hypothesis will be rejected for a p-value less than 0.05, and the distributions can be said not to be consistent at that confidence level. The Δχ² distribution is particularly useful here, because each sample in the distribution corresponds to the Δχ² values for corresponding noisy realisations.

In both binning configurations, the full calculation produces a distribution with mean ${\overline{\chi }}^2/\mathrm{DOF}\sim 1$, indicating that the simulated ‘data’ are reasonable realisations of the true result. In the six equi-populated bin configuration, the mean of the true CLASS distribution is at ${\overline{\chi }}^2/{\mathrm{DOF}}_{\mathrm{DR}1}\sim 0.9964$, while for the 13 bin configuration the mean of the true CLASS distribution is at ${\overline{\chi }}^2/{\mathrm{DOF}}_{\mathrm{DR}1}\sim 1.000$.

In Fig. 7 we show the distribution for the 1000 values of Δχ_X² ≔ χ_X² − χ_true². For the DR1 case, the result with the Limber approximation includes the dotted zero line, peaking around Δχ² ∼ 12, whilst the flat-sky approximation only incurs a Δχ² ∼ 0.2. The p-values for these cases are p_Limber = 0.038 and p_{flat − sky} = 0.393, respectively. This means that for the Limber approximation, the null hypothesis is rejected at the 95% confidence level. The flat-sky approximation result is not sufficient to reject the null hypothesis that the Δχ² distribution is consistent with zero.

Fig. 7. The Δχ2 frequency distributions of the flat-sky (green) and Limber (purple) approximations compared to the full result, with the spread provided by the Gaussian covariance for each survey configuration. The DR1 bin configuration is shown in lighter colours, while the eventual DR3 bin result is shown in darker colours. The improvement of the flat-sky over the Limber approximation is seen in the closer proximity of the former’s distribution to the dotted vertical black line at Δχ2 = 0, most notable in the DR3 bin case. The histograms are each normalised to their own peak values and the x-axis is log-scaled above Δχ2 = 50.

Comparatively, for the case of the DR3 data, the distinction between the Limber and flat-sky approximations is significantly larger. In terms of the Δχ², Fig. 7 shows that the flat-sky approximation only reduces the overall recovery of the full calculation result by Δχ² ∼ 5 on average, with p_{flat − sky} = 0.131, indicating that the null hypothesis that the distribution remains compatible with Δχ² = 0 still cannot be rejected at 95% confidence. In contrast, the Limber approximation results in a degradation with a distribution which now peaks around Δχ² ∼ 410, with a p-value at the limit of our number of simulations (p_Limber < 0.001), rejecting the null hypothesis at 99%.

5. Discussion

In this paper we compared the performance of the Limber and flat-sky approximations in the determination of the angular power spectra of galaxy number counts from the photometric survey of Euclid. We performed this study in order to assess whether the Limber approximation that enters present standard analysis packages for Euclid, like CLOE (Euclid Collaboration: Cardone et al. 2025; Euclid Collaboration: Joudaki et al. 2025) performs well enough to use for Euclid data. We showed that for the six-bin analysis foreseen with DR1, the Limber approximation is sufficient to model the angular power spectrum. While leading to a slight degradation of χ² as compared to the full or the flat-sky analysis, it is still able to recover the true result to within the 1σ error, given the survey configuration. This is no longer the case for the DR3, which will be analysed in 13 bins. These bins are sufficiently narrow to render the Limber approximation too crude, with a p-value below 0.1%. In this case the inaccuracies exceed the 1σ errors in the individual equal redshift correlations, resulting in a biased Δχ² distribution, far-removed from what can reasonably be expected from the covariance in the case of the full calculation. In contrast, the flat-sky approximation is able to maintain accurate recovery of the power spectra, with a p-value of 13.1% at worst, indicating that the overall error introduced by the approximation is small enough to avoid biasing the Δχ² distribution.

An accurate approximation is only useful if it can also be used to gain computational speed, compared to the full analysis. Since this work does not promote a particular code, we did not conduct an in-depth analysis of the numerical efficiency here. Nevertheless, it is useful to briefly relay the relevant computational timings encountered, to give an idea of the current state of the problem. The 10-fold improvement in performance compared to the full calculation, as reported in Gao et al. (2024), degrades when dealing with wide photometric windows, which require finer meshgrids to achieve accurate results. This involves increasing the number of discretisations in the comoving distance and its derivative, which are controlled by Nchi and Ndchi in the GZCphysics code, and consequently influence the computation time. In order to achieve the level of accuracy described in this work, we found a speed-up of at least two times for ℓ ≳ 130, for our modified version of the flat-sky code with respect to the full calculation. This is still much slower than the Limber approximation, which is, on average, 200 times faster than the full calculation in this case. Furthermore, for ℓ ≲ 100, the current implementation of flat-sky approximation does not save time compared to CLASS. Optimisations to the code in particular are beyond the scope of this work, and left for future analysis.

As the flat-sky calculation still involves a triple integral, a speed up along the lines of the BLAST code (Chiarenza et al. 2024) will be required. However, further optimisations to the flat-sky code are still possible, for example, by optimising the sampling of the radial window functions (Gao et al. 2024). As demonstrated here, the largest deviations of the Limber approximation occur for ℓ ≲ 20 so, in practice, some combination of the two approximations could be used, depending on the relevant scales of the application. Furthermore, for the D-L and L-L contribution, the Limber approximation remains sufficiently accurate, i.e. with deviation less than a few % of cosmic variance from the full numerical result, see Matthewson & Durrer (2021), so that it can certainly be used for well-separated bins where these terms dominate, though the total signal in these bins is small. Indeed, most of the unequal power spectra amplitudes have a contribution that is less than 30% of the corresponding equal redshift amplitude, and never exceeds 50%, which it attains only for the case of neighbouring bin correlation 7 × 8. Finally, the ℓ-dependence of the D-D and RSD-RSD equal redshift terms can be cast as a simple dependence in the lower boundary of the k-integration. Even though the spectrum calculation requires a triple integral, as explained in detail in Gao et al. (2024), the inner k-integration can be pre-computed, which means that the computation time in the case of MCMC parameter estimation would be reduced by of the order of a power of 2/3 for each point in the explored parameter space.

Finally, we address the relevance of this finding for the spectroscopic survey. The spectroscopic number counts are usually analysed by determining the 3D power spectrum which does not employ the Limber approximation. However, a 6 × 2pt analysis, including correlations of the spectroscopic number counts and the shear, is also planned. When using Limber, either one should use wide redshift bins for the spectroscopic survey also, in which case the 6 × 2 analysis does not add anything to the constraints (see Euclid Collaboration: Paganin et al. 2024), or one should choose slim redshift bins. For the redshift bins considered for DR3, we confirm from the results presented in this work that the Limber approximation will be insufficient, and that the flat-sky approximation provides a promising alternative.

Acknowledgments

RD and WM are partially supported by the Swiss National Science Foundation. SC acknowledges support from the Italian Ministry of University and Research (MUR), PRIN 2022 ‘EXSKALIBUR – Euclid-Cross-SKA: Likelihood Inference Building for Universe’s Research’, Grant No. 20222BBYB9, CUP D53D2300252 0006, from the Italian Ministry of Foreign Affairs and International Cooperation (MAECI), Grant No. ZA23GR03, and from the European Union – Next Generation EU. IT has been supported by the Ramon y Cajal fellowship (RYC2023-045531-I) funded by the State Research Agency of the Spanish Ministerio de Ciencia, Innovación y Universidades, MICIU/AEI/10.13039/501100011033/, and Social European Funds plus (FSE+). IT also acknowledges support from the same ministry, via projects PID2019-11317GB, PID2022-141079NB, PID2022-138896NB; the European Research Executive Agency HORIZON-MSCA-2021-SE-01 Research and Innovation programme under the Marie Skłodowska-Curie grant agreement number 101086388 (LACEGAL) and the programme Unidad de Excelencia María de Maeztu, project CEX2020-001058-M. This work has made use of CosmoHub, developed by PIC (maintained by IFAE and CIEMAT) in collaboration with ICE-CSIC. CosmoHub received funding from the Spanish government (MCIN/AEI/10.13039/501100011033), the EU NextGeneration/PRTR (PRTR-C17.I1), and the Generalitat de Catalunya. The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid, in particular the Agenzia Spaziale Italiana, the Austrian Forschungsförderungsgesellschaft funded through BMK, the Belgian Science Policy, the Canadian Euclid Consortium, the Deutsches Zentrum für Luft- und Raumfahrt, the DTU Space and the Niels Bohr Institute in Denmark, the French Centre National d’Etudes Spatiales, the Fundação para a Ciência e a Tecnologia, the Hungarian Academy of Sciences, the Ministerio de Ciencia, Innovación y Universidades, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Research Council of Finland, the Romanian Space Agency, the State Secretariat for Education, Research, and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site (www.euclid-ec.org).

References

All Tables

All Figures

	Fig. 1. Number density bins simulated for the photometric surveys of Euclid, in 6-bin and 13-bin equi-populated configurations. The numbering convention we used for the ith bin is shown in the coloured blocks above the maxima of each bin.
In the text

Fig. 2. Upper panel: Equal redshift angular power spectrum for the lowest redshift bin of Euclid, using six equi-populated bins. We show the full calculation result from CLASS as a black solid line, and the Limber approximation result as a purple dashed line. The green dashed line corresponds to the recalibrated version of the flat-sky approximation. Lower panel: Relative error (in %) of each approximation from the full calculation result.

In the text

Fig. 3. Upper panel: Equal redshift angular power spectrum for the lowest redshift bin of Euclid, using 13 equi-populated bins. We show the full calculation result from CLASS as a black solid line, and the Limber approximation result as a purple dashed line. The green dashed line corresponds to the recalibrated version of the flat-sky approximation. Lower panel: Relative error (in %) of each approximation from the full calculation result.

In the text

Fig. 4. Equal redshift angular power spectrum for the fifth redshift bin of Euclid (DR1). We show the full calculation result from CLASS as a black solid line, and the Limber approximation result as a purple dashed line. The green dashed line corresponds to the recalibrated version of the flat-sky approximation. The various grey lines represent the contributions, from observational effects, to the exact solution (calculated using CLASS). The L-L contribution in this correlation is too insignificant to be visible on these axes.

In the text

Fig. 5. Relative errors in the equal redshift angular power spectrum for redshift bin 1, 4, 8, and 13 of Euclid, using 13 equi-populated bins (DR3). We show the full calculation result from CLASS as a black solid line, and the Limber approximation result as a purple dashed line. The green dashed line corresponds to the recalibrated version of the flat-sky approximation. The relative error associated with the Gaussian covariance calculated at 1σ for the given survey configuration is shown as the grey contour. We note that the y-axis ranges are [ − 50, 50], in contrast to Figs. 2 and 3.

In the text

	Fig. 6. Same as in Fig. 5, but for un-equal redshift bin correlations 1 × 2, 4 × 5, 8 × 9, and 11 × 12 from DR3. The near vertical excursions in the 11 × 12 correlation is due to a zero crossing of the true correlation spectrum.
In the text

Fig. 7. The Δχ2 frequency distributions of the flat-sky (green) and Limber (purple) approximations compared to the full result, with the spread provided by the Gaussian covariance for each survey configuration. The DR1 bin configuration is shown in lighter colours, while the eventual DR3 bin result is shown in darker colours. The improvement of the flat-sky over the Limber approximation is seen in the closer proximity of the former’s distribution to the dotted vertical black line at Δχ2 = 0, most notable in the DR3 bin case. The histograms are each normalised to their own peak values and the x-axis is log-scaled above Δχ2 = 50.

In the text