# Precision calculations of the cosmic shear power spectrum projection

Martin Kilbinger, Catherine Heymans, Marika Asgari, Shahab Joudaki, Peter Schneider, Patrick Simon, Ludovic van Waerbeke, Joachim Harnois-Déraps, Hendrik Hildebrandt, Fabian Köhlinger, et al.

► To cite this version:

Martin Kilbinger, Catherine Heymans, Marika Asgari, Shahab Joudaki, Peter Schneider, et al. Precision calculations of the cosmic shear power spectrum projection. Monthly Notices of the Royal Astronomical Society, 2017, 472 (2), pp.2126-2141. 10.1093/mnras/stx2082. hal-01645619

HAL Id: hal-01645619

https://hal.science/hal-01645619v1

Submitted on 4 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Precision calculations of the cosmic shear power spectrum projection

Martin Kilbinger,<sup>1,2★</sup> Catherine Heymans,<sup>3</sup> Marika Asgari,<sup>3</sup> Shahab Joudaki,<sup>4,5</sup>

Peter Schneider,<sup>6</sup> Patrick Simon,<sup>6</sup> Ludovic Van Waerbeke,<sup>7</sup> Joachim Harnois-Déraps,<sup>3</sup>

Hendrik Hildebrandt,<sup>6</sup> Fabian Köhlinger,<sup>8</sup> Konrad Kuijken<sup>9</sup> and Massimo Viola<sup>9</sup>

Accepted 2017 August 10. Received 2017 August 9; in original form 2017 February 20

## ABSTRACT

We compute the spherical-sky weak-lensing power spectrum of the shear and convergence. We discuss various approximations, such as flat-sky, and first- and second-order Limber equations for the projection. We find that the impact of adopting these approximations is negligible when constraining cosmological parameters from current weak-lensing surveys. This is demonstrated using data from the Canada-France-Hawaii Telescope Lensing Survey. We find that the reported tension with *Planck* cosmic microwave background temperature anisotropy results cannot be alleviated. For future large-scale surveys with unprecedented precision, we show that the spherical second-order Limber approximation will provide sufficient accuracy. In this case, the cosmic-shear power spectrum is shown to be in agreement with the full projection at the sub-percent level for $\ell > 3$, with the corresponding errors an order of magnitude below cosmic variance for all $\ell$. When computing the two-point shear correlation function, we show that the flat-sky fast Hankel transformation results in errors below two percent compared to the full spherical transformation. In the spirit of reproducible research, our numerical implementation of all approximations and the full projection are publicly available within the package NICAEA at http://www.cosmostat.org/software/nicaea.

**Key words:** methods: statistical – cosmological parameters.

## 1 INTRODUCTION

The measurement of weak gravitational lensing by large-scale structures provides a powerful cosmological probe of dark matter, dark energy and modifications to gravity. As such, it is the primary science goal of several current (KiDS, DES, HSC<fnref n="1" />) and future (Euclid, LSST, WFIRST<fnref n="2" />) surveys. Interest in the results from these surveys is high as statistically significant deviations have been found between the cosmological parameter constraints from the cosmic microwave background (CMB) *Planck* experiment (Planck Collaboration XIII 2016) in comparison to weak-lensing constraints from both the Kilo-Degree Survey (KiDS; Hildebrandt et al. 2017) and the Canada-France-Hawaii Telescope Lensing Survey (CFHTLenS; Joudaki et al. 2017). If the source of this tension is not a result of so-far unconsidered sources of systematic errors in one or all experiments, extensions to the standard flat ΛCDM cosmological models need to be considered. Joudaki et al. (2016) have shown, for example, that the tension can be resolved with an evolving dark energy model.

<footnote n="1"><sup>1</sup>CEA/Irfu/SAp Saclay, Laboratoire AIM, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="2"><sup>2</sup>Institut d'Astrophysique de Paris, UMR7095 CNRS, Université Pierre and Marie Curie, 98 bis boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="3"><sup>3</sup>Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="4"><sup>4</sup>Centre for Astrophysics and Supercomputing, Swinburne University of Technology, PO Box 218, Hawthorn, VIC 3122, Australia</footnote>

<footnote n="5"><sup>5</sup>ARC Centre of Excellence for All-sky Astrophysics (CAASTRO)</footnote>

<footnote n="6"><sup>6</sup>Argelander-Institut für Astronomie, Auf dem Hügel 71, D-53121 Bonn, Germany</footnote>

<footnote n="7"><sup>7</sup>Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada</footnote>

<footnote n="8"><sup>8</sup>Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan</footnote>

<footnote n="9"><sup>9</sup>Leiden Observatory, Leiden University, Niels Bohrweg 2, NL-2333 CA Leiden, the Netherlands</footnote>

<footnote>* E-mail: martin.kilbinger@cea.fr</footnote>

<footnote n="1">1 KiDS: http://kids.strw.leidenuniv.nl, DES: http://www.darkenergysurvey.org, HSC: http://hsc.mtk.nao.ac.jp/ssp</footnote>

<footnote n="2">2 Euclid: http://sci.esa.int/euclid, LSST: http://www.lsst.org, WFIRST: http://wfirst.gsfc.nasa.gov</footnote>
In the era of the upcoming large-scale surveys that will provide measurements of cosmic shear with unprecedented precision, one needs to revisit the theoretical predictions of the observables to ensure that the accuracy of the models meet the accuracy of the observations. In this paper, we examine the widely used Limber approximation for the projected weak-lensing power spectrum. We consider spherical coordinates and the flat-sky approximation, and compute the full projection of the lensing power spectrum. The first-order extended Limber approximation provides sub-percent accuracy for $\ell > 10$ and is sufficient for present surveys. The associated errors are sub-dominant even for future large surveys.

We further show that the second-order Limber approximation is an accurate representation of the full projection, with better than percent level precision for scales $\ell > 3$. Since this approximation involves only 1D integrals over the matter power spectrum, it is very fast to calculate numerically and can readily be employed in Monte Carlo sampling methods to obtain precision constraints on cosmological parameters. We also compute the shear correlation function using a spherical transformation, and compare this to the flat-sky approximated, commonly used fast Hankel transformation.

This paper is organized as follows. In Section 2, we provide a pedagogical introduction to weak gravitational lensing theory, projections and power spectra for the flat-sky and spherical cases, following the seminal work by Hu (2000; see also Castro, Heavens & Kitching 2005). In Section 3, we derive weak-lensing observables using a second-order Limber approximation first introduced by Loverde & Afshordi (2008). We compare the shear power spectrum and the commonly used two-point shear correlation function for the full solution to a range of different approximations in Section 4, providing cosmological constraints for each case using CFHTLenS data from Kilbinger et al. (2013). This paper draws from several sources of literature that have previously discussed the full projection, or reviewed the accuracy of the Limber approximation for weak-lensing studies, namely Schmidt (2008); Bernardeau et al. (2012); Giannantonio et al. (2012) and Kitching et al. (2016); see also Lemos, Challinor & Efstathiou (2017) for a more recent work. We present a discussion and comparison of our results to these papers in Appendix B.

## 2 WEAK-LENSING PROJECTIONS AND POWER SPECTRA

In this section, we review the basic weak-lensing projection expressions, and compute lensing power spectra for a spherical case, and in the flat-sky approximation. For completeness, we provide a derivation of the weak-lensing power spectra in Appendix A.

### 2.1 The lensing potential

The lensing potential $\psi$ at a position on the sky $(\theta, \varphi)$ in the Born approximation is defined as the projected 3D metric potential $\Phi$ along the line of sight of a flat universe (Kaiser 1998; Bartelmann & Schneider 2001),

$$ \psi(\theta, \varphi) = \frac{2}{c^2} \int_0^\infty \frac{d\chi}{\chi} \Phi(\chi, \chi\theta, \chi\varphi; \chi) q(\chi), \quad (1) $$

where the lensing efficiency $q$ is given by

$$ q(\chi) = \int_{\chi}^{\chi_h} d\chi' n(\chi') \frac{\chi' - \chi}{\chi'}, \quad (2) $$

corresponding to a population of lensed galaxies with a normalized source redshift distribution $n_z(z) dz = n(\chi) d\chi$, with the limit being the comoving distance to the horizon $\chi_h$.<fnref n="3" /> Here, $c$ is the speed of light and the projection is carried out over comoving distances $\chi$. The last argument of the potential $\Phi$ is not to be understood as coordinate, but rather as a substitute of cosmic time, $t(\chi)$, to express the time at which the potential is evaluated. This is true in the following for all fields and functions thereof that dynamically change with cosmic epoch.

The form of the lensing efficiency $q$ in equation (2) assumes a homogeneous galaxy distribution without clustering, so that the redshift distribution in this approximation does not depend on the direction on the sky. Accounting for this position-dependence leads to corrections of weak-lensing quantities due to clustering of source galaxies with other sources (Schneider, Van Waerbeke & Mellier 2002), and with galaxies associated with lens structures (Bernardeau 1998; Hamana et al. 2002).

The 3D potential is related to the density contrast $\delta$ via the Poisson equation. Assuming General Relativity, this relation is written in Fourier space as

$$ \hat{\Phi}(\mathbf{k}; \chi) = -\frac{3}{2} \Omega_m H_0^2 k^{-2} a^{-1}(\chi) \hat{\delta}(\mathbf{k}; \chi), \quad (3) $$

where $\Omega_m$ is the matter density parameter, $H_0$ the Hubble constant, $\mathbf{k}$ a 3D Fourier wave vector with modulus $k$ being the comoving wavenumber and $a$ the scalefactor with $a = 1$ today. The Fourier transform of the potential and its inverse are defined as

$$ \hat{\Phi}(\mathbf{k}; \chi) = \int d^3\mathbf{r} \, \Phi(\mathbf{r}; \chi) e^{i\mathbf{k}\cdot\mathbf{r}}; \quad (4) $$

$$ \Phi(\mathbf{r}; \chi) = \int \frac{d^3\mathbf{k}}{(2\pi)^3} \hat{\Phi}(\mathbf{k}; \chi) e^{-i\mathbf{k}\cdot\mathbf{r}}, \quad (5) $$

<footnote n="3">³ In (1), we have replaced this limit without loss of generality with $\infty$, since $q(\chi > \chi_h) = 0$.</footnote>
where the integration range for both integrals is $\mathbb{R}^3$.

## 2.2 Lensing power spectra in the spherical case

### 2.2.1 Lensing potential power spectrum

Following Hu (2000), we decompose the lensing potential (equation 1) into spherical harmonics, in analogy to the CMB temperature, both of which are scalar functions on the sphere. This decomposition and its inverse are

$$ \psi(\theta, \varphi) = \sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{\ell} \psi_{\ell m} Y_{\ell m}(\theta, \varphi); \quad (6) $$

$$ \psi_{\ell m} = \int_{S^2} d\Omega \, \psi(\theta, \varphi) Y_{\ell m}^*(\theta, \varphi). \quad (7) $$

Complex conjugation is denoted with the superscript *. To specify tomographic redshift bins $i = 1...N_z$, we introduce a family of lensing efficiency functions $q_i$ defined by a corresponding family of redshift distributions $n_i$ via equation (2). The resulting lensing potential is denoted by $\psi_{\ell m, i}$. The tomographic power spectrum of the lensing potential between two redshift bins $i$ and $j$, $C_{ij}^{\psi}$ (Peebles 1980) is then defined as

$$ \langle \psi_{\ell m, i} \psi_{\ell' m', j}^* \rangle = \delta_{\ell \ell'} \delta_{mm'} C_{ij}^{\psi}(\ell). \quad (8) $$

Note that the argument $\ell$ is a discrete integer variable, and is often written as index, $C_\ell$. Using the properties of the spherical harmonics (see Appendix A for details), the power spectrum can be written as

$$ C_{ij}^{\psi}(\ell) = \frac{8}{c^4 \pi} \int_0^{\infty} \frac{d\chi}{\chi} q_i(\chi) \int_0^{\infty} \frac{d\chi'}{\chi'} q_j(\chi') \int dk k^2 j_\ell(k\chi) j_{\ell'}(k\chi') P_\Phi(k; \chi, \chi') \quad (9) $$

$$ = \frac{8}{\pi} A^2 \int_0^{\infty} \frac{d\chi}{\chi} \frac{q_i(\chi)}{a(\chi)} \int_0^{\infty} \frac{d\chi'}{\chi'} \frac{q_j(\chi')}{a(\chi')} \int \frac{dk}{k^2} j_\ell(k\chi) j_{\ell'}(k\chi') P_m(k; \chi, \chi'), \quad (10) $$

where $j_\ell$ is the spherical Bessel function of order $\ell$. For convenience, we introduce the normalization constant $A$ as

$$ A = \frac{3}{2} \Omega_m \left( \frac{H_0}{c} \right)^2. \quad (11) $$

The first line expresses $C_{ij}^{\psi}$ in terms of the 3D potential power spectrum $P_\Phi$, defined as

$$ \langle \hat{\Phi}(\mathbf{k}; \chi) \hat{\Phi}^*(\mathbf{k}'; \chi') \rangle = (2\pi)^3 \delta_D(\mathbf{k} - \mathbf{k}') P_\Phi(\mathbf{k}; \chi, \chi'). \quad (12) $$

The second line uses the 3D matter density power spectrum $P_m$, which is defined analogously as

$$ \langle \hat{\delta}(\mathbf{k}; \chi) \hat{\delta}^*(\mathbf{k}'; \chi') \rangle = (2\pi)^3 \delta_D(\mathbf{k} - \mathbf{k}') P_m(\mathbf{k}; \chi, \chi'), \quad (13) $$

and is related to $P_\Phi$ via the absolute square of the Poisson equation (3).

This type of cross-power spectrum between different cosmological epochs $\chi$ and $\chi'$ was introduced in Castro et al. (2005). In Sections 3.1.1 and 3.1.2, this unequal-time cross-spectrum (Kitching & Heavens 2016) will be further evaluated and simplified in the context of the Limber approximation. The oscillating Bessel functions in equation (10) ensure that only relatively close epochs contribute to the lensing potential correlation. This make sense since observed light rays from two galaxies at different positions on the sky that necessarily converge at the observer today pick up the density fluctuations at similar times while propagating through the large-scale structure. A similar argument has been made in Bartelmann & Schneider (2001): since the matter power spectrum scales with $k$ for $k \to 0$, there is decreasing power towards larger and larger scales. In particular, the correlation of cosmic fields decreases strongly above a coherence scale $|\chi - \chi'| \gtrsim L_{\text{coh}}$, which is significantly smaller than the horizon scale $\chi_h$.

In the following section, we will discuss the relations between shear and convergence to the lensing potential on the sphere, and derive the power spectrum of the former two fields.

### 2.2.2 Shear power spectrum

The shear $\gamma = \gamma_1 + i\gamma_2$ is related to the potential at linear order by the trace-free part of the Jacobi matrix. The involved differential operator on the sphere is called the edth derivative, $\delta$ (see Castro et al. (2005) for an in-depth mathematical discussion of this concept). The edth operator $\delta$ ($\delta^*$) raises (lowers) the spin $s$ of an object. Twice applying this operator to the scalar (spin-0) potential creates the spin-2 shear:

$$ \begin{aligned} \gamma(\theta, \varphi) &= \frac{1}{2} \bar{\partial} \bar{\partial} \psi(\theta, \varphi); \\ \gamma^*(\theta, \varphi) &= \frac{1}{2} \bar{\partial}^* \bar{\partial}^* \psi(\theta, \varphi). \end{aligned} \quad (14) $$
<table><caption>**Table 1.** The shear power spectrum approximations studied in this paper. ‘ID’ is the label used in the text and figures. The sixth (seventh) column indicates the mode $\ell_x$ such that for $\ell \le \ell_x$, the approximated power spectrum is more accurate than $x$, with $x = 0.1$ (0.01).</caption><thead><tr><th>Case</th><th>ID</th><th>equation</th><th>p(ℓ)</th><th>ν(ℓ)</th><th>ℓ<sub>0.1</sub></th><th>ℓ<sub>0.01</sub></th><th>Comment</th></tr></thead><tbody><tr><td>1<sup>st</sup>-order standard Limber, flat</td><td>L1Fl</td><td>(43)+(i)</td><td>1</td><td>ℓ</td><td>4</td><td>60</td><td>Pre-2014 CFHTLenS and DLS</td></tr><tr><td>1<sup>st</sup>-order extended Limber, flat</td><td>ExtL1Fl</td><td>(43)+(ii)</td><td>(ℓ<sup>4</sup>/(ℓ+1/2)<sup>3</sup>)</td><td>ℓ + 1/2</td><td>20</td><td>200</td><td>Converges only with O(ℓ<sup>-1</sup>)</td></tr><tr><td>1<sup>st</sup>-order extended Limber, hybrid, flat</td><td>ExtL1FlHyb</td><td>(43)+(iii)</td><td>1</td><td>ℓ + 1/2</td><td>4</td><td>10</td><td>Post-2014 CFHTLenS, DES-SV and KiDS</td></tr><tr><td>1<sup>st</sup>-order extended Limber, spherical</td><td>ExtL1Sph</td><td>(38)</td><td>(ℓ<sup>2</sup>/(ℓ+1/2)<sup>4</sup>)</td><td>ℓ + 1/2</td><td>3</td><td>12</td><td></td></tr><tr><td>2<sup>nd</sup>-order extended Limber, flat</td><td>ExtL2Fl</td><td>(43)+(44)+(iii)</td><td>(ℓ<sup>4</sup>/(ℓ+1/2)<sup>4</sup>)</td><td>ℓ + 1/2</td><td>19</td><td>200</td><td>Converges only with O(ℓ<sup>-1</sup>)</td></tr><tr><td>2<sup>nd</sup>-order extended Limber, hybrid, flat</td><td>ExtL2FlHyb</td><td>(43)+(44)+(iii)</td><td>1</td><td>ℓ + 1/2</td><td>5</td><td>16</td><td>Best flat-sky approximation</td></tr><tr><td>2<sup>nd</sup>-order extended Limber, spherical</td><td>ExtL2Sph</td><td>(38)+(42)</td><td>(ℓ<sup>2</sup>/(ℓ+1/2)<sup>4</sup>)</td><td>ℓ + 1/2</td><td>2</td><td>5</td><td>Best approximation</td></tr><tr><td>Full spherical</td><td>FullSph</td><td>(21)</td><td>—</td><td>—</td><td>—</td><td>—</td><td>Correct projection</td></tr></tbody></table>

To write the shear on the sphere in terms of the lensing potential $\psi$, we insert the harmonics expansion of the potential (6). This requires the calculation of second derivatives of the spherical harmonic functions. This operation defines a new object, the spin-weighted spherical harmonic $s Y_{\ell m}$. The shear can be written on the sphere in terms of these functions as a spherical harmonics multipole expansion with coefficients $\pm 2\gamma_{\ell m}$. This expansion together with its inverse is

$$ (\gamma_1 \pm i\gamma_2)(\theta, \varphi) = \sum_{\ell m} \pm 2\gamma_{\ell m} Y_{\ell m}(\theta, \varphi); \quad (15) $$

$$ 2\gamma_{\ell m} = \int_{S^2} d\Omega \, \gamma(\theta, \varphi) Y_{\ell m}^*(\theta, \varphi); \quad -2\gamma_{\ell m} = \int_{S^2} d\Omega \, \gamma^*(\theta, \varphi) Y_{\ell m}^*(\theta, \varphi). \quad (16) $$

The spin-weighted spherical harmonics $s Y_{\ell m}$ that are the basis function in the expansion of the shear (equation 15) can be calculated via the relations

$$ t(\ell, s) Y_{\ell m}(\theta, \varphi) = \bar{\partial}^s Y_{\ell m}(\theta, \varphi); \quad t(\ell, s) Y_{-\ell m}^*(\theta, \varphi) = (-1)^s (\bar{\partial}^*)^s Y_{\ell m}(\theta, \varphi), \quad (17) $$

for $0 \le s \le \ell$, with the spin pre-factor (Bernardeau et al. 2012)

$$ t(\ell, s) = \sqrt{\frac{(\ell + s)!}{(\ell - s)!}}. \quad (18) $$

Inserting the lensing potential expansion (equation 6) into the expression for the shear (equation 14), and using equation (17) to compute the derivatives, we find for the shear expansion coefficients (Hu 2000; Taylor 2001)

$$ \pm 2\gamma_{\ell m} = \frac{1}{2} t(\ell, 2) \psi_{\ell m}. \quad (19) $$

The two coefficients $\pm 2\gamma_{\ell m}$ and $-2\gamma_{\ell m}$ are identical since the potential $\psi$ is a real function.

The tomographic shear power spectrum, in analogy to equation (10), is defined by

$$ \langle {}_2\gamma_{\ell m, i} {}_2\gamma_{\ell' m', j}^* \rangle = \delta_{\ell \ell'} \delta_{mm'} C_{ij}^{\gamma}(\ell). \quad (20) $$

This is given by

$$ C_{ij}^{\gamma}(\ell) = \frac{1}{4} t^2(\ell, 2) C_{ij}^{\psi}(\ell) = \frac{2}{\pi} t^2(\ell, 2) A^2 \int_0^{\infty} \frac{d\chi}{\chi} \frac{q_i(\chi)}{a(\chi)} \int_0^{\infty} \frac{d\chi'}{\chi'} \frac{q_j(\chi')}{a(\chi')} \int_0^{\infty} \frac{dk}{k^2} P_m(k, \chi, \chi') j_{\ell}(k\chi) j_{\ell'}(k\chi'), \quad (\text{FullSph}) \quad (21) $$

where we use the label ‘FullSph’ (see Table 1 for a list of cases discussed in this work). The spherical spin pre-factor for the full projection shear power spectrum is $t^2(\ell, 2)$, which will be modified under flat-sky and Limber approximations below.

### 2.2.3 Convergence power spectrum

The convergence is related to the lensing potential on the sphere via the product of spin-raising and spin-lowering edth operators, which are identical to the spherical Laplacian differential operator.

$$ \kappa(\theta, \varphi) = \frac{1}{2} \bar{\partial} \bar{\partial}^* \psi(\theta, \varphi) = \frac{1}{2} \nabla^2 \psi(\theta, \varphi). \quad (22) $$

The spherical harmonics are eigenfunctions of the Laplacian,

$$ \nabla^2 Y_{\ell m}(\theta, \varphi) = -\ell(\ell + 1) Y_{\ell m}(\theta, \varphi) = -t^2(\ell, 1) Y_{\ell m}(\theta, \varphi). \quad (23) $$
The convergence power spectrum is then similar to the shear power spectrum (equation 20) with a different spherical pre-factor (Hu 2000; Joudaki & Kaplinghat 2012):

$$C_{ij}^{\kappa}(\ell) = \frac{1}{4} t^4(\ell, 1) C_{ij}^{\psi}(\ell) = \frac{\ell(\ell + 1)}{(\ell - 1)(\ell + 2)} C_{ij}^{\gamma}(\ell). \quad (24)$$

The convergence power spectrum is thus larger than the shear power spectrum, by 10 per cent for $\ell = 4$, 1 per cent for $\ell = 14$ and less than 0.1 per cent for $\ell > 45$.

## 2.3 Flat-sky approximation

The majority of cosmic shear analyses have used the predicted lensing power spectrum approximated on a flat sky, neglecting the sky curvature. This is a valid approach for past and current survey areas with an extent less than 10°. To account for the sky curvature of the observed data, the shear correlation functions from observed galaxy ellipticities are now routinely computed using spherical coordinates, since projecting to a Cartesian plane has been shown to cause significant biases of the two-point function on large scales (Fu et al. 2008), and lead to spurious B-modes (Asgari et al. 2017). Here, we examine the effect of sky curvature on the theoretical models of the power spectrum, and the effect on cosmological parameter inference (see Section 4.3).

For a flat-sky, the spherical harmonic expansions are approximated by Fourier transforms. The flat-sky equivalents of equations (6) and (7) are

$$\psi(\vartheta) = \int \frac{d^2\ell}{(2\pi)^2} e^{-i\ell\cdot\vartheta} \hat{\psi}(\ell), \quad (25)$$

$$\hat{\psi}(\ell) = \int d^2\vartheta e^{i\ell\cdot\vartheta} \psi(\vartheta), \quad (26)$$

where $\vartheta = (\theta, \varphi)$ is the vector describing a 2D angle on the sky. Instead of a harmonics coefficients $\psi_{\ell m}$, the Fourier representation of the potential $\hat{\psi}$ now depends on the vector $\ell \in \mathbb{R}^2$.

The flat-sky power spectrum, i.e. the flat-sky analogue of equation (8), is defined by

$$\langle \hat{\psi}_i(\ell) \hat{\psi}_j^*(\ell') \rangle = (2\pi)^2 \delta_D(\ell - \ell') P_{ij}^{\psi}(\ell). \quad (27)$$

Hu (2000) has shown that for small angles the harmonics expansion (equation 6) can be approximated by the Fourier representation in equation (25). They also demonstrated that the power spectra are approximately equal, $P^{\psi} \approx C^{\psi}$.

For a spin-2 field, Hu (2000) approximates the edth operator by Cartesian derivatives, and approximates equation (17) as

$$\ell^2 \pm 2Y_{\ell m}(\theta, \varphi) \approx e^{\mp 2i\psi_{\ell}} (\partial_1 \pm i\partial_2)^2 Y_{\ell m}(\theta, \varphi). \quad (28)$$

The spin pre-factor $t(\ell, 2) = \sqrt{(\ell - 1)\ell(\ell + 1)(\ell + 2)}$ is replaced by $\ell^2$ in flat coordinates, an approximation that holds for large $\ell$, since sky curvature can be neglected for small angular scales. We find for the flat-sky shear power spectrum

$$P_{ij}^{\gamma}(\ell) = \frac{2}{\pi} \ell^4 A^2 \int_0^{\infty} \frac{d\chi}{\chi} \frac{q_i(\chi)}{a(\chi)} \int_0^{\infty} \frac{d\chi'}{\chi' a(\chi')} \int_0^{\infty} \frac{dk}{k^2} P_m(k, \chi, \chi') j_{\ell}(k\chi) j_{\ell}(k\chi'), \quad (29)$$

with flat-sky pre-factor $\ell^4$. See Appendices B4 and B5 for discussions of alternative expressions for the flat-sky power spectrum.

# 3 SECOND-ORDER LIMBER APPROXIMATION FOR WEAK-LENSING

## 3.1 Spherical case

We follow Loverde & Afshordi (2008), who derive the second-order Limber expansion for general projections from 3D to 2D scalar fields in the spherical, all-sky case. We apply their general derivation to the case of weak lensing, and in contrast to Loverde & Afshordi (2008) account for a time-dependent power spectrum using two approaches presented in Sections 3.1.1 and 3.1.2.

First, we use the identity of Bessel functions

$$j_{\ell}(x) = \sqrt{\frac{\pi}{2x}} J_{\ell+1/2}(x) \quad (30)$$

in equation (21), where $J_\ell$ is the Bessel function of the first kind and order $\nu$. Next, Loverde & Afshordi (2008) solve integrals of the form

$$\int_0^\infty d\chi f(\chi) J_\nu(k\chi) = \int_0^\infty dx k^{-1} f(x/k) J_\nu(x) \quad (31)$$

by performing a Taylor expansion of an arbitrary differentiable function $f$ around $x = k\chi = \nu = \ell + 1/2$, where the Bessel function has its approximate maximum.
### 3.1.1 Geometric mean cross-correlation power spectrum

To separate the $k$ and $\chi$, $\chi'$-terms in equation (21), we first approximate the matter power cross-spectrum between two distances by the geometric mean of the two involved distances (Castro et al. 2005; Kitching & Heavens 2016):

$$P_m(k; \chi, \chi') = \sqrt{P_m(k; \chi) P_m(k; \chi')}. \quad (32)$$

This form is justified when considering the linear power spectrum, and follows when inserting the linearly evolving density contrast $\hat{\delta}(k; \chi) = D_+(\chi)\hat{\delta}_0(k)$ into equation (13), where $\delta_0$ is the present-day linearly extrapolated density contrast, and $D_+$ the linear growth factor with $D_+(0) = 1$. This is a good approximation also in the non-linear case as shown in Kitching & Heavens (2016).

With this, equation (21) is written as

$$C_{ij}^{\gamma}(\ell) \approx t^2(\ell, 2) A^2 \int_0^{\infty} \frac{dk}{k^3} \int_0^{\infty} \frac{d\chi}{\chi^{3/2}} \sqrt{P_m(k; \chi)} \frac{q_i(\chi)}{a(\chi)} J_{\ell+1/2}(k\chi) \int_0^{\infty} \frac{d\chi'}{\chi'^{3/2}} \sqrt{P_m(k; \chi')} \frac{q_j(\chi')}{a(\chi')} J_{\ell+1/2}(k\chi'). \quad (33)$$

Note that this equation has a pre-factor $t^2(\ell, 2)$ corresponding to a spin-2 field, in contrast to Loverde & Afshordi (2008), who show calculations for a scalar field.

Following Loverde & Afshordi (2008), we expand to third order

$$\lim_{\epsilon \to 0} \int_0^\infty dx \, e^{-\epsilon(x-\nu)} g(x) J_\nu(x) \approx g(\nu) - \frac{1}{2} \left. \frac{d^2 g}{dx^2} \right|_{x=\nu} - \frac{\nu}{6} \left. \frac{d^3 g}{dx^3} \right|_{x=\nu}, \quad (34)$$

with $g(x) = k^{-1}f(k, \chi)$, $\chi = x/k$ and its derivatives $g^{(n)}(x) = k^{-1-n}f^{(n)}(k, \chi)$ for a given $k$, where the derivatives of $f$ are with respect to the second argument $\chi$. In this series expansion, the replacement of the integral with the evaluation of $g$ and its derivatives at the maximum of the Bessel function is a good approximation of the integral if $g$ is varying more slowly than the oscillating Bessel function. As we will show below, $f$ is a slowly varying function of the comoving distance. In our case, the projection kernel is

$$f(k, \chi) = \sqrt{P_m(k; \chi)} a^{-1}(\chi) \chi^{-3/2} q(\chi). \quad (35)$$

In the tomographic case, the index $i$ is added to $q$ and $f$. Replacing both distance integrals in equation (33) by their Taylor-expansions around the maxima $\nu(\ell) = \ell + 1/2$ of the two Bessel functions, which are $k\chi$ and $k\chi'$, respectively, yields

$$C_{ij}^{\gamma}(\ell) \approx t^2(\ell, 2) A^2 \int_0^{\infty} \frac{dk}{k^3} k^{-2} \left[ f_i(k, \chi) - \frac{1}{2k^2} f_i''(k, \chi) - \frac{\nu(\ell)}{6k^3} f_i'''(k, \chi) + \dots \right] \left[ f_j(k, \chi) - \frac{1}{2k^2} f_j''(k, \chi) - \frac{\nu(\ell)}{6k^3} f_j'''(k, \chi) + \dots \right]. \quad (36)$$

Changing the integration to $\chi = \nu/k$ and collecting terms according to their $\nu$-dependence:

$$C_{ij}^{\gamma}(\ell) \approx C_{L1,ij}^{\gamma}(\ell) + C_{L2,ij}^{\gamma}(\ell) = \frac{t^2(\ell, 2)}{\nu^4(\ell)} A^2 \int_0^{\infty} d\chi \chi^3 \left\{ (f_i f_j)(\nu(\ell)/\chi, \chi) - \frac{1}{\nu^2(\ell)} \left[ \frac{\chi^2}{2} (f_i f_j'' + f_i'' f_j) (\nu(\ell)/\chi, \chi) + \frac{\chi^3}{6} (f_i f_j''' + f_i'' f_j') (\nu(\ell)/\chi, \chi) \right] + \mathcal{O}(\nu^{-4}) \right\}. \quad (37)$$

The first term corresponds to the well-known first-order Limber approximation (Limber 1953; Kaiser 1992), which is widely used in weak gravitational lensing. We retrieve the (spherical) standard expression by inserting back the projection kernel (equation 35),

$$C_{L1,ij}^{\gamma}(\ell) = \frac{t^2(\ell, 2)}{\nu^4(\ell)} A^2 \int d\chi \frac{q_i(\chi)q_j(\chi)}{a^2(\chi)} P_m\left(\frac{\nu(\ell)}{\chi}; \chi\right). \quad (\text{ExtL1Sph}) \quad (38)$$

In the Limber approximation, modes between structures at different epochs do not contribute to the single line-of-sight integration.

The second-order Limber term in equation (37) has an additional $\nu^{-2}$-dependence, and is therefore strongly suppressed for large $\ell$,

$$C_{L2,ij}^{\gamma}(\ell) = - \frac{1}{\nu^2(\ell)} \frac{t^2(\ell, 2)}{\nu^4(\ell)} \frac{A^2}{2} \int d\chi \chi^{7/2} a^{-1} P_m^{1/2}\left(\frac{\nu(\ell)}{\chi}; \chi\right) \left[ q_i f_j'' + f_i'' q_j + \frac{\chi}{3} (q_i f_j''' + f_i'' f_j) \right] (\nu(\ell)/\chi, \chi). \quad (39)$$

The higher order derivatives of the filter functions have to be computed numerically in the general case. These suffer from numerical noise and are sensitive to the set-up, for example the step size. The tabulation and interpolation of those derivatives is time-consuming since they depend on two arguments, $\nu$ and $\chi$. In the following section, we will separate the $k$- and $\chi$-dependent parts of the power spectrum to make the numerical derivatives faster and more smooth.

### 3.1.2 Approximation for small $\ell$

To further develop equation (32), we divide out the growth factor of the power spectrum,

$$P_m(k, \chi) := D_+^2(\chi) \tilde{P}_m(k, \chi) \approx D_+^2(\chi) \tilde{P}_m(k). \quad (40)$$

The function $\tilde{P}_m$ is, in general, not time-independent, except in the case of a linear matter power spectrum in the absence of massive neutrinos, and in General Relativity. However, the second-order Limber terms are expected to be important only for small $\ell$, since for large $\ell$, the
first-order Limber equation (38) is dominant. Equation (40) becomes a good approximation for small $ℓ$, since that means either small $k$, where the evolution is linear, or small $\chi$, where the lensing efficiency is small. The accuracy of a very nearby tomographic bin with low mean redshift should be further examined, but this is not the case for CFHTLenS.

With the definition (40), we factor out the function $\tilde{P}_m(k, \chi)$ from equation (35), and we define the separated kernel function $f_s$ as

$$f_s(\chi) = D_+(\chi)a^{-1}(\chi)\chi^{-3/2}q(\chi). \qquad (41)$$

The Limber equation up to second order of the shear power spectrum can then be approximated as

$$C_{ij}^{\gamma}(\ell) = C_{L1,ij}^{\gamma}(\ell) - \frac{1}{\nu^2(\ell)} \frac{\ell^2(\ell, 2)}{\nu^4(\ell)} \frac{A^2}{2} \int d\chi \chi^{7/2} a^{-1}(\chi) D_{+}^{-1}(\chi) P_m \left( \frac{\nu(\ell)}{\chi}; \chi \right) \left[ q_i f_{sj}'' + f_{si}'' q_j + \frac{\chi}{3} (q_i f_{sj}''' + f_{si}''' q_j) \right] (\chi). \quad (\text{ExtL2Sph}) \quad (42)$$

We compute the numerical higher derivatives as follows. The functions $f_s(\chi)$ are fitted as power laws with index $\approx -1.5$, which is expected if the growth suppression factor $D_+(a)/a$ varies slowly with $\chi$, and the lensing efficiency $q \approx 1$ for small and medium $\chi$, given the CFHTLenS redshift range. The fit is carried out between $\chi_{min} = 0.001$ Mpc h$^{-1}$ and $\chi_{max} = 500$ Mpc h$^{-1}$. At larger comoving distances, the kernel decreases faster than a power law, so we exclude this range from the fit. Even though on those scales the derivatives are larger due to the steeper decline, the associated errors are very small as these scales are down-weighed by the kernel function $f_s$ itself. At $\chi = 500$ Mpc h$^{-1}$, the filter function is less than $10^{-4}$ of its value at 1 Mpc h$^{-1}$.

## 3.2 Flat-sky

The extended flat-sky Limber approximation is readily derived from the spherical case equations (38, 42), by replacing the pre-factor $\ell^2(\ell, 2)$ with $\ell^4$,

$$P_{L1,ij}^{\gamma}(\ell) = p(\ell) A^2 \int d\chi \frac{q_i(\chi)q_j(\chi)}{a^2(\chi)} P_m \left( \frac{\nu(\ell)}{\chi}; \chi \right); \qquad (43)$$

$$P_{L2,ij}^{\gamma}(\ell) = -\frac{1}{\nu^2(\ell)} p(\ell) \frac{A^2}{2} \int d\chi \chi^{7/2} a^{-1}(\chi) D_{+}^{-1}(\chi) P_m \left( \frac{\nu(\ell)}{\chi}; \chi \right) \left[ q_i f_{sj}'' + f_{si}'' q_j + \frac{\chi}{3} (q_i f_{sj}''' + f_{si}''' q_j) \right] (\chi) \quad (44)$$

Further approximations can be made for the pre-factor $p(\ell) = \ell^4/\nu^4(\ell)$ and $\nu(\ell)$:

(i) $p(\ell) = 1$. This corresponds to $\nu(\ell) = \ell$, which is the standard Limber approximation (L1Fl). Until recently, i.e. for all pre-2014 CFHTLenS results and DLS (Jee et al. 2013) analyses, this was the approximation of choice. Note that we do not discuss the second-order Limber approximation with $p(\ell) = 1$.

(ii) $p(\ell) = \ell^4/(\ell + 1/2)^4$. This corresponds to the extended Limber approximation (ExtL1Fl, ExtL2Fl) with $\nu(\ell) = \ell + 1/2$; however the following case is typically employed.

(iii) $p(\ell) = 1$, but keeping as argument of the power spectrum $\nu(\ell) = \ell + 1/2$. This is a hybrid between standard and extended Limber approximation (ExtL1FlHyb, ExtL2FlHyb), and the first-order case was used in Hildebrandt et al. (2017); Joudaki et al. (2017, 2016) and Abbott et al. (2016). As is shown below, this is a better approximation to the full projection than case (ii). In equation (44), the second-order suppression factor is also left to be $\nu^{-2}(\ell) = (\ell + 1/2)^{-2}$, providing a slightly more accurate approximation compared to $\nu^{-2}(\ell) = \ell^{-2}$.

# 4 RESULTS

## 4.1 Comparison of the approximations for the lensing power spectrum

In Fig. 1, we present the full spherical projection of the shear power spectrum in comparison to shear power spectra derived assuming the range of different approximations listed in Table 1. The adopted redshift distribution corresponds to CFHTLenS (Kilbinger et al. 2013) and we assume their best-fitting flat $\Lambda$CDM cosmology with $\Omega_m = 0.279$, $\Omega_b = 0.046$, $\sigma_8 = 0.79$, $h = 0.701$, $n_s = 0.96$. For $\ell > 100$, we find that all shear power spectra predictions agree with the full spherical solution to better than one percent, with the majority of the approximations tested accurate to better than 0.1 per cent.

Considering first the flat-sky cases, the standard first-order Limber approximation, L1Fl, that was adopted for all pre-2014 CFHTLenS analyses and DLS analyses, we find it to be accurate to better than 10 per cent for $\ell > 3$, converging slowly to the true projection with percent level precision at $\ell > 100$. For the extended Limber approximations ‘hybrid’ cases (ExtL1FlHyb and ExtL2FlHyb), despite decreased accuracy for $\ell < 8$ in comparison to the standard first-order Limber case, the errors with respect to the true power spectrum decrease much faster, as $\ell^{-2}$, such that percent-level precision is reached at $\ell > 15$. The first-order extended Limber approximation ‘hybrid’ cases (ExtL1FlHyb) was adopted by Joudaki et al. (2017), Joudaki et al. (2016), DES-SV (Abbott et al. 2016) and Hildebrandt et al. (2017).<fnref n="4" />

<footnote n="4">⁴ We confirm that there is a typographical error in equation (4) of Hildebrandt et al. (2017) and in equation (2) of Abbott et al. (2016) (private communication with Joe Zuntz), which should include the extra term of ‘+0.5’ in $\nu(\ell)$ that was incorporated in both cosmological analyses.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The shear power spectrum for different approximations listed in Table 1. Limber to first order: standard with flat-sky (L1Fl), extended for flat sky (ExtL1Fl), extended hybrid for flat sky (ExtL1FlHyb) and extended in the spherical expansion (ExtL1Sph); second-order Limber approximations: extended flat sky (ExtL2Fl), extended hybrid flat sky (ExtL2FlHyb) and extended spherical expansion (ExtL2Sph); full (exact) spherical projection (FullSph). The left-hand panel shows the total shear power spectrum. The right-hand panel shows the fractional difference resulting from each approximation, relative to the full spherical projection of the shear power spectrum. The two light grey curves on the top show the cosmic variance for KiDS- and Euclid-like surveys with areas of 1 500 and 15 000 deg$^2$, respectively.</figcaption></figure>

The outlier in the flat-sky cases is the extended Limber approximation (ExtL1Fl), which performs relatively poorly, and reaches 10 per cent precision only at $\ell > 100$. The same slow convergence can be observed for the corresponding second-order flat case, ExtL2Fl. To our knowledge, this form of the flat-sky approximation has not been used in any cosmic shear studies to date, and should not be used in any future studies given these results. The poor behaviour of this case, in comparison with the ‘hybrid’ case (e.g. ExtL1FlHyb) can be understood by considering Taylor expansions of the different pre-factors. The spherical pre-factor,

$$p(\ell) = \frac{(\ell + 2)(\ell + 1)\ell(\ell - 1)}{(\ell + 0.5)^4} = 1 - \frac{5}{2\ell^2} + \mathcal{O}(\ell^{-3}), \qquad (45)$$

can be compared to the flat-sky extended Limber pre-factor

$$p(\ell) = \ell^4/(\ell + 0.5)^4 = 1 - \frac{2}{\ell} + \frac{5}{2\ell^2} + \mathcal{O}(\ell^{-3}), \qquad (46)$$

showing it to be more deviant from the full spherical solution, than the ‘hybrid’ $p(\ell) = 1$ case.

Considering now the spherical-sky cases, we find that including the spherical pre-factor decreases the difference between the extended Limber approximated cases (ExtL1Sph and ExtL2Sph) and the full spherical solution by a factor of a few for $\ell < 5$. We find that using the spherical-sky second-order extended Limber approximation (ExtL2Sph) yields percent accuracy down to $\ell > 3$. The numerical calculation of the second-order extended Limber approximation is a factor of 15 times faster than the calculation of the full spherical solution (averaged over the first 18 $\ell$-modes). We note that the sub-0.1 per cent-fluctuations seen in the right-hand panel of Fig. 1 is due to numerical noise arising from numerical integration errors in the calculation of the full spherical solution when $\ell > 20$.

We note that in all cases the second-order Limber expansion adds power to the first-order term. In the flat hybrid case, which overestimates the full spherical solution, this results in the second-order expansion being less precise compared to first-order.

Compared to the statistical power of current and future surveys, all approximations discussed here are subdominant to the cosmic variance, $\Delta C(\ell)/C(\ell) = [f_{sky}(2\ell + 1)]^{-1/2}$ (Kaiser 1992), where $f_{sky}$ is the fraction of sky observed by the survey. The uncertainties from the Limber approximation in the case of ExtL2Sph is an order of magnitude below the cosmic variance of a Euclid-like survey (sky area 15 000 deg$^2$) for all $\ell$.

## 4.2 Effects on the shear correlation function

The majority of cosmic shear analyses to date have adopted real-space correlation statistics, since these can be measured directly from an observed galaxy shape catalogue. The baseline quantity is the two-point correlation function (Miralda-Escude 1991; Kaiser 1992; Bartelmann & Schneider 2001), given in the flat-sky approximation by

$$\xi_{+}(\theta) = \langle \gamma \gamma^{*} \rangle (\theta) = \frac{1}{2\pi} \int d\ell \, \ell \, J_{0}(\ell\theta) P^{\gamma}(\ell); \quad \xi_{-}(\theta) = \langle \gamma \gamma \rangle (\theta) = \frac{1}{2\pi} \int d\ell \, \ell \, J_{4}(\ell\theta) P^{\gamma}(\ell). \qquad (47)$$
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The difference of the two-point shear correlation functions $\xi_+$ (left) and $\xi_-$ (right) of the ExtL2Sph projection relative to the full spherical case (FullSph). Two cases of the shear correlation function transformation for ExtL2Spha are shown, the full spherical case (equation 48, green solid lines), and the flat-sky Hankel transform (equation 47, red dashed curves).</figcaption></figure>

The flat-sky shear power spectrum $P^\gamma$ can be related to the underlying matter power spectrum through equation (43) when adopting a first-order extended Limber approximation, or equation (44) when adopting a second-order extended Limber approximation.

On a sphere, correlation functions formally cannot be related to the power spectrum by the Hankel transform in equation (47), and should be replaced by the spherical transform (Ng & Liu 1999; Chon et al. 2004)

$$ \xi_+(\theta) = \frac{1}{4\pi} \sum_{\ell=2}^{\infty} (2\ell + 1)C^\gamma(\ell) d_{2,\ell}^+(\theta); \quad \xi_-(\theta) = \frac{1}{4\pi} \sum_{\ell=2}^{\infty} (2\ell + 1)C^\gamma(\ell) d_{-2,\ell}^-(\theta), \qquad (48) $$

where $d_{m,n}^\ell$ are the reduced Wigner D-matrices (see Appendix C for details on their numerical calculation).

The spherical power spectrum is formally not defined for non-integer $\ell$ (see Castro et al. 2005, for alternative spherical-sky formulae for the two-point correlation function), as functions defined on the sphere are necessarily periodic. As we have shown in Section 4.1, however, the spherical second-order extended Limber approximation provides a percent-level precision representation of the full spherical projection for $\ell > 3$. The spherical pre-factor $l(l, 2)$ (equation 18) can be generalized to non-integer arguments and is positive for $\ell \ge 2$. We can thus use the spherical power spectrum with the Hankel transformation in equation (47) to compute the two-point correlation functions. This has the advantage being able to employ fast FFT numerical implementations of the Hankel transforms (Hamilton 2000), when Monte Carlo sampling.

We compare predictions for the two-point shear correlation function using the Hankel transformation and the full spherical transformation in Fig. 2. We show the full projection and the best approximation, ExtL2Sph. Note that for the ‘FullSph’ case we replace the full projection with ExtL2Sph for $\ell > 200$ to reduce computation time. We find that the Hankel transform (equation 47) is accurate to better than 5 (0.2) percent for $\xi_+$ ($\xi_-$). The difference between the second-order Limber and full projection solution using the spherical transform (equation 48) is well below one (0.03) percent on scales of $\vartheta < 300$ arcmin for $\xi_+$ ($\xi_-$). The red lines in Fig. 2 present the limit of precision that we can achieve with the current fast Hankel transform implementation for the correlation function. The green lines show the limit of the second-order Limber approximation on the correlation function.

Fig. 3 shows the two-point correlation functions $\xi_+$ (left) and $\xi_-$ (right) using the different cases for the shear power spectrum listed in Table 1. The component $\xi_+$ is calculated with the appropriate transformation, i.e. Hankel for the flat cases, and spherical involving Legendre polynomials for the spherical cases. For $\xi_-$, which is less dominated by large scales and thus Limber and flat-sky approximations, in comparison to $\xi_+$, we use the Hankel transform in all cases. In addition, for $\xi_-$, the approximation ‘ExtL2Sph’ is our reference case. The adopted CFHTLenS redshift distribution and fiducial cosmological model are the same as in Fig. 1. As is clear by the red dotted curve, using the Planck cosmology (Planck Collaboration XIII 2016) induces a significant change in the amplitude of the shear correlation function in comparison to the different projection methods.

## 4.3 Application to CFHTLenS data

The CFHTLenS represented a major step forward for the field of weak gravitational lensing, in terms of improved accuracy in data reduction (Erben et al. 2013), the implementation of PSF-Gaussianized matched multiband photometry (Hildebrandt et al. 2012), cross-correlation clustering analysis between photometric redshift slices to verify tomographic redshift distributions (Benjamin et al. 2013), accurate calibrated
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The two-point shear correlation functions $\xi_+$ (left) and $\xi_-$ (right). In the spherical cases (ExtL1Sph, ExtL2Sph, FullSph), $\xi_+$ and $\xi_-$ have been computed using the spherical transform (equation 48). For the flat cases, the Hankel transform in equation (47) was used. The upper panels show the total shear correlation functions for the range of cases listed in Table 1. The lower panels show the relative differences to the spherical sky second-order extended Limber approximation, ExtL2Sph. The theoretical models correspond to the CFHTLenS best-fitting cosmological parameters with $\Omega_m = 0.279$, $h = 0.701$ and $\sigma_8 = 0.79$ (Kilbinger et al. 2013). For comparison, we also show, in the upper panels, the spherical sky second-order extended Limber approximation model for the Planck-best fitting cosmology with $\Omega_m = 0.3$, $h = 0.67$ and $\sigma_8 = 0.83$ (Planck Collaboration XIII 2016).</figcaption></figure>

shape measurements (Miller et al. 2013) and a full suite of informative systematic tests to select a clean data set (Heymans et al. 2012). Since the public release of this survey in 2013, the community has continued to scrutinize and advance our understanding of CFHTLenS by identifying a number of areas where analyses could improve:

(i) Choi et al. (2016) identified significant biases in the tomographic photometric redshift distributions using a more effective clustering analysis, in comparison to Benjamin et al. (2013), by incorporating newly overlapping spectroscopic data from the Sloan Digital Sky Survey. The conclusion of this work was that any re-analysis of CFHTLenS should include systematic error terms to account for bias and scatter, with a prediction that accounting for these biases would *reduce* the recovered amplitude of $\sigma_8$ by ~4 per cent. Additional new techniques to calibrate the redshift distribution of tomographic bins were introduced recently in Hildebrandt et al. (2017).

(ii) The CFHTLenS tomographic cosmological analysis was then revisited by Joudaki et al. (2017) in order to include a full redshift error analysis based on the results from Choi et al. (2016). The impact of correcting for these biases, including their associated errors, served to reduce the overall constraining power of the survey and hence also the tension between CFHTLenS and CMB constraints.

(iii) Asgari et al. (2017) used the stringent COSEBI statistic (Schneider, Eifler & Krause 2010) to identify significant non-lensing B-mode distortions when the CFHTLenS data was split into tomographic slices.
<table><caption>**Table 2.** Mean and 68 per cent credible interval for $\sigma_8(\Omega_m/0.27)^{0.6}$ and $\sigma_8(\Omega_m/0.3)^{0.6}$ for various approximations to the lensing power spectrum projections listed in Table 1.</caption><thead><tr><th>ID</th><th>σ<sub>8</sub>(Ω<sub>m</sub>/0.27)<sup>0.6</sup></th><th>σ<sub>8</sub>(Ω<sub>m</sub>/0.3)<sup>0.6</sup></th></tr></thead><tbody><tr><td>L1Fl</td><td>0.787<sup>+0.031</sup><sub>-0.033</sub></td><td>0.739<sup>+0.029</sup><sub>-0.031</sub></td></tr><tr><td>ExtL1Fl</td><td>0.792 ± 0.032</td><td>0.744 ± 0.030</td></tr><tr><td>ExtL1FlHyb</td><td>0.788<sup>+0.031</sup><sub>-0.033</sub></td><td>0.740<sup>+0.029</sup><sub>-0.031</sub></td></tr><tr><td>ExtL2FlHyb</td><td>0.788<sup>+0.031</sup><sub>-0.033</sub></td><td>0.740<sup>+0.029</sup><sub>-0.031</sub></td></tr><tr><td>ExtL2Sph(Hankel)</td><td>0.789<sup>+0.031</sup><sub>-0.032</sub></td><td>0.740<sup>+0.029</sup><sub>-0.030</sub></td></tr></tbody></table>

(iv) Kuijken et al. (2015) showed that the CFHTLenS shear calibration corrections derived in Miller et al. (2013) were underestimated as a result of an imperfect match between the galaxy populations in the data and image simulations.

(v) Fenech Conti et al. (2017) demonstrated that the CFHTLenS data would have been subject to a weight bias that favours galaxies that are more intrinsically oriented with the point-spread function. They also showed that the impact of calibration selection biases, which were not considered in Miller et al. (2013), would have lead to the overcorrection of multiplicative shear bias in the CFHTLenS analyses, by a few percent.

(vi) Joudaki et al. (2017) updated the CFHTLenS covariance matrices using larger box numerical simulations that were less subject to the lack of power on large scales. A complementary accurate estimate of the covariance matrix using analytical methods will be published soon (Joachimi et al., in preparation)

(vii) Takahashi et al. (2012) provided a more accurate non-linear power spectrum correction than that used in the original CFHTLenS analyses, and the halo model from Mead et al. (2015) allowed for simultaneous modelling of baryonic modifications to the non-linear power spectrum.

All these advances in our understanding were incorporated and accounted for in the recent KiDS cosmic shear analysis (Hildebrandt et al. 2017), which reports a 2.3σ tension with *Planck*. Efforts are now underway to fully re-analyse CFHTLenS using the advanced KiDS analysis pipeline with revised shape measurements and calibrations for the shear and photometric redshifts. Until this analysis is complete, we note that these known shortcomings with the original CFHTLenS results impact in different ways the cosmological conclusions that one can draw. As CFHTLenS has similar statistical power to current weak-lensing surveys, however, it nevertheless provides a very useful testbed with which to demonstrate the impact of adopting different approximations when constraining cosmological parameters.

In this work, we focus on the weak-lensing power spectrum projection, and assess the impact of various approximations on cosmological constraints from CFHTLenS. For consistency with the original analysis (Kilbinger et al. 2013), we adopt the same priors and non-linear power spectrum corrections from Smith et al. (2003).

We re-analyse the 2D CFHTLenS measurement of the two-point shear correlation function $\xi_{\pm}(\theta)$ from Kilbinger et al. (2013), defined in equation (47). As in Kilbinger et al. (2013), we fit both components $\xi_+$ and $\xi_-$ between angular scales $\theta = 0.8$ and 350 arc min, and use a N-body simulation estimate of the non-Gaussian covariance including the cross-covariance between both components. Bayesian Population Monte Carlo parameter sampling is performed using the publicly available software COSMOPMC<fnref n="5" /> (Wraith et al. 2009; Kilbinger et al. 2010). The cosmological modelling part includes the various lensing projections, calculated using the software library NICAEA.<fnref n="6" />

For a first-order standard Limber flat-sky approximation (L1Fl), we find $\sigma_8(\Omega_m/0.27)^{0.6} = 0.787^{+0.031}_{-0.033}$, the same result that was published in Kilbinger et al. (2013). Using the second-order extended Limber flat-sky hybrid approximation (ExtL2FlHyb) results in $\sigma_8(\Omega_m/0.27)^{0.6} = 0.788 \pm 0.032$, a negligible change of the amplitude that is well within the Monte Carlo sampling noise. The largest difference is measured with the deprecated case ExtL1Fl, for which the recovered amplitude is larger by 16 per cent of the statistical error. These negligible changes to the error bars were to be expected owing to the high level of statistical noise and cosmic variance in comparison to the low-level impact of the various approximations shown in Fig. 1.

Table 2 lists the mean and 68 per cent credible interval for $\sigma_8\Omega_m^{0.6}$ for the various approximations to the lensing power spectrum projections listed in Table 1. Note again that these values do not represent the state-of-the-art cosmological results, since many of the above listed analysis advancements made since 2013 have not been taken into account. As an example of a significant effect, when using the revised non-linear power spectrum of Takahashi et al. (2012) in place of Smith et al. (2003), there is a decrease of $0.6\sigma$ with $\sigma_8(\Omega_m/0.27)^{0.6} = 0.768^{+0.029}_{-0.031}$, using L1Fl.

Considering the cosmological constraints from tomographic Kilo-Degree Survey (KiDS), we conclude that these are robust to flat-sky and Limber approximations. The case ExtL1FlHyb that was used for the analysis of KiDS data in Hildebrandt et al. (2017) and Joudaki et al. (2016) introduces errors that are more than an order of magnitude lower than the cosmic variance for that survey, and thus this approximation has a negligible impact on the cosmological parameters.

<footnote n="5">⁵ http://www.cosmostat.org/software/cosmopmc</footnote>

<footnote n="6">⁶ http://www.cosmostat.org/software/nicaea</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Integrand of $\xi_+$ (upper), $\xi_-$ (upper middle), $E_1$ (lower middle, E-COSEBIs) and $\langle M_{ap} \rangle^2$ (lower panel). All integrands are of the form $\ell F(\ell)P(\ell)$, where $F(\ell)$ is the corresponding weight-function for each statistic and $P(\ell)$ is the E-mode convergence power spectrum, with the exception of $\xi_\pm$, for which $P(\ell)$ is equal to the sum of the E and B-mode power spectra. Two cases are shown for each statistic as listed in each caption. For the aperture mass statistic, $\theta_{\max} = 2\theta$ is shown. Note that higher order COSEBIs modes generally probe larger $\ell$-modes, hence here we show only the lowest mode $E_1$. All values are normalized with respect to their maximum value. This figure illustrates how different two-point cosmic shear statistics have different dependences between the angular scales sampled and the $\ell$-range probed.</figcaption></figure>

## 4.4 Alternative two-point shear statistics; the mass aperture statistic and COSEBIs

The two-point shear correlation function $\xi_\pm$ represents the current baseline observable for cosmic shear measurements. As shown in Fig. 3, however, using the standard first-order extended Limber flat-sky approximation (ExtL1FlHyb) can result in errors exceeding 10 per cent, on angular scales $\theta > 300$ arcmin. This is a result of the weight given to low $\ell$ modes in the $\xi_+$ statistic, as illustrated in Fig. 4, which shows the integrand of $\xi_+$ and $\xi_-$ (upper two panels) for two cases ($\theta = 100$ and $350$ arcmin), normalized to their maximum value. This error does not impact CFHTLenS analyses, given the low signal to noise of the measurements on these scales. It will however become increasingly important for upcoming wider field surveys that will accurately probe these scales.

In this paper, we provide a solution in the form of the second-order extended Limber approximation, but another option to consider is the use of alternative two-point shear statistics that are less sensitive to accuracy in shear power spectrum measurement at low $\ell$. Both the aperture-mass dispersion, $\langle M_{ap}^2 \rangle$ (Schneider et al. 1998), and the Complete Orthogonal Sets of E/B-mode Integrals (COSEBIs), $E_n$ (Schneider et al. 2010) statistics satisfy this requirement and are linearly related to the shear power spectrum in the flat-sky approximation via integrals of the form

$$ \langle M_{ap}^2 \rangle(\theta) = \frac{1}{2\pi} \int_0^\infty d\ell \, \ell \, \hat{U}^2(\theta\ell) P^\gamma(\ell), \qquad (49) $$

$$ E_n = \frac{1}{2\pi} \int_0^\infty d\ell \, \ell \, W_n(\ell) P^\gamma(\ell), \qquad (50) $$

where the Fourier-space filter functions $\hat{U}$ and $W_n$ are defined in Schneider et al. (1998) and Schneider et al. (2010), respectively. Fig. 4 shows the integrands of these statistics, again normalized to their maximum value, where the integrands are of the form $\ell F(\ell)P^\gamma(\ell)$. The lower middle panel in Fig. 4 shows the COSEBIs integrands for two angular ranges, $[1', 100']$ and $[0.8', 350']$, where we show only the integrands
for the lowest COSEBIs mode, $E_1$, as the higher modes generally probe larger $\ell$-modes. The lowest panel shows the integrands of aperture mass dispersion statistics, for the same two maximum angular ranges.

Note that the development of the aperture-mass dispersion statistic, $(M_{ap}^2)$, was initially motivated to enable the separation of the measured signal into an E-mode (cosmological signal) and B-mode (systematics). This statistic is, however, a lossy conversion and is biased by small angular separations, where blending of galaxies makes shear measurement challenging (Kilbinger, Schneider & Eifler 2006). The COSEBIs statistic tackles both these shortcomings. Kilbinger et al. (2013) present a detailed comparison of cosmological constraints obtained from this range of different two-point shear statistics finding consistent results.

## 5 CONCLUSIONS

In this paper, we evaluate precision theoretical calculations for cosmic shear observables, bringing together different sources from the literature to provide a pedagogical review of the impact of adopting flat-sky and Limber approximations. We demonstrate that for current surveys, such as CFHTLenS and KiDS, these approximations have a negligible impact on cosmological parameter constraints.

For future surveys, the decrease in statistical errors places higher requirements on the accuracy of the theoretical modelling. There is also, however, the need to be able to rapidly sample the multidimensional cosmological parameter likelihoods. This requirement for computational speed is incompatible with a theoretical analysis that calculates a full spherical solution for the shear power spectrum, without adopting any approximation. We therefore present alternative solutions, revisiting the work of Bernardeau et al. (2012), who showed that adopting the second-order extended Limber approximation of Loverde & Afshordi (2008) provides a representation of the full spherical solution for the shear power spectrum that is accurate at the sub-percent level for $\ell > 3$. We have verified this result and provide to the community our fast numerical implementation of all the approximations studied in this analysis, and the slow calculation of the full projection within the publicly available package NICAEA at http://www.cosmostat.org/software/nicaea.

Finally, we propose that future surveys seek to optimize the statistical analyses of their cosmic shear data. For example, moving from the standard two-point shear correlation function statistic to the more stringent 'COSEBI' statistic (Schneider et al. 2010) renders the cosmic shear measurement insensitive to the low-$\ell$ scales where the Limber and flat-sky approximations have an impact on the precision of the theoretical modelling.

We have considered a flat universe throughout this paper. To generalize the calculations to non-flat models, one needs to modify the comoving angular diameter distance to account for the spatial curvature $K \neq 0$. In addition, the spherical Bessel functions are replaced with hyperspherical Bessel functions (Abbott & Schaefer 1986). In a universe with positive curvature $K > 0$, the 3D wave modes become discrete integer variables. For non-flat models, we do however not expect qualitative differences from our results.

## ACKNOWLEDGEMENTS

The authors thank Tommaso Giannantonio, Sarah Bridle, Ami Choi, Donnacha Kirk, Lance Miller, Benjamin Joachimi, Chris Blake, Joe Zuntz, Joanne Cohn, Alex Hall and Adam Amara for very helpful discussions. This work was financially supported by the Deutsche Forschungsgemeinschaft (DFG) (Emmy Noether grant Hi 1495/2-1; SI 1769/1-1; TR33 'The Dark Universe'), the Alexander von Humboldt Foundation, the European Research Council (ERC) (grants 279396, 647112), the Seventh Framework Programme of the European Commission (Marie Sklodwoska Curie Fellowship grant 656869), the Netherlands Organization for Scientific Research (NWO) (grants 614.001.103), Natural Sciences and Engineering Research Council of Canada (NSERC), Canadian Institute for Advanced Research (CIFAR) and the World Premier International Research Center Initiative (WPI), Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. Parts of this research were conducted by the Australian Research Council Centre of Excellence for All-sky Astrophysics (CAASTRO), through project number CE110001020i, and by the German Federal Ministry for Economic Affairs and Energy (BMWi) under project no. 50QE1103.

## REFERENCES

* Abbott L. F., Schaefer R. K., 1986, ApJ, 308, 546

* Abbott T. et al., 2016, Phys. Rev. D, 94, 022001

* Asgari M., Heymans C., Blake C., Harnois-Deraps J., Schneider P., Van Waerbeke L., 2017, MNRAS, 464, 1676

* Bartelmann M., Schneider P., 2001, Phys. Rep., 340, 291

* Benjamin J. et al., 2013, MNRAS, 431, 1547

* Bernardeau F., 1998, A&A, 338, 375

* Bernardeau F., Bonvin C., Van de Rijt N., Vernizzi F., 2012, Phys. Rev. D, 86, 023001

* Blanco M. A., Flórez M., Bermejo M., 1997, J. Mol. Struct., 419, 19

* Castro P. G., Heavens A. F., Kitching T. D., 2005, Phys. Rev. D, 72, 023516

* Choi A. et al., 2016, MNRAS, 463, 3737

* Chon G., Challinor A., Prunet S., Hivon E., Szapudi I., 2004, MNRAS, 350, 914

* Erben T. et al., 2013, MNRAS, 433, 2545

* Fenech Conti I., Herbonnet R., Hoekstra H., Merten J., Miller L., Viola M., 2017, MNRAS, 467, 1627

* Fu L. et al., 2008, A&A, 479, 9

* Giannantonio T., Porciani C., Carron J., Amara A., Pillepich A., 2012, MNRAS, 422, 2854
Hamana T., Colombi S. T., Thion A., Devriendt J. E. G. T., Mellier Y., Bernardeau F., 2002, MNRAS, 330, 365

Hamilton A. J. S., 2000, MNRAS, 312, 257

Heymans C. et al., 2012, MNRAS, 427, 146

Hildebrandt H. et al., 2012, MNRAS, 421, 2355

Hildebrandt H. et al., 2017, MNRAS, 465, 1454

Hu W., 2000, Phys. Rev. D, 62, 043007

Jee M. J., Tyson J. A., Schneider M. D., Wittman D., Schmidt S., Hilbert S., 2013, ApJ, 765, 74

Joudaki S., Kaplinghat M., 2012, Phys. Rev. D, 86, 023526

Joudaki S. et al., 2016, MNRAS, 471, 1259

Joudaki S. et al., 2017, MNRAS, 465, 2033

Kaiser N., 1992, ApJ, 388, 272

Kaiser N., 1998, ApJ, 498, 26

Kilbinger M., Schneider P., Eifler T., 2006, A&A, 457, 15

Kilbinger M. et al., 2010, MNRAS, 405, 2381

Kilbinger M. et al., 2013, MNRAS, 430, 2200

Kitching T. D., Heavens A. F., 2016, Phys. Rev., 95, 063522

Kitching T. D., Alsing J., Heavens A. F., Jimenez R., McEwen J. D., Verde L., 2016, MNRAS, 469, 2737

Kuijken K. et al., 2015, MNRAS, 454, 3500

Lemos P., Challinor A., Efstathiou G., 2017, J. Cosmol. Astropart. Phys., 5, 014

Limber D. N., 1953, ApJ, 117, 134

Loverde M., Afshordi N., 2008, Phys. Rev. D, 78, 123506

Mead A. J., Peacock J. A., Heymans C., Joudaki S., Heavens A. F., 2015, MNRAS, 454, 1958

Miller L. et al., 2013, MNRAS, 429, 2858

Miralda-Escude J., 1991, ApJ, 370, 1

Ng K.-W., Liu G.-C., 1999, Int. J. Mod. Phys. D, 8, 61

Peebles P. J. E., 1980, The Large-Scale Structure of the Universe. Princeton Univ. Press

Planck Collaboration XIII, 2016, A&A, 594, A13

Schmidt F., 2008, Phys. Rev. D, 78, 043002

Schneider P., Van Waerbeke L., Jain B., Kruse G., 1998, MNRAS, 296, 873

Schneider P., Van Waerbeke L., Mellier Y., 2002, A&A, 389, 729

Schneider P., Eifler T., Krause E., 2010, A&A, 520, A116

Smith R. E. et al., 2003, MNRAS, 341, 1311

Takahashi R., Sato M., Nishimichi T., Taruya A., Oguri M., 2012, ApJ, 761, 152

Taylor A. N., 2001, Phys. Rev. Lett. submitted in press

Van De Rijt N., 2012, PhD thesis, École Polytechnique

Wraith D., Kilbinger M., Benabed K., Cappé O., Cardoso J.-F., Fort G., Prunet S., Robert C. P., 2009, Phys. Rev. D, 80, 023507

# APPENDIX A: DERIVATION OF THE WEAK-LENSING POWER SPECTRA

The following derivations are detailed in Hu (2000) and Castro et al. (2005), and are provided here for completeness.

## A1 Spherical case

### A1.1 Lensing potential power spectrum

To obtain the power spectrum of the lensing potential, we insert the lensing projection (equation 1) into the inverse harmonics expansion (equation 7) and write the 3D potential as its Fourier transform (equation 5) to get

$$ \psi_{\ell m} = \frac{2}{c^2} \int d\Omega Y_{\ell m}^*(\theta, \varphi) \int_0^\infty \frac{d\chi}{\chi} q(\chi) \int \frac{d^3 k}{(2\pi)^3} \hat{\Phi}(\mathbf{k}; \chi) e^{-i\mathbf{k}\cdot\mathbf{r}}. \quad (A1) $$

The 3D position vector $\mathbf{r}$ is a 3D position vector with polar coordinate $\mathbf{r} = \chi$ and polar angles $(\theta, \varphi)$. Similarly, we denote with $\theta_k, \varphi_k$ the polar angles of the 3D Fourier vector $\mathbf{k}$. We insert the expansion of a plane wave into spherical harmonics,

$$ e^{i\mathbf{k}\cdot\mathbf{r}} = 4\pi \sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{\ell} i^{\ell} j_{\ell}(k\chi) Y_{\ell m}(\theta, \varphi) Y_{\ell m}(\theta_k, \varphi_k). \quad (A2) $$

Making use of the orthogonality of the spherical harmonics

$$ \int d\Omega Y_{\ell m}(\theta, \varphi) Y_{\ell'm'}^*(\theta, \varphi) = \delta_{\ell\ell'} \delta_{mm'}, \quad (A3) $$

the expression for $\psi_{\ell m}$ simplifies to

$$ \psi_{\ell m} = \frac{i^\ell}{c^2 \pi^2} \int_0^\infty \frac{d\chi}{\chi} q(\chi) \int d^3 k \, \hat{\Phi}(\mathbf{k}; \chi) j_\ell(k\chi) Y_{\ell m}(\theta_k, \varphi_k). \quad (A4) $$

To obtain the potential power spectrum, we take the absolute square of the last equation and use the definition of the 3D potential power spectrum (equation 12). The delta-function resolves one 3D Fourier integral. We split the second integration into radial and spherical
coordinates, $d^3k = dk^2d\Omega_k$ and use once again the orthogonality of the spherical harmonics to resolve the spherical integral. This leads to the potential power spectrum in equation (9).

## APPENDIX B: DISCUSSION AND COMPARISON TO PREVIOUSLY PUBLISHED WORK

In this section, we briefly discuss previously published work on the full projection and second-order Limber equation of the lensing power spectrum, cross-checking and comparing their results with our independent findings.

### B1 Kitching et al. 2016 (version 1)

Kitching et al. (2016) compute the full projection of the weak-lensing power spectrum, which they present as spherical-radial representation of the 3D shear field. Our results in equation (21) correspond to their equations (7) and (8) assuming a flat universe and the case of perfect photometric redshifts, with $p(z|z_p) = \delta_D(z - z_p)$, and for a bin function that selects the redshift bin of z, $W^{SR}(z, z_p)$ is unity if $z_p$ is in the redshift bin denoted by z, and zero otherwise. We find that equation (7) in Kitching et al. (2016) is missing a factor $2/\pi$.

Kitching et al. (2016) derive the spherical and extended Limber approximation starting from the full spherical projection in their Appendix A. We find that the filter function $q$ defined in their equation (31) has an additional factor of comoving distance $r$, and an additional factor of $\pi/2$.

As shown in this paper, we are unable to reproduce the differences that Kitching et al. (2016) report, between the full spherical solution and the different approximations, neither for the power spectrum nor for the shear correlation function.

### B2 Bernardeau, Bonvin, Van de Rijt, Vernizzi (2012); Van de Rijt (2012)

Bernardeau et al. (2012) present the non-tomographic full projection $C(\ell)$ in the approximation of the 3D potential power spectrum $P_\Phi$ separating into $k$- and $\chi$-dependent functions in their equation (44). Their expression holds for a single source redshift.

The PhD thesis of Van De Rijt (2012) presents an explicit calculation of the second-order Limber approximation. They carry out the derivatives of the kernels $f$ under the assumption of a constant growth-suppression factor $D_+(a)/a$. Then, for a constant source comoving distance $\chi_S$, the lensing efficiency in equation (2) is $q(\chi) = (\chi_S - \chi)\chi_S^{-1}$, and the derivative of the separated kernel function in equation (41) can be calculated analytically,

$$f''_s(\chi) + \frac{\chi}{3} f'''_s(\chi) = -\frac{1}{8} \chi^{-5/2} \left( \frac{5}{\chi} + \frac{1}{\chi_S} \right) \frac{D_+(\chi)}{a(\chi)} \quad (B1)$$

This result confirms the power-law behaviour for $\chi \ll \chi_S$ of the filter function, which we exploited earlier to fit this function.

Inserting equation (B1) into the second-order Limber power spectrum in equation (42) and using the inverse Poisson equation to replace the matter with the potential power spectrum, we obtain the same expression as Van De Rijt (2012) (their equation 7.19).

In Fig. B1, we reproduce fig. 7.3 from Van De Rijt (2012) using a similar set up, a flat $\Lambda$ Universe with $\Omega_m = 0.3$, $h = 0.65$, $\Omega_b = 0.0461$, $\sigma_8 = 0.8$, $n_s = 0.96$. All source galaxies are at redshift $z_S = 1$. The non-linear 3D matter power spectrum from Takahashi et al. (2012) is used. The ratio of the first- and second-order Limber approximated power spectra to the full projection shows excellent agreement at the sub-percent level.

### B3 LoVerde & Afshordi (2009)

This paper introduces the extended Limber approximation to second order that we apply in this work. Although they present the specific case of 2D galaxy clustering, their calculations are general enough to apply to a weak-lensing context. Their equation (5) is a spherical cross-power spectrum of two scalar fields A and B, projected from 3D to 2D via projection kernels $F_A$ and $F_B$ defined in their equation (4). Comparing their expressions with the weak-lensing potential (1), we set $F_A(\chi) = F_B(\chi) = 2c^{-2}D_+(\chi)q(\chi)\chi^{-1}$. With that, their equation (5) is identical to equation (9).

The second-order Limber approximation in LoVerde & Afshordi (2009) is presented in equation (12). This is consistent with our first-order (equation 38) and second-order (equation 39) Limber approximation terms, when accounting for the difference between lensing potential and shear 2D power spectrum, and 3D potential and matter power spectrum.

### B4 Schmidt (2009)

Schmidt (2008) derive the lensing power spectrum in the flat-sky limit (see their equation 9). Inserting the Poisson equation and the growth function, and writing the redshift filter function $W_\kappa$ (eq. 10) in terms of the lensing efficiency and comoving distances, $W_\kappa[z(\chi)] = H(z)^{-1}\chi q(\chi)$, their expression, using our notations, reads

$$P_{ij}^{\gamma}(\ell) = \frac{2}{\pi} A^2 \int_{0}^{\infty} d\chi \chi \frac{q_i(\chi)}{a(\chi)} \int_{0}^{\infty} d\chi' \chi' \frac{q_j(\chi')}{a(\chi')} \int_{0}^{\infty} dk k^2 P_m(k; \chi, \chi') j_{\ell}(k\chi) j_{\ell}(k\chi'), \quad (B2)$$
<figure><img src="image_7.png" /><figcaption>**Figure B1.** The relative differences in percentage of the spherical first- and second-order Limber shear power spectra with respect to the full projection as function of wave mode $\ell$ (see Table 1). In this figure, the redshift distribution is chosen to be at a single source plane at $z_S = 1$, with the cosmological parameters (see text) are chosen to match Van De Rijt (2012, see their fig. 7.3 for comparison).</figcaption></figure>

This is consistent with our equation (29) under the additional assumption that mainly modes with $k\chi \approx k\chi' \approx \ell$ contribute to the integral, that is modes around the maxima of the Bessel functions. Then, we can draw out of the integral the factor $\ell^4 \approx k^4\chi^2\chi'^2$, to recover equation (29). In fact, using the approximation $k\chi = \ell$ seems to go too far, this is already halfway the Limber approximation.

## B5 Giannantonio et al. (2012)

Giannantonio et al. (2012) derive the flat-sky lensing power spectrum in their equations (25) and (26). Their window function $W^{\epsilon_i}$ defined in their equation (25) for a flat universe equals $Aq(\chi)/a(\chi)$, since $dz(dN(z)/dz) = d\chi n(\chi)$; however, due to a typo there is a factor $r_K[r(z)]$ missing in the window function, which translates into a missing $r_K[r]r_K[r']$ in the full projection integral.<fnref n="7" /> This also leads to an erroneous $r_K^{-2}(r)$ in their Limber equation (27). With these factors accounted for, and making the additional approximation $k\chi \approx k\chi' \approx \ell$ (see their appendix B4), we reproduce the expressions of Giannantonio et al. (2012).

## APPENDIX C: FAST EVALUATION OF THE SHEAR CORRELATION FUNCTION ON THE SPHERE

The calculation of the shear correlation function on the sphere requires the estimation of the reduced Wigner D-matrices $d_{2,2}^\ell$ and $d_{2,2-\ell}$. The general calculation of $d_{m,n}^\ell$ is cumbersome, but there are quick and numerically stable recurrence relations if we are interested in only a subset of these matrices. In particular, following Blanco, Flórez & Bermejo (1997), we can show that

$$d_{m,n}^{\ell} = \frac{\ell(2\ell - 1)}{\sqrt{[\ell^2 - m^2][\ell^2 - n^2]}} \left[ \left( d_{0,0}^1 - \frac{m}{\ell(\ell - 1)} \right) d_{m,n}^{\ell-1} - \frac{\sqrt{[(\ell - 1)^2 - m^2][(\ell - 1)^2 - n^2]}}{(\ell - 1)(2\ell - 1)} d_{m,n}^{\ell-2} \right], \quad (C1)$$

which allows us to calculate all required reduced Wigner D-matrices from the first two elements.

<footnote n="7">⁷ Typo confirmed by T. G., priv. comm.</footnote>