The unequal-time matter power spectrum: impact on weak lensing observables

We investigate the impact of a common approximation on weak lensing power spectra: the use of single-epoch matter power spectra in integrals over redshift. We disentangle this from the closely connected Limber's approximation. We derive the unequal-time matter power spectrum at one-loop in standard perturbation theory and effective field theory to deal with non-linear physics. We compare these formalisms and conclude that the unequal-time power spectrum using effective field theory breaks for larger scales. As an alternative, we introduce the midpoint approximation. We also provide, for the first time, a fitting function for the time evolution of the effective field theory counterterms based on the Quijote simulations. Then we compute the angular power spectrum using a range of approaches: the Limber's approximation, and the geometric and midpoint approximations. We compare our results with the exact calculation at all angular scales using the unequal-time power spectrum. We use DES Y1 and LSST-like redshift distributions for our analysis. We find that the use of the Limber's approximation in weak lensing diverges from the exact calculation of the angular power spectrum on large-angle separations, $\ell<10$. Even though this deviation is of order $2\%$ maximum for cosmic lensing, we find the biggest effect for galaxy clustering and galaxy-galaxy lensing. We show that not only is this true for upcoming galaxy surveys, but also for current data such as DES Y1. Finally, we make our pipeline and analysis publicly available as a Python package called unequalpy.


Introduction
Cosmology is living a golden era of high-precision observations with upcoming galaxy surveys probing the Universe on unprecedented smaller scales, such as DESI 1 , Euclid 2 , the Rubin Observatory 3 and the Nancy Grace Roman Space Telescope 4 . This upcoming galaxy data will have implications in the prediction of cosmic lensing observables.
When we look at the sky, we observe a two-dimensional projection of different source points at different cosmological distances. In doing so, we look back in time, capturing all the light from our past light cone, see Figure 1, observing objects at different time (or redshift) slices 5 . Objects that seem close to each other are actually farther apart and might belong to different time slices. In addition, if the distance between two objects is not large, their look-back time can be thought to be equal. This is the widely-used thin-shell approximation, whose validity should be tested in this era of high-precision cosmology. To do so, we explore the angular correlation functions and the angular power spectrum for different quantities in weak lensing. For two general fields, A and B, the angular correlation function reads ω (i,j) AB (θ) = dr 1 dr 2 f i A (r 1 )f j B (r 2 )ξ AB (r 12 ; r 1 , r 2 ) (1.1) Figure 1: Observations within our past light-cone. This diagram represents two different time slices within our past light-cone, t 1 and t 2 where t 1 > t 2 . It shows two source points, A and B, at different cosmological distances, r 1 and r 2 (or redshift, z 1 and z 2 ). Note that structure has evolved between these two time slices.
with f k C (r) the weight functions for redshift bins k. ξ AB (r 12 ; r 1 , r 2 ) is the spatial unequaltime correlation function for spatial separations r 12 , c.f. Fig. 1. This quantity measures the correlation between the fields A and B at two different time slices and its Fourier transform is the unequal-time power spectrum, P AB (k; r 1 , r 2 ). Likewise, after expanding in spherical harmonics, the angular power spectrum yields C (i,j) AB ( ) = 2 π dkk 2 dr 1 dr 2 f i A (r 1 )f j B (r 2 )j (kr 1 )j (kr 2 )P AB (k; r 1 , r 2 ) .
(1. 2) with j (kr) the spherical Bessel function of order . For a detailed derivation, the reader can refer to [1], for example. Solving these integrals can be non-trivial, mostly when working in harmonic space where the spherical Bessel functions present a highly-oscillatory behaviour. Through history, many numerical and analytical approaches have been tested. For a summary on the numerical approaches, one could refer to our accompanying paper [2]. On the analytical side, one could focus on approximations to the unequal-time correlators, or make assumptions about the filters or approximate the spherical Bessel functions to the amplitude of their first peak.
At the level of the power spectrum, the most widely-used approximation to the unequaltime correlators is the geometric approximation. For this one computes the geometric mean of the two equal-time power spectra at comoving distances or redshifts z 1 and z 2 , respectively P (k; z 1 , z 2 ) ≈ P (k; z 1 )P (k; z 2 ) .
(1. 3) This means that instead of obtaining the full correlation between two time slices, Fig. 1, one simply computes the power spectrum at each time slice . The validity of this approximation is restricted to very large scales where linear physics is exact. That being said, there exist some work on the computation of unequal-time correlators. Kitching et al. [1] were among the first ones in deriving such description using one-loop standard perturbation theory. They address the accuracy issue for the equal-time power spectrum approximation, claiming that even the use of the geometric mean power spectrum (1.3) may result in biased predictions of the cosmic shear power spectrum for Euclid-or LSST-like weak lensing experiments. They stress that in order to compute unequal-time correlators to sufficient accuracy, advance in perturbation theory on non-linear scales is required. This is part of the main focus of this work. Conversely, Chisari and Pontzen in [3] show that the Zel'dovich approximation in Lagrangian perturbation theory provides a much more accurate (< 10 %) analytical description to model unequal-time correlators and validate their results against N-body simulations.
Developing an accurate description of the unequal-time power spectra is one of the main goals of our work. For this we re-derive the unequal-time power spectrum in standard perturbation theory [1]. We deal with non-linear physics using effective field theory and obtain a new approximation to the unequal-time prescription: the midpoint approximation.
At the level of the angular power spectra, Limber made the first attempt to compute the angular correlation functions. In [4], they develop a method to analyse the counts of the extra-galactic "nebulae", i.e. galaxies, in terms of a fluctuating density field for the "nebulae" in space. For this, they assume that the filters are smooth and the correlators fall off fast enough so that they can compute all these quantities at the mean redshift. Then they apply their methods to the counts obtained by Shane and Wirtanen at the Lick Observatory [5]. Twenty years later, Peebles [6] presents a number of theoretical results for the analysis of data distributed on a sphere, applicable when the survey covers only a portion of the sky, employing a discrete version of the Limber's equation. In the late 90s the approximation starts taking its familiar shape when Kaiser [7] re-derives the Limber's equation within the flat-sky approximation in Fourier space for a homogeneous and isotropic universe with spatial curvature. A similar derivation for the spherical case can be found in Lemos et al. [8].
More recently, Simon [9] revisits the Limber's equation. They distinguish between two different regimes: small-angle and large-angle separations, showing that Limber's and the so-called thin-layer equations are approximations for these two extremes. They also study Limber's accuracy for a power-law spatial correlation and claims that Limber's approximation diverges when galaxies are closely distributed. This implies that Limber's equation possibly over-estimates the angular correlation to some degree. For historical reasons, when employing Limber's equation, the small-angle approximation is automatically used. Simon explains that this is accurate to about 10% for angles smaller than θ 40 • . This type of analysis has led to the idea that such small-angle approximations could contribute significantly to the tension between the CMB measurements and weak lensing data -c.f. [10]. However, Lemos et al. [8] conclude that the impact of small-angle approximations on cosmological parameter estimation is negligible for current data.
In addition, some efforts were devoted to extend the Limber's approximation in its Dirac-delta version. These are the so-called post-Limber approximations. One example is Lo Verde et al. [11] where they develop a systematic derivation for the Limber approximation to the angular cross-power spectrum of two random fields as a series expansion in 1/( + 1/2). Nontheless, it is not clear how this would alleviate the divergence of the series expansion at small .
Finally, Fang et al. [12] present a new method based on a generalised FFTLog algorithm for the efficient computation of angular power spectra beyond the Limber approximation. This simplifies the calculation and improves the numerical speed and stability. They implement their method for galaxy clustering and galaxy-galaxy lensing power spectra, and find that the Limber's approximation for galaxy clustering in future analyses like LSST Year 1 and DES Year 6 may cause significant biases in cosmological parameters, indicating that going beyond the Limber approximation is necessary for these analyses. One of their key points in their method is the separation of the integrals of the angular power spectrum into large scales and small scales. By doing so, they drop the Limber's approximation in the linear contribution and they never have to compute the unequal-time power spectrum for the non-linear term. However, for the latter, they employ the Limber's approximation that transforms the unequal-time correlation into the usual equal-time power spectrum. Even though this splitting seems natural and logical, we wonder whether they might be losing accuracy on the small-scale contribution and whether the non-linearities are well accounted for. For this reason, we understand that efforts should be focused on implementing an all-angle method to compute the exact unequal-time calculation.
In the following, we develop our analysis on the unequal-time matter power spectrum and its impact on weak lensing. In Section 2, we derive the unequal-time matter power spectrum in standard perturbation theory and effective field theory, and we present the new midpoint approximation. For those more interested in weak lensing analysis, we recommend to skip Section 2 and read Section 3. In this section 3, we compute the main observables in weak lensing using a range of approximations and the all-angle approach presented in our accompanying paper [2]. We then analyse the relevance of the unequal-time description using DES Y1, and LSST-like redshift distributions. A brief summary and the main conclusions can be found in Section 4, followed by the appendices A, B and C. Our final product is the publicly available python package unequalpy 6 [13] with functionality to reproduce all the results and analysis presented in this paper.

Unequal-time power spectrum
The main goal of this section is to develop an accurate description of the unequal-time power spectra dealing with non-linear physics. To do so, we re-derive the unequal-time matter power spectrum up to third order in standard perturbation theory [1]. The equal-time power spectrum is defined as whereas the unequal-time power spectrum reads Note that the definition above (2.2) ensures an isotropic power spectrum, with no privileged direction in the local coordinate system. Therefore, the power spectrum only depends on the length of the wavenumber k and not on its direction. In addition, we deal with non-linear physics using effective field theory. Some authors [14,15] showed the capability of this framework to encode small-scale physics and to provide highly accurate predictions on increasingly smaller scales in the context of equal-time power spectra. To show whether such improvement propagates for unequal-time correlators, we also derive the one-loop unequal-time matter power spectrum using effective field theory of large scale structures. At the end of this section, we introduce a new approximation to the unequal-time prescription, the midpoint approximation, as an alternative to model non linearities.

Standard Perturbation Theory
To compute the matter power spectrum up to third order in perturbation theory, we need to solve perturbatively the equation of motion for the matter density contrast, δ = (ρ − ρ 0 )/ρ 0 (with ρ the matter density field and ρ 0 the background density) up to third order. Then we compute the two-point correlation functions using the Einstein-de Sitter approximation, EdS [16], and assuming the initial density perturbation δ * k to be a Gaussian random field. For the curious reader, we also derive the full-time dependence in appendix A. However, for the purpose of our analysis, we stick to the EdS approximation. Then there are three contributions labelled P 11 , P 22 and P 13 , where k is the common magnitude of the wavevectors k 1 and k 2 . The linear contribution P 11 is described as the tree-level power spectrum, and the sum P 22 + 2P 13 is the one-loop contribution (c.f. [15,16] for a detailed expression of these terms). Then the one-loop unequal-time matter power spectrum reads where we factored out the time dependence. D(z) is the normalised linear growth function, which we compute using SkyPy 0.3 functions [17]. Finally, the one-loop equal-time power spectrum in standard perturbation theory is retrieved by setting z 1 = z 2 = z in the above definitions (2.3). This yields P SPT (k, z) = D 2 (z)P 11 (k) + D 4 (z)P 22 (k) + 2D 4 (z)P 13 (k) . (2.5)

Effective Field Theory
In order to compute the unequal-time matter power spectrum using effective field theory, we follow the same renormalisation programme described in [15]. We split the 13-and 22type loop integrals into the linear and the non-linear regimes. The linear contribution is calculated by using the standard perturbation results, where we know this is exact. For the non-linear contribution, we Taylor expand the integrands in terms of k/k N L , dropping order four contributions. Thus all the microscopic physics is encoded in a redshift-dependent free-fitting parameter, the so-called counterterm, c(z) ≡ c s (z)/k N L . Then, the unequal-time power spectrum in effective field theory reads Finally, the equal-time power spectrum is retrieved when z 1 = z 2 = z in equation (2.6) In order to fit the counterterms, we use data from the Quijote simulations [18]. These simulations are very well documented, offer a great number of cosmologies, provide the matter power spectrum, and, most importantly, offer enough data for different mean redshifts.

Time evolution of the counterterms
The Quijote simulations [18] are a set of 43100 full N-body simulations spanning more than 7000 cosmological models, providing full snapshots of the simulations at redshifts 0, 0.5, 1, 2 and 3. In the following, we explain how we use this set of simulations to fit our counterterms and how we obtain their time evolution.
We use the Quijote fiducial model to perform a Bayesian analysis to find the best value of the counterterms (2.7) at every redshift available. We use a flat prior for the free fitting parameter and emcee3.0.2 7 . We use matter power spectrum data up to k = 0.4h/Mpc, since two-loop contributions start to dominate on smaller scales. The results are shown in Figure 2. The computed values are well parametrised by a fitting function of the form: with the best fit values m = 2.564 +0.008 −0.007 Mpc 2 /h 2 , n = 0.036 +0.001 −0.001 Mpc 2 /h 2 and a = 1.961 +0.008 −0.008 . This is already a new result in cosmology. For example, in [15] authors used a different approach to fit for the counterterms in real space, using an estimator for the fit at every single redshift and using the halofit matter power spectrum given by CAMB 8 .
In the next section, we analyse which of the above prescriptions provides the best framework to analyse unequal-time correlators, and weak lensing observables.

Distinguishing between non-linear and unequal-time effects
In this section, we assess the performance of the unequal-time matter power spectrum and we explore the effects coming from non-linear physics and unequal-time correlators, c.f. Table  1. To do so we apply the geometric approximation (1.3) to the unequal-time matter power spectra. In general, where ∆P (k; z 1 , z 2 ) represents the error when using the geometric approximation. This error vanishes when z 1 = z 2 and for linear theory. Equation (2.9) seems the best description to analyse effects coming from correlations between different time slices. Then we can disentangle non-linear effects by studying the effective field theory counterparts. Table 1: Grid of theoretical descriptions to analyse. SPT stands for "standard perturbation theory", and EFT for "effective field theory". Comparing items top to bottom (same column) would describe unequal-time effects, whereas the left-to-right comparison (same row) would describe the effect of the counterterms and non-linear physics. Comparing in diagonal would be a mix of both effects.
Non-linear effects: Left-to-right comparison in Table 1. We compare standard perturbation theory and effective field theory at the level of the geometric approximation and the full unequal-time correlator. Here we compute the error committed when using standard perturbation theory, instead of effective field theory, both for the geometric approximation and the unequal-time calculation. For both cases the power spectrum of reference is the one from the effective field theory: We calculate this error for a given mean redshift z m = 0.5 and different redshift separations: ∆z = 0 (same time slice) and ∆z = 0.2. In Figure 3, we observe that the error increases on the linear regime and there is little difference between using the geometric and unequal-time calculations. We observe a breaking scale where the prediction diverges which is shifted to lower values of k for higher redshift separations. This breaking scale indicates that higher-order loop corrections are needed (extra counterterms and stochastic parameters), therefore the prediction can no longer be trusted. When we look at the same time slice, ∆z = 0, the absolute error does not vanish, it becomes proportional to the counterterms multiplied by the power spectrum, ∝ 2c 2 (z m )P SPT (k, z m ). For increasingly larger redshift separations, the relative error tends to shift to the left and become larger, preserving the shape.   Time effects Figure 3: Non-linear effects (left) and unequal-time effects (right) when using the geometric approximation instead of the unequal-time description for both standard perturbation theory and effective field theory. On the left, we plot the fractional difference between standard perturbation theory and effective field theory using either the geometric approximation (red lines), equation (2.10), or the full unequal-time power spectra (black lines), equation (2.11). On the right, we plot the error when using the geometric approximation instead of the unequal-time description (equation 2.12) for both standard perturbation theory (black lines) and effective field theory (red lines). This is done for a mean redshift z m = 0.5 and different redshift separations, ∆z. The divergences in these plots show where the theory breaks down and higher-order corrections become dominant, therefore the prediction cannot be trusted on those scales. We also observe how this phenomena occurs at lower k values for effective field theory.
Unequal-time effects: Top-to-bottom comparison in Table 1. We compute the error committed when using the geometric approximation instead of the unequal-time calculation, for both standard perturbation theory and effective field theory. The power spectrum of reference is the unequal-time counterpart with theory = {SP T, EF T }. Again we analyse the error for z m = 0.5 and widths ∆z = {0, 0.1, 0.2}. In Figure 3, we observe that the effect of using the geometric approximation instead of using the unequal-time calculation yields an error that increases monotonically. This shows that the geometric approximation is very good on very large scales, where linear theory is exact, but not accurate enough to deal with non-linear physics. In principle, the effective field theory framework should improve such predictions, but it breaks before the SPT counterpart. Again, this effect is due to higher-order effects becoming dominant. We observe a similar behaviour on large scales for both formalisms, with a larger error for larger redshift separations. For measurements at the same time slice, ∆z = 0, the error vanishes. This is true because the unequal-time correlator equals the geometric approximation when In conclusion, Figure 3 shows the relative comparison between different perturbation formalisms and unequal-time descriptions, c.f. Table 1. However, answering the question "which description is best?" is no trivial task. The effective field theory prediction breaks on intermediate scales and we believe that the current prescription of unequal-time prediction within effective field theory is more sensitive to homogeneity breaking along the line of sight. Therefore we stick with the standard perturbation formalism for the rest of our analysis.
For the purpose of our weak lensing analysis, we work within the standard framework and develop a new approximation to the unequal-time prescription which alleviates some divergences on the mild non-linear regime. This is the so-called midpoint approximation.

The midpoint approximation
As a result of our incapability to determine the best prescription of deal with non-linear physics along the radial direction, we present an alternative and derive the new midpoint approximation, defined as with z m = (z 1 +z 2 )/2. This means that instead of using the continuous information along the line of sight, we only use the power spectrum at the mean redshift, regardless the width ∆z.
There are other choices to define the mean redshift, e.g. the more natural mean on comoving distance, some weighted mean redshift, or choosing the redshift such that you can retrieve the exact growth on large scales. Nonetheless, we make this particular choice for simplicity and the exact definition of the midpoint would have an insignificant effect compared to the use of Limber's approximation or the exact projection method. Equation (2.13) is different from the geometric mean approximation (1.3) which combines the information from the two ends of the redshift shell, z 1 and z 2 . In principle, equation (2.13) accepts any perturbation theory. We know that in the context of equal-time correlators, effective field theory predictions are several percent levels more accurate than the standard formalism. However, when using effective field theory, the error is larger. Then, until we understand which formalism truly reflects the reality of our universe, we will work within the standard framework.
We now show how this new approximation improves our predictions on smaller scales. We compute the error when using one of these approximations instead of the unequal-time power spectrum in standard perturbation theory. For the geometric approximation, this is given by equation (2.12). For the new midpoint approximation, the error reads (2.14) We show the results in Figure 4. In the left panel, we fix the midpoint value of redshift z m = 0.5 and consider different redshift separations. Again, we observe how the geometric approximation is exact on very large scales but the error grows on small scales, and again the approximation is worse for larger separations of redshift. The midpoint approximation shows a higher error on larger scales, but the prediction improves on mild non-linear scales. On even smaller scales, the midpoint approximation tends to meet the geometric approximation curve.
The new approximation shows a particular feature on a single scale around k ≈ 0.1h/Mpc. This is where the approximation equals the unequal-time value and is characteristic of the midpoint approximation. The geometric approximation only equals the unequal-time value when we look at the same time slice. Beyond such scale, the midpoint approximation underpredicts the actual power spectrum whereas the geometric approximation overpredicts it.
When we fix the redshift separation ∆z = 0.1 and vary the mean value z m (right panel in Figure 4), we observe that this turnover happens for increasingly smaller scales as we increase the mean redshift. This means that the unequal-time calculation and the midpoint approximation are equal at increasingly smaller scales for redshift bins that are farther away from the observer.   .14), and for the geometric approximation (red lines), equation (2.12), for a fixed mean redshift, z m = 0.5. It shows the effect of increasing redshift separation between the two epochs and the particular feature for the midpoint approximation at the scale where it equals the prediction from the unequal-time power spectrum. On the right, we show the same error but for a fixed redshift width, ∆z = 0.1. It shows the impact of increasing the mean redshift: lower error for larger mean redshift and the midpoint feature shifting to smaller scales.
In conclusion, the geometric approximation is better on very large scales. We introduced the midpoint approximation hoping that it would predict mild non-linear physics more accurately. However, none of these approximations are completely satisfactory at the level of the matter power spectrum. In the next section, we analyse whether these features propagate when computing weak lensing observables. For the curious reader, we recommend to read Appendix C where we show how these differences on large scales between the midpoint and the geometric approximations are small when computing the angular correlations.

Impact on weak lensing observables
In this section we compute the angular power spectrum in cosmic lensing. First, we calculate the angular correlation function w(θ) of the cosmic convergence field, galaxy clustering and galaxy-galaxy lensing, c.f. [19]. As we mentioned above in Sec. 1, they present a general form w where (i, j) corresponds to the redshift bins, a and b refer to the fields, x is comoving distance, f k c (x) are the filters and ξ(r 12 ; t(x 1 ), t(x 2 )) is the unequal-time matter correlation function for the separation at cosmic times (or redshift) corresponding to x 1 and x 2 -c.f. Fig. 1. The filters f k c (x) for the k-th redshift bin depend on the selection function of the corresponding galaxy survey and the field, c. The main fields we work with in weak lensing are • Cosmic convergence: given the three-dimensional matter density contrast field δ( r; t) in comoving coordinates r and cosmic time t, the convergence κ( Θ) in direction Θ is the integral with the filter defined as where H 0 and Ω m are the cosmological parameters, c is the speed of light, x is comoving distance, a(x) is the scale factor corresponding to x, and q i (x) is the lensing efficiency given by for the distribution n i (x) of observed sources within the i-th redshift bin.
• Galaxy number density: given the three-dimensional matter density contrast field δ( r; t), the galaxy density field δ g ( Θ; t) is the integral with the filter defined as where b(x) is the bias parameter and H(x) the Hubble parameter at a redshift corresponding to a comoving distance x.
In the following we introduce the main observables to be analysed. For simplicity, in this work we ignore intrinsic alignments, redshift-space distortions and lensing magnification.
• Convergence and cosmic shear. The convergence κ is the cosmic lensing quantity most directly related to the matter field, of which it is the projection along the line of sight. The angular correlation function of cosmic convergence reads using the filters (3.4).
In practice, it is not cosmic convergence but cosmic shear that is observable. The two-point statistics are related through their respective angular power spectra, C κκ l and C γγ l , with The results we obtain for cosmic convergence are therefore readily applied to cosmic shear.
• Galaxy clustering. The angular correlation function quantifies correlations between galaxy number density fields: using the filters (3.7).
• Galaxy-galaxy lensing. The angular correlation function quantifies the correlation between the shape of background (or source) and foreground (lens) galaxy number density. In the weak lensing regime, the observed galaxy shape is the sum of an intrinsic (unlensed component) and a shear due to gravitational lensing. For simplicity, we only consider the shear component: using the filters (3.4) and (3.7).

Computation of the angular power spectrum
There are different approaches to compute the angular power spectra in cosmic lensing (1.2): 1. Brute-force calculation: This involves integrals of highly-oscillatory spherical Bessel functions, making the computation time-consuming with potential sources of inaccuracies. This is not due to numerical issues but due to a mathematical behaviour that we cannot control using sophisticated numerical methods.
2. Approximations to the unequal-time matter power spectrum, which do not reduce the number of integrals.
3. Assumptions about the filters: the Limber's approximation.
The curious reader can refer to Appendix B to understand why the Limber's and the geometric approximation seem synonyms and why their combination is completely unnecessary. In addition, note that the Dirac-delta approximation reduces the number of integrals to one, whereas the geometric approximation involves a triple integral. The Limber's equation reads • In angular space (original derivation [4]): one can write the angular correlation function (3.1) as the integral over the mean radial distance x = (x 1 + x 2 )/2 and the radial separation where the distance between the points in terms of x and r 12 is now Limber [4] introduced the approximation for the integral (3.17) using the assumption i) that the filters and correlation function change slowly and can be approximated by their midpoint values, (3.19) ii) that the angle separation θ between the points is small, so that the distance r 12 can be approximated as r 12 ≈ x 2 θ 2 + R 2 ; (3.20) iii) that the integral over R can be extended over the entire real line, assuming the spatial correlation function falls off fast enough. Limber's approximation for the correlation function (3.17) is thus where ξ L is Limber's matter correlation function, defined as 4. In this work, we employ a novel approach by computing the angular correlation function in real space, which eliminates the issue of integrating over highly oscillatory functions. We use the inverse Fourier transform of the unequal-time matter power spectrum (2.4) to obtain the unequal-time correlation function in configuration space, sin kr kr k 2 dk , (3.24) which can be evaluated efficiently over a logarithmic range of r values using the FFTLog algorithm [22]. To obtain the angular power spectrum C l from the angular correlation function (3.23), we use the general relation between the latter and the former, By evaluating the angular correlation function over a grid of θ values, and truncating this series at a suitable l max , the modes C l can be recovered by a least squares fit.
We have implemented all of the above steps in the corfu 9 package for Python [23].

DES and LSST surveys
In this section we show the results for cosmic lensing, galaxy clustering and galaxy-galaxy lensing using DES Y1 data [24], and LSST-like Y10 data [25]. We compare the results for a range of approaches (the Limber's (3.21), the geometric (3.13) and the midpoint approximations (3.14)) against the exact calculation (3.12). For this we employ the numerical methods derived in our accompanying paper [2].

LSST redshift distributions
We reproduce the redshift distribution of galaxies, following the LSST Dark Energy Science Collaboration document [25] using 10 redshift bins for the Y10 forecast.
• Lens sample. Ten photometric redshift bins z ph ∈ (0.2, 1.2) with ∆z = 0.1 convolving each bin with a Gaussian photo-z scatter The parametric redshift distribution reads with parameters (z 0 , α) = (0.28, 0.90). Then the true distribution of galaxies n i (z) that fall in the i-th photo-z bin is defined with a probability distribution p(z ph |z) in z ph at a given redshift, z. Instead of solving this integral, we use the error function given in equation (6) by Ma et al. [26].
Finally, the linear bias parameters are defined: is the normalised linear growth function, which we compute using SkyPy 0.3 functions [17].
• Source sample. Ten photometric bins z ph defined with equal numbers of source galaxies per bin. This is done using the true redshift distribution, and then the bins are convolved with the photo-z error distribution to make the photo-z distributions: The parameters in (3.27) are now (z 0 , α) = (0.11, 0.68).

Results
We show the results for cosmic lensing, galaxy-galaxy lensing and galaxy clusters in figures 6, 7, 8, 9 and 10. We perform the analysis for DES Y1, and LSST-like redshift distributions. However, to avoid cluttering we only choose to show the results for DES Y1. The LSST-like data generates similar results to DES Y1, even with narrower redshift bins. Simon [9] and Fang et al. [12] already warned us that the Limber's approximation could lead to biased predictions for upcoming galaxy surveys like DES Y6, LSST or Euclid. Here we show that this is true even for DES Y1 data. Likewise, findings by Chisari and Pontzen [3] and Kitching and Heavens [1] suggested that non-Limber predictions are sufficient using the equal-time correlators because they are insensitive to the small-scale physics. We agree with this statement, our results draw similar conclusions regardless the unequal-time prescription (one-loop standard perturbation theory, the midpoint or the geometric approximation). For larger-angle separations, c.f. Fig. 6 and 7, a full unequal-time description seems to provide a more accurate prediction. Those figures show how the Limber's curve deviates on larger angles and this effect propagates and causes the divergence of the angular power spectrum on small . The geometric and midpoint curves also deviates from the unequal-time prediction on large-angle separations. Figures 8, 9 and 10 show the comparison between the angular power spectrum when using one of these approximations (Limber, geometric or midpoint) with respect to the exact calculation at all angular scales using the unequal-time matter power spectrum. The most striking conclusion is the fact that the results from the geometric and midpoint approximations seem to mimic perfectly all the features from the exact calculation. This means that all the features that we observed at the level of the matter power spectrum, Figure 4, get washed out by the integrals over the line of sight (see Appendix C for details). Nonetheless, there is a tiny difference for very large-angle separations for galaxy-galaxy lensing and for the farthest redshift bins, for example (4,4). We also observe is that the Limber's approximation deviates from the exact calculation on large angle separations, < 10. This is something that we already knew, this effect was predicted by Simon [9]. This deviation is of order ∼ 2% maximum for cosmic lensing. For galaxy-galaxy lensing the minimum deviation is ∼ 1% and it can be very large for some redshift bins. The biggest effect is found for galaxy clustering.

Conclusions
In this paper we have advanced in the field of high precision cosmology, working on the theory side to match the demanding accuracy of upcoming galaxy surveys. These upcoming high-precision observations will have a considerable impact on weak lensing measurements. Therefore, we have devoted our efforts to compute the angular power spectrum exactly at all angular scales using higher-order unequal-time matter power spectrum in perturbation theory. Remember that we define equal-time correlators as the correlation between fields at the same time slice. Likewise, unequal-time correlators measure the correlation between fields at different times slices or redshift. Until now the most successful approach to compute the angular power spectrum was the Limber's approximation. Its success resides on the reduction of a triple integration to a single integral over the k range. In doing so, it also reduces the complexity at the level of the matter power spectrum as only the equal-time correlator is needed. In Fourier space, the exact calculation involves products of spherical Bessel functions that are highly oscillatory. This makes the computation numerically expensive and highly difficult. Nonetheless, many authors have already shown that such approximation will lead to biased predictions of the cosmological parameters with the upcoming galaxy surveys. Our analysis not only did support this statement, but also concluded that this is true for current data, such as DES Y1.   We computed the one-loop unequal time matter power spectrum using standard perturbation theory and effective field theory to deal with non-linear physics. We have compared with the traditional geometric approximation, concluding that effective field theory breaks at larger scales and tends to give a higher error than the prediction from standard pertur-bation theory. For these reasons, we continued the rest of our analysis within the standard formalism. We also presented a new unequal-time prescription, the midpoint approximation. We arbitrarily chose our definition of the mean redshift as the average between two different redshift slices. We explained that one can make many other definitions of the mean redshift, with the most natural choice being the mean on comoving distance. However, the same conclusions of this analysis apply regardless of such definition. We also conclude that, at the level of the matter power spectrum, the geometric approximation is much better on very large scales, whereas the midpoint approximation seems to equal the unequal-time prediction at some scale that depends on the mean redshift.
In our final section, we showed our results for cosmic lensing, galaxy clustering and galaxy-galaxy lensing. We have used data from DES Y1 and LSST-like Y10 for a whole range of approaches (the Limber's approximation, and the geometric and the midpoint approximations at the level of the power spectrum), comparing against the exact calculation. To compute these quantities we have used the numerical methods derived in our accompanying paper [2] where instead of computing the results in Fourier space, we compute the angular correlations in real space.
In the following, we present a list summarising the main results and conclusions of our work. At the level of the matter power spectrum: • We use the Quijote simulations to find the best value of the counterterms at every redshift available and we parametrise their time evolution. This was never done in the literature before.
• We assess the performance of the unequal-time matter power spectrum and we explore the effects coming from non-linear physics and unequal-time correlators. We conclude that the effective field theory framework does not provide an improved description of the non-linear unequal-time power spectrum and we introduce the new midpoint approximation.
• The geometric approximation is better on very large scales and the midpoint approximation show some features that might predict more accurately mild non-linear physics. However, these features do not propagate when computing weak lensing observables.
We showed how neither the effective field theory formalism nor the standard framework provide a completely satisfactory prediction of the matter power spectrum along the line of sight. One reason for this might be the assumption of a homogeneous and isotropic Universe when we define two-point statistics. Observing within our past light-cone breaks homogeneity along the line of sight, and only the spherical symmetry on the two-dimensional sky is preserved. This is worth investigating in the future. One promising alternative could be kinetic field theory [27]. This theory allows the straightforward computation of highly-accurate non-linear unequal-time correlators on small scales without free fitting parameters, infinite loop corrections or N-body simulations. We also hope that our work motivates simulation groups to investigate alternative methods to include continuous information along the radial direction, so that we can draw strong conclusions on the best formalism.
Finally, the implications on weak lensing are: • The results from the geometric and midpoint approximations mimic perfectly all the features from the exact calculation, washing out any difference at the level of the matter power spectrum.
• The particular choice of the midpoint redshift has negligible impact compared to the projection method. In other words, no matter your prescription of the unequal-time matter power spectrum (as long as the correlation functions have unequal-time information), once the integrals are computed over the entire angular range, the angular correlations would be very similar. To compute these integrals, one would need to drop Limber and use our corfu method [23].
• We also observe that the result from the Limber's approximation deviates from the exact calculation on large-angle separations, < 10. This deviation is of order ∼ 2% maximum for cosmic lensing. For galaxy-galaxy lensing the minimum deviation is ∼ 1%, being very large for some redshift bins. The biggest effect is found for galaxy clustering.
Fang et al. [12], Chisari and Pontzen [19] and Kitching and Heavens [1] already explained that full-sky observables need to be modelled beyond Limber, e.g. intrinsic galaxy alignments and galaxy shapes. Here we showed that this is true not only for upcoming data, but also for current DES Y1 data. It would be interesting to see the impact on the prediction of cosmological parameters and check whether it alleviates some of the tensions between galaxy surveys. Although this is out of the scope of our present work, it would be worth investigating in the future. Our final product is the publicly available python package unequalpy [13] with functionality to reproduce all the results and analysis presented in this paper.
for the equal-time standard perturbation and effective field theory one-loop matter power spectra can be found in that reference. In this appendix, the reader can find the results for their unequal-time counter-parts.
• Standard Perturbation Theory. Within the full-time approach, the linear and the one-loop contributions now read P SPT (k; z 1 , z 2 ) = P 11 (k; z 1 , z 2 ) + P 22 (k; z 1 , z 2 ) + 2P 13 (k; z 1 , z 2 ) (A.1) with P 11 (k; z 1 , z 2 ) = D(z 1 )D(z 2 )P 11 (k) (A.2a) (A.2c) The quantities P i appearing in these expressions can be found in de la Bella et al. [15]. In order to compare with the geometric approximation, one would need to square the expression above.
• Effective Field Theory. Within the full-time approach, the effective field theory unequal-time contribution reads as equation (2.6) but using equation (A.1) instead of the Einstein-de Sitter version. The same applies for the expressions for the squared equal-time power spectrum and the squared of the unequal-time counterpart. Note that using the full time dependence of the non-linear growth functions, the best fit values of the counterterms might vary.

B On why the Limber's and the geometric approximations seem synonyms
It is common to find in the literature the use of the geometric approximation combined with the Limber's approximation. Some authors decide to apply it before using Limber in its Dirac-delta version (3.15), c.f [8]. However, this is totally unnecessary. When applying the Limber's approximation, the use of the geometric approximation is unnecessary. Let us write the angular power spectrum of the weak lensing potential, φ, for two different redshift bins, i and j: dk νk 2 f i φ (ν/k)f j φ (ν/k)P (k; ν/k, ν/k) (B.2) and from (2.9), we are safe to say that P (k; ν/k, ν/k) = P (k; ν/k)P (k; ν/k) is exact. Therefore, equation (B. . It simply comes from measuring the correlation between fields at the same time slice, and therefore it is exact. It naturally comes out from using the Limber's approximation. Therefore, this result is subjected to all its assumptions and there is no need to impose any extra approximation on the unequal-time correlator.
C On why the midpoint and the geometric approximations have similar angular power spectra In this section, we show how the features and differences between the geometric and midpoint approximation get washed out after integrating along the line of sight when computing the angular power spectrum. For simplicity, we work with a particular toy model where the distributions of galaxies are given by Dirac deltas. For this we compute the effective growth functions with ≡ ∂/∂z. Remember, z m = (z 1 + z 2 )/2 and ∆z = (z 1 − z 2 )/2. Figure 11 shows that the difference between these approximations (C.3) only depends on the redshift width for a fixed mean redshift. This is reasonable, given that the midpoint approximation does not use information from the redshift separations. If we now consider a toy model, c.f. 12 where the  which is a very small quantity, since both factors are small themselves. This is shown in Figure 13, where the oscillations come from the product of the spherical Bessel functions. We can conclude that any big differences between approximations at the level of the matter power spectrum, seems to have little impact when integrating along the line of sight.  Figure 13: Difference between the angular power spectra using the geometric and the midpoint approximations, equation (C.7). We see the difference is compatible with zero and oscillations due to the spherical Bessel functions.