# Full sky Integrated Sachs-Wolfe maps for the MICE Grand Challenge lightcone
simulation

Krishna Naidoo,1,2 Pablo Fosalba,3,4 Lorne Whiteway2 and Ofer Lahav2

1Center for Theoretical Physics, Polish Academy of Sciences, Al. Lotników
32/46, 02-668 Warsaw, Poland
2Department of Physics & Astronomy, University College London, Gower Street,
London WC1E 6BT, UK
3Institut d’Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain
4Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans,
s/n, 08193 Barcelona, Spain
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present full-sky maps of the Integrated Sachs-Wolfe effect (ISW) for the
MICE Grand Challenge lightcone simulation up to redshift $1.4$. The maps are
constructed in the linear regime using spherical Bessel transforms. We compare
and contrast this procedure against analytical approximations found in the
literature. By computing the ISW in the linear regime we remove the
substantial computing and storage resources required to calculate the non-
linear Rees-Sciama effect. Since the linear ISW is $\sim 10^{2}\,\times$
larger in $\Delta T/T$ this has a neglegible impact on the maps produced and
only becomes relevant on scales which are dominated by cosmic microwave
background (CMB) anisotropies. The MICE simulation products have been
extensively used for studies involving current and future galaxy surveys. The
availability of these maps will allow MICE to be used for future galaxy and
CMB cross-correlation studies, ISW reconstruction studies and ISW void-
stacking studies probed by galaxy surveys such as DES, DESI, _Euclid_ and
Rubin LSST. The pipeline developed in this study is provided as a public
Python package pyGenISW. This could be used in future studies for constructing
the ISW from existing and future simulation suites probing vast sets of
cosmological parameters and models.

###### keywords:

cosmic background radiation – large-scale structure of Universe – methods:
numerical

††pubyear: 2020††pagerange: Full sky Integrated Sachs-Wolfe maps for the MICE
Grand Challenge lightcone simulation–Full sky Integrated Sachs-Wolfe maps for
the MICE Grand Challenge lightcone simulation

##  1 Introduction

The Integrated Sachs-Wolfe effect (ISW; Sachs & Wolfe, 1967), caused by the
evolution of gravitational potentials in large scale structure (LSS), imprints
features onto the cosmic microwave background (CMB). The strength of these
features is sensitive to the underlying cosmological model, in particular the
quantity and nature of dark energy (Crittenden & Turok, 1996). However, the
CMB is dominated by primordial anisotropies meaning the ISW is only detectable
by cross-correlating the CMB with tracers of LSS. This has been performed on
several galaxy surveys to constrain the standard cosmological model $\Lambda$
Cold Dark Matter ($\Lambda$CDM) and to test expanded/alternative models
(Peiris & Spergel, 2000; Boughn & Crittenden, 2002; Fosalba et al., 2003;
Scranton et al., 2003; Nolta et al., 2004; Afshordi et al., 2004; Corasaniti
et al., 2005; Padmanabhan et al., 2005; Giannantonio et al., 2006; Vielva et
al., 2006; McEwen et al., 2007; Giannantonio et al., 2008; Ho et al., 2008;
Xia et al., 2009; Dupé et al., 2011; Giannantonio et al., 2012; Goto et al.,
2012; Giannantonio et al., 2014; Barreiro et al., 2013; Planck Collaboration
et al., 2014; Cabass et al., 2015; Ferraro et al., 2015; Nicola et al., 2016;
Planck Collaboration et al., 2016; Stölzner et al., 2018; Hang et al., 2020).

Future galaxy redshift surveys will generate enormous catalogues of the
position and redshift of galaxies providing a large observational data set
with which we may conduct cross-correlation studies of LSS with the CMB. In
order to understand the observational sensitivities and systematics of the ISW
to cosmological parameters we need to be able to construct ISW maps and
corresponding galaxy mocks for a wide range of cosmological parameters and
models. At present this is computationally expensive as accurate ISW maps
require regular snapshots of the gravitational potential. This calculation
therefore requires some foresight and cannot be done ad-hoc after the
simulations have been run. For this reason only a handful of such simulations
exist (see Cai et al., 2010; Watson et al., 2014; Carbone et al., 2016; Adamek
et al., 2020) using either the best fit Planck cosmology or exploring a
limited set of cosmological parameters/models.

This limitation is most striking for ISW void-stacking studies (Granett et
al., 2008; Pápai et al., 2011; Nadathur et al., 2012; Flender et al., 2013;
Hernández-Monteagudo & Smith, 2013; Ilić et al., 2013; Granett et al., 2015;
Kovács & García-Bellido, 2016; Nadathur & Crittenden, 2016; Cai et al., 2017;
Kovács, 2018; Kovács et al., 2019) which, to compare to theoretical
predictions, require the simulation of realistic galaxy catalogues with
corresponding ISW maps. Some of these studies have measured a $\gtrsim
2\sigma$ excess in the ISW from large voids, the source of which remains
unclear, but due to the computational cost of running ISW simulations have
been limited to comparisons to the fiducial $\Lambda$CDM model (frequently
comparing to the Jubilee ISW maps; Watson et al., 2014).

The discovery of a void along the line-of-sight (LOS) of the CMB Cold Spot
(CS) anomaly (Szapudi et al., 2015) led Kovács & García-Bellido (2016) to
speculate whether the anomaly was actually caused by the same anomalous excess
found for the ISW of other large voids. Naidoo et al. (2016) considered a
$\Lambda$CDM solution to the CS anomaly – multiple voids along the LOS (later
discovered by Mackenzie et al., 2017). However, even in the most extreme of
scenarios multiple voids were unable to explain the full CS profile. Nadathur
et al. (2014) showed that the significance of the CS is low ($\sim 2\sigma$)
and could be explained as a tail-end Gaussian fluctuation. The significance of
the CS was placed into further question when the effect of masking was
considered and found to enhance the CS significance. By removing this effect
the CS was found to be significant at only $\sim 1.9\sigma$ (Naidoo et al.,
2017). If the CS is evidence of the same anomalous excess in the void ISW
signal, it will be difficult to establish from measurements of the CS alone.
Instead, studies will need to further explore the signals from voids in data,
using new and larger data sets, and by comparing to ISW simulations.

To explore the ISW for a larger set of cosmological parameters and models,
simulated ISW maps need to be easier to produce (so that these maps can be
produced ad-hoc for a large set of existing and future $N$-body simulations,
without pre-planning). To do this we construct the ISW in the linear regime,
removing the costly requirement for regular snapshots of the gravitational
potential across cosmic time. This comes at a small cost: we lose the non-
linear ISW (known as the Rees-Sciama effect; Rees & Sciama, 1968), but since
the ISW is $\sim 10^{2}$ times larger in temperature (Nadathur et al., 2014)
this effect is neglegible for the scales of interest in observational
cosmology (which, due to the low signal-to-noise of the measurement, are
generally limited to spherical harmonic modes of $\ell<100$; Hang et al.
2020). For $\ell\gtrsim 500$ the Rees-Sciama effect starts to dominate (Cai et
al., 2009) and for $\ell\gtrsim 5000$ it is the dominant source of CMB
anisotropies (Seljak, 1996).

The ISW is constructed for the MICE Grand Challenge lightcone simulation
(Fosalba et al., 2015b) using a spherical Bessel transform method (Shapiro et
al., 2012) conducted in spherical polar coordinates – a natural coordinate
system for lightcone simulations provided in spherical shells (Fosalba et al.,
2008).

The paper is organised as follows: in Section 2 we describe the MICE
simulation data used and the ISW construction methods implemented in this
study; in Section 3 we compare the ISW maps and their statistics to each other
and theoretical expectations; lastly, in Section 4 we discuss the results and
the relevance of the data and pipeline produced in this study for future work.

##  2 Method

###  2.1 Data: The MICE Density Field

The MICE Grand Challenge lightcone simulation (Fosalba et al., 2015b; Crocce
et al., 2015; Fosalba et al., 2015a; Carretero et al., 2015; Hoffmann et al.,
2015) is a large $N$-body simulation constructed with Gadget-2 (Springel,
2005). The simulation was run with a comoving box of length $3072h^{-1}{\rm
Mpc}$, $4096^{3}$ dark matter particles, with a particle-mesh grid of $4096$
used to calculate large-scale forces computed with Fast Fourier Transforms.
The simulation used a flat $\Lambda$CDM cosmological model consistent with the
best-fit WMAP 5-year data (Dunkley et al., 2009) – i.e. $\Omega_{\rm
m,0}=0.25$, $\Omega_{\Lambda,0}=0.75$, $\Omega_{\rm b,0}=0.044$, $n_{\rm
s}=0.95$, $\sigma_{8}=0.8$ and $h=0.7$.

The simulation’s density contrast field is provided in an ‘onion’
configuration (Fosalba et al., 2008), where particles from the MICE lightcone
have been binned onto 400 HEALPix maps (Górski et al., 2005). The redshift
slices are thinnest at low redshift and the thicknesses of the slices are well
below the smallest scales of interest ($k_{\max}=0.1\,h{\rm Mpc^{-1}}$). The
simulation is large enough that no box repetitions are required along the LOS
up to redshift $z=1.4$. For this reason the analysis in this paper is limited
to $z\leq 1.4$. Corresponding galaxy mock catalogues can be obtained from the
online database CosmoHub111https://cosmohub.pic.es/home (Carretero et al.,
2017; Tallada et al., 2020).

###  2.2 Theory: The Integrated Sachs-Wolfe Effect

We work in the Newtonian gauge with no anisotropic stress, so that metric
perturbations can be parameterised by a single perturbation variable $\Phi$,
the gravitational potential. The ISW is the imprint on the CMB of the
evolution of $\Phi$ in large scale structure (LSS). The effect alters the CMB
temperature $T$ in the LOS direction $\boldsymbol{\hat{n}}$ by

$\frac{\Delta
T_{\textnormal{ISW}(\boldsymbol{\hat{n}})}}{T}=\frac{2}{c^{2}}\int_{\eta_{\textnormal{LS}}}^{\eta_{0}}\Phi^{\prime}(r\boldsymbol{\hat{n}},\eta)\
d\eta,$ (1)

where the intergral is over a photon path from LS (= last scattering time) to
$0$ (= the observer’s time), $r=c(\eta_{0}-\eta)$ is the comoving distance on
this path corresponding to conformal time $\eta$, $c$ is the speed of light,
and $\Phi^{\prime}$ is the partial derivative of $\Phi$ with respect to $\eta$
(holding comoving position fixed). Changing the integration variable to
comoving distance $r$ yields

$\frac{\Delta
T_{\textnormal{ISW}(\boldsymbol{\hat{n}})}}{T}=\frac{2}{c^{3}}\int_{0}^{r_{\textnormal{LS}}}\dot{\Phi}(r\boldsymbol{\hat{n}},\eta(r))\
a(r)\ dr.$ (2)

Here $\eta(r)=\eta_{0}-r/c$ and we have additionally changed the integrand to
be the partial derivative $\dot{\Phi}$ of $\Phi$ with respect to time (holding
comoving position fixed); this introduces the scale factor $a=1/(1+z)$ where
$z=z(r)$ is the redshift. The relationship between redshift and comoving
distance can be approximated (during and after matter domination in
$\Lambda$CDM) by

$r(z)\simeq 3000\,h^{-1}{\rm Mpc}\int_{0}^{z}\frac{1}{\Omega_{\rm
m,0}(1+z^{\prime})^{3}+\Omega_{\rm\Lambda,0}}dz^{\prime},$ (3)

where $\Omega_{\rm m,0}$ is the current matter density and
$\Omega_{\Lambda,0}$ the current dark energy density (relative to the critical
density).

The non-linear contributions to the ISW (the Rees-Sciama effect) have been
shown to be subdominant ($\sim 10^{-2}\times$ the linear ISW in temperature)
in $\Lambda$CDM (Seljak, 1996; Cai et al., 2010; Nadathur et al., 2014).
Therefore the ISW can be approximated to within $\sim 1\%$ using linear
perturbation theory, which we discuss below.

####  2.2.1 Linear Theory Approximation

We seek to calculate Eq. 2 in the linear regime. We start with the Poisson
equation:

$\nabla^{2}\Phi(\boldsymbol{x},t)=\frac{3}{2}H_{0}^{2}\Omega_{\rm
m,0}\frac{\delta(\boldsymbol{x},t)}{a}.$ (4)

In what follows we work in the linear regime. Here, density perturbations can
be separated:

$\delta(\boldsymbol{x},t)=D(t)\delta(\boldsymbol{x}),$ (5)

where the linear growth factor $D(z)$ is defined to be

$D(z)\propto
H(z)\int_{0}^{a}\frac{1}{(a^{\prime}H(a^{\prime}))^{3}}da^{\prime}$ (6)

with $D(0)=1$ and $H(z)$ the Hubble expansion rate at redshift $z$. The latter
may be approximated during and after matter domination as

$H(z)\simeq H_{0}\left(\Omega_{\rm
m,0}(1+z)^{3}+\Omega_{\rm\Lambda,0}\right)^{1/2},$ (7)

where $H_{0}$, the Hubble constant, is the present expansion rate. Since we
assume flat curvature, the energy density for dark energy is given by
$\Omega_{\rm\Lambda}=1-\Omega_{\rm m}$ at any epoch.

Let $\nabla^{-2}$ denote the formal inverse of the Laplace operator. Combining
Eqs. 4 and 5 yields

$\Phi(\boldsymbol{x},t)=\frac{3}{2}H_{0}^{2}\Omega_{\rm
m,0}\frac{D(t)}{a}\nabla^{-2}\delta(\boldsymbol{x}),$ (8)

whence

$\dot{\Phi}(\boldsymbol{x},t)=\frac{3}{2}H_{0}^{2}\Omega_{\rm m,0}\
\frac{\partial}{\partial t}\left(\frac{D(t)}{a}\right)\
\nabla^{-2}\delta(\boldsymbol{x}).$ (9)

But

$\frac{\partial}{\partial
t}\left(\frac{D(t)}{a}\right)=\frac{D(t)}{a}H(t)\left[f(t)-1\right];$ (10)

here

$\displaystyle f$ $\displaystyle\equiv{\rm d}\ln D/{\rm d}\ln a$
$\displaystyle\approx\Omega_{\rm m}(z)^{0.55}$ (11)

(Peebles, 1980; Lahav et al., 1991), and

$\Omega_{\rm m}(z)=\frac{\Omega_{\rm m,0}(1+z)^{3}}{\Omega_{\rm
m,0}(1+z)^{3}+\Omega_{\rm\Lambda,0}}.$ (12)

Combining Eqs. 8, 9, and 10 yields

$\dot{\Phi}(\boldsymbol{x},t)=H(t)\,\left[f(t)-1\right]\,\Phi(\boldsymbol{x},t),$
(13)

which coupled with Eq. 2 gives the linear theory approximation for the ISW:

$\frac{\Delta T_{\rm
ISW}(\boldsymbol{\hat{n}})}{T}=\frac{2}{c^{3}}\int^{r_{\rm
LS}}_{0}H(t)\,[f(t)-1]\,\Phi(r\boldsymbol{\hat{n}},\eta(r))\,a(r)\,dr.$ (14)

Note when $f\simeq 1$, a condition which is true in an Einstein-de Sitter
Universe (where $\Omega_{\rm m}=1$ and $\Omega_{\Lambda}=0$) and during matter
domination in $\Lambda$CDM, $\Delta T_{\rm ISW}\simeq 0$. In this scenario the
ISW is dominated by the non-linear Rees-Sciama effect. For this reason the
most significant contributions to the ISW are at low redshift during $\Lambda$
domination.

####  2.2.2 Theoretical Angular Power Spectra

The angular power spectra $C_{\ell}^{XY}$ for sources $X$ and $Y$ is

$C_{\ell}^{XY}=\frac{2}{\pi}\int
k^{2}\,P(k)\,I_{\ell}^{X}(k)\,I_{\ell}^{Y}(k)\,dk,$ (15)

where

$I_{\ell}^{X}(k)=\int_{0}^{\infty}D(z)\,W^{X}(r,k)\,j_{\ell}(kr)\,dr,$ (16)

$j_{\ell}$ is the spherical Bessel function and $W^{X}(r,k)$ is a source-
specific window function. For the ISW the window function is given by

$W^{I}(r,k)=\frac{3\Omega_{m,0}H_{0}^{2}}{c^{3}k^{2}}\,H(z)\,\left[1-f(z)\right]\,{\rm
Rect}(r;r_{\min},r_{\max}),$ (17)

where ‘$I$’ is used as a shorthand for ISW and $\rm Rect$ is the rectangular
step function defined as

${\rm Rect}(r;r_{\min},r_{\max})={\rm H_{step}}(r-r_{\min})-{\rm
H_{step}}(r-r_{\max}),$ (18)

$\rm H_{step}$ is the heavyside step function, $r_{\min}=r(z_{\min})$ and
$r_{\max}=r(z_{\max})$. This is used to specify the redshift interval
$z_{\min}\leq z\leq z_{\max}$ for the contribution to the ISW. The galaxy
window function is given by

$W^{G}(r)=b(r)\,\Theta(r),$ (19)

where $b(r)$ is a function describing the galaxy bias,

$\Theta(r)=\frac{r^{2}\,n(r)}{\int x^{2}\,n(x)\,dx},$ (20)

and $n(r)$ is the galaxy redshift distribution. However, in this study we have
direct access to the true underlying density contrast field meaning (unlike
data from a galaxy survey) the redshift distribution and bias modelling can be
simply set to $n(r)=1$ and $b(r)=1$ in the interval $z_{\min}\leq z\leq
z_{\max}$ and $0$ otherwise. The window function reduces to

$W^{\delta}(r)=\begin{dcases}\frac{3r^{2}}{r^{3}_{\max}-r^{3}_{\min}},&\text{for
}z_{\min}\leq z\leq z_{\max},\\\ 0,&\text{otherwise,}\end{dcases}$ (21)

where ‘$\delta$’ is used as a shorthand for the density contrast field in the
interval $z_{\min}\leq z\leq z_{\max}$. It is common to calculate this using
the Limber approximation (Limber, 1954; Afshordi et al., 2004), replacing the
spherical Bessel function with a Dirac delta function $\delta_{\rm D}$ of the
form

$j_{\ell}(x)\simeq\sqrt{\frac{\pi}{2\ell+1}}\,\delta_{\rm
D}\left(\ell+\frac{1}{2}-x\right),$ (22)

giving the following approximation for the angular power spectra,

$C_{\ell}^{XY}\simeq\int\frac{D^{2}(z)}{r^{2}}\left[P(k)\,W^{X}(r,k)\,W^{Y}(r,k)\right]_{k=k_{\ell}}dr,$
(23)

where $k_{\ell}=(\ell+1/2)/r$.

Figure 1: The ISW effect for the MICE simulation constructed for contributions
in the range $0\leq z\leq 1.4$ using a SBT with normal boundary conditions.
The map is shown with a blue-to-red diverging colormap indicating cold and hot
features given in units $\mu{\rm K}$. The features of the map are quite broad
in comparison to CMB maps which typically have features on the scale of $\sim
1\,{\rm deg}$ and are smaller in amplitude typically of the order of $\sim
10\,\mu{\rm K}$ in the ISW compared to $\sim 100\,\mu{\rm K}$ for the CMB.

###  2.3 Techniques for Constructing Integrated Sachs-Wolfe Maps

In this section we describe the numerical techniques used to construct ISW
maps. The focus of this study is the spherical Bessel transform (SBT)
technique (Shapiro et al., 2012), described in Section 2.3.1. We also compare
this to a Spherical Harmonic scaling (which we call SHS) relation, described
in Section 2.3.2, and the Francis & Peacock (2010) technique, described in
Section 2.3.3.

####  2.3.1 Spherical Bessel Transform Method

A three dimensional field $\zeta$ described in spherical polar coordinates $r$
(radial axis), $\theta$ and $\phi$ (latitude and longitude respectively) can
be represented by its SBT coefficients $\zeta_{\ell mn}$ as

$\zeta(r,\theta,\phi)=\sum^{\infty}_{\ell=0}\sum^{\ell}_{m=-\ell}\sum_{n=1}^{\infty}\zeta_{\ell
mn}\,R_{\ell n}(r)\,Y_{\ell m}(\theta,\phi),$ (24)

where

$R_{\ell n}(r)=\frac{1}{\sqrt{N_{\ell n}}}j_{\ell}(k_{\ell n}r),$ (25)

$Y_{\ell m}$ are spherical harmonics, $k_{\ell n}=q_{\ell n}/r_{\max}$ and
$q_{\ell n}$ is the locations of the $n^{\rm th}$ zero of $j_{\ell}(x)$ for
normal boundary conditions (i.e. $j_{\ell}(k_{\ell n}r_{\max})=0$) or
$\partial_{x}j_{\ell}(x)$ for derivative boundary conditions (i.e.
$\partial_{x}j_{\ell}(k_{\ell n}r_{\max})=0$). Lastly, $N_{\ell n}$ is a
normalisation constant defined as

$N_{\ell n}=\begin{dcases}\frac{r_{\max}^{3}}{2}j_{\ell+1}^{2}(k_{\ell
n}r_{\max}),&\text{normal boundary,}\\\
\frac{r_{\max}^{3}}{2}\left(1-\frac{\ell(\ell+1)}{k_{\ell
n}^{2}r_{\max}^{2}}\right)j_{\ell}^{2}(k_{\ell n}r_{\max}),&\text{derivative
boundary}.\end{dcases}$ (26)

See Wang et al. (2009) for an overview of SBT and related transforms. The SBT
basis functions are eigenfunctions of the Laplace operator on a ball-shaped
domain. This will be crucial, as it is the Laplace operator in the Poisson
equation that links the gravitational potential (which drives the ISW) to the
overdensities (which are the observable quantities within the simulation).

The SBT coefficients are calculated from

$\zeta_{\ell
mn}=\int^{r_{\max}}_{0}\int^{\pi}_{0}\int^{2\pi}_{0}\zeta(r,\theta,\phi)\,R_{\ell
n}(r)\,Y^{*}_{\ell m}(\theta,\phi)\,r^{2}\,\sin(\theta)\,d\phi\,d\theta\,dr,$
(27)

where we assume the field is defined up to a maximum radius $r_{\max}$. We
will be concerned with $\delta_{\ell mn}$, the SBT coefficients of the density
contrast field $\delta$ $=\delta(\boldsymbol{x})$ defined on three-dimensional
space at $z=0$, and with $\Phi_{\ell mn}$, the SBT coefficients of the
gravitational potential defined on the observer’s lightcone.

Following Leistedt et al. (2012) we calculate $\delta_{\ell mn}$ in two
stages. The first step derives the spherical harmonics of each slice
$\delta^{i}$:

$\delta_{\ell
m}^{i}=\int_{0}^{\pi}\int_{0}^{2\pi}\frac{\delta^{i}(\theta,\phi)}{D(z_{\rm
eff})}\,Y_{\ell m}^{*}(\theta,\phi)\,\sin(\theta)\,d\phi\,d\theta.$ (28)

Here we have corrected for the linear evolution of $\delta$ by dividing by the
linear growth function at $z_{\rm eff}$ (the effective redshift of the slice).
This ensures that the SBT coefficients of $\delta$ are taken at $z=0$. The
second step calculates the SBT coefficients from

$\delta_{\ell mn}=\sum_{i}\delta_{\ell
m}^{i}\int^{r^{i}_{\max}}_{r^{i}_{\min}}\,R_{\ell n}(r)\,r^{2}\,dr,$ (29)

where $r^{i}_{\min}$ and $r^{i}_{\max}$ are the minimum and maximum radius for
the $i^{\rm th}$ slice.

Let $a^{I}_{\ell m}$ be the spherical harmonic coefficients for the ISW:

$\frac{\Delta
T_{\textnormal{ISW}}}{T}=\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell}a^{I}_{\ell
m}Y_{\ell m}(\theta,\phi).$ (30)

Equating spherical harmonic coefficients in Eqs. 14 and 30 gives

$a^{I}_{\ell m}=\frac{2}{c^{3}}\int_{0}^{r_{\rm
LS}}H(r)[f(r)-1]\sum_{n=1}^{\infty}\Phi_{\ell mn}R_{\ell n}(r)a(r)dr.$ (31)

The SBT basis function $R_{\ell n}Y_{\ell m}$ is an eigenfunction of the
Laplace operator with eigenvalue $-{k_{\ell n}}^{2}$. Combining this with a
SBT representation of the Poisson Eqn. 4 and applying Eqn. 5 yields

$\Phi_{\ell
mn}=-\frac{3}{2}H_{0}^{2}\Omega_{m,0}\frac{D(t)}{a}\frac{1}{{k_{\ell
n}}^{2}}\delta_{\ell mn};$ (32)

combining the last two results gives our target expression for the ISW (due to
Shapiro et al. (2012)):

$a_{\ell
m}^{I}=\frac{3H_{0}^{2}\Omega_{m,0}}{c^{3}}\sum_{n=1}^{n_{\max}}\frac{\delta_{\ell
mn}}{{k_{\ell n}}^{2}}\int_{0}^{r_{\max}}D(r)H(r)\left[1-f(r)\right]R_{\ell
n}(r)dr.$ (33)

In this study we compute only the SBT coefficients that correspond to Fourier
modes in the range $k_{\rm F}\leq k\leq k_{\max}$ where $k_{\max}=0.1\,h{\rm
Mpc^{-1}}$ and $k_{\rm F}$ is the fundamental frequency of the simulation
($k_{\rm F}=2\pi/L_{\rm box}$ where $L_{\rm box}$ is the length of the
simulation box). This means in our analysis the SBT coefficients are computed
for $\ell=2$ to $\ell_{\max}$ (ignoring the monopole and dipole components)
and $n=1$ to $n_{\max}$, where $\ell_{\max}=\lfloor r_{\max}k_{\max}\rfloor$
and $n_{\max}=\lfloor\ell_{\max}/\pi\rfloor$ (the floor function ${\rm
Floor}(x)=\lfloor x\rfloor$ ensures $l_{\max}$ and $n_{\max}$ are integers).

Figure 2: The ISW constructed using the SBT with normal boundary conditions
shown for contributions in redshift slices of $\Delta z=0.2$ in the range
$0\leq z\leq 1.2$. The maps are divided by their standard deviation to
highlight the scale of the features in these different redshift slices. This
change of resolution originates from the size of structures and the scale of
homogeneity, at closer redshift the features are broad due to the size of
nearby cosmic structures but at high redshift features are smaller as these
cosmic structures are further away.

####  2.3.2 Spherical Harmonic Scaling

A common alternative approach to constructing ISW maps from the density field
is to use the SHS relation (Manzotti & Dodelson, 2014; Muir & Huterer, 2016),

$a^{I}_{\ell
m}=\frac{C_{\ell}^{I\delta}}{C_{\ell}^{\delta\delta}}a^{\delta}_{\ell m},$
(34)

where $C_{\ell}^{\delta\delta}$ is the auto-angular power spectrum for the
density contrast, $C_{\ell}^{I\delta}$ the cross-angular power spectrum for
the ISW-density contrast and $a^{\delta}_{\ell m}$ the spherical harmonic
coefficients for the density field integrated over the redshift region of
interest.

Figure 3: The ISW maps are shown for four methods for contributions from the
redshift range $0\leq z\leq 1.4$. In the top left panel is the SBT (normal)
map, in the top right panel is the SBT (derivative) map, in the bottom left
panel is the SHS map and in the bottom right panel is the FP map. The maps are
in fairly good agreement but the approximations on the bottom panels show
variations in the amplitude of certain features while the SBT methods are
almost identical. Figure 4: Differences between the SBT (normal) map and maps
generated by the three other methods (normalised by the standard deviation of
the SBT (normal) map) measured in the redshift range $0\leq z\leq 1.4$. From
left to right are the difference for the SBT (derivative) map, for the SHS
map, and for the FP map (all with respect to the SBT (normal) map). This shows
only a small difference between the SBT methods (which differ only in their
boundary condition) and shows a significant difference between the SHS and FP
maps (which unlike the SBT methods are computed without any LOS information).

####  2.3.3 Francis & Peacock Approximation

Francis & Peacock (2010) derive an approximation for the ISW (now referred to
as the FP approximation),

$a_{\ell m}^{I}\simeq\frac{3\Omega_{\rm
m,0}H_{0}^{2}}{\ell(\ell+1)c^{3}}\,H(z_{\rm eff})\left[1-f(z_{\rm
eff})\right]\,\Delta r\,r_{\rm eff}^{2}\,a_{\ell m}^{\delta},$ (35)

where $r_{\rm eff}$ is the effective comoving radius of the density contrast
field,

$r_{\rm
eff}=\frac{3}{4}\left(\frac{r_{\max}^{4}-r_{\min}^{4}}{r_{\max}^{3}-r_{\min}^{3}}\right),$
(36)

$z_{\rm eff}$ is the redshift corresponding to $r_{\rm eff}$, $\Delta
r=r_{\max}-r_{\min}$, $r_{\min}$ is the minimum comoving radius and $r_{\max}$
is the maximum comoving radius. This approximation is best suited for cases
where $\Delta r$ is thin; for this reason with this approximation we use thin
slices and then combine the maps.

####  2.3.4 Software Pipeline

The Python (Van Rossum & Drake, 2009) package pyGenISW is made publicly
available222https://github.com/knaidoo29/pyGenISW and can be used to construct
ISW maps from data provided in spherical shells (given in HEALPix format)
using the SBT method (which is the focus of this paper) as well as the
alternative SHS approach and the FP approximation.

The package depends on TheoryCL333https://github.com/knaidoo29/TheoryCL which
computes the linear growth functions and angular power spectra for the ISW and
density contrast sources; camb444https://camb.readthedocs.io/ (Lewis et al.,
2000) to compute the linear power spectrum,
healpy555https://healpy.readthedocs.io/ (Górski et al., 2005; Zonca et al.,
2019) for computing and manipulating maps and carrying out spherical harmonic
operations; SciPy used for spherical Bessel related functions, integration
functions and interpolation functions and NumPy666https://numpy.org/ (Harris
et al., 2020). Note, pyGenISW uses a mixture of SciPy777https://www.scipy.org/
(Virtanen et al., 2020) and its own iterative spherical Bessel root finding
function to determine $q_{\ell n}$; this is because the current implementation
in SciPy is unstable for large $\ell$ and $n$.

Figure 5: The auto- and cross-angular power spectra for the ISW and density
field for contributions in the redshift range $0\leq z\leq 1.4$. In the top
two panels we show the ISW angular power spectra $C_{\ell}^{II}$ for the SBT
methods (left) and the SHS and FP approximations (middle). In the bottom
panels we show the cross-angular power spectra for the ISW and density field
for the SBT method (left) and the SHS and FP approximations (middle) and show
the auto-angular power spectra for the MICE density field (right). In the
subplots we show the significance of deviations from the full theoretical
$C_{\ell}$ (i.e. with no cuts in $k$-range shown with full black lines) where
$\sigma(C_{\ell})=(C_{\ell}^{\rm Measured}-C_{\ell}^{\rm Theory})/\Delta
C_{\ell}$. The grey bands indicate $1\sigma$ confidence regions. We see that
$C_{\ell}^{II}$ and $C_{\ell}^{I\delta}$ follow the trend of $C_{\ell}(k_{\rm
F}\leq k\leq k_{\max})$ (shown with a dashed black line), i.e. beginning to
fall in amplitude at high-$\ell$ with respect to the full $C_{\ell}$ mainly
due to non-linear contributions at low redshift and drop at low-$\ell$ due to
the lack of modes larger than the simulation box. The SHS approximation is
calculated using the very $C_{\ell}$ that we compare to, so this comparison is
rather redundant as we expect a reasonably good fit by construction. However
in the case of the SBT and FP approximation no $C_{\ell}$ are given and this
comparison represents a true test of the accuracy of these methods. The SBT
ISW maps show the closest match to the theoretical $C_{\ell}$, while the
approximations appear to be slightly too high in amplitude for the FP
approximation and too low for the SHS approximation. Furthermore, we see the
SBT methods deviate only at very low $\ell$.

##  3 Results: MICE Integrated Sachs-Wolfe Maps

###  3.1 Using Spherical Bessel Transforms

The ISW maps for MICE are constructed using the SBT with both normal and
derivative boundary conditions. We compute the SBT coefficients up to an
$r_{\max}$ that corresponds to $z=2$, exceeding the maximum redshift of
$z=1.4$ for which the ISW is computed. This is to ensure that no artifacts are
measured near the boundary. The choice of $r_{\max}=r(z=2)$ and
$k_{\max}=0.1\,h{\rm Mpc^{-1}}$ means we only need to compute the SBT up to
$\ell_{\max}=380$ and $n_{\max}=121$. The full-sky MICE ISW map for
contributions in the range $0\leq z\leq 1.4$ is shown in Fig. 1 using the SBT
with normal boundary conditions. The ISW map constructed using the SBT with
derivative boundary conditions is not highlighted here, as the map is almost
identical to the one shown in Fig. 1; it will be discussed in Section 3.2.

To highlight the scales of the features that contribute at different
redshifts, we plot in Fig. 2 the ISW constructed using the SBT with normal
boundary conditions for redshift slices of $\Delta z=0.2$ between $z=0$ and
$z=1.2$. Features in the map at low redshift are large in angular scale and
become smaller at higher redshift.

###  3.2 Comparison to SHS and Francis & Peacock Approximation

The ISW maps for MICE have been constructed using four methods: the first two
use the SBT with normal and derivative boundary conditions, while the latter
two use the SHS and the FP approximation. The latter two assume no cuts in the
Fourier modes used, while for the SBT the modes have been explicitly limited
to $k_{\rm F}\leq k\leq k_{\max}$. To ensure the SHS and the FP approximations
are comparable to those of the SBT, and do not amplify spurious Fourier modes,
we multiply the spherical harmonic coefficients by $C^{II}_{\ell}(k_{\rm
F}\leq k\leq k_{\max})/C_{\ell}^{II}$. Furthermore in the case of the FP
approximation, since this method works best for thin shells we compute the ISW
contributions in slices of $\Delta z=0.2$ which are then combined for
contributions in the range $0\leq z\leq 1.4$.

In Fig. 3 the ISW maps for the four methods are shown for contributions in the
range $0\leq z\leq 1.4$. In Fig. 4 the difference with respect to the SBT with
normal boundary conditions is shown. These Figures show strong agreement
between SBT with normal and derivative boundary conditions, but find the
difference between the SHS approximation and the FP approximation to be fairly
significant (deviations of the order of $\sim 2\sigma$) driven by an absence
of LOS density information for the latter two methods.

Figure 6: The standard deviation of the ISW maps for contributions in the
range $0\leq z\leq 1.4$ with $\Delta z=0.2$. The four methods are compared to
Gaussian realisations, where the $1\sigma$ confidence region is shown in grey.
The theoretical standard deviation calculated from Eq. 37 is shown with grey
dotted lines. In the top panel we consider the full ISW maps with no
alterations and in the bottom panel we limit the maps to spherical harmonic
coefficients with $\ell>10$. The Figure shows that the FP approximation is
consistently too large and the SHS slightly too low when compared to Gaussian
realisations. The SBT methods are the most consistent with Gaussian
realisation with the exception of some ranges in $z$ where they appear
slightly lower for the full ISW map.

###  3.3 Angular Power Spectra: Comparisons to Theory

All the methods in this study express the ISW in terms of its spherical
harmonics. For this reason the auto- and cross-angular power spectra of the
ISW and density field maps are calculated using healpy’s alm2cl function. In
Fig. 5 we plot the auto- and cross-angular power spectra for the four ISW maps
and for the density field. This is compared to theoretical angular power
spectra with and without cuts in the Fourier modes considered. For the full
$C_{\ell}$ the integration in $k$ space is carried out for a $k$-range of
$10^{-4}\leq k\leq 1\,h{\rm Mpc^{-1}}$, while the cut-$C_{\ell}$ uses a
$k$-range of $k_{\rm F}\leq k\leq 0.1\,h{\rm Mpc^{-1}}$. The cut-$C_{\ell}$
are used to determine the ranges in $\ell$ for which the maps are valid.
Furthermore, they indicate on which $\ell$-ranges the non-linear scales become
relevant (where ‘non-linear’ refers to small scales with $k>0.1\,h{\rm
Mpc^{-1}}$). For the theoretical $C_{\ell}$ we calculate it fully without
using the Limber approximation. This is because comparisons of the ISW for
thin shells with $\Delta z=0.2$ showed an offset in the amplitude at all
$\ell$ when compared to the measured $C_{\ell}$. The $C_{\ell}$ are binned
into $\Delta\ell=5$ and are shown in Fig. 5 to agree very well (for all
methods) to the theory. The $C_{\ell}$ also reveal that the SBT ISW maps are
in closest agreement to theoretical expectations, the SHS map appears to be
low in power while FP appears to be slightly higher in power than would be
expected (however these deviations are consistent with cosmic variance $\Delta
C_{\ell}=\sqrt{2C_{\ell}^{2}/(2\ell+1)}$).

###  3.4 Dependence on Redshift and Comparison to Gaussian Realisations

In this section we measure the standard deviation of the ISW maps for
contributions in the range $0\leq z\leq 1.4$ with $\Delta z=0.2$. These are
compared to Gaussian realisations computed using healpy’s synalm function from
the theoretical $C_{\ell}(k_{\rm F}\leq k\leq k_{\max})$. The results are
shown in Fig. 6 where in the top panel we consider the standard deviation of
the ISW maps with no alteration and in the bottom panel we limit to only
spherical harmonic coefficients $\ell>10$. The theoretical standard deviations
(shown with dotted grey lines) are calculated from

$\sigma=\left(\sum_{\ell}\frac{2\ell+1}{4\pi}C_{\ell}\right)^{1/2},$ (37)

(Tegmark, 1997). The comparisons show that the standard deviation in the ISW
for the FP approximation is consistently too high and for the SHS is slightly
low with respect to Gaussian realisations. In comparison the SBT methods are
consistent with Gaussian realisations for $\ell>10$ but are low for certain
redshifts when considering the full map with no alterations. The plot also
indicates where the biggest contributions to the ISW occur, peaking around
$z\simeq 0.3$ for all scales and at $z\simeq 0.5$ for features with $\ell>10$.
This demonstrates the benefit of using the SBT methods over the SHS or FP
approximation as they are better able to reproduce the theoretically expected
standard deviation for the ISW temperature maps.

###  3.5 Maps and Ancillary Data

The ISW maps constructed in this study are
public.888https://doi.org/10.5281/zenodo.4088697 They are constructed for the
full redshift range $0\leq z\leq 1.4$ and for spherical shells with redshift
width $\Delta z=0.2$. The maps themselves are provided in fits format as
HEALPix maps with nside=256. The spherical harmonics are also provided so that
the maps can be generated to the desired nside. Theoretical auto- and cross-
angular power spectra are provided for the ISW and density field with and
without using the Limber approximation.

##  4 Discussion

In this paper we construct the ISW map for the MICE lightcone simulation and
develop a pipeline for quickly constructing the ISW for future simulations;
the SBT transform (the bottleneck of this computation) was calculated in $\sim
10$ minutes on a single core on commercially available hardware for
$\ell_{\max}=380$ and $n_{\max}=121$. A significant computational limitation
of constructing simulated ISW maps is the calculation of the time derivative
of the gravitational potential $\dot{\Phi}$. Such a calculation requires the
regular output of snapshots, which is memory intensive, and cannot be
performed ad-hoc. To remove this obstacle we calculate the ISW maps fully in
the linear regime, sacrificing the construction of the non-linear ISW or Rees-
Sciama effect. Since the Rees-Sciama effect is comparatively very small ($\sim
1\%$ of the ISW; Seljak, 1996; Nadathur et al., 2014) this sacrifice has a
neglegible impact on the maps created. This approximation enables the ISW to
be computed from the density field at $z=0$, where we insert the time
evolution of the gravitational potential analytically using linear growth
functions. For lightcone simulations such as MICE we divide the density field
at redshift $z$ by the linear growth $D(z)$ to approximate the density field
at $z=0$.

The Python module developed in this study, pyGenISW, is publicly
available999https://github.com/knaidoo29/pyGenISW and can be used to construct
the ISW using the SBT method (the focus of this paper), as well as using the
SHS and the FP approximation. The ISW maps and ancillary data (spherical
harmonic coefficients and theoretical angular power spectra) are made
public101010https://doi.org/10.5281/zenodo.4088697 to facilitate future LSS-
CMB cross-correlation studies using MICE.

In Fig. 1 we highlight one of the main data products of this study – the ISW
map for contributions in the redshift range $0\leq z\leq 1.4$ for the SBT
method with normal boundary conditions. In Fig. 2 the maps are constructed for
redshift intervals with $\Delta z=0.2$. These maps show the redshift evolution
of the scales of features and will enable isolated studies of the ISW at
different redshift. We compare the SBT methods, the SHS and the FP
approximation in Fig. 3 and subtract the SBT (normal) map from the other three
in Fig. 4 demonstrating that the two SBT methods are virtually identical
(since they differ only in their boundary conditions) and that there are
significant large angular scale differences between the SBT methods and the
SHS and FP approximation (driven by the absense of LOS density information).
In Fig. 5 the auto- and cross-angular power spectra for the ISW and density
field are shown. In these plots we see that the SBT method is in greater
agreement with theoretical expectations while SHS appears to be slightly too
low and the FP slightly too large at all scales, demonstrating the benefit of
using the SBT methods over the SHS or FP approximation. This effect is
illustrated again in Fig. 6 where we compare the standard deviation of the ISW
maps to those from Gaussian realisations. Once again the SBT methods are in
greater agreement with theoretical expectations showing that these particular
maps are the best maps to use for future MICE LSS-CMB cross-correlation
studies.

These maps will be of particular relevance to large-area galaxy surveys, such
as Dark Energy Survey (DES)111111http://www.darkenergysurvey.org, Dark Energy
Spectroscopic Instrument (DESI)121212http://desi.lbl.gov/, _Euclid_
131313http://www.euclid-ec.org/, and the Rubin Legacy Survey of Space and Time
(LSST)141414https://www.lsst.org/. Cross-correlation studies have been
performed on several galaxy surveys to constrain the standard cosmological
model $\Lambda$CDM and to test extensions/alternative models (see e.g. Hang et
al., 2020). The availability of the MICE ISW maps will enable future studies
to test their cross-correlation pipelines and to test the results against the
predictions of the fiducial $\Lambda$CDM model. Furthermore, ISW
reconstruction methods (Barreiro et al., 2008; Granett et al., 2009; Manzotti
& Dodelson, 2014; Muir & Huterer, 2016) will be able to use these maps as a
ground truth.

ISW void-stacking measurements are in slight tension with predictions from
$\Lambda$CDM (a tension of $\gtrsim 2\sigma$ Granett et al., 2008; Pápai et
al., 2011; Nadathur et al., 2012; Flender et al., 2013; Hernández-Monteagudo &
Smith, 2013; Ilić et al., 2013; Granett et al., 2015; Kovács & García-Bellido,
2016; Nadathur & Crittenden, 2016; Cai et al., 2017; Kovács, 2018; Kovács et
al., 2019). An explanation for the source of this tension remains unclear.
Recently, Kovács et al. (2020) showed that the signal could be explained if
modifications to the growth history were applied but Hang et al. (2020) show
this solution is incompatible with cross-correlation studies. To determine
whether variations in the value of cosmological parameters or models could
explain this excess signal we will need to construct the ISW map for a large
range of simulations with realistic galaxy catalogues to construct observable-
like void catalogues. The pipeline developed in this study enables future work
to construct the ISW from existing and future simulation suites, allowing us
to understand the parameter and model dependence of the ISW void-stacking
measurement. Accurate ISW maps from simulations such as MICE will enable
cosmologists to fully exploit future galaxy surveys by combining probes from
the early universe in the form of the CMB and the late universe in form of
LSS. This will provide further tests of the standard cosmological model
$\Lambda$CDM and may be crucial in establishing the validity of extended and
alternative models. In future work we plan to provide MICE ISW maps for
contributions at higher redshift and to higher $\ell_{\max}$ (essentially
pushing our scale limits from the linear to the quasi-linear regime) as well
as generating ISW maps for existing simulations suites probing different
cosmological models and parameters.

## ACKNOWLEDGEMENTS

We thank Andras Kovacs, Peter Coles and Andrew Pontzen for providing useful
comments and discussions. KN acknowledges support from the Science and
Technology Facilities Council grant ST/N50449X and from the (Polish) National
Science Centre grant #2018/31/G/ST9/03388. PF acknowledges support from MINECO
through grant ESP2017-89838-C3-1-R, the H2020 European Union grants LACEGAL
734374 and EWC 776247 with ERDF funds, and Generalitat de Catalunya through
CERCA to grant 2017-SGR-885 and funding to IEEC. OL acknowledges support from
an STFC Consolidated Grant ST/R000476/1.

The analysis in this paper was constructed using the following Python (Van
Rossum & Drake, 2009) modules: NumPy (Harris et al., 2020), SciPy (Virtanen et
al., 2020), healpy (Górski et al., 2005; Zonca et al., 2019) and camb (Lewis
et al., 2000); figures were constructed using Matplotlib (Hunter, 2007) and
Cartopy (Met Office, 2015).

## Data Availability

The MICE density field maps used in this study can be provided upon reasonable
request to the main author. The ISW maps and ancillary data are public and can
be obtained from https://doi.org/10.5281/zenodo.4088697.

## References

  * Adamek et al. (2020) Adamek J., Rasera Y., Corasaniti P. S., Alimi J.-M., 2020, Phys. Rev. D, 101, 023512
  * Afshordi et al. (2004) Afshordi N., Loh Y.-S., Strauss M. A., 2004, Phys. Rev. D, 69, 083524
  * Barreiro et al. (2008) Barreiro R. B., Vielva P., Hernandez-Monteagudo C., Martinez-Gonzalez E., 2008, IEEE Journal of Selected Topics in Signal Processing, 2, 747
  * Barreiro et al. (2013) Barreiro R. B., Vielva P., Marcos-Caballero A., Martínez-González E., 2013, MNRAS, 430, 259
  * Boughn & Crittenden (2002) Boughn S. P., Crittenden R. G., 2002, Phys. Rev. Lett., 88, 021302
  * Cabass et al. (2015) Cabass G., Gerbino M., Giusarma E., Melchiorri A., Pagano L., Salvati L., 2015, Phys. Rev. D, 92, 063534
  * Cai et al. (2009) Cai Y.-C., Cole S., Jenkins A., Frenk C., 2009, MNRAS, 396, 772
  * Cai et al. (2010) Cai Y.-C., Cole S., Jenkins A., Frenk C. S., 2010, MNRAS, 407, 201
  * Cai et al. (2017) Cai Y.-C., Neyrinck M., Mao Q., Peacock J. A., Szapudi I., Berlind A. A., 2017, MNRAS, 466, 3364
  * Carbone et al. (2016) Carbone C., Petkova M., Dolag K., 2016, J. Cosmology Astropart. Phys., 2016, 034
  * Carretero et al. (2015) Carretero J., Castander F. J., Gaztañaga E., Crocce M., Fosalba P., 2015, MNRAS, 447, 646
  * Carretero et al. (2017) Carretero J., et al., 2017, in Proceedings of the European Physical Society Conference on High Energy Physics. 5-12 July. p. 488
  * Corasaniti et al. (2005) Corasaniti P.-S., Giannantonio T., Melchiorri A., 2005, Phys. Rev. D, 71, 123521
  * Crittenden & Turok (1996) Crittenden R. G., Turok N., 1996, Phys. Rev. Lett., 76, 575
  * Crocce et al. (2015) Crocce M., Castander F. J., Gaztañaga E., Fosalba P., Carretero J., 2015, MNRAS, 453, 1513
  * Dunkley et al. (2009) Dunkley J., et al., 2009, ApJS, 180, 306
  * Dupé et al. (2011) Dupé F. X., Rassat A., Starck J. L., Fadili M. J., 2011, A&A, 534, A51
  * Ferraro et al. (2015) Ferraro S., Sherwin B. D., Spergel D. N., 2015, Phys. Rev. D, 91, 083533
  * Flender et al. (2013) Flender S., Hotchkiss S., Nadathur S., 2013, J. Cosmology Astropart. Phys., 2013, 013
  * Fosalba et al. (2003) Fosalba P., Gaztañaga E., Castander F. J., 2003, ApJ, 597, L89
  * Fosalba et al. (2008) Fosalba P., Gaztañaga E., Castander F. J., Manera M., 2008, MNRAS, 391, 435
  * Fosalba et al. (2015a) Fosalba P., Gaztañaga E., Castander F. J., Crocce M., 2015a, MNRAS, 447, 1319
  * Fosalba et al. (2015b) Fosalba P., Crocce M., Gaztañaga E., Castander F. J., 2015b, MNRAS, 448, 2987
  * Francis & Peacock (2010) Francis C. L., Peacock J. A., 2010, MNRAS, 406, 2
  * Giannantonio et al. (2006) Giannantonio T., et al., 2006, Phys. Rev. D, 74, 063520
  * Giannantonio et al. (2008) Giannantonio T., Scranton R., Crittenden R. G., Nichol R. C., Boughn S. P., Myers A. D., Richards G. T., 2008, Phys. Rev. D, 77, 123520
  * Giannantonio et al. (2012) Giannantonio T., Crittenden R., Nichol R., Ross A. J., 2012, MNRAS, 426, 2581
  * Giannantonio et al. (2014) Giannantonio T., Ross A. J., Percival W. J., Crittenden R., Bacher D., Kilbinger M., Nichol R., Weller J., 2014, Phys. Rev. D, 89, 023511
  * Górski et al. (2005) Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759
  * Goto et al. (2012) Goto T., Szapudi I., Granett B. R., 2012, MNRAS, 422, L77
  * Granett et al. (2008) Granett B. R., Neyrinck M. C., Szapudi I., 2008, ApJ, 683, L99
  * Granett et al. (2009) Granett B. R., Neyrinck M. C., Szapudi I., 2009, ApJ, 701, 414
  * Granett et al. (2015) Granett B. R., Kovács A., Hawken A. J., 2015, MNRAS, 454, 2804
  * Hang et al. (2020) Hang Q., Alam S., Peacock J. A., Cai Y.-C., 2020, arXiv e-prints, p. arXiv:2010.00466
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hernández-Monteagudo & Smith (2013) Hernández-Monteagudo C., Smith R. E., 2013, MNRAS, 435, 1094
  * Ho et al. (2008) Ho S., Hirata C., Padmanabhan N., Seljak U., Bahcall N., 2008, Phys. Rev. D, 78, 043519
  * Hoffmann et al. (2015) Hoffmann K., Bel J., Gaztañaga E., Crocce M., Fosalba P., Castander F. J., 2015, MNRAS, 447, 1724
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Ilić et al. (2013) Ilić S., Langer M., Douspis M., 2013, A&A, 556, A51
  * Kovács (2018) Kovács A., 2018, MNRAS, 475, 1777
  * Kovács & García-Bellido (2016) Kovács A., García-Bellido J., 2016, MNRAS, 462, 1882
  * Kovács et al. (2019) Kovács A., et al., 2019, MNRAS, 484, 5267
  * Kovács et al. (2020) Kovács A., Beck R., Szapudi I., Csabai I., Rácz G., Dobos L., 2020, MNRAS, 499, 320
  * Lahav et al. (1991) Lahav O., Lilje P. B., Primack J. R., Rees M. J., 1991, Monthly Notices of the Royal Astronomical Society, 251, 128
  * Leistedt et al. (2012) Leistedt B., Rassat A., Réfrégier A., Starck J. L., 2012, A&A, 540, A60
  * Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473
  * Limber (1954) Limber D. N., 1954, ApJ, 119, 655
  * Mackenzie et al. (2017) Mackenzie R., Shanks T., Bremer M. N., Cai Y.-C., Gunawardhana M. L. P., Kovács A., Norberg P., Szapudi I., 2017, MNRAS, 470, 2328
  * Manzotti & Dodelson (2014) Manzotti A., Dodelson S., 2014, Phys. Rev. D, 90, 123009
  * McEwen et al. (2007) McEwen J. D., Vielva P., Hobson M. P., Martínez-González E., Lasenby A. N., 2007, MNRAS, 376, 1211
  * Met Office (2015) Met Office 2010 - 2015, Cartopy: a cartographic python library with a Matplotlib interface. Exeter, Devon, https://scitools.org.uk/cartopy
  * Muir & Huterer (2016) Muir J., Huterer D., 2016, Phys. Rev. D, 94, 043503
  * Nadathur & Crittenden (2016) Nadathur S., Crittenden R., 2016, ApJ, 830, L19
  * Nadathur et al. (2012) Nadathur S., Hotchkiss S., Sarkar S., 2012, J. Cosmology Astropart. Phys., 2012, 042
  * Nadathur et al. (2014) Nadathur S., Lavinto M., Hotchkiss S., Räsänen S., 2014, Phys. Rev. D, 90, 103510
  * Naidoo et al. (2016) Naidoo K., Benoit-Lévy A., Lahav O., 2016, MNRAS, 459, L71
  * Naidoo et al. (2017) Naidoo K., Benoit-Lévy A., Lahav O., 2017, MNRAS, 472, L65
  * Nicola et al. (2016) Nicola A., Refregier A., Amara A., 2016, Phys. Rev. D, 94, 083517
  * Nolta et al. (2004) Nolta M. R., et al., 2004, ApJ, 608, 10
  * Padmanabhan et al. (2005) Padmanabhan N., Hirata C. M., Seljak U., Schlegel D. J., Brinkmann J., Schneider D. P., 2005, Phys. Rev. D, 72, 043525
  * Pápai et al. (2011) Pápai P., Szapudi I., Granett B. R., 2011, ApJ, 732, 27
  * Peebles (1980) Peebles P. J. E., 1980, The large-scale structure of the universe. Princeton University Press Princeton, N.J
  * Peiris & Spergel (2000) Peiris H. V., Spergel D. N., 2000, ApJ, 540, 605
  * Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571, A19
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A21
  * Rees & Sciama (1968) Rees M. J., Sciama D. W., 1968, Nature, 217, 511
  * Sachs & Wolfe (1967) Sachs R. K., Wolfe A. M., 1967, ApJ, 147, 73
  * Scranton et al. (2003) Scranton R., et al., 2003, arXiv e-prints, pp astro–ph/0307335
  * Seljak (1996) Seljak U., 1996, ApJ, 460, 549
  * Shapiro et al. (2012) Shapiro C., Crittenden R. G., Percival W. J., 2012, MNRAS, 422, 2341
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Stölzner et al. (2018) Stölzner B., Cuoco A., Lesgourgues J., Bilicki M., 2018, Phys. Rev. D, 97, 063506
  * Szapudi et al. (2015) Szapudi I., et al., 2015, MNRAS, 450, 288
  * Tallada et al. (2020) Tallada P., et al., 2020, Astronomy and Computing, 32, 100391
  * Tegmark (1997) Tegmark M., 1997, Phys. Rev. D, 56, 4514
  * Van Rossum & Drake (2009) Van Rossum G., Drake F. L., 2009, Python 3 Reference Manual. CreateSpace, Scotts Valley, CA
  * Vielva et al. (2006) Vielva P., Martínez-González E., Tucci M., 2006, MNRAS, 365, 891
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Wang et al. (2009) Wang Q., Ronneberger O., Burkhardt H., 2009, IEEE Transactions on Pattern Analysis and Machine Intelligence, 31, 1715
  * Watson et al. (2014) Watson W. A., et al., 2014, MNRAS, 438, 412
  * Xia et al. (2009) Xia J.-Q., Viel M., Baccigalupi C., Matarrese S., 2009, J. Cosmology Astropart. Phys., 2009, 003
  * Zonca et al. (2019) Zonca A., Singer L., Lenz D., Reinecke M., Rosset C., Hivon E., Gorski K., 2019, Journal of Open Source Software, 4, 1298