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

Lucia F. de la Bella  Nicolas Tessore  and Sarah Bridle

###### Abstract

We investigate the impact of a common approximation of 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 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.

##  1 Introduction

Cosmology is living a golden era of high-precision observations with upcoming
galaxy surveys probing the Universe on unprecedented smaller scales, such as
DESI111https://www.desi.lbl.gov , Euclid222https://www.euclid-ec.org, the
Rubin Observatory 333Formerly known as LSST, https://www.lsst.org and the
Nancy Grace Roman Space Telescope444Formerly known as WFIRST,
https://roman.gsfc.nasa.gov. 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) slices555Note that the
observed large-scale structure does not display full spatial symmetry because
all observations are done within our past light-cone, breaking homogeneity
along the line of sight, and preserving the spherical symmetry on the two-
dimensional sky.. 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

$\omega^{(i,j)}_{AB}(\theta)=\iint
dr_{1}dr_{2}f_{A}^{i}(r_{1})f_{B}^{j}(r_{2})\xi_{AB}(r_{12};r_{1},r_{2})$
(1.1)

with $f^{k}_{C}(r)$ the weight functions for redshift bins $k$.
$\xi_{AB}(r_{12};r_{1},r_{2})$ is the spatial unequal-time 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})$.

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.

Likewise, after expanding in spherical harmonics, the angular power spectrum
yields

$C^{(i,j)}_{AB}(\ell)=\frac{2}{\pi}\int dkk^{2}\iint
dr_{1}dr_{2}f^{i}_{A}(r_{1})f^{j}_{B}(r_{2})j_{\ell}(kr_{1})j_{\ell}(kr_{2})P_{AB}(k;r_{1},r_{2})\;.$
(1.2)

with $j_{\ell}(kr)$ the spherical Bessel function of order $\ell$. 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
unequal-time 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})\approx\sqrt{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 $\theta\lesssim 40^{\circ}$. 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/(\ell+1/2)$. Nontheless, it is not clear
how this would alleviate the divergence of the series expansion at small
$\ell$.

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`666https://github.com/Lucia-Fonseca/unequalpy.git [13] with
functionality to reproduce all the results and analysis presented in this
paper.

##  2 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

$\langle\delta(\vec{k},z)\delta(\vec{k^{\prime}},z)\rangle\equiv(2\pi)^{3}\delta_{\text{D}}(\vec{k}+\vec{k}^{\prime})P(k;z)\;,$
(2.1)

whereas the unequal-time power spectrum reads

$\langle\delta(\vec{k},z_{1})\delta(\vec{k^{\prime}},z_{2})\rangle\equiv(2\pi)^{3}\delta_{\text{D}}(\vec{k}+\vec{k}^{\prime})P(k;z_{1},z_{2})\;.$
(2.2)

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.

###  2.1 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, $\delta=(\rho-\rho_{0})/\rho_{0}$ (with $\rho$ the matter density
field and $\rho_{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
$\delta^{\ast}_{\bm{\mathrm{{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}$,

$\displaystyle\langle\delta_{\bm{\mathrm{{k}}}_{1},1}(z_{1})\delta_{\bm{\mathrm{{k}}}_{2},1}(z_{2})\rangle$
$\displaystyle=(2\pi)^{3}\delta(\bm{\mathrm{{k}}}_{1}+\bm{\mathrm{{k}}}_{2})P_{11}(k,z_{1},z_{2})$
(2.3a)
$\displaystyle\langle\delta_{\bm{\mathrm{{k}}}_{1},2}(z_{1})\delta_{\bm{\mathrm{{k}}}_{2},2}(z_{2})\rangle$
$\displaystyle=(2\pi)^{3}\delta(\bm{\mathrm{{k}}}_{1}+\bm{\mathrm{{k}}}_{2})P_{22}(k,z_{1},z_{2})$
(2.3b)
$\displaystyle\langle\delta_{\bm{\mathrm{{k}}}_{1},1}(z_{1})\delta_{\bm{\mathrm{{k}}}_{2},3}(z_{2})+\delta_{\bm{\mathrm{{k}}}_{1},3}(z_{1})\delta_{\bm{\mathrm{{k}}}_{2},1}(z_{2})\rangle$
$\displaystyle=(2\pi)^{3}\delta(\bm{\mathrm{{k}}}_{1}+\bm{\mathrm{{k}}}_{2})2P_{13}(k,z_{1},z_{2}),$
(2.3c)

where $k$ is the common magnitude of the wavevectors $\bm{\mathrm{{k}}}_{1}$
and $\bm{\mathrm{{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. [16, 15] for a detailed expression of these terms). Then
the one-loop unequal-time matter power spectrum reads

$\begin{split}P_{\mathrm{\tiny
SPT}}(k;z_{1},z_{2})&=D(z_{1})D(z_{2})P_{11}(k)\\\
&+D^{2}(z_{1})D^{2}(z_{2})P_{22}(k)+D(z_{1})D(z_{2})\left[D^{2}(z_{1})+D^{2}(z_{2})\right]P_{13}(k)\end{split}$
(2.4)

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_{\mathrm{\tiny
SPT}}(k,z)=D^{2}(z)P_{11}(k)+D^{4}(z)P_{22}(k)+2D^{4}(z)P_{13}(k)\;.$ (2.5)

###  2.2 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 22-type 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_{NL}$, dropping
order four contributions. Thus all the microscopic physics is encoded in a
redshift-dependent free-fitting parameter, the so-called counterterm,
$c(z)\equiv c_{s}(z)/k_{NL}$. Then, the unequal-time power spectrum in
effective field theory reads

$P_{\mathrm{\tiny EFT}}(k;z_{1},z_{2})=P_{\mathrm{\tiny
SPT}}(k;z_{1},z_{2})-\frac{1}{2}\left(c^{2}(z_{1})+c^{2}(z_{2})\right)D(z_{1})D(z_{2})k^{2}P_{11}(k)$
(2.6)

Finally, the equal-time power spectrum is retrieved when $z_{1}=z_{2}=z$ in
equation (2.6)

$P_{\mathrm{\tiny EFT}}(k,z)=P_{\mathrm{\tiny
SPT}}(k,z)-c^{2}(z)D^{2}(z)k^{2}P_{11}(k)\;.$ (2.7)

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.

####  2.2.1 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`777https://emcee.readthedocs.io. We use matter power spectrum data
up to $k=0.4h/\text{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:

$c^{2}(z)=m\textrm{e}^{-az}+n$ (2.8)

with the best fit values $m=2.564^{+0.008}_{-0.007}\text{Mpc}^{2}/h^{2}$,
$n=0.036^{+0.001}_{-0.001}\text{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 CAMB888https://camb.info.

Figure 2: The left plot shows the parametrization of the counterterms as a
function of redshift. The plot on the right represents the correlation between
the parameters involved in the time evolution of the counterterms. In this
case, the slope and intercept in equation (2.8) show the weakest correlation.

In the next section, we analyse which of the above prescriptions provides the
best framework to analyse unequal-time correlators, and weak lensing
observables.

###  2.3 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,

$P(k;z_{1},z_{2})^{2}=P(k,z_{1})P(k,z_{2})+\Delta P(k;z_{1},z_{2})\;.$ (2.9)

where $\Delta 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.

| Equal-time
---
SPT
Equation (2.5)
| Equal-time
---
EFT
Equation (2.7)
| Unequal-time
---
SPT
Equation (2.4)
| Unequal-time
---
EFT
Equation (2.6)
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:

$\epsilon_{geom}=\frac{|P_{\mathrm{\tiny SPT}}(k,z_{1})P_{\mathrm{\tiny
SPT}}(k,z_{2})-P_{\mathrm{\tiny EFT}}(k,z_{1})P_{\mathrm{\tiny
EFT}}(k,z_{2})|}{P_{\mathrm{\tiny EFT}}(k,z_{1})P_{\mathrm{\tiny
EFT}}(k,z_{2})}$ (2.10) $\epsilon_{uet}=\frac{|P_{\mathrm{\tiny
SPT}}(k,z_{1},z_{2})^{2}-P_{\mathrm{\tiny
EFT}}(k,z_{1},z_{2})^{2}|}{P_{\mathrm{\tiny EFT}}(k,z_{1},z_{2})^{2}}\;.$
(2.11)

We calculate this error for a given mean redshift $z_{m}=0.5$ and different
redshift separations: $\Delta z=0$ (same time slice) and $\Delta 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, $\Delta z=0$, the
absolute error does not vanish, it becomes proportional to the counterterms
multiplied by the power spectrum, $\propto 2c^{2}(z_{m})P_{\mathrm{\tiny
SPT}}(k,z_{m})$. For increasingly larger redshift separations, the relative
error tends to shift to the left and become larger, preserving the shape.

#####  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

$\epsilon_{theory}=\frac{|P_{\mathrm{\tiny theory}}(k,z_{1})P_{\mathrm{\tiny
theory}}(k,z_{2})-P_{\mathrm{\tiny
theory}}(k,z_{1},z_{2})^{2}|}{P_{\mathrm{\tiny theory}}(k,z_{1},z_{2})^{2}}$
(2.12)

with $theory=\\{SPT,EFT\\}$. Again we analyse the error for $z_{m}=0.5$ and
widths $\Delta 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, $\Delta z=0$, the error
vanishes. This is true because the unequal-time correlator equals the
geometric approximation when $z_{1}=z_{2}$.

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, $\Delta 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.

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.

###  2.4 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

$P(k;z_{1},z_{2})\approx P(k,z_{m})$ (2.13)

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 $\Delta 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

$\epsilon_{mid}^{spt}=\frac{|P_{\mathrm{\tiny
SPT}}(k,z_{m})^{2}-P_{\mathrm{\tiny
SPT}}(k;z_{1},z_{2})^{2}|}{P_{\mathrm{\tiny SPT}}(k;z_{1},z_{2})^{2}}\;.$
(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\approx 0.1h/\text{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 $\Delta 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.

Figure 4: These plots represent the relative error committed when using the
traditional geometric approximation (red lines), and the mean redshift
approximation in standard perturbation theory (black lines) instead of the
unequal-time description, equation. On the left, we plot the error in the
squared power spectrum for the midpoint approximation (black lines), equation
(2.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, $\Delta 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.

##  3 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(\theta)$ 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^{(i,j)}_{ab}(\theta)=\iint_{0}^{\infty}dx_{1}\,dx_{2}\,f^{i}_{a}(x_{1})f^{j}_{b}(x_{2})\,\xi(r_{12};t(x_{1}),t(x_{2}))\;,$
(3.1)

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
$\xi(r_{12};t(x_{1}),t(x_{2}))$ is the unequal-time matter correlation
function for the separation

$r_{12}=\sqrt{x_{1}^{2}+x_{2}^{2}-2x_{1}x_{2}\cos\theta}$ (3.2)

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
$\delta(\vec{r};t)$ in comoving coordinates $\vec{r}$ and cosmic time $t$, the
convergence $\kappa(\vec{\Theta})$ in direction $\vec{\Theta}$ is the integral

$\kappa^{i}(\vec{\Theta})=\int_{0}^{\infty}\\!f^{i}_{\kappa}(x)\,\delta(x\vec{\Theta};t(x))\,dx$
(3.3)

with the filter defined as

$f^{i}_{\kappa}(x)=\frac{3H_{0}^{2}\Omega_{m}}{2c^{2}}\frac{q_{i}(x)\,x}{a(x)}$
(3.4)

where $H_{0}$ and $\Omega_{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

$q_{i}(x)=\int_{x}^{\infty}\\!\frac{x^{\prime}-x}{x^{\prime}}\,n_{i}(x^{\prime})\,dx^{\prime}$
(3.5)

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 $\delta(\vec{r};t)$, the galaxy density field
$\delta_{g}(\vec{\Theta};t)$ is the integral

$\delta_{g}^{i}(\vec{\Theta};t)=\int_{0}^{\infty}\\!f^{i}_{g}(x)\,\delta(x\vec{\Theta};t(x))\,dx$
(3.6)

with the filter defined as

$f^{i}_{g}(x)=\frac{1}{c}b(x)n_{i}(x)H(x)$ (3.7)

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 $\kappa$ 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

$w^{(i,j)}_{\kappa\kappa}(\theta)=\iint_{0}^{\infty}dx_{1}\,dx_{2}\,f^{i}_{\kappa}(x_{1})f^{j}_{\kappa}(x_{2})\,\xi(r_{12};t_{1},t_{2})\;,$
(3.8)

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}^{\kappa\kappa}$ and $C_{l}^{\gamma\gamma}$, with

$C^{(i,j)}_{\gamma\gamma}(\ell)=\frac{(l-1)\,(l+2)}{l\,(l+1)}\,C^{(i,j)}_{\kappa\kappa}(\ell)\;.$
(3.9)

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:

$w^{(i,j)}_{gg}(\theta)=\iint_{0}^{\infty}dx_{1}\,dx_{2}\,f^{i}_{g}(x_{1})f^{j}_{g}(x_{2})\,\xi(r_{12};t_{1},t_{2})\;,$
(3.10)

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:

$w^{(i,j)}_{\kappa
g}(\theta)=\iint_{0}^{\infty}dx_{1}\,dx_{2}\,f^{i}_{\kappa}(x_{1})f^{j}_{g}(x_{2})\,\xi(r_{12};t_{1},t_{2})\;,$
(3.11)

using the filters (3.4) and (3.7).

###  3.1 Computation of the angular power spectrum

There are different approaches to compute the angular power spectra in cosmic
lensing (1.2):

  1. 1.

Brute-force calculation:

$\begin{split}C^{(i,j)}_{ab}(\ell)=\frac{2}{\pi}\int_{0}^{\infty}dkk^{2}\iint_{0}^{\infty}dx_{1}dx_{2}\,f^{i}_{a}(x_{1})f^{j}_{b}(x_{2})j_{\ell}(kx_{1})j_{\ell}(kx_{2})P\left(k;x_{1},x_{2}\right)\;.\end{split}$
(3.12)

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. 2.

Approximations to the unequal-time matter power spectrum, which do not reduce
the number of integrals.

     * •

The geometric approximation (1.3):

$\begin{split}C^{(i,j)}_{ab}(\ell)\approx\frac{2}{\pi}\int_{0}^{\infty}dkk^{2}&\int_{0}^{\infty}dx_{1}\,f^{i}_{a}(x_{1})j_{\ell}(kx_{1})\sqrt{P(k;x_{1})}\\\
\times&\int_{0}^{\infty}dx_{2}f^{j}_{b}(x_{2})j_{\ell}(kx_{2})\sqrt{P(k;x_{2})}\;.\end{split}$
(3.13)

     * •

The new midpoint approximation (2.13):

$\begin{split}C^{(i,j)}_{ab}(\ell)\approx\frac{2}{\pi}\int_{0}^{\infty}dkk^{2}\iint_{0}^{\infty}dx_{1}dx_{2}\,f^{i}_{a}(x_{1})f^{j}_{b}(x_{2})j_{\ell}(kx_{1})j_{\ell}(kx_{2})P\left(k;x_{m}\right)\;.\end{split}$
(3.14)

where $x_{m}$ corresponds to the radial distance at the midpoint redshift
$z_{m}=\frac{z_{1}+z_{2}}{2}$.

  3. 3.

Assumptions about the filters: the Limber’s approximation.

     * •

In harmonic space, one can use the Dirac-delta version of the Limber’s
approximation where the spherical Bessel functions are approximated to the
amplitude of their first peak

$j_{\ell}(kx)\simeq\sqrt{\pi/(2\nu)}\delta(\nu-kx)$ (3.15)

with $\nu\equiv\ell+1/2$, c.f. [20, 21]. Some authors also apply the geometric
approximation before using equation (3.15) and call that “Limber’s
approximation”, c.f. [8]. 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

$C^{(i,j)}_{ab}(\ell)\approx\frac{1}{\nu}\int_{0}^{\infty}\\!dkk^{2}f^{i}_{a}(\nu/k)f^{j}_{b}(\nu/k)P(k,\nu/k)\;.$
(3.16)

     * •

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 $R=x_{1}-x_{2}$

$w(\theta)=\int_{0}^{\infty}\\!\int_{-2x}^{2x}\\!f_{1}(x+R/2)\,f_{2}(x-R/2)\,\xi(r_{12};t_{1},t_{2})\,dR\,dx\;,$
(3.17)

where the distance between the points in terms of $x$ and $r_{12}$ is now

$r_{12}=\sqrt{2x^{2}\,(1-\cos\theta)+R^{2}\,(1+\cos\theta)/2}\;.$ (3.18)

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,

$f_{1}(x+R/2)\,f_{2}(x-R/2)\,\xi(r_{12};t_{1},t_{2})\approx
f_{1}(x)\,f_{2}(x)\,\xi(r_{12};t)\;,$ (3.19)

_ii)_ that the angle separation $\theta$ between the points is small, so that
the distance $r_{12}$ can be approximated as

$r_{12}\approx\sqrt{x^{2}\theta^{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

$w_{\rm L}(\theta)=\int_{0}^{\infty}\\!f_{1}(x)\,f_{2}(x)\,\xi_{\rm
L}(x\theta;t)\,x\theta\,dx\;,$ (3.21)

where $\xi_{\rm L}$ is Limber’s matter correlation function, defined as

$\xi_{\rm
L}(r;t)=\frac{1}{r}\,\int_{-\infty}^{\infty}\\!\xi(\sqrt{r^{2}+R^{2}};t)\,dR\;.$
(3.22)

  4. 4.

In this work, we employ a novel approach by computing the angular correlation
function in real space,

$w(\theta)=\iint_{0}^{\infty}\\!f_{1}(x_{1})\,f_{2}(x_{2})\,\xi(r_{12};t_{1},t_{2})\,dx_{1}\,dx_{2}\;,$
(3.23)

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,

$\xi(r;t_{1},t_{2})=\frac{1}{2\pi^{2}}\int_{0}^{\infty}\\!P(k;t_{1},t_{2})\,\frac{\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,

$w(\theta)=\sum_{l}\frac{2l+1}{4\pi}\,C_{l}\,\mathcal{P}_{l}(\cos\theta)\;.$
(3.25)

By evaluating the angular correlation function over a grid of $\theta$ 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
corfu999https://github.com/ntessore/corfu package for Python [23].

###  3.2 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].

#### DES Y1 redshift distributions

We obtain all DES Y1 quantities by running the `des-y1-test` in
`CosmoSIS`101010https://bitbucket.org/joezuntz/cosmosis, without intrinsic
alignment, c.f Figure 5. Specifications:

  * •

Cosmology:

$\\{\Omega_{m}=0.2678,H_{0}=67.5,\Omega_{b}=0.0483,n_{s}=0.965,A_{s}=2.1\times
10^{-9},w=-1\\}.$

  * •

Convergence filter (3.4): we use the source redshift distribution of galaxies
for four different redshift bins.

  * •

Galaxy filter (3.7): we use the lens sample of redshift distribution of
galaxies and bias for five different redshift bins. The bias parameter at each
redshift:

$b_{DES}(z)=\\{1.45,1.55,1.65,1.8,2.0\\}\;.$

Figure 5: Filters for the source and lens samples for DES-Y1 redshift
distributions.

#### 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}\in(0.2,1.2)$ with $\Delta
z=0.1$ convolving each bin with a Gaussian photo-z scatter

$\sigma_{z}=0.03(1+z)\;.$ (3.26)

The parametric redshift distribution reads

$n(z)\propto z^{2}e^{-(z/z_{0})^{\alpha}}$ (3.27)

with parameters $(z_{0},\alpha)=(0.28,0.90)$. Then the true distribution of
galaxies $n_{i}(z)$ that fall in the $i$-th photo-z bin is defined

$n_{i}(z)=\int_{z_{ph}^{i}}^{z_{ph}^{i+1}}dz_{ph}n(z)p(z_{ph}|z)$ (3.28)

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:

$b_{LSST}(z_{ph})=\frac{b}{D(z_{ph})}$ (3.29)

with $b=0.95$. $D(z)$ 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:

$\sigma_{z}=0.05(1+z)\;.$ (3.30)

The parameters in (3.27) are now $(z_{0},\alpha)=(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.

Figure 6: Cosmic lensing. Angular correlation for the whole range of
approximations and the fuller calculation at all angular scales. For the
unequal-time (black), the geometric (red) and the midpoint (green)
approximations, positive values are represented by solid lines and negative
values by dashed lines. The Limber’s prediction is represented by blue dots
(filled for positive and open for negative values). Figure 7: Cosmic lensing.
Comparison between the angular correlation function when using an
approximation and the fuller calculation at all angular scales for one of the
redshift bins, e.g. $(1,1)$. The Limber’s approximation does not provide a
good prediction on large-angle separation. The midpoint and geometric
approximation seems to deal equally well with large-angle correlations,
although they deviate from the full unequal-time description. This effect
propagates to the angular power spectrum, showing as deviations on small
multipoles, $\ell$. The feature at $\theta\approx 9\times 10^{2}$ (arcmin)
shows the change of sign in the angular correlation function. Figure 8: Cosmic
lensing. Comparison between the angular power spectrum when using an
approximation and the fuller calculation at all angular scales. Figure 9:
Galaxy-galaxy lensing. Otherwise the same as Figure 5. Figure 10: Galaxy
clustering. Otherwise the same as Figure 5.

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 $\ell$. 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, $\ell<10$. This is
something that we already knew, this effect was predicted by Simon [9]. This
deviation is of order $\sim 2\%$ maximum for cosmic lensing. For galaxy-galaxy
lensing the minimum deviation is $\sim 1\%$ and it can be very large for some
redshift bins. The biggest effect is found for galaxy clustering.

##  4 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 perturbation 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, $\ell<10$. This deviation is
of order $\sim 2\%$ maximum for cosmic lensing. For galaxy-galaxy lensing the
minimum deviation is $\sim 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.

## Acknowledgments

We would like to acknowledge all of the insightful comments from our
colleagues at the University of Manchester and the SkyPy Collaboration. We
specially thank I. Harrison and J. P. Cordero for their insights on Bayesian
analysis and redshift distributions. We also thank N. MacCrann for his help
with `CosmoSIS`. We also acknowledge A. Heavens, E. Chisari and A. Pontzen,
and D. Thomas for their insightful comments and feedback.

The preparation of this manuscript was made possible by a number of software
packages: `NumPy`, `SciPy` [28], `Astropy` [29], `Matplotlib` [30],
`IPython/Jupyter` [31], `emcee` [32]. This research made use of `SkyPy`, a
Python package for forward modeling astronomical surveys [17], and `CosmoSIS`
[33]. This work produced and made use of `corfu` [23] and `unequalPy` [13].

The work reported in this paper has been supported by the European Research
Council in the form of a Consolidator Grant with Grant Agreement No. 681431
(LFdlB, NT, SB).

LFdlB would also like to thank all NHS and key workers around the world for
their unstoppable work during the global pandemic. LFdlB would also like to
add her support to the “Black Lives Matter”111111https://blacklivesmatter.com
movement and their fight for human rights against racism inside and outside
academia.

##  Appendix A Full-time matter power spectrum

De la Bella et al. [15] proved the Einstein-De Sitter approximation to lose
percent-level accuracy with respect to the full-time computation of the growth
functions. The calculations 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_{\mathrm{\tiny
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

$\displaystyle P_{11}(k;z_{1},z_{2})$
$\displaystyle=D(z_{1})D(z_{2})P_{11}(k)$ (A.2a)
$\displaystyle\begin{split}P_{22}(k;z_{1},z_{2})&=D_{A}(z_{1})D_{A}(z_{2})P_{AA}(k)+\Big{(}D_{A}(z_{1})D_{B}(z_{2})+D_{A}(z_{2})D_{B}(z_{1})\Big{)}P_{AB}(k)\\\
&+D_{B}(z_{1})D_{B}(z_{2})P_{BB}(k)\end{split}{}$ (A.2b)
$\begin{split}P_{13}(k;z_{1},z_{2})=\frac{1}{2}P^{\ast}(k)\Bigg{\\{}&\Big{[}D(z_{1})\Big{(}D_{D}(z_{2})-D_{J}(z_{2})\Big{)}+D(z_{2})\Big{(}D_{D}(z_{1})-D_{J}(z_{1})\Big{)}\Big{]}P_{D}(k)\\\
+&\Big{[}D(z_{1})D_{E}(z_{2})+D(z2)D_{E}(z_{1})\Big{]}P_{E}(k)\\\
+&\Big{[}D(z_{1})\Big{(}D_{F}(z2)+D_{J}(z_{2})\Big{)}+D(z_{2})\Big{(}D_{F}(z_{1})+D_{J}(z_{1})\Big{)}\Big{]}P_{F}(k)\\\
+&\Big{[}D(z_{1})D_{G}(z_{2})+D(z_{2})D_{G}(z)\Big{]}P_{G}(k)\\\
+&\frac{1}{2}\Big{[}D(z_{1})D_{J}(z_{2})+D(z2)D_{J}(z_{1})\Big{]}\Big{[}P_{J2}(k)-2P_{J1}(k)\Big{]}\Bigg{\\}}.\end{split}$
(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.

##  Appendix 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, $\phi$,
for two different redshift bins, $i$ and $j$:

$C^{(i,j)}_{\phi\phi}(\ell)=\frac{2}{\pi}\int_{0}^{\infty}\frac{dk}{k^{2}}\,\iint_{0}^{\infty}dx_{1}\,dx_{2}\,f^{i}_{\phi}(x_{1})f^{j}_{\phi}(x_{2})\,j_{\ell}(kx_{1})j_{\ell}(kx_{2})\,P(k;x_{1},x_{2})\;.$
(B.1)

By not using the geometric approximation (1.3) and simply applying the
Limber’s approximation in the form (3.15), the equation above yields

$C^{(i,j)}_{\phi\phi}(\ell)\approx\int_{0}^{\infty}\\!\frac{dk}{\nu
k^{2}}f^{i}_{\phi}(\nu/k)f^{j}_{\phi}(\nu/k)P(k;\nu/k,\nu/k)$ (B.2)

and from (2.9), we are safe to say that

$\begin{split}P(k;\nu/k,\nu/k)=&\sqrt{P(k;\nu/k)P(k;\nu/k)}\\\
=&P(k;\nu/k)\end{split}$ (B.3)

is exact. Therefore, equation (B.2) reads

$C^{(i,j)}_{\phi\phi}(\ell)\approx\int_{0}^{\infty}\\!\frac{dk}{\nu
k^{2}}f^{i}_{\phi}(\nu/k)f^{j}_{\phi}(\nu/k)P(k;\nu/k)\;.$ (B.4)

Note here, equation (B.3) is not the geometric approximation (1.3). 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.

##  Appendix 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

  * •

Equal-time version

$D_{eff}(k,z)^{2}=\frac{P(k,z)}{P(k,0)}$ (C.1)

  * •

Unequal-time version

$D_{eff}(k,z_{1},z_{2})^{2}=\frac{P(k,z_{1},z_{2})}{P(k,0)}\;.$ (C.2)

Then the difference between the midpoint and geometric approximations at the
level of the power spectrum becomes

$|P_{geom}-P_{mid}|=|D_{eff}(k,z_{1})D_{eff}(k,z_{2})-D_{eff}^{2}(k,z_{m})|P(k,0)\;.$
(C.3)

We now perform a Taylor expansion of the effective growth functions around the
mean redshift value, $z_{m}$:

$D_{eff}(k,z)\approx D_{eff}(k,z_{m})+(z-z_{m})\frac{\partial
D_{eff}(k,z)}{\partial
z}\bigg{|}_{z_{m}}+(z-z_{m})^{2}\frac{\partial^{2}D_{eff}(k,z)}{\partial
z^{2}}\bigg{|}_{z_{m}}$ (C.4)

Therefore, at leading order

$|D_{eff}(k,z_{1})D_{eff}(k,z_{2})-D_{eff}^{2}(k,z_{m})|\approx\left(\Delta
zD_{eff}^{\prime}(k,z_{m})\right)^{2}$ (C.5)

with ${}^{\prime}\equiv\partial/\partial z$. Remember, $z_{m}=(z_{1}+z_{2})/2$
and $\Delta 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.

Figure 11: The the difference between the geometric and midpoint
approximations depends on the redshift width for a fixed mean redshift,
equation (C.5).

If we now consider a toy model, c.f. 12 where the redshift distribution of
galaxies follow this law

$n(z)=\delta_{D}(z-z_{a})+\delta_{D}(z-z_{b})$ (C.6)

in a way that $z_{m}=(z_{a}+z_{b})/2$ and $\Delta z=(z_{a}-z_{b})/2$.

Figure 12: The redshift distribution of galaxies for our toy model.

When we compute the angular power spectra for both approximations, the
difference reads

$|C_{geom}(\ell)-C_{mid}(\ell)|\approx\left(\Delta
zD_{eff}^{\prime}(k,z_{m})\right)^{2}\int\frac{dk}{k^{2}}j_{\ell}(kz_{a})j_{\ell}(kz_{b})P(k,0)$
(C.7)

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.

## References

  * [1] T. D. Kitching and A. F. Heavens, _Unequal-Time Correlators for Cosmology_ , _Phys. Rev._ D95 (2017) 063522 [1612.00770].
  * [2] N. Tessore and L. F. de la Bella, _Cosmic lensing two-point functions with exact projection_ , .
  * [3] N. E. Chisari and A. Pontzen, _Unequal time correlators and the Zel’dovich approximation_ , _Phys. Rev._ D100 (2019) 023543 [1905.02078].
  * [4] D. N. Limber, _The Analysis of Counts of the Extragalactic Nebulae in Terms of a Fluctuating Density Field._ , _APJ_ 117 (1953) 134.
  * [5] D. N. Limber, _The Analysis of Counts of the Extragalactic Nebulae in Terms of a Fluctuating Density Field. II._ , _”APJ”_ 119 (1954) 655.
  * [6] P. J. E. Peebles, _Statistical Analysis of Catalogs of Extragalactic Objects. I. Theory_ , _”APJ”_ 185 (1973) 413.
  * [7] N. Kaiser, _Weak lensing and cosmology_ , _The Astrophysical Journal_ 498 (1998) 26–42.
  * [8] P. Lemos, A. Challinor and G. Efstathiou, _The effect of Limber and flat-sky approximations on galaxy weak lensing_ , _JCAP_ 05 (2017) 014 [1704.01054].
  * [9] P. Simon, _How accurate is Limber’s equation?_ , _Astron. Astrophys._ 473 (2007) 711 [astro-ph/0609165].
  * [10] T. D. Kitching, J. Alsing, A. F. Heavens, R. Jimenez, J. D. McEwen and L. Verde, _The limits of cosmic shear_ , _Monthly Notices of the Royal Astronomical Society_ 469 (2017) 2737–2749.
  * [11] M. LoVerde and N. Afshordi, _Extended Limber Approximation_ , _Phys. Rev. D_ 78 (2008) 123506 [0809.5112].
  * [12] X. Fang, E. Krause, T. Eifler and N. MacCrann, _Beyond Limber: Efficient computation of angular power spectra for galaxy clustering and weak lensing_ , _JCAP_ 05 (2020) 010 [1911.11947].
  * [13] L. F. de la Bella, _unequalpy_ , Nov., 2020.  10.5281/zenodo.4268314.
  * [14] J. J. M. Carrasco, M. P. Hertzberg and L. Senatore, _The Effective Field Theory of Cosmological Large Scale Structures_ , _JHEP_ 09 (2012) 082 [1206.2926].
  * [15] L. F. de la Bella, D. Regan, D. Seery and S. Hotchkiss, _The matter power spectrum in redshift space using effective field theory_ , _JCAP_ 1711 (2017) 039 [1704.05309].
  * [16] F. Bernardeau, S. Colombi, E. Gaztanaga and R. Scoccimarro, _Large scale structure of the universe and cosmological perturbation theory_ , _Phys. Rept._ 367 (2002) 1 [astro-ph/0112551].
  * [17] S. Collaboration, _Skypy_ , Apr., 2020.  10.5281/zenodo.4071945.
  * [18] F. Villaescusa-Navarro et al., _The Quijote simulations_ , 1909.05273.
  * [19] N. E. Chisari, D. Alonso, E. Krause, C. D. Leonard, P. Bull, J. Neveu et al., _Core cosmology library: Precision cosmological predictions for lsst_ , _The Astrophysical Journal Supplement Series_ 242 (2019) 2.
  * [20] N. Afshordi, Y.-S. Loh and M. A. Strauss, _Cross - correlation of the Cosmic Microwave Background with the 2MASS galaxy survey: Signatures of dark energy, hot gas, and point sources_ , _Phys. Rev. D_ 69 (2004) 083524 [astro-ph/0308260].
  * [21] LSST Dark Energy Science collaboration, _Core Cosmology Library: Precision Cosmological Predictions for LSST_ , _Astrophys. J. Suppl._ 242 (2019) 2 [1812.05995].
  * [22] A. J. S. Hamilton, _Uncorrelated modes of the non-linear power spectrum_ , _”MNRAS”_ 312 (2000) 257 [astro-ph/9905191].
  * [23] N. Tessore and L. F. de la Bella, _corfu_ , Nov., 2020.  10.5281/zenodo.4268486.
  * [24] T. Abbott, F. B. Abdalla, A. Alarcon, J. Aleksić, S. Allam, S. Allen et al., _Dark energy survey year 1 results: Cosmological constraints from galaxy clustering and weak lensing_ , _Physical Review D_ 98 (2018) .
  * [25] T. L. D. E. S. Collaboration, R. Mandelbaum, T. Eifler, R. Hložek, T. Collett, E. Gawiser et al., _The lsst dark energy science collaboration (desc) science requirements document_ , 2018.
  * [26] Z. Ma, W. Hu and D. Huterer, _Effects of photometric redshift uncertainties on weak‐lensing tomography_ , _The Astrophysical Journal_ 636 (2006) 21–29.
  * [27] M. Bartelmann, E. Kozlikin, R. Lilow, C. Littek, F. Fabis, I. Kostyuk et al., _Cosmic structure formation with kinetic field theory_ , _Annalen der Physik_ 531 (2019) 1800446.
  * [28] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau et al., _SciPy 1.0: fundamental algorithms for scientific computing in Python_ , _Nature Methods_ 17 (2020) 261 [1907.10121].
  * [29] A. M. Price-Whelan, B. M. Sipőcz, H. M. Günther, P. L. Lim, S. M. Crawford, S. Conseil et al., _The astropy project: Building an open-science project and status of the v2.0 core package_ , _The Astronomical Journal_ 156 (2018) 123.
  * [30] J. D. Hunter, _Matplotlib: A 2D Graphics Environment_ , _Computing in Science and Engineering_ 9 (2007) 90.
  * [31] F. Perez and B. E. Granger, _IPython: A System for Interactive Scientific Computing_ , _Computing in Science and Engineering_ 9 (2007) 21.
  * [32] D. Foreman-Mackey, D. W. Hogg, D. Lang and J. Goodman, _emcee: The MCMC Hammer_ , _PASP_ 125 (2013) 306 [1202.3665].
  * [33] J. Zuntz, M. Paterno, E. Jennings, D. Rudd, A. Manzotti, S. Dodelson et al., _Cosmosis: Modular cosmological parameter estimation_ , _Astronomy and Computing_ 12 (2015) 45–59.