Flat-sky angular power spectra revisited

We revisit the flat-sky approximation for evaluating the angular power spectra of projected random fields by retaining information about the correlations along the line of sight. For the case of projections with broad, overlapping radial window functions, these line-of-sight correlations are suppressed and are ignored in the commonly adopted Limber approximation. However, retaining the correlations is important for narrow window functions or unequal-time spectra but introduces significant computational difficulties due to the highly oscillatory nature of the integrands involved. We deal with the integral over line-of-sight wave-modes in the flat-sky approximation analytically, using the FFTlog expansion of the 3D power spectrum. This results in an efficient computational method, which is a substantial improvement compared to any full-sky approaches. We apply our results to galaxy clustering (with and without redshift-space distortions), CMB lensing and galaxy lensing observables in a flat ΛCDM universe. In the case of galaxy clustering, we find excellent agreement with the full-sky results on large (percent-level agreement) and intermediate or small (subpercent agreement) scales, dramatically out-performing the Limber approximation for both wide and narrow window functions, and in equal- and unequal-time cases. In the cases of lensing, we show on the full-sky that the angular power spectrum of the lensing convergence can be very well approximated by projecting the 3D Laplacian (rather than the correct angular Laplacian) of the gravitational potential, even on large scales. Combining this approximation with our flat-sky techniques provides an efficient and accurate evaluation of the CMB lensing angular power spectrum on all scales. We further analyse the clustering and lensing angular power spectra by isolating the projection effects due to the observable- and survey-specific window functions, separating them from the effects due to integration along the line of sight and unequal-time mixing in the 3D power spectrum. All of the angular power spectrum results presented in this paper are obtained using a Python code implementation, which we make publicly available.

A Full-sky angular power spectrum in real and redshift space 31 B Approximating the CMB lensing convergence 32 1 Introduction Next-generation galaxy surveys, such as Euclid [1], DESI [2], Rubin [3], Roman [4], SPHEREx [5,6] and others, aim to address a range of cosmological questions, from uncovering the nature of dark energy and tests of general relativity on large scales [7][8][9], to constraining the properties of the initial conditions of the universe by measuring signals of primordial non-Gaussianity [10][11][12][13][14].
To succeed in these tasks, reliable, accurate and efficient measurements of the galaxy-overdensity statistics are paramount, amongst which the two-point functions (e.g., angular spectra and 3D power spectra) take a leading role.The angular power spectrum has been an observable of choice in many surveys, especially for the study of weak gravitational lensing and the cosmic microwave background (CMB), where the extensive comparisons of theoretical predictions and observations are used with the goal of constraining cosmological parameters.Even for galaxy clustering surveys, where the 3D power spectrum is the most commonly used statistic, the angular power spectrum has certain advantages, e.g., it is defined in terms of variables -redshifts and angular coordinates -which are cosmology independent.Moreover, on large angular scales, the standard 3D power spectrum analysis starts to exhibit effects related to the fixed line of sight.On the other hand, angular power spectra are naturally defined on the full sky.However, the evaluation of full-sky angular power spectra is computationally demanding (see, e.g., [15][16][17][18][19] for discussion), as it involves integration over two spherical Bessel functions, making the integrands highly oscillatory.Furthermore, for spectroscopic surveys many angular power spectra are required between observables projected in narrow radial window functions to avoid loss of information.There have been several attempts to speed up the full-sky evaluation based on the FFTLog algorithm [20][21][22][23]: by expanding the 3D power spectrum in (complex) powers of the wavenumber k, the integration over spherical Bessel functions can be performed analytically in terms of special functions.This provides significant improvements relative to direct computations, but the evaluation still poses substantial computational challenges and thus motivates the search for alternative approaches.One commonly used approach relies on a set of approximations resulting in the 'Limber approximation' solution [24,25], generically accurate on small scales (large multipoles ) and appropriate only for broad and overlapping radial window functions.However, relying on the Limber approximation for galaxy clustering in future surveys could lead to a biased cosmological analysis (e.g., [23]), particularly if the focus is on large or intermediate scales to avoid other modelling challenges (e.g., scale-dependent bias, non-linear clustering and baryonic effects).For these reasons, a middle-ground between the full-sky treatment and Limber-like approximations may be useful.Furthermore, next-generation surveys will exhibit improved photometric-redshift accuracy in clustering measurements with narrower radial window functions, allowing efficient cross-correlation of galaxy fields.This is precisely the regime where the Limber approximation is known not to be valid and thus cannot be relied upon.
An intermediate regime between the full-sky results and Limber approximations has received much less attention in the literature, with the existing works varying in the degree of approximations and corrections they consider [13,16,[26][27][28][29][30].The difference between Limber and these flat-sky results lies in the treatment of the correlations along the line-of-sight, which Limber neglects.Recently, ref. [31] explored the accuracy of the flat-sky approximation by comparing it to the full-sky results in various setups, finding that it can reach subpercent agreement in galaxy number counts and galaxy lensing scenarios.We highlight that the precise form of these various flat-sky approximations, presented in the literature, can differ in certain details that can be traced to the different choices in the geometric and other approximations.Most of the approaches choose the 3D two-point correlation function as the starting point.Moreover, these differences tend to become starker in the unequal-time case (when cross-correlating galaxy sources from different redshift bins), and expressions often get increasingly more complex compared to the equal-time case.
In this work, we revisit the flat-sky approximation enhancing the treatment in several ways: (i) we consider the case of unequal-time projections and discuss how to define a single effective radius for connecting the scale of transverse spatial fluctations to angular fluctuations; (ii) we optimize the flat-sky calculation using the FFTLog algorithm, which can greatly speed up the calculation; and (iii) we compare our results to the full-sky results in the case of galaxy clustering and CMB lensing, providing a detailed analysis of where the flat-sky approximation deviates from the full-sky.We investigate the dependence on 3D scales k of the angular power spectrum for each of these observables, analysing the importance of contributions from different integration regions before and after projection over the radial window functions.
The paper is organized as follows.In section 2, after reviewing the full-sky angular power Symbol Definition δ K ij Kronecker symbol Dirac delta function in n dimensions W (χ) Radial window function; related to the specific observable and survey δ(x) 3D over-density field of matter or biased tracer δ(θ, φ) 2D projected field in angular coordinates on the sky P (k; χ, χ ) Unequal-time power spectrum of the 3D density field C (χ, χ ) Unequal-time angular power spectrum (in the narrow window function limit) C Projected angular power spectrum (with finite width window functions) Table 1: Notation used for the most important quantities throughout this paper.
spectrum and Limber approximation, we derive the corresponding result within the flat-sky approximation retaining integration over the line-of-sight wave-modes.We present two independent derivations of these flat-sky results, starting in Fourier space or real space, respectively.The latter provides and motivates the commonly used geometric recalibration of scales in the flat-sky approximation, which we adopt for our numerical results.In the same section, we describe our implementation of the FFTLog algorithm for a fast and optimized flat-sky calculation.In section 3, we apply our results to galaxy clustering (including also redshift-space distortions) and CMB lensing spectra.Section 4 provides a detailed discussion and comparison of the full-and flat-sky results.We also perform a simple asymptotic analysis of our flat-sky angular power spectrum clarifying its functional dependencies.We conclude in section 5.In appendix A we give a short review of the full-sky angular power spectrum results, while appendix B provides further discussion of approximations we adopt for lensing power spectra.Throughout this paper, we work with the Planck best-fit spatially flat ΛCDM cosmology [32], with CDM density Ω c h 2 = 0.11933, baryon density Ω b h 2 = 0.02242, Hubble constant H 0 = 100h km s −1 Mpc −1 with h = 0.6766, scalar spectral index n s = 0.9665, and fluctuation amplitude σ 8 = 0.8103.In table 1 we summarise our notation for the key quantities that feature throughout the paper.All results presented here are derived with a Python code that we have developed and that we make publicly available on the GitHub repository. 1 It is built using Python3 [33], while all the basic cosmological functions are directly imported from the Boltzmann codes CAMB [34] (for calculation of the matter power spectrum and comparisons with the full-sky angular power spectrum of the CMB lensing potential) and CLASS [35] (for calculation of the linear growth factor and logarithmic growth rate).

Angular power spectra of projected 3D fields
In this section, we study the angular power spectrum for a general observable obtained by projecting a 3D random field.We start with a short review of the full-sky results and the commonly used Limber approximation, also stating the assumptions that the latter is based on.In the rest of the section, we then derive the expression for the flat-sky approximation for the angular power spectrum that, crucially, relies on incorporating the wave-modes along the line of sight.We also discuss the origin of the geometrical recalibration of multipoles , which is often made in the Limber approximation.Finally, we present our efficient numerical algorithm, based on FFTlog, for evaluation of angular power spectra in the flat-sky approximation.

Full-sky angular power spectrum
We consider the projection δ(θ, φ) of some 3D field δ, where θ and φ are spherical coordinates.It is convenient to use the spherical-harmonic expansion of δ, where the coefficients δ ,m are easily obtained from the orthogonality of the spherical harmonics: We consider the projected field δ(θ, φ) in terms of the 3D field δ(χ, θ, φ; η) in real space, where W (χ) is a radial window function describing the observing characteristics of a certain tracer, as well as the specifics of the survey geometry, and η 0 − χ is the conformal look-back time.
The quantity δ(k) is the Fourier transform2 of the 3D over-density field, and is related to the 3D power spectrum as For compactness, we use χ to label the look-back time η 0 − χ in the power spectrum.To maintain generality, we are interested in the angular cross-power spectrum of projected fields δ(θ, φ) and δ (θ, φ), where δ is constructed with a window function W (χ); this is given by where δ K denotes the Kronecker delta symbol.Note how the projected fields are statistically isotropic, with no correlations between different multipoles, a property inherited from the statistical homogeneity and isotropy of the 3D field being projected.As usual, we define the angular power spectrum C as and we have the relation between the angular power spectrum and the unequal-time 3D power spectrum where we have introduced the unequal-time angular power spectrum C (χ, χ ), given by The unequal-time angular power spectrum C can also be viewed as the usual projected angular power spectrum in the limit of narrow window functions, W (χ) → δ D (χ − χ * ).Equation (2.6) cleanly separates the survey-specific radial selection, as encoded in the window functions, and the cosmology dependence, which is all contained in C through its dependence on the unequal-time 3D power spectrum and the mixing and projecting of 3D wave-modes and radial distances.We will focus on these purely geometric and survey-independent projection effects in the next sections when we develop the flat-sky approximation.
In the next sections, we will consider the unequal-time and projected angular power spectra in flat-sky and Limber approximations.In order to distinguish between the different quantities and their approximations, we introduce the label to the above spectra so that, in the full-sky case, we have C full (χ, χ ) and C full correspondingly for the unequal-time and projected angular power spectra.

Limber approximation
In studying the statistical properties of extragalactic nebulae [24], Limber first introduced several approximations to simplify the evaluation of the two-point angular correlation function (the inverse Legendre transform of the above angular power spectrum).Following this result, refs.[25,36,37] derived the Limber approximation directly in Fourier space.Due to its simplicity, involving a single integral for the Fourier-space version, the Limber approximation has remained one of the most commonly used means to evaluate the angular power spectrum.Its validity has been investigated in [38][39][40].
There are two basic assumptions underlying the Limber approximation in Fourier space: (i) the sky is considered as flat; and (ii) the radial window functions are so broad compared to scales of interest for the inhomogeneities that modes with k = 0 are suppressed to negligible levels by the radial integration.The second assumption is the one that we will subsequently relax by including the physical effects of the modes along the line of sight.We briefly re-derive the Limber approximation here, from these assumptions.Let θ be angular position in the plane of the sky around some line of sight n, and δ(θ) the projected field at θ. Starting from eq. (2.2), the projected two-point correlation function in real space becomes where the 3D position x = χ( n + θ) and similarly for x .Using the second assumption, i.e., neglecting the k contributions in the 3D power spectrum P (k; χ, χ ) ≈ P (|k ⊥ |; χ, χ ), we effectively constrain χ = χ , and we can write (2.9) The angular power spectrum in the Limber approximation is then easily obtained by performing the 2D Fourier transform: We shall return to the Limber approximation in section 3, where we further discuss its validity in relation to exact, full-sky results and our new flat-sky approximation.

Flat-sky angular power spectrum
We now revisit the flat-sky angular power spectrum, but retaining the integration over modes along the line of sight.By Fourier transforming the flat-sky projected field δ(θ), we obtain Taking the two-point correlation of this quantity gives us where, in the first line, k = k n + /χ and k = k n + /χ .Note that in the 3D power spectrum, we have set k ⊥ = √ / √ χχ using the 2D Dirac delta function.This is convenient since it preserves the symmetry of the integral on the right under simultaneous exchange of l and l and W (χ) and W (χ ), even without reference to the delta function, which will be important later when we approximate the delta function further.Finally, changing the variables as χ = χ + δχ/2 and χ = χ − δχ/2, we obtain Notice that in eq.(2.13) the 2D Dirac delta function depends on the χ and χ variables and thus cannot just be taken out of the radial integrals.The interpretation of this observation is that the unequal-time two-point correlation function, in the flat-sky approximation, is not diagonal in Fourier space, and there are also non-vanishing contributions when + = 0.However, we saw that this does not happen in the full-sky case in eq.(2.4), and is thus an artefact of performing the flat-sky projections at two different χs (unequal-time).We can understand this behaviour as follows.At radial distance χ, translating θ by a in the flat-sky approximation produces a transverse 3D displacement of χa.The unequal-time two-point correlation function in real space will therefore not be a function of |θ − θ |, but rather |χθ − χ θ | since this quantity is invariant under 3D transverse translations.The 2D delta function δ 2D ( /χ + /χ ) in eq.(2.13) is required to generate this dependence for the real-space angular correlation function.
Nonetheless, we are most interested in scenarios when the two variables χ and χ and close to the mean, i.e., χ ≈ χ ≈ χ δχ.As argued in [19,41], we can formally expand the delta function around δ 2D ( + ) as where δ = δχ/(2 χ) and ∆ = − .The two-point angular correlation function can thus be expanded around its diagonal component as where the components of the projected angular spectra are given by and we have introduced a flat-sky version of the unequal-time angular power spectrum This can be compared to its full-sky counterpart given in eq.(2.7).Recently, ref. [19] showed that this flat-sky result can be obtained as a formal asymptotic limit of the full-sky case when is large and |δχ/ χ| ∝ 1/ .In section 3, we will compare the two and show their close agreement in various cosmological scenarios.
Regarding the higher components of the projected angular spectrum C (n) ( ), we can see that they are suppressed by the powers of δχ/ χ, which for suitable and compact choices of window functions W (χ) is an excellent expansion.In the rest of the paper, we thus consider only the diagonal piece C flat ( ) ≡ C (0) ( ) as our flat-sky approximation of the projected angular power spectrum.Note that for this term, we may replace the argument √ of C (0) in eq.(2.15) with since it is multiplied by the (undifferentiated) delta function.One might be tempted to consider the higher components C (n) ( ) as corrections to the C flat ( ) with the goal of providing better agreement with the full-sky result C full .This is, however, not the correct interpretation.Indeed, as discussed earlier, from the set-up of our flat-sky result we have sacrificed the isotropy property (manifested as the translational invariance in the 2D plane) in the unequal-time case.The correct interpretation of the off-diagonal terms is thus as the estimate of the error on the validity of our flat-sky C flat ( ) result as an asymptotic approximation of the full-sky result C full , where the statistical isotropy is exactly realised.
At this point it is instructive to revisit the Limber approximation.As we have mentioned in the previous subsection, in addition to the flat-sky approximation, the Limber approximation also neglects the dependence of the 3D power spectrum on the wave-modes along the line of sight.
Adopting this additional approximation, our result in equations (2.16) and (2.17) immediately reduces to in agreement with the Limber approximation (eq.2.10).
We emphasise that our flat-sky results and the Limber approximation, therefore, differ only in that we keep the non-zero k information, i.e., we retain the correlations along the line-of-sight.As we shall see, keeping these wave-modes is crucial to capture features observed in the full-sky projected angular power spectrum, especially for projected clustering at unequal times (see [41] for the related discussion).
Before moving on, let us consider in more detail under which physical conditions we can expect the Limber approximation to be applicable, i.e., when the wave-modes along the line of sight can be neglected (see also [38]).We start from the expression for C flat (i.e., C (0) ) given in eq.(2.16).Neglecting the residual dependence on δχ in C flat [all except the term exp(ik δχ)], we see that the wavenumber dependence of the 3D power spectrum takes the form of k 2 = k 2 + ( / χ) 2 .For the Limber assumption that P (k) ≈ P ( / χ) to hold, the effective k domain has to be suppressed so that / χ k .This suppression can arise from the integral over δχ in the case of broad and overlapping radial window functions.If this integral is performed first, neglecting the sub-leading δχ contributions in the 3D power spectrum, it reduces to taking the Fourier transform of the δχ-dependent part of the product of the two window functions. 3Generically, if these δχ window function contributions are broad (we can associate to them some scale σ) their Fourier transform will be narrow and constrained to the range |k | 1/σ.Finally, we conclude that for broad enough window functions (and when the δχ integral bounds modes along the line of sight to the narrow range |k | 1/σ), for angular modes with / χ 1/σ we can use P (k) ≈ P ( / χ), and consequently the Limber approximation.

Geometric recalibration
Our new derivation of the flat-sky approximate result was obtained by working in Fourier space.It is instructive also to tackle the problem starting from the two-point correlation function in real space.This approach has recently been taken in refs.[31,42].As we shall see below, this real-space analysis provides us with some insights into a geometric recalibration that we include in our flat-sky approximation to improve its accuracy. 3For concreteness, let us consider a simple example of these δχ dependencies for equal Gaussian window functions, centred on χ * and of width σ.We have and the δχ integral in the expression for C (0) gives  We aim to compute the unequal-time angular power spectrum C (χ, χ ) from a Legendre transform of the angular correlation function: (2.21) Here, ξ(θ, χ, χ ) is simply the 3D two-point correlation function ξ 3D of the density contrast δ evaluated at radii χ and χ and angular separation θ.It is convenient to choose these two points A and B as in the left-hand panel of figure 1.The 3D distance between the points satisfies This is the same distance as between points A and B , which lie in planes at perpendicular distances χ and χ from the observer, respectively, and have transverse separation √ χχ ω.Here, ω ≡ 2 sin(θ/2) plays the role of a flat-sky angular coordinate and the effective radial coordinate χ g ≡ √ χχ (the geometric mean of χ and χ ) is used to connect ω with transverse distances.This geometry is illustrated in the right-hand panel of figure 1.We note that the mapping from the full-sky θ to the flat-sky ω is area preserving, d cos θ = ωdω.Furthermore, ω naturally arises in the asymptotic expansion of the Legenedre polynomials, as we discuss below.Expressing the 3D distance in terms of δχ and χ g ω, we have The 3D correlation function is the Fourier transform of the (unequal-time) 3D power spectrum, so we can write where ω is a 2D vector in the plane of the sky with magnitude 2 sin(θ/2).We note that we have not made any approximations so far.However, to make further progress in evaluating the Legendre transform in eq. ( 2.21), we switch the integration variable to ω and extend the domain of integration from 0 ≤ ω ≤ 2 to 0 ≤ ω ≤ ∞.Moreover, we approximate the Legendre polynomial by the asymptotic expansion P (cos θ) ≈ J 0 ( ( + 1)ω), where J 0 is the zeroth-order Bessel function of the first kind, valid for large and small θ.This gives The integral over ω can be evaluated using the integral representation of J 0 : where φ ω is the 2D polar angle of ω and L is a 2D vector with magnitude ( + 1).Finally, evaluating the integral over k ⊥ in eq.(2.25), we find This can be integrated over χ and χ [or δχ and χ = (χ + χ )/2] to obtain where ) and (2.28) can be compared to our previous flat-sky results, C flat ( ) in eq.(2.17) and C flat ( ) = C (0) ( ) in eq.(2.16).The only differences are the replacement → ( + 1) and the presence of χ 2 g = χ2 (1 − δ 2 ), rather than simply χ2 , in the prefactor in eq.(2.27).The replacement → ( + 1) ≈ + 1/2 is widely used in applications of the Limber approximation.We refer to it as geometric recalibration.As we show in section 3, we find universally better agreement between our flat-sky results and their full-sky counterparts when including this geometric recalibration.For the results in this paper, we use the prefactor 1/ χ2 rather than 1/χ 2 g as both give very similar results.

Evaluation of the flat-sky angular power spectrum
In this subsection we discuss our strategy for numerical evaluation of the flat-sky angular power spectrum given in, for example, eqs.(2.16) and (2.17).The integral over k is a one-dimensional Fourier transform, which makes the evaluation fairly straightforward.Nonetheless, direct evaluation using the discrete Fourier transform can still be somewhat time consuming.This is especially so when the residual δχ dependencies in the 3D power spectrum are kept since then one cannot use a single Fourier transform to obtain a grid of corresponding δχs, but rather a single transform for each δχ point is required.
In this paper we take a different approach, relying on the fact that the 3D power spectra (be it linear or nonlinear) can be well represented by the discrete Mellin transform.The pioneering application of this method in cosmology was fast-Fourier-transforming the 3D power spectrum and correlation function in ref. [20]; the algorithm has been named FFTLog.Since then, FFTLog has also been used in the computation of nonlinear corrections to cosmological correlators [43][44][45][46][47], as well as in the evaluation of the full-sky angular power spectrum [21][22][23]48].We will follow a similar route to ref. [21], but with the difference of applying the FFTLog algorithm to compute the flat-sky angular power spectrum.The advantage, as we shall see, is that the resulting flatsky expressions are significantly simpler than their full-sky counterparts, yielding a significant computationally speed-up.
Our starting point is thus to represent the 3D power spectrum P (k) in terms of a sum of (complex) powers of the wavenumber k, i.e., P (k) i α i k ν i .In this work, we constrain our analysis to the linear version of the 3D power spectrum for which the time dependence is separable: where p(k) is the 3D linear power spectrum at z = 0 and D(χ) is the linear growth factor normalised to unity at z = 0. Working with the linear form of the power spectrum provides certain simplifications, however, it is important to stress that in no significant way does this represent a limitation of the method, and a similar procedure can be adopted when using any nonlinear P (k) results.As discussed earlier, the relation between multipoles and 3D wavenumbers k is, after geometric recalibration, where we introduced the shorthand notation ˜ ≡ ( + 1)/( χ√ 1 − δ 2 ).We note that ˜ has dimensions of inverse length.We use the FFTLog algorithm to obtain the coefficients and powers of the expansion p(k) = i α i k ν i , where we note that the frequencies ν i are complex numbers with fixed negative real part (known as the bias).One can always change the real bias term to the frequencies, with an associated change in the coefficients α i , by multiplying p(k) by the appropriate power of k before performing the transform.In our implementation, we have found that using (ν i ) = −0.5 gives excellent recovery of the 3D power spectrum, however, one can always change the biasing term in our code, if desired.
The flat-sky result for the unequal-time angular power spectrum given in eq.(2.17) can thus be written as ) where K ν (z) is the modified Bessel function of the second kind, and where we absorbed the D ( χ + δχ/2) D ( χ − δχ/2) factor in the α i coefficients.It is useful to introduce a function defined as which enables us to rewrite the unequal-time angular power spectrum C flat ( ) in the following form: We note that M (2) , which is a property inherited from the fact that C flat ( , χ, δχ) = C flat ( , χ, −δχ).The representation of our results in terms of the M (2) ν (x) functions provides us with several benefits.First, we note that these functions do not carry any information on the cosmological parameters, i.e., in any cosmological analysis (and once the grid of ν i has been fixed) M (2) ν i (x) can be pre-computed and interpolated in the x variable.This allows for almost instantaneous evaluation of these functions for any choice of x = ˜ δχ.Note that the argument ˜ δχ depends on cosmology via the definition of the radial distances; however, we can ensure that M (2) ν i (x) is pre-computed over a sufficiently wide range to encompass any reasonable variation in cosmology.The rest of the cosmological information is, of course, carried by the α i coefficients obtained by the single FFTLog decomposition of the linear power spectrum.The second beneficial property is that we require (pre-evaluation) of only N ν special functions M (2) ν i (x) functions, where N ν is the number of frequencies used in the FFTLog expansion of p(k).This is to be contrasted with the evaluation of the full-sky C full with the FFTLog expansion, as put forward in refs.[21,22] and which we summarise in Appendix A. The equivalent full-sky representation of the angular power spectrum is given in eq.(A.4); it depends on the functions I (ν, t) (defined in eq.A.3 and where t = χ /χ).For a given ν i , separate evaluations of special functions (in this case, involving the hypergeometric function 2 F 1 ) are required over a grid of and t values as the arguments do not combine into the single x = ˜ δχ as in the flat-sky case.The full-sky calculation thus requires a much larger number of special-function evaluations.Moreover, in each evaluation the hypergeometric function 2 F 1 is required, which is more involved and thus computationally more costly than the M (2) ν (x).Being able to pre-compute the functions M (2) and then interpolate the argument makes the computational performance advantage of the flat-sky approximation even starker.
In our default setup, we use N ν = 201 modes for the expansion of p(k), which is sufficient to reproduce precisely the linear power spectrum at z = 0. We typically choose k min = 10 −8 h/Mpc and k max = 52.0h/Mpc when taking the transform of p(k).With these choices, the fractional error in the reconstructed p(k) is below 0.5% for 5 × 10 −4 h/Mpc ≤ k ≤ 8 h/Mpc.One can further optimise the choice of frequencies if required.As mentioned, we then pre-calculate the set of M (2) ν i (x) functions choosing 2000 sampling points for their x arguments.For projected galaxy clustering, if we sample χ and δχ at 50 and 100 points, respectively (which we found to provide a converged result), it takes about 12 s to calculate 150 multipoles in Python3 on a personal laptop.Note that these timings do not include generation of the arrays of window functions and growth factors, and depend on the choice of the length of sampling arrays.The algorithm could certainly be further optimised for application in likelihood analysis, where computational time is critical.The details of our calculation and code documentation can be found in the GitHub repository. 5

Results for galaxy clustering and CMB lensing
In this section, we first use our results to compute the projected angular power spectrum in the case of galaxy clustering.We then extend this analysis by including redshift-space distortions.Lastly, we turn to the case of CMB lensing.We compare our results to the full-sky FFTLog-based calculation following ref.[21], which we have also implemented in our code.In each case, we also show comparisons with the Limber approximation.

Galaxy clustering
We focus first on projected galaxy clustering, where galaxies are selected within a given redshift range according to the radial window function W (χ). We assume a simple case where the window function is Gaussian, specified by a central distance χ * and width σ χ , i.e., where, for narrow window functions, σ χ is related to the standard deviation of redshift via σ χ = cσ z /H(z) for a 0 = 1.We highlight that our computation of the projected angular power spectrum C , as given in eq.(2.28), is not particularly sensitive to characteristics of W (χ) such as its smoothness.In particular, we do not require W (χ) to be differentiable (a property used, for example, in [21]) or to have a smooth and well-behaved Fourier transform (utilised in [23]).This may prove helpful in analyses of survey data, where the window functions may be complicated due to photometric redshift uncertainties, for example.Rather, we focus on developing fast computational methods for the unequal-time angular power spectrum C , which is then integrated against the pair of window functions, without needing to put additional constraints on these, to obtain the projected spectrum C .
In figure 2, we present results for the angular power spectrum of galaxy clustering projected with equal window functions centred on redshifts z = 1.0 and 2.0, and with widths σ z = 0.05 and σ z = 0.3, respectively.We find excellent overall agreement of the full-and flat-sky results, on all scales, while the Limber approximation deviates significantly from these on large scales.For multipoles 50, the flat-sky approximation (before geometric recalibration) only deviates from  In both panels, we compare the full-sky expression (black dash-dotted line), the Limber approximation (blue dashed line), and the flat-sky approximation (red solid line).The green solid line and cyan dashed line are geometrically recalibrated flat-sky and Limber approximations, as described in section 2.4.In the bottom panels we show fractional residuals compared to the full-sky results.
the target full-sky angular power spectrum by a few percent at most.These small differences between the full-and flat-sky results are removed very effectively when the additional geometric recalibration → ( + 1) is implemented in the flat-sky result (see the discussion in section 2.4).Geometric recalibration in the Limber approximation also improves its agreement with the fullsky result somewhat, but starting from much larger errors than our flat-sky results.We note that for the relatively narrow window functions (left panel of figure 2), the fractional error in the Limber approximation remains at the percent level or higher up to multipoles of a few hundred.The Limber approximation performs better for the higher-redshift and broader window functions (the right panel of figure 2), as expected from the discussion in section 2.2.
We now consider the cross-correlation between the signal projected with unequal-time window functions.The magnitude of the angular power spectrum is shown in figure 3 for two setups: first with central redshifts z = 1.0, z = 1.25 and equal widths σ z = 0.05; and second with z = 2.0, z = 3.5 and σ z = 0.3.We again find very good overall agreement of the full-and flat-sky results, especially once the geometric recalibration is taken into account.In both scenarios considered, the angular power spectrum is negative (anti-correlated) on large scales while remaining positive on smaller scales.This anti-correlation feature is also captured very well by the flat-sky result, providing us again with sub-percent agreement with the full-sky result on all scales (away from the zero crossing) after geometric recalibration.We highlight that even without the recalibration, the flat-sky result correctly captures the sign change and the shape of the full-sky results.The effect that the → ( + 1) recalibration provides is a slight horizontal shift to lower multipoles  on large scales.In contrast, the Limber approximation fails to capture the sign change on large scales, being strictly positive on all scales, and, in general, performs poorly for cross-correlations with little overlap between the radial window functions.
As mentioned in section 2, the key difference between the Limber approximation and our flat-sky result is that the latter retains information about the k wave-modes, which the Limber approximation explicitly disregards.As a result, the flat-sky approximation can model both equal-time and unequal-time correlations.The inclusion of these modes along the line-of-sight is essential to capture the correct behaviour of the cross-power spectrum for unequal-time window functions on large scales, where they lead to anti-correlations (figure 3).Keeping the integration over k is thus the primary source of the improvements we observe in the flat-sky results over those from the Limber approximation.Both of these approximations, of course, fail to account for the effects of sky curvature beyond what is already captured by the geometric recalibration → ( + 1).However, our results suggest that these are only relatively minor corrections for any realistic clustering survey setup and geometry, even on the largest scales ( < 10).

Galaxy clustering with redshift-space distortions
In addition to the Hubble expansion, galaxies experience an additional, large-scale peculiar motion due to gravitational infall .This motion, via the Doppler effect, causes the apparent displacement of the real-space galaxy distribution along the line-of-sight direction when observed in redshift.In linear theory (assuming the flat-sky and distant-observer approximation), the effect on the redshift-space over-density of sources (e.g., galaxies) can be expressed as where f ≡ d ln D/d ln a is the logarithmic growth rate, and b is the linear galaxy bias, which we set to unity throughout this paper.Equation (3.2) is the well-known Kaiser effect [49,50] describing redshift-space distortions (RSD).The unequal-time 3D linear power spectrum of δ s is then simply given as where P L (k; χ, χ ) is the 3D linear matter power spectrum at lookback times χ and χ .Note that P s (k) depends separately on k and k ⊥ as the fixed line-of-sight breaks statistical isotropy.Using the 3D RSD power spectrum in the expression for the unequal-time angular power spectrum (eq.2.17), and keeping track separately of k and k ⊥ , we have where k = k 2 + ˜ 2 .We decompose the linear power spectrum as a sum of power laws in the same way as in section 2.5.To handle the additional k 2 and k 4 terms, we take derivatives of the expression given in eq. ( 2.31) with respect to δχ (treating ˜ and χ as parameters that are not differentiated) and simultaneously shift the ν i index down.We obtain the following: where the functions M (2) ν (x) are defined in eq.(2.32).Evaluating the derivatives using relations for derivatives of modified Bessel functions, we have 8) The treatment and numerical implementation of these functions is similar to our discussion of M (2) ν (x) in section 2.5.Using the FFTLog expansion, eq.(3.4) can now be written as (where the linear growth factors are again absorbed into the α i coefficients)  This is the expression we have implemented in our code.Implementation requires two additional sets of pre-calculated special functions compared to the non-RSD case.As a result, the total precalculation time increases approximately three-fold (taking around 15 min on a personal laptop, using 2000 sampling points for ˜ |δχ|).For a given choice of window functions and their sampling, the computational time for C flat ( ) including RSD also increases five-fold to approximately 0.45 sec per multipole.Before presenting numerical results, let us briefly discuss the Limber approximation in the presence of RSD.In its strictest form, RSD do not contribute in the Limber approximation since they require non-zero k .However, by relaxing the Limber assumptions to a certain level, ref. [51] provides a way of extending the Limber approximation to include RSD. 2 We have implemented this model in order to compare and contrast the performance of our results.
Figures 4 and 5 show the angular power spectra for projected galaxy clustering, including RSD, radial window function W (χ), and then using the following effective weight in the usual Limber approximation (with the 3D power spectrum evaluated at k = ( + 1/2)/χ):  for the same equal-and unequal-time setups as in section 3.1.Our flat-sky approach is in very good agreement with the full-sky results.Similar to the case without RSDs, before geometric recalibration, the flat-sky approximation deviates from the full-sky result by less than 4% (for the equal-time case) and less than 6% (for the unequal-time case) on all scales of interest.Our results again correctly capture the large-scale anti-correlation feature in the unequal-time case.When the → ( + 1) geometric recalibration is taken into account, our flat-sky results further improve, reducing the relative error to less than 2% for both equal-time and unequal-time cases, even on the largest scales.
In comparison, the extended Limber approximation of ref. [51] exhibits similar behaviour as without RSD for the case of equal-time window functions (figure 4).As expected, the approximation is only accurate at smaller scales and performs better with broader window functions, although the RSD are suppressed for such window functions.In the case of unequal-time window functions (figure 5), with RSD included, the extended Limber approximation does capture the large-scale anti-correlation feature in the angular power spectrum, although the shape is not accurately reproduced.We note that while the usual Limber approximation without RSD is strictly positive on all scales, the RSD contributions to G Limber (χ) of the extended approximation may be negative on large scales for some range of χ, allowing the power to be negative too.Since RSD make a significant contribution on large scales in figure 5 (compare with figure 3), this may account for the rough agreement there between the extended Limber approximation and the full-sky spectrum.

CMB lensing
Moving on from the galaxy clustering case, in this subsection we apply the flat-sky approximation to CMB lensing (see ref. [52] for a review).The variable of interest is the CMB lensing potential φ, which is related to the line-of-sight integral of the Newtonian potential Ψ (or, more carefully, the Weyl potential) in the Born approximation: where χ * is the comoving distance of the CMB last-scattering surface (note that in all the calculations, we assume a flat universe).Two differences emerge, therefore, when comparing the galaxy clustering and CMB lensing calculations, as follows.
(i) The choice of window function: (ii) The field being considered is the Newtonian potential Ψ rather than the matter over-density.They are related by the Poisson equation (in Fourier space and assuming anisotropic stress is negligible) where ρm (z) is the background matter density.
It is straightforward to relate the 3D power spectra of Ψ and δ according to where Ω m (χ) is the fraction of the energy density in matter, and H(χ) is the Hubble parameter.Defining N Ψ (χ) ≡ 3Ω m (χ)a 2 (χ)H 2 (χ)/2, our flat-sky approximation to the angular power spectrum of the lensing potential is dk 2π e ik δχ 1 k 4 P δ (k; χ, δχ) , (3.14) where the total wavenumber k is again given by eq.(2.30).
For the 3D power spectrum of the matter over-density, P δ , we again take the unequal-time linear theory prediction given in eq.(2.29).After performing the FFTLog expansion of p(k), all we need is to shift every ν i in eqs.(2.31) and (2.33) to ν i − 4, while keeping the α i unchanged.The angular power spectrum of the lensing potential can then be written as where the appropriate unequal-time angular power spectrum for lensing is In the CMB lensing literature, when employing the Limber approximation, the angular power spectrum of the lensing convergence κ is generally calculated rather than the lensing potential φ.The convergence is related to the angular Laplacian of the lensing potential: κ = −∇ 2 θ φ/2, so that the spherical power spectra are related by C κκ = [ ( + 1)] 2 C φφ /4.In the flat-sky approximation, we have The field being projected is, therefore, the transverse Laplacian of Ψ, which in Fourier space is If we calculate the angular power spectrum of κ directly with our flat-sky approximation, we have where ˜ = ( + 1)/χ g with χ g = χ√ 1 − δ 2 = √ χχ as usual (after geometric recalibration).Here, the factor χ4 1 − δ 2 2 = (χχ ) 2 comes from the conversions between angular and transverse Laplacians at radii χ and χ , and ˜ 4 from k 4 ⊥ .Their product is [ ( + 1)] 2 and so, comparing with eq.(3.14), we recover in agreement with the full-sky relation.As we shall discuss shortly, the flat-sky expression (3.14) for the angular power spectrum of the CMB lensing potential is less accurate on large scales than the clustering results presented in the previous section.This is because of the different scale dependencies of the 3D power spectra P (k) (gravitational potential versus over-density), with the lensing case having most power on large scales, and the window functions, with lensing being very broad and extending to χ = 0 where W Ψ ∝ 1/χ.We propose an alternative flat-sky approximation whereby we ignore the distinction between the transverse and full Laplacian in the lensing convergence, approximating In this case, the scale dependence of the field being projected is the same as for clustering.The approximation wrongly includes radial derivatives of the gravitational potential.We test the accuracy of this approximation firstly in the full-sky case.In figure 6 (left panel), we compare the full-sky angular power spectrum of the approximate convergence (eq.3.21) with the full-sky [ ( + 1)] 2 C φφ /4; the maximum difference is below 0.5% (and much less on smaller scales).The effect of the radial derivatives is very small, even on large scales, because of the integration over the broad lensing window function.We provide further insight into this approximation in appendix B. We therefore proceed to compute the angular power spectrum of the approximate convergence in eq.(3.21) using the flat-sky approximation.This is simply given by eq. ( 6: Angular power spectrum of the CMB lensing potential.In the left panel, the green solid line is our flat-sky approximation computing the power spectrum of the lensing potential directly (eq.3.14), while the red solid line is our alternative flat-sky approximation, computing the convergence directly with the full 3D Laplacian rather than just the transverse Laplacian (eq.3.22).The cyan dashed line is the Limber approximation.This, and the flat-sky results in this panel, include geometric recalibration.The black dot-dashed line is the exact full-sky angular power spectrum, while the blue solid line is from the full-sky convergence power spectrum but computed with the full 3D Laplacian, rather than the transverse Laplacian, acting on the gravitational potential.The latter wrongly includes radial derivatives of the potential, but introduces only a very small error.In the right panel, we compare the exact full-sky result to the flat-sky (working from the lensing potential) and Limber approximations without geometric recalibration.Comparing with the same colour lines in the left panel illustrates the importance of recalibration (i.e.replacing with ( + 1)) on large scales for lensing.
in the integral over k replaced with 1/ ˜ 4 , so that dk 2π e ik δχ P δ k 2 + ˜ 2 ; χ, δχ .(3.22)We shall also compare to the Limber approximation.For CMB lensing this is while with geometric recalibration we replace with ( + 1). Figure 6 compares several approximations to the CMB lensing angular power spectrum with the full-sky result.The flat-sky result in eq.(3.14), which starts from the lensing potential φ, and the Limber approximation in eq.(3.23) 6, but for the lens source plane at z = 2 rather than the CMB last-scattering surface.
after geometric recalibration.The importance of recalibration is shown in the right panel of the figure; without it, the fractional errors in the flat-sky and Limber approximations are significantly larger at low multipoles and not until 100 are they at the percent level.For 10, the recalibrated flat-sky and Limber approximations are less accurate.As noted above, these larger errors arise from the shape of P Ψ (k) and the form of the window function for CMB lensing.However, we find that if instead we compute the convergence power spectrum approximating the transverse Laplacian with the full 3D Laplacian (eq.3.22), the errors on large scales are reduced to below one percent and are much below this on smaller scales.This provides a uniformly accurate approximation, and outperforms the Limber approximation on all scales.
Curiously, the Limber approximation after geometric recalibration is rather more accurate on large scales than the flat-sky approximation that takes the lensing potential as its starting point (compare the cyan and green lines in figure 6).This is only the case after geometric recalibration, which appears to overly suppress the CMB lensing angular power spectrum.

Galaxy lensing
We find similar results as for CMB lensing when the source plane is at lower redshift, as appropriate for lensing of galaxies.In figure 7, we show the angular power spectrum of the lensing potential for sources at z = 2 (compared to z ≈ 1080 for the CMB last-scattering surface).In figure 8, we consider the cross-correlation between the lensing convergence field for these sources and projected clustering centred at z = 0.5 with σ z = 0.2.In both cases, replacing the transverse Laplacian by the full 3D Laplacian in the spherical convergence is still a very good approximation on all scales.Furthermore, the (recalibrated) flat-sky approximation with the 3D Laplacian very accurately reproduces its full-sky equivalent.4 Unequal-time angular power spectra, C In the previous section, we compared the full-sky and flat-sky-approximated projected angular power spectra for galaxy clustering and CMB lensing.In both cases, the results are given as integrals of the unequal-time angular power spectrum C over the observable-and survey-specific window function.The expression for the full-sky case is given in eq.(2.6), while for the corresponding flat-sky case we refer to eq. (2.16) (or eq.2.28).Any errors in the flat-sky approximation for the projected angular power spectra must therefore6 arise from errors in the unequal-time spectra C .Moreover, C is entirely independent of the specific survey and functional forms of the chosen window functions and thus solely determined by cosmology.In this section, we analyse these unequal-time spectra in more detail, comparing the full-sky C , given in eq.(2.7), and the flat-sky C flat ( ), given in eq.(2.17).For CMB lensing, we consider the calculation of the lensing potential, wherein we project the gravitational potential, to highlight the effect of the very different scale dependence of the 3D power spectrum.The relevant flat-sky C flat ( ) is given in eq.(3.16).Throughout this section, we adopt the geometric recalibration of the flat-sky spectra, replacing with ( + 1).Recall that the key differences between the calculations for matter clustering and the CMB lensing potential are: (i) the scale dependence of the 3D power spectra P (k) (over-density versus gravitational potential), with the lensing case having most power on large scales; and (ii) the   window functions, with lensing being very broad and extending to χ = 0 where W Ψ ∝ 1/χ.Both differences would be expected to exacerbate errors in the flat-sky approximation at low multipoles, as seen in figure 6.Only the difference in the 3D power spectra impacts the C directly, but the window functions control what ranges of χ and χ (or χ and δχ) contribute significantly to the projected spectra C .In figures 9 (for clustering) and 10 (lensing), we plot in the top two rows C full and C flat ( ) as a function of χ and χ for multipoles = 2, 10 and 100.Comparing the two figures, we see clearly how the dependence of C on distances is determined by the shape of the 3D power spectrum.For clustering, the support of C is strongly clustered around χ = χ (δχ = 0) and so correlations between different radii fall sharply with their separation.For window functions broad compared to the off-diagonal width of C , the Limber approximation will be very accurate.In constrast, for the lensing case there are significant correlations between widely separated radii, with the distributions becoming broader at low multipoles.These differences arise from the very different k dependence of the 3D spectra P ( k 2 + ˜ 2 ) when projecting δ or Ψ, being much more concentrated around k = 0 for the lensing case.Since the δχ dependence of C is determined mostly by the Fourier transform of the 3D spectrum with respect to k , the narrower 3D spectra for lensing give broader correlations in δχ.The bottom two rows in figures 9 and 10 multiply the C by the appropriate pair of window functions; the integral of the resulting quantity over χ and χ gives the observable projected spectra.For clustering, the window functions are Gaussians centred on z = 1 with widths σ z = 0.05.At low multipoles, the lensing angular power spectrum receives significant contributions from nearby structures, a point we shall return to below.For both clustering and lensing, we see good agreement between the full-and flat-sky results, and indeed at the level of figures 9 and 10 it is difficult to see any differences.
In order to highlight specific aspects of the the unequal-time angular power spectra, we plot slices through C (χ, χ ) in figures 11 (clustering) and 12 (lensing).The top rows show the dependence on χ ≡ (χ + χ )/2 at δχ ≡ χ − χ = 0.For the clustering case, C monotonically decreases with χ, while for lensing it monotonically increases.In section 4.1 we explain how the asymptotic behaviours of the C at small and large χ, for δχ = 0, follow from the shape of the specific 3D power spectrum being projected.The bottom rows of the figures show the changes in the compactness of C considered as a function of δχ for a fixed χ = 2381 Mpc/h (corresponding to z = 1.05).In both clustering and lensing cases, the peak of C over δχ becomes narrower as increases so that correlations between different radii are suppressed.As also seen in the 2D plots in figures 9 and 10, the peak is much broader and flat-topped for lensing than clustering.It is also noteworthy that the clustering C can be negative for low multipoles, giving anti-correlations between different radii, but the lensing C cannot.This behaviour arises from the 3D clustering power spectrum, P δ ( k 2 + ˜ 2 ), having a local minimum at k = 0 for ˜ k eq , where k eq is the matter-radiation equality scale.
Let us focus next on the comparison of the full-sky and flat-sky C in the sectional plots in figures 11 (clustering) and 12 (CMB lensing).At low multipoles, differences start to appear in the lensing case for small radii χ and χ or, equivalently, at low χ for δχ = 0 (as in the top row of figure 12), or for |δχ| ≈ 2 χ (bottom row of the figure).These differences are exacerbated by multiplication by the lensing radial window functions, which, recall, behave as W Ψ ∝ 1/χ for χ χ * , as shown in the second row of the figure.At the lowest multipoles, significant errors in C flat ( ) from lenses at χ 200 Mpc/h lead to the relatively large errors in the angular power spectrum of the CMB lensing potential seen in figure 6.However, we emphasise again that a simple and accurate work-around is to calculate the flat-sky convergence power spectrum, replacing the transverse Laplacian by the full 3D Laplacian.
Finally, let us look in more detail at the range of distances that contribute to the CMB lensing angular power spectrum.We show in figure 13 the cumulative contribution to C φφ from χ ≤ χup (integrating over all allowed δχ).We compare the full-sky result, the Limber approximation and the flat-sky approximation (proceeding via the lensing potential).As expected from the   9 but for the CMB lensing case, where the relevant 3D power spectrum is δ (k)/k 4 rather than P δ (k).In the bottom two rows, the lensing kernels (eq.3.11) are included.
projection, CMB lensing picks up more contributions from larger distances at higher multipoles.At = 2, around 10 % of the angular power comes from χ < 200 Mpc/h, where figure 12 shows there are significant errors in the flat-sky approximation.

analysis of the flat-sky C
In this subsection, we discuss the asymptotic behaviour of the unequal-time angular power spectrum for small and large χ at δχ = 0, corresponding to the top rows in figures 11 and 12.In the earlier sections, we have found good agreement between the full-sky and the flat-sky behaviour; for simplicity, we therefore focus on analyzing the flat-sky case only.For this purpose, we also reduce the 3D power spectrum to a very simple form, roughly capturing the behaviour of the ΛCDM matter power spectrum in the IR/UV regime.We thus have where k eq is the power spectrum equality scale and 0 < < 1.
Setting δχ = 0 in eq.(2.17), C flat ( ) simply reduces to integrating the matter power spectrum over k .We consider first the behaviour as χ → 0 so that ˜ k eq for all .For k ˜ , we approximately consider the power spectrum as independent of k , i.e., k 3 eq p(k) ∼ ( ˜ /k eq ) −3+ , while for k ˜ , we can ignore the contribution of ˜ in k, so that k 3 eq p(k) ∼ (k /k eq ) −3+ .For clustering, therefore, we have where we have used the asymptotic property of the linear growth factor at low χ, D( χ) → const.The only change in the lensing case is that we divide P (k) by an additional factor of k 4 , which gives the asymptotic χ4− behaviour.These asymptotic behaviours can be observed in the top rows of figures 11 and 12 for projected clustering and lensing, respectively.
Similarly, for large χ, we can take ˜ to be much smaller than k eq .We have to consider three regions in the matter power spectrum.For k ˜ , we still have a constant k 3 eq p(k) ∼ ˜ /k eq .For ˜ k k eq , we take k 3 eq p(k) ∼ k /k eq , while for k k eq , we take k 3 eq p(k) ∼ (k /k eq ) −3+ .Collecting the three regimes gives us For projected clustering and k eq χ 1, this implies that C flat ∼ D 2 ( χ) χ−2 , consistent with the steeper fall-off in the top row of figure 11 at large χ.Following the same calculation for lensing (and dividing the power spectrum by k 4 ), we have where we have already used the property that D( χ)N Ψ ( χ) is approximately constant at large χ.This means that C flat ( ) approaches to a constant for large χ, as we observe in the top row of figure 12.

Conclusion
We revisited the performance of the flat-sky approximation in evaluating the angular power spectrum of projected fields.Unlike the commonly used Limber approximation, we retained the contribution of wave-modes along the line of sight (labelled by k ).We showed that with this inclusion, very accurate spectra can be obtained with the flat-sky approximation including cases where the projection is with narrow radial window functions or window functions with limited overlap.Moreover, the approximation generally remains accurate at low multipoles, where curved-sky effects are expected to be significant, with a simple geometric recalibration of the flat-sky results.Our flat-sky approximation also provides a self-consistent way of including redshift-space distortion effects, given that we explicitly retain the dependence on the k modes.Our results were developed in a general form, independent of the specifics of a given observable or survey.For this purpose, we separated the computation of the projected angular power spectrum into computation of the unequal-time angular power spectrum, C (χ, χ ), and integration over survey-and observable-specific radial window functions, W (χ) and W (χ ).The C (χ, χ ) includes only the projections of the relevant unequal-time 3D power spectrum onto radii χ and χ , and takes the form of a Fourier integral over k with conjugate variable δχ = χ − χ .Given the specifics of each observable and survey, the treatment of the radial integrals over the window functions can then be further optimised.
A key step in our numerical evaluation of the angular power spectrum in the flat-sky approximation is accelerating the evaluation of C (χ, χ ) based on a discreet Mellin transform of the 3D power spectrum, commonly called the FFTLog algorithm.We utilise FFTLog to expand the 3D power spectrum as a sum of power-laws in k with complex frequencies ν i .The Fourier integral required for C (χ, χ ) can then be expressed as a sum of modified Bessel functions of the second kind K ν i (x).In comparison, an equivalent procedure in the calculation of the fullsky angular power spectrum results in a sum of terms involving the hypergeometric functions 2 F 1 (a, b; c; z) [21], increasing the computational complexity and evaluation time.Moreover, these modified Bessel functions depend only on the variable x = δχ ( + 1)/ √ χχ .We can therefore easily pre-compute and store the K ν i (x) functions and interpolate them as required to evaluate the angular power spectrum at different points in parameter space when sampling over cosmological parameters during inference.For simplicity, we presented results using only the linear-theory 3D power spectrum.However, the method is general and can be applied to any choice of the 3D power spectrum (not necessarily given in perturbative form).We note that further numerical optimisations are possible, e.g., implementing an optimised set of the complex frequencies ν i and improved sampling for the integration over the radial window functions.We are exploring the former, and emulation of the associated Mellin transform, in ongoing work.Even so, and irrespective of these further optimisations, our flat-sky approach already provides a significant computational improvement compared to the full-sky FFTLog-based algorithm.
We presented results for projected clustering (e.g., galaxy surveys) with and without redshiftspace distortions, gravitational lensing and their cross-correlation.In general, we found excellent (percent-level or better) agreement of the flat-sky and full-sky results on all scales (including also the lowest s), with the flat-sky results outperforming the Limber approximation.However, we found rather larger errors from the flat-sky approximation when evaluating the angular power spectrum of the CMB lensing potential at low multipoles.These errors arise since the power spectrum of the gravitational potential, whose projection determines the CMB lensing potential, is very red with most power on large scales.Furthermore, the radial window function rises as 1/χ at small radial distances, so at low multipoles non-negligible contributions to the lensing power spectrum come from nearby lenses.We presented a simple work-around that restores the accuracy of the flat-sky approximation on large scales.This involves working with an approximate form of the lensing convergence, which approximates the transverse part of the Laplacian of the gravitational potential (coming from the angular Laplacian in the convergence) by the full 3D Laplacian including radial derivatives.In this manner, the approximate convergence involves the projection of the matter over-density, which has a power spectrum that peaks at the equality scale k eq .As an aside, we verified that using an equivalent approximation in the spherical convergence is very accurate at all multipoles.
To conclude, in this paper, we investigated the performance of the flat-sky approximation taking into account the correlations along the line of sight.This gives an accurate and efficient approximation of the full-sky angular power spectrum (at sub-percent level) for galaxy clustering (including redshift-space distortions) and gravitational lensing.We have developed an efficient Python implementation of our flat-sky method, based on the FFTLog expansion, which we make publicly available. 7With this method and further optimisation, one can compute the projected angular power spectrum at speeds comparable to the Limber approximation.This opens an alternative route for efficient computation of angular power-spectrum observables in parameter searches.Such an efficient framework is especially important in light of many upcoming surveys, taking data on large and intermediate scales and at high signal-to-noise ratios.Moreover, effects due to the photon geodesic projections (i.e., relativistic corrections) will be an important consideration for these surveys.The addition of these effects goes beyond our current work; however, incorporating them into our framework should be fairly straightforward.Moreover, and in line with the aforementioned surveys, the cross-correlations of unequal-time narrow window functions are expected to play an increasingly important role in the future cosmological analysis.Thus, our computational framework offers a valuable and efficient means to analyse these upcoming data sets.
where W = W (χ) and W = W (χ ).We use this expression to compute our reference full-sky results presented in section 3.
If we include RSD, the unequal-time angular power spectrum on the full sky is modified to C χ, χ = 4π k 2 dk 2π 2 P (k; χ, χ ) b j (kχ) − f j (2) (kχ) b j (kχ ) − f j (2) (kχ ) , (A.5) where j (2) (kχ) is the second derivative of the spherical Bessel function evaluated at kχ.We have included the galaxy bias, which depends on k and χ generally, although we have set it to unity in the rest of the paper.Recall, also, that f ≡ d ln D/d ln a is the logarithmic growth rate.To deal with the derivatives of the spherical Bessel functions, ref. [21] made use of integration by parts to reduce all integrals to the form in eq.(A.2).However, this moves the derivatives onto the window functions (and growth factors and rates), which may introduce complexities in the case of irregular window functions, which is often the case for realistic surveys.We proceed, instead, by expressing the second derivative of the spherical Bessel function as follows: Using this expression in eq.(A.5) gives k 2 χ 2 j +1 (kχ)j (kχ ) + 4f f k 4 χχ j +1 (kχ)j +1 (kχ ) , (A.7) which, after applying the decomposition of P (k; χ, χ ) as done in eq.(A.1), is fully determined by the integrals given in eq.(A.3).

B Approximating the CMB lensing convergence
In section 3.3, we introduced the approximation of using the full 3D Laplacian rather than its transverse part when computing the CMB lensing convergence.Here, we aim to provide further insight into why this is very accurate, even on large scales.We consider the full-sky case, so that eq.(3.17) becomes for the true convergence If instead we replace the transverse (i.e., angular) part of the 3D Laplacian, ∇ 2 ⊥ , with the full Laplacian, ∇ 2 , we introduce an error, ∆κ, involving the radial derivatives: It is instructive to consider the case of an Einstein-de Sitter universe, where the potential Ψ does not evolve in time.In that case, the radial derivatives in ∆κ can be replaced by total derivatives along the line of sight.Integrating by parts and using the explicit form for the lensing window function, eq.where the approximation assumes a scale-invariant power spectrum of primordial curvature perturbations with amplitude A s .At = 2, C ∆κ∆κ 2 ≈ 7.5 × 10 −10 compared to C κκ 2 ≈ 7.9 × 10 −8 , so the error from the power spectrum of ∆κ is sub-percent.The error term −2C κ∆κ from the cross-correlation between ∆κ and κ is also expected to be small, even on large scales, since the potential fluctuations on the last-scattering surface are weakly correlated with the fluctuations at lower χ that dominate the convergence.

. 20 )
which for large σ clearly constrains the support of the integral over k to the region |k | 1/σ.

Figure 1 :
Figure 1: Geometric layout of the full-sky (left panel) and flat-sky (right panel) unequal-time, two-point correlation function for points A and B separated by an angle θ and at radial distances χ and χ, respectively.The observer is at O. The line-of-sight axis n is the angle bisector in the full-sky geometry.Points A and B are corresponding points in the flat-sky set-up having the same 3D separation as A and B in the spherical case.The flat-sky angular coordinate is ω ≡ 2 sin(θ/2) and transverse separations are determined from ω through the effective distance χ g ≡ √ χχ .

Figure 2 :
Figure2: Angular power spectra of projected galaxy clustering with equal-time window functions.The left panel is for central redshifts z = z = 1.0 and width σ z = 0.05, while the right panel is for z = z = 2.0, with σ z = 0.3.In both panels, we compare the full-sky expression (black dash-dotted line), the Limber approximation (blue dashed line), and the flat-sky approximation (red solid line).The green solid line and cyan dashed line are geometrically recalibrated flat-sky and Limber approximations, as described in section 2.4.In the bottom panels we show fractional residuals compared to the full-sky results.

Figure 3 :
Figure 3: As figure 2 but for unequal-time window functions.The left panel is for central redshifts z = 1.0, z = 1.25, with σ z = 0.05 in both cases, while the right panel is for z = 2.0, z = 3.5, with σ z = 0.3.The true angular power spectrum is negative on large scales (the absolute value is plotted).

Figure 4 :
Figure4: Angular power spectra of projected galaxy clustering, including RSD, with equal-time window functions.The left panel is for central redshifts z = z = 1.0 and width σ z = 0.05, while the right panel is for z = z = 2.0, with σ z = 0.3.In both panels, we compare the full-sky calculation (black dash-dotted line), the extended Limber approximation (blue dashed line), and the flat-sky approximation without geometric recalibration (red solid line) and with recalibration (green solid line).

Figure 5 :
Figure 5: As figure 4 but for unequal-time window functions.The left panel is for central redshifts z = 1.0, z = 1.25, with σ z = 0.05, while the right panel is for z = 2.0, z = 3.5, with σ z = 0.3.

Figure 8 :
Figure 8: Angular power spectrum for the cross-correlation between the lensing convergence for a source plane at z = 2 and projected clustering with central redshift z = 0.5 and width σ z = 0.2.

2 )Figure 9 :
Figure 9: 2D plots of the unequal-time angular power spectrum C (χ, χ ) for projected clustering.The columns correspond to multipoles = 2, 10 and 100.The first row shows the flat-sky approximation C flat ( ), while the second row is the full-sky C full .The third and fourth rows are the flat-sky and full-sky C (χ, χ ), respectively, multiplied by radial window functions that are Gaussians centred at z = z = 1.0, with widths σ z = σ z = 0.05.Note the colour scales are logarithmic and the dynamic range is greatly expanded in the bottom two rows to probe the tails of the window functions.

14 )Figure 12 :
Figure 12: Slices through the unequal-time angular power spectrum C (χ, χ ) appropriate to the CMB lensing potential.The columns correspond to multipoles = 2, 10 and 100.The first row compares the flat-sky approximation C flat ( ) and the full-sky C full as a function of χ at δχ = 0.The second row multiplies these by the lensing window functions [so the quantity plotted has dimensions of (h/Mpc) 2 ].The third row shows slices (without multiplication by the window functions) at χ = 2381 Mpc/h (around z = 1) as a function of δχ.

Figure 13 :
Figure 13: Cumulative contribution to the angular power spectrum of the CMB lensing potential, C φφ , from χ ≤ χup .Results are shown for the full-sky (blue solid lines), Limber approximation (green dashed lines) and flat-sky approximation (eq.3.15; orange dot-dashed lines).Each panel represents a different multipole .