Multiple Lensing of the Cosmic Microwave Background anisotropies
algorithm.tex
%!TEX root = jcap_article_MLcode.tex
% !TEX spellcheck = en-EN

\section{The Algorithm}
\label{sec:algorithm}
%
In the previous Section we have described the basic formalism and equations to evaluate the weak lensing effects on the full sky. In this Section we proceed outlining the basic steps of the algorithm used to lens the CMB photons throughout:
\begin{itemize}
\item[(i)] starting from an N-Body simulation, we create 3D
  simulated matter distribution around a chosen observer;
\item[(ii)] taking advantage of the proper sampling in redshift of the
  simulation, we select different shells of matter at different times to
  reconstruct our past lightcone and mimic cosmological evolution;
\item[(iii)] we project all the matter in a given shell over a single 2D
  spherical map which acts as lensing plane;
\item[(iv)] we solve the full-sky Poisson equation in the harmonic
  domain and compute the lensing potential map for the single lens plane
  and for the integrated potential of Eq.~\ref{eq:born_phi_dis};
\item[(v)] we use this lensing potential map to lens the
  CMB source plane; %applying the \lstinline!LensPix! algorithm;
\item[(vi)] we repeat step (iii)-(v) for all the selected shells, thus
  following the evolution of the source plane
\end{itemize}
\noindent
In our analysis we have used a N-Body simulation of cosmic
structure formation in a flat $\Lambda$CDM universe
with an underlying cosmology described by the following set of cosmological parameters: 
%% \begin{equation*}
%% \{ \Omega_{dm},  \Omega_{b}, \Omega_{\Lambda}, n_s, \sigma_8, H_0
%%   \} = \{ 0.2260, 0.045, 0.729, 0.96, 0.77, 72 \: \rm{Km / s /
%%     Mpc} \}. 
%% \end{equation*}
% \begin{equation*}
% \{ \Omega_{dm},  \Omega_{b}, \Omega_{\Lambda}, n_s, \sigma_8, H_0
%   \} = \{ 0.27, 0.05, 0.68, 0.96, 0.83, 67 \: \rm{Km / s /
%     Mpc} \}. 
% \end{equation*}
\begin{equation*} %CoDECS
\{ \Omega_{dm},  \Omega_{b}, \Omega_{\Lambda}, n_s, \sigma_8, H_0
  \} = \{ 0.226, 0.0451, 0.729, 0.966, 0.809, 70.3 \: \rm{Km / s /
    Mpc} \}. 
\end{equation*}
The simulation follows the evolution of the matter distribution in a
 cubic (comoving) volume $($1000 $h^{-1} \text{Mpc})^3$ from redshift
 $z = 10$ to present time using a modified TreePM version of
 \lstinline!GADGET!\footnote{\url{http://www.mpa-garching.mpg.de/gadget}},
 specifically developed to include all the additional physical effects
 that characterize different dark energy models (see \cite{Baldi12} for a
 detailed description of the code). The whole numerical project goes
 under the name of COupled Dark Energy Cosmological Simulation, or
 CoDECS\footnote{\url{www.marcobaldi.it/CoDECS}}. At present, they
 include two distinct set of publicly available runs, the L-CoDECS and the H-CoDECS. 
The L-CoDECS simulations consist in $1024^3$ CDM and $1024^3$ baryon
particles, both treated with collisionless dynamics only, which means that
baryonic particles are not considered as gas particles but just as a
different family of collisionless particles distinguished from CDM. This is done in order to account for
the effect of the uncoupled baryon fraction in cDE models which would not be correctly represented
by CDM-only simulations. The mass resolution at $z = 0$ for this set
of simulations is $M_{CDM} = 5.84 \times 10^{10} M_{\odot} / h$ for CDM and
$M_b = 1.17 \times 10^{10} M_{\odot} /h$ for baryons, while the gravitational softening
is set to $\epsilon_s = 20$ comoving $kpc/h$, corresponding to 0.04 times the mean linear inter-particle
separation. The H-CoDECS simulations are instead adiabatic
hydrodynamical simulations on much smaller scales, which we do not consider in the present work.
In this paper we will analyse only the $\Lambda$CDM simulation of the L-CoDECS suite, while the analysis on different DE models will be discussed in a future paper. 

\subsection{Constructing mass shells}
\label{sec:stacking}
N-Body simulations are usually stored as a series of snapshots, each
representing the simulation box at a certain stage of its evolution. 
%The algorithm uses these boxes, cutting a series of mass shells to reconstruct the backwards lightcone: a whole universe from a single, box-size-limited numerical simulation. 
As a first step, we fix the observer. We consider the last
snapshot, at redshift $z=0$, and compute the center of mass of the whole simulation. 
This centre represents the origin of our new system
of reference, which sees all the CMB
radiation around it. As we explore the universe around the new origin, the further
we move in space, the more we look back in time and see how structures
develop and grow, until we reach a volume large enough to carry out
the integration for CMB lensing. 
One of the difficulties in this approach is that, even though the size of the simulation box is limited, we need to use the box to reconstruct the full backwards lightcone. Therefore,
we need to replicate the box volume several times in space, so that the
entire observable volume is covered. 
In particular, as described in \cite{Carbone08}, the simulation volume needs to
be repeated  along both the positive and
negative directions of the three principal Cartesian axes x, y,
and z, keeping the origin centered on the observer. 

To construct the all-sky past light cone we exploit the simulation
outputs at different times which are equally spaced in the logarithm
of the scale factor, $log_{10}(a)$, corresponding to an average spacing of
150 Mpc$h^{-1}$ comoving.
The need to repeat the simulation volume due to its finite size
immediately means that, without augmenting large-scale structures, the
maps will suffer from a deficit of lensing power on large angular
scales, due to the finite box size. More importantly, a scheme is
required to avoid the repetition of the same structures along the line
of sight. Previous studies that constructed simulated lightcone maps
for small patches of the sky typically simply randomized each
of the repeated boxes along the past lightcone by applying independent random translations and reflections 
(e.g. \cite{Springel01}). However, in the present application
this procedure would produce artifacts like ripples in the
simulated deflection-angle field, because the gravitational
field would become discontinuous at box boundaries, leading
to jumps in the deflection angle. It is therefore mandatory
that the simulated lensing potential of our all sky maps is
everywhere continuous on the sky, which requires that the
3D tessellation of the peculiar gravitational potential is
continuous transverse to every line of sight.

Following \cite{Carbone08, Carbone09}, our solution is to divide up the volume out to the redshift covered
by the simulation $z_{max}$ into larger spherical shells, each of
thickness 1 Gpc$h^{-1}$ comoving (as the box size)
All the simulation boxes falling into the same larger shell are made to undergo
the same, coherent randomization process, i.e. they are all
translated and rotated with the same random vectors generating a
homogeneous coordinate transformation throughout the shell. But this
randomization changes from shell to shell.
% Once it is copied everywhere, we exploit the time resolution and
% dense sampling in redshifts of our simulation to proceed. Each output has a
% spacing in time logarithmic in the scale factor, $\log a$,
% corresponding
% to an average distance of $\sim 150 \: h^{-1}$ Mpc (comoving)
% on the past light-cone, i.e., data corresponding
% to a particular output in time are assigned to a spherical shell having the same thickness around the observer. The existing periodic boundary conditions in the simulation also guarantee a smooth
% matter distribution between the box edges.
% More importantly, a scheme is required to avoid the repetition
% of the same structures along the line of sight, 
% and we must be sure that the simulated lensing potential is
% continuous in the direction transverse to the line of sight.
%  For this purpose, every time a new set of period replicas is required,
%   i.e. when we reach a comoving distance greater than the box size of
%   the simulation, 
% all the replicated boxes undergo
% the same, coherent randomisation process. This requires that all the boxes are
% translated and rotated with respect to the same random vectors, generating thus
% an homogeneous coordinate transformation. When the integration along
% the comoving distance is
% limited within the boundaries of the simulation box, the coordinate transformation
% is kept fixed while it changes (generating a new set of random
% vectors) when we need to fill the volume of the
% lightcone by stacking again the simulation box around the observer.
% All the replicated simulation
% boxes falling into the same spherical shell undergo
% the same, coherent randomization process, i.e. they are all
% translated and rotated with the same random vectors generating
% an homogeneous coordinate transformation throughout
% the replicated boxes.
%%Figure~\ref{fig:stacking} shows a schematic sketch of this stacking
%%technique.
\noindent
See Figure~1 of \cite{Carbone08} for a schematic sketch of this stacking technique.
%Figure \ref{fig:stacking} shows a schematic sketch of this process,
%which was also exploited in earlier works (Carbone {\it et al.}, 2008;
%Carbone {\it et al.}, 2009 \cite{Carbone08,Carbone09}).
%Note that this  eliminates any preferred direction in the simulated
%all-sky maps (for instance, used by Hilbert {\it et al.}, 2009 \cite{Hilbert09}).
As already mentioned, the need to repeat the simulation volume due to its finite
size implies that the maps will
suffer from a non proper description of the large angular scales. We
note however that if the box size is sufficiently big like in
e.g. \cite{Watson13} this procedure is no longer necessary, at least
up to the redshift covered by the simulation size.
The final results of the whole process is a series of concentric shells that substitutes our
snapshots. For our specific input N-body simulation, we get 25 matter
shells, building a lightcone spanning from $z=0$ to $z_{max} = 10$.

% %
% \begin{figure}
% \centering
% \includegraphics[scale=0.5]{stacking.ps}
% \caption{
% Sketch of the adopted stacking process. The crossing
% of CMB photons through the dark matter distribution
% of the Universe is followed by stacking and replicating the boxes of the
% N-Body simulation to fill the entire shell volume. Each shell has a
% thickness of $\sim 150 h^{-1}$ Mpc (comoving).
%  }\label{fig:stacking}
% \end{figure}
% %
\subsection{From shells to maps}
\label{sec:maps2shell}
Following the scheme proposed in \cite{DasBode}, we convert the position of a
particle distributed within a 3D matter shell into its angular position on the 2D
spherical map of the (projected) surface matter density. Note
  that among all the particles in the simulation, only the ones falling
  within the radii of the spherical shell of width 150 Mpc$h^{-1}$ are projected onto the 2D
  spherical map. 
We then assign each particle to a specific sky pixel in the \lstinline!HEALPix! \footnote{\url{http://healpix.sourceforge.net}} pixelisation scheme starting from its spherical
coordinates ($\theta$,$\phi$)  and using the \lstinline!ang2pix!
routine of the \lstinline!HEALPix! suite. We place the particle
mass into the pixel so that each sky
pixel reads $\Sigma^{\theta}_p = m_p / \Omega_{pix}$, where
$\Omega_{pix}$ is the area of a pixel in steradians. If $n$
particles fall inside the beam defined by a pixel, the pixel will have a surface mass density of $n\Sigma^{\theta}_p$. In
\lstinline!HEALPix!, the resolution is controlled by the parameter
\lstinline!NSIDE! which determines the number $N_{pix}$ of pixels of equal area
into which the entire sphere is divided through the relation
$N_{pix} = 12 \times $\lstinline!NSIDE!$^2$, so that each pixel has
an area of $\Omega_{pix} = 4 \pi/N_{pix}\: \textrm{sterad}$.
The angular resolution is often expressed through the
number $\theta_{res} = \sqrt{\Omega_{pix}}$.
For a value of \lstinline!NSIDE! set to 2048 (4096), the corresponding
an angular resolution is $1^{\prime}.717$ ($0^{\prime}.858$). \\*
%
The real interesting quantity in our lensing calculation, in addition to
the surface mass density, is the convergence map $K^{(k)}({\bf
  \hat{n}})$ on the lens-plane $k$. To get this quantity we first
  compute the overdensity maps ($\Delta \Sigma^{\theta}$) using the average surface
mass density of the 2D map. Then we multiply by this map by its geometrical weight $(1+z_k)/D_A(\eta_k)$,
depending on the lens plane distance from the observer and its
redshift, assumed to be an average between the time at the beginning
and the end of the shell (see Eq.~\ref{eq:key}). As a final step, we produce a convergence map
from each shell of the  lightcone which will become our
lensing planes to lens the CMB signal.\\*
%
%A small technical details: to refine our calculation, we decide to cut
%another set of smaller slices inside the bigger shell of universe made up with the stacking
%technique. The idea is to decompose the latter into
%smaller, $30 h^{-1}$ Mpc slice of matter; the purpose is to mimic the integration of
%the lensing potential by having more convergence maps with different
%weights in redshift. In the end, about 300 such shells where produces from the simulation (50
%snapshots, about 6 $30 h^{-1}$ Mpc shell per snapshot), 300 convergence maps
%alias lensing plane ready to be used for our calculation. 
%
From the convergence map $K^{(k)}(\hat{{\bf n}})$ we then extract the
gravitational potential $\psi^{(k)}$ following
Eq.~(\ref{eq:kappa_to_phi}), using the \lstinline!HEALPix! spherical
harmonics transform (SHT) routines to decompose $K^{(k)}(\hat{{\bf
    n}})$ into its harmonic coefficients $K_{\ell m}$. Note that we
correct the smoothing of the true underlying continuous field on the
pixel scale directly in the harmonic domain when we solve for the
Poisson equation (see Appendix \ref{appendix:ps} for more details).

\subsection{Lensing the CMB}
To propagate the CMB photons through the different shells we adopt a
pixel-based approach first presented in \cite{Lewis05}. Starting from
the $\psi_{\ell m}$ coefficients, we compute the deflection field for
each shell $\boldsymbol\alpha^{(k)}$ evaluating Eq.~(\ref{eq:gradphi})
in the harmonic domain. Being the deflection field purely a gradient
field (i.e. a spin-1 curl-free vector field), it can can be easily evaluated
with a spin-1 SHT. The E and B decomposition of the field reads:
\be
_{1}\alpha^{E}_{\ell m}{}^{(k)}=\sqrt{\ell(\ell+1)}\psi^{(k)}_{\ell m} \qquad _{1}\alpha^{B}_{\ell m}{}^{(k)}=0.
\ee
Once the deflection field is give, each pixel-based method remaps the CMB field
as a function of the position on the sky assuming the lensed signal
observed along a direction ${\bf \hat{n}}$ equal to the signal
coming from another direction ${\bf \hat{n}^{\prime}}$, 
%We take advantage again of the analysis and synthesis operation of the
%\lstinline!HEALPix! implementation: given $\Phi^{(k)}(\hat{{\bf n}})$, and its
%expansion set $\Phi^{(k)}_{\ell m}$, we can easily construct maps
%of the gradient from the harmonic components of the lensing field
%itself \cite{Lewis05}. Once we obtained the deflection field $\boldsymbol\alpha^{(k)}$, we can
%proceed to re-map the anisotropy in the CMB as functions of the new
%position ${\bf \hat{n}^{\prime}}{}^{(k)}$, with:
\be
{\bf \hat{n}^{\prime}}{}^{(k)} = \hat{{\bf n}}^{(k-1)} + \boldsymbol\alpha^{(k)},
\ee
where $\hat{{\bf n}}^{(k-1)}$ represents the unlensed position of the
CMB photons at the previous step. To our level of approximation $\| \nabla \Phi \|$ is
assumed to be constant between ${\bf \hat{n}}^{(k-1)}$ and ${\bf
  \hat{n}^{\prime}}{}^{(k)}$, consistent with working out the lensing
potential in the Born approximation between two consecutive
shells. %Pixel-based methods, though in general very efficient, are subjected to several numerical problems which are discussed in detail in Sect. \ref{sect:numerical-pixel-based}.
%The first one is related to the bandwidths of the lensed CMB fields
%generated as a result of the calculation, which are naturally not
%band-limited because lensing acts as a convolution in the harmonic
%domain \cite{Hu2000}. The other challenge arises from the fact that
%the displaced direction at each iteration ${\bf
%  \hat{n}^{\prime}}{}^{(k)}$ does not correspond in general to pixel
%centers of any iso-latitudinal grid on the sphere. The values of the
%CMB signal at those locations cannot be computed with the aid of fast
%SHT algorithm and a more elaborated approach is needed.
In this work we adopted the publicly available code
\lstinline!LensPix! to propagate the CMB signal through all the
lensing shells. \lstinline!LensPix! implements a pixel-based lensing
method using a bi-cubic polynomial interpolation scheme to evaluate the
source plane along the displaced direction. This method has been shown
to be accurate at the sub-percent level to produce temperature and
E-modes signal. However, the recovery of the B-modes of the CMB polarization
is more difficult because B-modes are more sensitive to numerical
effects like the involved resolution and the choice
of the band-limit (i.e. the power cut-off $\ell_{max}$) in the calculation. We will
discuss the impact on the relevant numerical effects in Section~\ref{sect:numerical-pixel-based}
and we refer the reader to \cite{Fabbian13} for a complete discussion
of the numerical problems and accuracy of pixel based methods.\\*
%, leaving to future studies the implementation of more advanced techniques for interpolation in the sub-arc minute regime.
%%% The correct remapping equations can be easily derived from identities of spherical
%%% triangles. 
%%% For completeness, we give the
%%% derivation in the Appendix A.
%%
%It is worth pointing out that the re-mapping is done in a
%much faster way by generating an un-lensed sky at higher resolution;
%\cite{Lewis05} found that using $12 \times 8192^2$ $\sim 10^9$ pixels
%for the high resolution map is sufficient to get a $12 \times
%1024^2$-pixel lensed sky with $C_{\ell}$  accurate to $0.5$\% at
%$\ell \sim 2000$. The fractional accuracy on the polarization is similar, except
%for the B-modes induced by lensing which are a more sensitive: $12
%\times 16384^2 \sim 3 \times 10^9$ pixels are required to get a lensed
%power spectrum $C^{BB}_{\ell}$ accurate at percent level for $\ell \sim
%2000$. The B-spectrum is quite sensitive to $\ell_{max}$ and the
%non-linear power spectrum, so that this accuracy level is only
%notional. Note that choosing $\ell_{max}$ too low generally
%underestimates the B-mode power, whereas pixelization errors lead to
%an overestimation; it is possible to get spuriously accurate power
%spectra without actually having computed the lensed sky accurately. 
%
%In this work we adopt the \lstinline!LensPix! interpolation scheme, leaving to
%future studies the implementation of more advanced techniques for
%interpolation in the sub-arch minute regime.
%
Finally, note that the simulated lightcone recovers the distribution of
matter in the Universe up to
$z_{max}=10$, and therefore the primordial CMB fields are lensed by LSS
only in this specific redshift interval. In other words, photons are ray-traced
in a Universe evolving from $z_{max}$ to $z=0$. The impact of
high-redshifts contributions is besides the goal of this algorithm,
which will be tested against analytical and semi-analytical
computations which we have modified accordingly to perform CMB lensing
only in this redshift range. 

conclusions.tex
%!TEX root = jcap_article_MLcode.tex
% !TEX spellcheck = en-EN

\section{Discussion and Conclusions}
\label{sec:concl}

We have developed and tested a new algorithm to study the gravitational lensing of the
CMB on the full sky. Starting from snapshots of an N-Body
simulation, we reconstructed the whole
lightcone around the observer between $z\in [0,10]$ overcoming the
finite size of N-Body box through the box stacking technique developed
in \cite{Carbone08}. We sliced the
lightcone into 25 different spherical maps onto which the matter
distribution was projected. The spherical shells were then used as
source planes to lens the incoming CMB photons either adopting an
effective method based on the Born approximation or using a
multiple-lens ray-tracing approach. This is, to our knowledge, the
first attempt to apply this kind of algorithm to the CMB polarization. 
%In such approach, the continuous deflections of photons is approximated by a series of multiple lens-planes that bend the
%light rays. In practice, this means that at each step the CMB fields 
%are re-mapped from their previous positions in the sky to a new set of
%coordinates, with the displacement given by the intensity of matter field the
%light is passing through.\\*
%
For this reason, we performed a detailed analysis of the numerical effects involved in the ML method coming both from the N-Body simulation and from the ray-tracing procedure itself. 
%multiple plane ray-tracing and the effective, Born approximation version of the algorithm 
The Born approximation, which has been widely tested in the literature, was used as benchmark to
highlight the multiple-lens range of validity and to asses its virtues in
reproducing non-linearities from the N-Body simulation at small scales. 
In particular, the projection of the N-Body matter distribution onto
concentric spherical maps allows to compress all the interesting
information from the N-Body simulation %(large scale structures, non-linearities, ...)
into a more manageable lightcone, mimicking a realistic distribution of large scale structure as
observed by present and future large galaxy surveys.\\* %carlo + giulio
%
%The first tests show the robustness of our map-making procedure, and
%the fact that our lensing quantities (both the convergence and the
%potential) fairly reproduce and overlap semi-analytical
%realizations. %carlo 
%At the level of maps, the algorithm permits to compress all the interesting
%information from the simulation (large scale structures, non-linearities, ...)
%into a more manageable lightcone, which mimics the real Universe as
%observed by present and future satellite missions. %carlo 
%We perform a statistical validation based on the two point correlation functions 
%in the harmonic domain, evaluated in each redshift shell across the lightcone, 
%by computing the angular power spectra and maps for lensing observables:  
%we are able to reconstruct the signal up to very small scales, grasping the intrinsic 
%non-linearities of the simulation all across the various stages of their evolution 
%through the lightcone. 
%
We validated the lensing planes reconstruction both on the map level
and on the statistical level using the 2-point correlation function in
the harmonic domain evaluated for all the spherical maps constructed
across 
the past lightcone. We found the latter to be reproduced fairly well
by semi-analytical approximations to the non-linear evolution
implemented in widely used Boltzmann codes, though deviations at the
percent level were clearly observed.
%
We also analysed the final, lensed CMB anisotropies in both temperature and
polarization for the effective as well as the multiple
plane approach, paying a particular attention to the B-modes of polarization. 
%As for the basic lensing quantities at various redshifts, 
%the resulting angular power spectra are in a very good agreement
%with semi-analytical realizations up to $\ell \sim 4100$; 
These are in fact the most sensitive quantities both to the overall lensing
process and to the numerical effects. In the latter case, we discussed
in detail how to minimize their impact.
% , and we have taken into
%account various numerical artifacts that can spoil this signal. In
%particular, we have computed - both analytically and using Monte Carlo
%techniques - the influence of the shot-noise coming from the finite
%resolution of the numerical simulation. 
We found however that these numerical effects are usually negligible for the temperature and E-modes polarization field and important only for B-modes. The B-modes signal was found to be lower than the one computed using the semi-analytical, following the general trend observed in the extracted lensing potential power spectrum.\\*
%negligible for the total intensity as well as
%E-mode polarization field, while - as expected - is important at
%percent level to reconstruct the BB signal. Moreover, we have adopted a tight
%constraint ($\ell_{max} =$ \lstinline!NSIDE!) to keep residual power aliasing under
%control while synthesize the maps in the harmonic domain. The final, lensed
%BB signal is found to be in excess over a wide range of multipoles with
%respect to the semi-analytical angular power spectrum, accordingly to
%what results in the lensing potential computed from the simulation. This is in fact welcome, since
%it underlines the ability of lensed BB signal to recover the correct
%evolution of structures in the Universe. 
%
%carlo + giulio
Finally, we have extended the control of the validity of the Born approximation to the limiting 
resolution of the present setup of our simulations. Our results indicate that, when checking the 
angular power spectra of lensing observables, including CMB lensed fields, the latter approximation 
describes well the ray-tracing performed by the ML approach. However we expect the latter to perform 
better for studies aiming at investigating the statistics of the signal at smaller angular scales, 
or in presence of distortion from isolated, sharp structures. \\*
For what concerns the behaviour of the signal at large
multipoles $l>3000$, in a recent paper \cite{Hagstotz14} the authors model corrections to the Born
approximation by using perturbation theory applied to the lensing magnification matrix. 
This kind of analysis was first presented in \cite{CoorayHu02} and \cite{ShapiroCooray06}, who 
computed the second-order corrections to the angular power spectra of 
the lensing observables (convergence, shear and rotation). 
By adopting the Peacocks \& Dodds  matter 
power spectrum \cite{PeacockDodds96}, they concluded that these corrections are not relevant for the galaxy 
weak-lensing case being two orders 
of magnitude lower than the first-order contribution. On the other 
hand, using the matter
 power spectra extracted from \lstinline!CAMB!, \cite{Hagstotz14} applies 
the same analytical framework of \cite{CoorayHu02} and
\cite{ShapiroCooray06} to the CMB-lensing 
case, and finds a large excess in the B-modes power spectrum with respect to 
the first-order contribution. For a cross-check of these results, following the same analytical setup, we have 
independently computed second order corrections to the CMB convergence power 
spectrum, using the \lstinline!CAMB! matter power spectra as input. In agreement with 
\cite{Hagstotz14}, we find that these corrections seem to affect the CMB
 lensing potential at very small scales and, consequently, the B-mode power
 spectrum at all the multipoles. We have also applied the same analytical second order computations to the galaxy weak-lensing
 case with sources at $z_s=1,2$, using again \lstinline!CAMB! matter power spectra 
as input. In this case, we find in the convergence power spectrum an
 excess of power at $\ell \approx 2000$ of the order of $10\%$ with 
respect to first order contributions, largely 
in tension with the numerical
analyses of multiple-lens ray-tracing from N-body simulations present in the literature \cite{Hilbert09,Becker12}, which conversely find no
evidence of important differences with respect to the Born 
approximation for $z_s=1,2$. \\*
Given this tension, it is necessary to investigate in more detail the 
validity of the analytical corrections. To this purpose, we are developing an 
improved numerical setup to properly simulate both the CMB lensing 
per se and the LSS evolution below the smallest angular scales
 considered in this work.  
We plan to address this issue in a future paper including also an 
extension of our formalism to propagate the whole lensing
 magnification matrix, in order to trace more accurately all second-order
corrections for all the lensing observables. 
% using numerical simulation and by \cite{} extending the perturbation
% theory to fourth order in the potential.
% With the current set-up of our machine we are not able to
% reach an accuracy high enough to be able to distinguish this signal,
% which is ``masked'' by the presence of the shot-noise from N-Body 
% simulation especially for small $\ell$ (see Figures??).  
%However we stress that the
%analytical calculation may require a further expansion of the series,
%such that higher perturbative order corrections may cancel out the
%contributions computed by \cite{CoorayHu02,ShapiroCooray06}.


%
%carlo 
%Our work sets the ground for high resolution CMB lensing analyses beyond the statistics of two point 
%correlation functions. Alternative and more efficient control of possible aliasing power introduced 
%by re-mapping due to the deflections of CMB photons passing though LSS will be tested in a forthcoming 
%work. 
% \begin{figure}[!htb]
% \centering
% \includegraphics[scale=0.55,angle=90]{Images/NoiseSub_BB_integral.ps}
% \caption{Angular power spectrum for the lensing B-modes. The signal is
%   computed at different redshifts as it evolves in time.}
% \label{fig:bb_lensed_int}
% \end{figure}
The spherical map matter projection, in both the Born
approximation and ML implementations, can be particularly useful in
cross-correlation
studies between the CMB with other tracers of mass 
and foreground sources, for characterizing the simulation of mock
catalogues of observables built from N-Body simulations. The
tomography of LSS, which is an intrinsic feature of the two lensing
approaches analysed in this work,
can be exploited to investigate different cosmological scenarios,
looking at the effects of different DE models on small scales as well
on the whole evolution of the matter in Universe. These feature will
be of great importance for upcoming projects such as the Euclid
satellite that can fully exploit the capabilities of cross-correlation
as cosmological probe. 

intro.tex
%!TEX root = jcap_article_MLcode.tex
% !TEX spellcheck = en-EN

\section{Introduction}
 
The recent measurements of the Planck satellite \cite{PlanckXVI} has unveiled a
Universe well described by a cosmological model known as
$\Lambda$CDM. In this framework, the Universe is expanding 
accelerated by a Dark Energy (DE)
component well described by a cosmological constant $\Lambda$. Cold Dark Matter
(CDM) is responsible for the matter halos around galaxies and galaxy
clusters, while leptons and baryons take a minor part in
the entire cosmic energy density budget.
%the primordial density spectrum is Gaussian within errors, with a
%power spectrum which shows a small departure from scale
%invariance.
%, $P(k) \propto k^{n_s}$ , where $n_s \sim 0.96$. 
%No direct detection of
%cosmological Gravitational Waves (GWs) is proved to date. 
The robustness of the $\Lambda$CDM model comes primarily from 
advanced and enduring studies on the 
anisotropies in the CMB, combined with
LSS and Type Ia Supernovae observations.
Cosmological perturbations, scalar
(such as density), vector and tensor ones can be
directly observed as primary anisotropies in the CMB sky, imprinted via Compton
scattering at the epoch of recombination. 
%At recombination, matter perturbations are imprinted and
%stored by Thomson scattering in photons perturbations, producing a total intensity
%(T) anisotropies as well as linear polarisation ones, described by the Stokes
%parameters Q and U. The latter can be further combined into
%gradient modes (or electric-type, E)
%and curl ones (or magnetic-type, B). At the linear level, while the E
%modes are sourced by all kinds
%of perturbations, B modes are excited only by vorticity and GWs (see e.g.
%\citep{Hu2000} for review). 
%The Planck satellite was able to show these CMB anisotropies in total
%intensity (T) over the
%whole sky down to a resolution of 5 arc-minutes, confirming theoretical
%predictions of $\Lambda$CDM.
%with a sensitivity of a few 
%$\mu K$ per resolution element. 
%(as a comparison, the Wilkinson
%Microwave Anisotropy Probe
%(WMAP) \cite{WMAP} reached an angular resolution of about 20
%arc-minutes). 
A step further in terms of constraining power on cosmological
parameters will be achieved with the high accuracy measurements of the
CMB polarization provided by the upcoming Planck polarization data and current and future high sensitivity ground-based and balloon-borne experiment like, e.g, EBEX \cite{EBEX}, and POLARBEAR \cite{Polarbear}, SPTpol \cite{SPTpol}, ACTpol \cite{ACTpol}, Spider \cite{Spider} and Keck array \cite{keck}\footnote{See \url{http://lambda.gsfc.nasa.gov/links/experimental\_sites.cfm} for a complete list of all the latest missions and upcoming experiments}. These are designed to measure in particular the so-called B-modes of polarization \cite{Zaldarriaga96, Kamionkowski96} which could provide a direct evidence for primordial gravitational waves (tensor perturbation) generated in many inflationary scenarios \cite{Seljak97} if such signal is detected on the degree scale. The BICEP2 collaboration recently announced a robust measurement of a signal compatible with inflationary B-modes \cite{bicep2} and upcoming dataset will be crucial to understand if the signal is genuinely primordial or contaminated by diffuse astrophysical polarised emission \cite{flauger}. \\*
%
%However, despite of the accuracy and precision of the Planck data, we still do not know
%much about the dark sector, neither its physics, nor its origin, nor
%its behaviour. Again CMB data can come to our rescue: 
In addition to primary
anisotropies previously described, the so-called secondary anisotropies \cite{Aghanim} are
imprinted in the CMB by the interaction of its photons 
with LSS along their
paths from the last scattering surface to the observer. 
%% To start, electrons
%% freed by cosmological reionisation at the beginning of structure
%% formation can re-scatter CMB photons: a boost in polarisation
%% anisotropies is detected on the angular scales subtended by the horizon at that
%% epoch. Secondly, photons can be affected by time-evolving
%% gravitational potentials when they travel across the LSS, the so-called
%% Integrated Sachs-Wolf (ISW) effect, which is now been observed by Planck \cite{PlanckXIX}.
%The subject of this work concerns the possibility 
%of LSS to act as gravitational lenses:
One of the most important sources of secondary anisotropies is the
gravitational lensing effect on CMB photons by on growing matter
inhomogeneities which bend and modify the geodesic path of the
light. 
The net-effect of those deviations is a reshuffling of power in the modes of the
primordial power spectrum of the total intensity
 and gradient-like (E-modes) component of the
CMB polarization. Moreover, lensing distorts the primordial polarization
patterns converting E-modes into B-modes pattern, generating power on
the sub-degree scale where we expect the primordial signal to be negligible (for a complete review, see \cite{LewisChallinor06}).
%We take advantage of the fact that the CMB 
%lensing cross section is non-zero at
%redshifts roughly between 0.1 and 10, peaking at
%$z \simeq 2-3$, rather independently of the particular cosmological
%model considered \cite{AcquavivaBacci06,Smith06}. Moreover,
%since the lensing efficiency is related to the behaviour of
%the expansion and perturbations
%at the corresponding epoch, the CMB lensing signal carries important
%information about the cosmic acceleration, a direct probe 
%of the DE behaviour. 
The progress
towards the detection of lensing in the CMB data has been slow, since measurements of the CMB
precise enough to enable a detection of weak lensing were
not available until recently. Moreover, picking out
non-Gaussian signatures - which arise from mode mixing on different
scales induced by lensing - in the measured CMB sky by
itself is extremely difficult, due to confusion from systematics,
foregrounds, and limited angular resolution.
The first robust detections of the lensing effect were
done using temperature data only by ACT \cite{ACT} and SPT \cite{SPT, vanEngelen}, and later
confirmed by Planck with a significance greater than 25$\sigma$. Only
recently, however, the evidence of lensing was detected for the first
time in polarization data by POLARBEAR \cite{pblens, pbbb} and SPTpol
\cite{hanson13}.\\* 
The interest on CMB lensing for cosmological application lies in the possibility of extracting information about the projected large scale structure potential, and thus to
constrain the late time evolution of the Universe and, e.g. the Dark
Energy and massive neutrinos properties \cite{AcquavivaBacci06,
  Smith06, hu-dark,LewisChallinor06}. A step forward to probe DE would be to cross correlate the CMB lensing with the
observations of the actual lenses in LSS surveys as seen by
independent tracers of the matter distribution. This option has
already been exploited to obtained astrophysical and cosmological
information by, e.g., SPT, ACT and, more recently, by POLARBEAR
collaborations  \cite{Bleem12, Sherwin12,act-cfht, pb-cib,hanson13},
but a major improvement is expected in about a decade with the
observations of the ESA-Euclid satellite. 
This will combine arc-second imaging of billions of galaxies with photometric
redshift accuracy corresponding to the percent level, 
between redshifts $0\lesssim z\lesssim 2$
\cite{Laureijs11,Amendola12}. 
%The mission primary
%target will be to constrain the DE behaviour in the
%corresponding redshift interval, correlating maps of galaxy clustering
%with weak lensing shear associated to the induced ellipticity in
%galaxy images. 
To fully exploit the capabilities of these cross correlation studies \cite{pearson-zahn}, accurate prediction for the common observables are crucial. \\*
%
The most accurate way to obtain predictions for observables of weak-lensing surveys
is to perform ray-tracing through large, high resolution N-Body numerical
simulations to study the full non-linear and hierarchical
growth of cosmic structures. Although this approaches are
computationally very demanding, they allow to check and balance 
the approximations and assumptions made in widely-used semi-analytic
models, adjusting and extending these models if necessary.
Numerous ray-tracing methods have been developed so far in the context
of both strong and weak gravitational lensing. %The main idea is to
                                %follow the path of a bunch of light
                                %rays drawn from a known source,
                                %through the matter distribution in
                                %the Universe. 
Though exact algorithm are available \cite{Killedar11}, they are not
suitable for application targeting observations of large fraction of
the sky for computational reasons. 
%One direct approach
%\cite{Killedar11} is to integrate numerically the null geodesic from the
%observer to source; a highly precise method but time-consuming for
%full-sky applications. 
A simpler and popular approach consists in using the matter
distribution in the N-body simulations to calculate 
lensing observables by photon ray-tracing along ``unperturbed'', i. e.
undeflected light paths in the so-called Born approximation
(e.g. \cite{Hilbert07a, Couchman99,Carbone09}). %In this contest \cite{Couchman99,Carbone09}, the lensing signal is obtained by integration of the three-dimensional gravitational potential. 
In particular, \cite{Carbone13}
applied this technique to study a set of N-Body simulation with
different cosmologies and DE dynamics, to investigate the variation of
the lensing pattern with respect to the standard $\Lambda$CDM model. \cite{antolini}, conversely, showed that the
correct integrated matter distribution used to lens incoming CMB
photons in the Born approximation can be properly
reconstructed using standard lensing reconstruction techniques
\cite{hu-okamoto}. %In \cite{Fosalba08} a lightcone is made of several matter shells from a N-Body simulation to create an all-sky lensing map.\\*
However, when facing a complex cosmological structure, we must take
into account that each light ray undergoes several
distortions due to matter
inhomogeneities, i.e. approximating the actual path of a photon instead of adopting a single effective deviation from the unperturbed, line-of-sight integral assumed in Born
approximation. The single effective lens must therefore be replaced by
a multiple-lens (ML) approach, where large volumes of matter are
projected onto a series of lens planes \cite{BlandfordNarayan, Jain00, Pace07, Hilbert09,Becker12} so that the
continuous deflection experienced by a light ray is approximated by finite deflections at each
of the planes. %Many simulations of weak lensing by clusters and LSSs
%\cite{Jain00,Pace07} employ this method to trace light
%rays. \cite{Hilbert09} use ray-tracing through the Millennium
%Simulation to calculate cosmic-shear and galaxy-galaxy-lensing
%statistics.
%, finding that first-order analytical prediction of convergence power
%spectrum agrees with ray-tracing results.
% \cite{Becker12} combines the
%previous techniques to compare the galaxy shear B-modes and rotation power
%spectra at a very high resolution. 
A ML algorithm full-sky CMB lensing application was sketched in 
%Recent simulations of CMB lensing use generalizations of the single- or multiple plane approximation, to take into account the curvature of the sky: 
\cite{DasBode} who simulated lensed CMB maps in temperature only with
an angular resolution of $0.9'$. However detailed comparison of the
effective Born approximation method and the ML approach was not
discussed and only the results derived in the Born approximation scenario were presented.\\*
%using an highly accurate interpolation method to enable sampling of
%the deflection field (and the un-lensed CMB) on an irregular grid.
%

In this work, we implemented a multiple plane ray-tracing
algorithm to lens CMB temperature and polarization fields combining the aforementioned
work by \cite{DasBode} and using the lightcone reconstruction from a single N-Body simulation of \cite{Fosalba08}.
The final rationale will be to investigate DE effects in different cosmologies at
arc-minutes scales, where are expected to be most noticeable. At these
scales, the Born approximation is likely expected not to trace with
high accuracy local deviations due to
small-scales inhomogeneities, and thus a more precise and realistic method is needed.
Moreover, in order to be successful, we need to be able to control and
discriminate between physical non-linearities of the
N-Body simulations from numerical issues connected to the various approximations in
both the lensing algorithm and the simulation itself. A detailed
analysis of these issues together with their impact on the lensed CMB
observables will thus be presented.
%In this first paper
%we go beyond the Born approximation, starting from the aforementioned
%work by \cite{DasBode} where they show results for the temperature field
%only computed within the
%Born approximation scenario. 
%%Here, we present and validate a multiple plane
%%ray-tracing algorithm to lens photon coming from the last scatter
%%surface, showing how the CMB is affected in both temperature and
%%polarization. 
%%In this way we can better trace local deviations due to
%%small-scales inhomogeneities, as well as second-order effects arising
%%from the intrinsic non-linearity of the N-Body simulation, which are
%%naturally washed out by first-order approximation. 
%For the first time
%in the literature, we analyse the lensed temperature and polarization fields at a very high
%resolution, carefully taking under control numerical
%artifacts arising from the interpolation scheme adopted, as well as
%modelling the shot-noise spectrum induced by the N-Body
%simulation. 
%%All the results are then tested against the
%%semi-analytical computations currently available, with the multiple
%%plane approximation always compared to the effective/Born case.
%The starting point of our analysis is a simulated box, which is used by
%our algorithm to produce a past light-cone composed by different shell
%of matter, which represents a set of spherical lens planes. Each plane will then lens
%for the photons coming from the last scattering surface, gradually
%building up the signal of a lensed CMB field.
%%Notice that we have to extend the standard formalism of flat planes to the sphere and
%%``curved'' planes, to account both for the curvature of the sky and to use
%%the spherical harmonics formalism, commonly used in CMB analysis. 
%The exact simulation of the lensed CMB sky requires the computation of
%spin-0 and spin-2 spherical harmonics on an irregular grid, which is caused by the
%deflections of the photons from their original positions on the CMB
%surface. For this work, we follow the approach by \cite{Lewis05}, 
%%This direct method is computationally expensive
%%and requires robust parallelization. \cite{Lewis05} suggested an
%%alternative 
%where we would re-sample the un-lensed CMB sky, generated with finer resolution,
%at these deflected positions. This approach was implemented in the
%publicly available
%\lstinline!LensPix!\footnote{\url{http://cosmologist.info/lenspix}} code.
%The advantage of this whole approach is the freedom to test
%whatever cosmological scenario we want, simply by changing the
%N-Body simulation we would like to analyse. 
The present paper is intended to be the
first one of a series investigating the response of CMB lensing upon
DE and/or modified gravity cosmologies or neutrino physics, as well as the
feasibility and the constraining power of cross-correlation studies involving Planck and
simulated Euclid data. %The tomography of the matter distribution -possible within the lightcone approach - would give us insights on how the primordial CMB gets gradually lensed; the final results can be tested against current de-lensing algorithms that permit to ``clean'' the observed, lensed CMB fields from the intrinsic lensing signal.
%
%In a forthcoming paper we
%will use this algorithm to test different DE scenario as well as the
%impact of neutrino physics on lensing.  
%
This paper is organized as follows: in the Section 2, we introduce
the theoretical background and notation used for our lensing
algorithm. In Section 3 we discuss our ray-tracing technique, which is
then tested and evaluated in Section 4. Section 5 shows results for the lensed
temperature and polarization fields of the CMB,
with particular emphasis on differences and similarities between the
Born approximation and the multiple-lens approach. The last Section
draws the conclusions.

jcap_article_mlcode.tex
% !TEX encoding = UTF-8
% !TEX spellcheck = en-EN

\documentclass[a4paper,11pt]{article}
%\pdfoutput=1 % if your are submitting a pdflatex (i.e. if you have
             % images in pdf, png or jpg format)

\usepackage{jcappub} % for details on the use of the package, please
                     % see the JCAP-author-manual

%\usepackage[T1]{fontenc} % if needed

\usepackage{graphicx}
%\usepackage{epstopdf}
%\epstopdfsetup{update}
\usepackage{booktabs}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{subfig}
\usepackage{url}
\usepackage{listings}
\lstloadlanguages{c++}\lstset{basicstyle=\small\ttfamily}

%\usepackage[round]{natbib}
\usepackage[]{natbib}

\renewcommand{\vec}[1]{\bmath{#1}}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\ba}{\begin{eqnarray}}
\newcommand{\ea}{\end{eqnarray}}
\newcommand{\brr}{\begin{array}}
\newcommand{\err}{\end{array}}
\newcommand{\bc}{\begin{center}}
\newcommand{\ec}{\end{center}}
\newcommand{\hG}{\,h^{-1}{\rm Gpc}}
\newcommand{\hm}{\,h^{-1}{\rm Mpc}}
\newcommand{\hk}{\,h^{-1}{\rm kpc}}
\newcommand{\msun}{\,h^{-1}{\rm M}_\odot}
\newcommand{\hMpc}{\mbox{h^{-1}{\rmn{Mpc}}~}}
\newcommand{\rvir}{\mbox{$R_{\rmn{vir}}$}}
\newcommand{\mvir}{\mbox{$M_{\rmn{vir}}$}}
\newcommand{\tmw}{\mbox{$T_{\rmn{mw}}$}}
\newcommand{\hMpcI}{\mbox{h\,{\rmn{Mpc}}^{-1}}}
\newcommand{\lb}{{\left<\right.}}
\newcommand{\rb}{{\left.\right>}}
\newcommand{\lum}{\,{\rm erg\,s^{-1}}}
\newcommand{\kms}{km s$^{-1}$}
\newcommand{\surf}{M$_\odot$ pc$^{-2}$}
\newcommand{\ssfr}{M$_\odot$ yr${-1}$ kpc$^{-2}$}
\newcommand{\msunyr}{M$_\odot$ yr$^{-1}$}
\newcommand{\tell}{\tilde{\ell}}
\newcommand{\parall}{\mathbin{\!/\mkern-5mu/\!}}

%\graphicspath{{../}}

\newcommand{\apj}{ApJ}
\newcommand{\apjl}{ApJ Letters}
\newcommand{\mnras}{MNRAS}
\newcommand{\aap}{A\&A}
\newcommand{\prd}{Phys.\ Rev.\ D}
\newcommand{\physrep}{Physics Reports}
\newcommand{\jcap}{JCAP}
\newcommand{\procspie}{Proc. of the SPIE}
\title{Multiple Lensing of the Cosmic Microwave Background anisotropies}

% more complex case: 4 authors, 3 institutions, 2 footnotes
\author[a]{M. Calabrese,}
\author[b,e]{C. Carbone,}
\author[a,f]{G. Fabbian,}
\author[c,d,e]{M. Baldi,}
\author[a,f]{C. Baccigalupi}

\affiliation[a]{SISSA, Via Bonomea 256, 34136, Trieste (TS), Italy}
\affiliation[b]{I.N.A.F, Osservatorio Astronomico di Brera, Via Bianchi 46, 23807, Merate (MI), Italy}
\affiliation[c]{Dipartimento di Fisica e Astronomia, Alma Mater
  Studiorum Universita' di Bologna, viale Berti Pichat 6/2, I-40127, Bologna (BO), Italy}
\affiliation[d]{I.N.A.F, Osservatorio Astronomico di Bologna, via Ranzani 1, I-40127
Bologna (BO), Italy}
\affiliation[e]{I.N.F.N. - Sezione di Bologna, viale Berti Pichat 6/2,
  I-40127, Bologna (BO), Italy}
\affiliation[f]{INFN, Sezione di Trieste, Via Valerio 2, I-34127 Trieste, Italy}

% e-mail addresses: one for each author, in the same order as the authors
\emailAdd{mcalabre@sissa.it}
\emailAdd{carmelita.carbone@brera.inaf.it}
\emailAdd{gfabbian@sissa.it}
\emailAdd{marco.baldi5@unibo.it}
\emailAdd{bacci@sissa.it}

\abstract{
%  Context: 
%The lensing of the Cosmic Microwave Background (CMB) is a
%  well-observed physical phenomena and it is now established as one of
%  the most promising and useful probes to understand how structures form and evolve in the Universe.
%  Aims: 
We study the gravitational lensing effect on the Cosmic
Microwave Background (CMB) anisotropies performing a ray-tracing of
the primordial CMB photons through intervening large-scale structures
(LSS) distribution predicted by N-Body numerical simulations with a
particular focus on the precise recovery of the lens-induced polarized
counterpart of the source plane.
%  Methods: 
%{\bf Our algorithm implements a slicing of the LSS distribution as a function of redshift and a projection of the simulated matter distribution into a series of concentric mass shells acting as lensing plane mimicking then the 3D distribution of structures in the Universe} (forse da togliere). 
We apply both a multiple plane ray-tracing and an effective deflection
approach based on the Born approximation to deflect the CMB photons
trajectories through the simulated lightcone. 
% Results: 
We discuss the results obtained with both these methods together with
the impact of LSS
non-linear evolution on the CMB temperature and polarization power
spectra.  We compare our results with semi-analytical approximations implemented in Boltzmann codes like, e.g., CAMB.
We show that, with our current N-body setup, the predicted lensing
power is recovered with good accuracy in a wide range of multipoles
while excess power with respect to semi-analytic prescriptions is observed in the lensing potential
on scales $\ell \gtrsim 3000$. We quantify the impact
of the numerical effects connected to the resolution in the N-Body
simulation together with the resolution and band-limit chosen to
synthesise the CMB source plane. We found these quantities to be
particularly important for the simulation of B-mode polarization power spectrum.
}
 
\keywords{CMB - gravitational lensing - cosmology: theory - methods:
  numerical }

\begin{document}
\maketitle
\flushbottom

\input{./intro.tex}
\input{./theory.tex}
\input{./algorithm.tex}
\input{./tests.tex}
\input{./results.tex}
\input{./conclusions.tex}

\acknowledgments 

We acknowledge the use of the publicly available Code for Anisotropies
in the Microwave Background (\lstinline!CAMB!), \lstinline!LensPix! and
the Hierarchical Equalised Latitude Pixel spherical pixelization
scheme (\lstinline!HEALPix!).
The L-CoDECs simulations were carried out on the Power6 cluster at the RZG computing centre in Garching and on the SP6 machine at the ``Centro Interuniversitario del Nord-Est per il
Calcolo Elettronico'' (CINECA, Bologna).
The ray-tracing computations have been performed on the IBM Fermi
cluster at CINECA, with CPU time assigned
under several CINECA class-C calls,  and at the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.
Part of this work was supported by the InDark INFN Grant. CC
acknowledges financial support to the ``INAF Fellowships Programme
2010''
and to the European Research Council through the Darklight Advanced
Research Grant (\# 291521).
MC acknowledges the
financial support through the INAF funds for the Euclid mission. MB is
supported by the  Marie Curie Intra European Fellowship
``SIDUN''  within the 7th Framework  Programme of the European Commission.

\appendix

%% \section{Remapping algorithm}
%% \begin{figure}
%% \centering
%% \includegraphics[scale=0.45]{Immagini/spherical_change.eps}
%% \caption{
%% Geometry illustrating how the point A is re-mapped in B.
%%  }\label{fig:remapping}
%% \end{figure}
%% In Fig.~\ref{fig:remapping}, let the initial and final position of the ray in question
%% be the points $A \equiv (\theta,\phi)$ and $B \equiv (\theta ', \phi +
%% \Delta \phi)$, respectively. The North pole of the sphere is indicated
%% as C, so that the dihedral angle at A is also the angle
%% between $\tilde{\boldsymbol\alpha}$ and $- {\bf e}_{\theta}$, so that
%% \be
%% \tilde{\boldsymbol\alpha} = -\tilde{\alpha} \cos \delta {\bf e}_{\theta} +
%% \tilde{\alpha} \sin \delta {\bf e}_{\phi}.
%% \ee
%% Now, applying the spherical cosine rule to the triangle
%% ABC, we have
%% \be
%% \cos \theta ' = \cos \theta \cos \tilde{\alpha} + \sin \theta \sin
%% \tilde{\alpha} \cos \delta,
%% \ee
%% and applying the sine rule
%% \be
%% \sin \Delta \phi = \sin \tilde{\alpha} \frac{\sin \delta}{\sin \theta'}
%% \ee
%% We use these equations to remap points on the CMB sky
%% and on the intermediate spherical shells in the multiple
%% plane case, as described below.
\section{Measuring angular power spectrum}\label{appendix:ps}
The lensing potential is extracted from N-body simulation through a
binned map-making procedure of particles contained in a sky pixel. For
this reason when we want to extract its underlying power spectrum we
have to correct for the pixel window function of the \lstinline!HEALPix!
grid. Assuming  an azimuthally symmetric patch, as it is the case for
the full sky, the pixel window function is azimuthally symmetric and
can be used to correct the pseudo power spectrum extracted from the
full sky map using a simple spherical harmonic analysis operation.
A continuous field sampled on the \lstinline!HEALPix!
sphere is a smoothed version of the true underlying field due to the finite pixel size, 
i.e. the value of the field in pixel $i$ is given by
\be
\Phi^{pix(i)} = \int {\rm d}^2 \hat{{\bf n}}w^{(i)}(\hat{{\bf n}})\Phi(\hat{{\bf n}}),
\ee
where $w^{(i)}$ is the window function of the $i$-th pixel as is
given by
\be
w^{(i)}(\hat{{\bf n}}) =
\begin{cases}
\Omega_{pix}^{-1}, & \text{inside pixel $i$} \\
0 & \text{elsewhere}.
\end{cases} 
\ee
Expanding the true field $\Phi$ in terms of spherical harmonics
as
\be
\Phi(\hat{{\bf n}}) = \sum_{\ell m} \Phi_{\ell m} Y_{\ell m}(\hat{{\bf n}}),
\ee
we have
\be
\Phi^{pix(i)} = \sum_{\ell m}w_{\ell m}^{(i)} \Phi_{\ell m};
\label{eq::phi_pix}
\ee
where
\be
w^{(i)}_{\ell m} = \int d^2 \hat{{\bf n}} w^{(i)}(\hat{{\bf n}})
Y_{\ell m}(\hat{{\bf n}})
\label{eq::pixel_w}
\ee
is the spherical harmonic transform of the pixel window
function. The computation of these coefficients for each and every
pixel is required for a complete analysis in the \lstinline!HEALPix! scheme; however this calculation
becomes computationally unfeasible, even for a moderate \lstinline!NSIDE!. Therefore, it is
advantageous to ignore the azimuthal variation and rewrite
Equation~(\ref{eq::pixel_w}) as
\be
w^{(i)}_{\ell m} = w^{(i)}_{\ell} Y_{\ell m}(\hat{{\bf n}}_i),
\label{eq::47}
\ee
defining an azimuthally averaged window function as
\be
w^{(i)}_{\ell} = \frac{4 \pi}{(2\ell+1)} \left[ \sum^{\ell}_{m=-\ell} \left
  | w_{\ell m}
  \right |^2 \right]^{1/2}.
\ee
It follows immediately from Equations~(\ref{eq::phi_pix}) and (\ref{eq::47})  that
the estimate of the power spectrum of the pixelated field
is given by
\be
C_{\ell}^{\Phi,pix} = w^2_{\ell} \left \langle \Phi_{\ell m} \Phi^{*}_{\ell m} \right \rangle,
\ee
where the pixel averaged window function is defined as
\be
w_{\ell} = \left( \frac{1}{N_{pix}} \sum^{N_{pix}-1}_{i=0} \left(w^{(i)}_{\ell} \right)^2
\right)^{1/2}.
\ee
This function is available for $\ell < 4\times$\lstinline!NSIDE! in the
\lstinline!HEALPix! distribution. As we divide the computed power spectrum
by the square of the above function, it is possible to correct the effect of the
pixel window; in our case we act directly 
on the spherical harmonics coefficient $\Phi_{\ell m}$ recovered from our spherical maps,
using the actual $w_{\ell}$.

\bibliographystyle{JHEP}
\bibliography{multilens}

\end{document}

results.tex
 %!TEX root = jcap_article_MLcode.tex
 % !TEX spellcheck = en-EN
 
\section{Results}

\subsection{Numerical effects of the pixel-based methods}\label{sect:numerical-pixel-based}
Pixel-based methods for CMB lensing, though in general very efficient, are subjected to several numerical problems.
The first one is related to the bandwidths of the lensed CMB fields
generated as a result of the calculation. %, which are naturally not
%band-limited because lensing acts as a convolution in the harmonic
%domain \cite{Hu2000}.
Because the lensing effect happens before the intervention of any
instrumental response, the synthesis and analysis of relevant sky
signals in the pixel-based lensing methods (CMB source plane and
lensing potential map) require using a resolution sufficient to
support the signal up to the intrinsic bandwidth $\ell_{max}$ set by the user specific application and its required accuracy.
However, since mathematically the lensing effects act as a convolution in
the harmonic domain, the bandwidth of the resulting lensed field 
is broader than the one of the unlensed CMB and of the lensing
potential. Therefore, given the bandwidth used to synthesize the CMB source plane ($\ell_{max}^{CMB}$), i.e before undergoing any deflection, and the one used to solve the Poisson equation and to create the deflection field for a given shell ($\ell_{max}^{\psi}$), the
resulting lensed CMB will have an approximate band-limit of
$\ell_{max}^{CMB}+\ell_{max}^{\psi}$\footnote{We note that in general this bandlimit is only approximate and the resulting function is strictly not band-limited unlike the input source plane.}.
Consequently, the lensed map should
have its resolution appropriately chosen to eliminate potential power
aliasing effects on these angular scales \citep{Fabbian13}. We note
moreover that these aliasing effects are even more important
in the case of the ML approach because the bandwidth extension induced by
lensing happens each time the lensed CMB is propagated through a
single shell.\\*
The second challenge arises from the fact that
the displaced direction at each iteration ${\bf \hat{n}^{\prime}}{}^{(k)}$ does not correspond in general to the pixel
centers of any iso-latitudinal grid on the sphere. The values of the
CMB signal at those locations thus cannot be computed with the aid of fast
SHT algorithm and a more elaborated approach is needed.
In the context of pixel-based simulation methods, interpolation is the
most popular workaround of the need to directly calculate values of
the unlensed fields for every displaced directions %, which typically
%will not correspond to grid points of any iso-latitudinal
%pixelization chosen for the problem. 
The exact solution, which consist in a direct resummation of the
spherical harmonics at the displaced position, is in fact unfeasible even for moderate
resolutions \cite{Lewis05, Fabbian13}.
Any interpolation in this context, however, is not without its dangers
because interpolations tend to smooth the underlying signals and - as a
consequence - to hide aliasing effects in the lensed maps. For this
reason we chose the bandwidth of our signal
($\ell_{max}^{CMB}=\ell_{max}^{\psi}=4096$) and the resolution of our grid
(\lstinline!NSIDE!=4096) following the recipe provided in
\cite{Fabbian13} to minimize all of these effects simultaneously. This choice however limits the multipole range where the lensing signal can be reproduced with high reliability especially in the case of B-modes of polarization to $\ell\lesssim 2500$ (see Sec~\ref{sec:angular_power_spectrum}).

\subsection{Shot-noise contribution}
\label{sec:shotnoise}
In this section we estimate the impact of the intrinsic
discretisation of the N-Body simulation on the final lensed CMB power
spectra. Since we expect changes in the power spectrum on the order of few
percent, it is mandatory to be able to
control numerical artifacts with the same level of precision. 
% A good way to quantify the latter is to compute the shot-noise contribution to the lensing potential spectrum, propagate it to the lensed-CMB fields. 
%The idea is to use an analytical approach we tried to understand better the role of
%numerical effects on the recovered signal and isolate them in the
%context of full-sky lensing
%simulation of CMB. 
For this study we use the analytical modelling of weak lensing in the harmonic
domain discussed in \citep{Hu2000}. This treats lensing as an
effective convolution in Fourier space between the unlensed CMB and
the lensing potential and it is based on a second order Taylor expansion
around the unperturbed photons direction. 
The formulae are accurate to better than the percent level on the
angular scales considered in this work, especially in the case of
B-modes, and allow to quantify more easily the impact of the choice of
the band-limit on the recovered result. In the specific case of
B-modes, the convolution reads:
%on intermediate scales where we expect the assumption of the
%analytical approximation (second order Taylor expansion for the
%computation of lensed power spectra) to be less accurate. In the range
%$2000\lesssim \ell \lesssim 4000$, in fact, the expected displacement
%from cosmological weak lensing has the same magnitude of the angular
%scales in exam (few arc-minutes) and therefore a perturbative approach is less accurate.
%The analytic approximation defines the lensing modes denoted with a
%tilde ($\tilde{C}_{\ell}$) as a
%convolution between un-lensed modes and lensing potential:
%\begin{equation}
%\tilde{C}_{\tell_{T}}^{TT}=\left( 1 - (\tell^2_{T} + \tell_{T}) R
%\right)C^{TT}_{\tell_T} + \frac{1}{2}\sum_{\ell_{T},\ell_{\psi}}\frac{| _0
%F_{\tell_{T}\ell_{T}\ell_{\psi}}|^{2}}{(2\tell_{T}+1)}C^{TT}_{\ell_{X}}C_{\ell_{\psi}}^{\psi\psi}, 
%\label{eq:analyticalt}
%\end{equation}
%for the temperature, and
\begin{equation}
%\tilde{C}_{\tell_{X}}^{X}=\left( 1 - (\tell^2_{X} + \tell_{X} - 2) R
%\right)C^X_{\tell_X} + \frac{1}{2}\sum_{\ell_{X},\ell_{\psi}}\frac{| _2
%F_{\tell_{X}\ell_{X}\ell_{\psi}}|^{2}}{(2\tell_{X}+1)}C^{EE}_{\ell_{X}}C_{\ell_{\psi}}^{\psi\psi}, 
%-----> R factor \left( 1 - (\tell^2_{B} + \tell_{B} - 2) R\right)
\tilde{C}_{\tell_{B}}^{BB}=O(\tell_{B},C_{\ell}^{\psi\psi})\cdot C^{BB}_{\tell_B}
+ \frac{1}{2}\sum_{\ell_{E},\ell_{\psi}}\frac{| _2
F_{\tell_{B}\ell_{E}\ell_{\psi}}|^{2}}{(2\tell_{B}+1)}C_{\ell_{\psi}}^{\psi\psi}\left[C^{EE}_{\ell_{E}}(1-(-1)^{L}) +C^{BB}_{\ell_{E}}(1+(-1)^{L})\right] , 
\label{eq:analyticalb}
\end{equation}
%for the polarization E or B.
%where $\alpha$ is an integer which is different for each CMB spectra,
%$\alpha=2$ for X=E, $\alpha=0$ for X=T and $\alpha=1$ for X=TE. 
where we denote with tilde a lensed quantity and
$L=\tell_B+\ell_{E}+\ell_{\psi}$. We refer the reader to
\citep{Fabbian13} for a detailed discussion of the properties of the
convolution kernels $F_{\tell_{B}\ell_{E}\ell_{\psi}}$ and to the
$O(\tell_{B},C_{\ell}^{\psi\psi})$ factor. Similar expression can also
be derived for the TT, TE and EE power spectra \citep{Hu2000}.
Since this formalism does not make any assumption on the explicit form
nor the origin of $C_{\ell}^{\psi\psi}$, we can plug in
Eq.~(\ref{eq:analyticalb}) the shot-noise power spectrum instead of the
lensing potential extracted from the N-body simulation, to estimate the
fraction of the recovered signal generated by the limited resolution
of our simulated data.
%Note that since we do not assume any primordial B-mode ($C^{BB}_{\ell}=0$),
%the BB power is generated only by the convolution of E-mode and
%the lensing potential, i.e. it is totally due to the lensing
%process. For the lensing potential $\psi$ alone, the $\tell$
%and standard $\ell$ mode are equivalent. 
%In the first term of Eqs~(\ref{eq:analyticalt},\ref{eq:analyticalb}), the unlensed field is multiplied by a factor $R$, which reads:
%\be
%R = \frac{1}{2}\sum \ell (\ell +1) \frac{2\ell +1}{4 \pi} C_{\ell}^{\psi\psi}.
%\ee
%The coupling kernels
%$_{\alpha} F_{\ell\tell_{X}\tell_{\psi}}$ are defined as \cite{Hu2000}
%\begin{eqnarray}
%_{\alpha} F_{\ell_{1}\ell_{2}\ell_{3}}&=&\frac{1}{2}[\ell_{2}(\ell_{2}+1) + \ell_{3}(\ell_{3}+1)-\ell_1(\ell_1+1)]\sqrt{\frac{(2\ell_1+1)(2\ell_{2}+1)(2\ell_{3}+1)}{4\pi}}\nonumber\\
%&\cdot&\begin{pmatrix}\ell_1&\ell_{2}&\ell_{3}\\\alpha&0&-\alpha\end{pmatrix}\cdot(1-(-1)^{L}) \qquad\qquad \textrm{where}\  L=\ell_1+\ell_{2}+\ell_{3}.
%\end{eqnarray}
For this purpose we truncated the sum of Eq.~(\ref{eq:analyticalb}) to
the same band-limit value used in the lensing simulation,
i.e. $\ell^{E}_{max}=\ell^{\psi}_{max}=4096$. \\*
We first evaluate the shot-noise contribution to the $\psi^{eff}$
lensing potential starting from Eq.~(\ref{eq:clkappa_noise}) and
assuming the average number density to be the one of the
$\psi^{eff}$ field, $\bar{n}^{eff}_{k}$. This is then used as an input
for the analytical formulae of Eq.~(\ref{eq:analyticalb}), assuming the
primordial B-modes $C_{\ell}^{BB}$ to be zero as it is the case for
the unlensed CMB realizations used in the following. To validate the
analytical shot-noise predictions we also produce 100 Monte Carlo
realizations of shot-noise for the effective lensing
potential. We use those maps to extract a deflection field which is
then given as input to \lstinline!LensPix! to lens a random
realization of the CMB signal. The final average of all the power
spectra of these set of lensed CMB maps contains thus only the
lensing effect due to the shot-noise acting on primordial
anisotropies.\\*
To evaluate the shot-noise contribution to the ML method we compute 
Eq.~(\ref{eq:clkappa_noise}) for each $k$-shell and then apply the analytical convolution %\ref{eq:analyticalb}
iteratively assuming as input CMB spectrum for the $k$-th shell the lensed one
emerging from lensing of the previous $(k-1)$-th iteration. \\*
As discussed in the
following section, if we assume a power cut-off for the incoming CMB
and the lensing potential to be $\ell_{max}^{CMB}$ and
$\ell_{max}^{\psi}$ respectively, the lensed CMB after each deflection shows
power up to a multipole $\tell \approx  \ell_{max}^{CMB}+\ell_{max}^{\psi}$,
due to the properties of the lensing convolution kernels in the harmonic domain
\citep{Fabbian13}. The evaluation of the lensing kernels
requires in fact a computationally-heavy summation of Wigner-$3j$
symbols up to high multipoles. We therefore have assumed that for the
iteration $k>0$ the incoming CMB has power at most up to
$\ell=8192$. This additional cutoff is high enough not to affect
significantly the results on the scale considered in this work. The
analytical formulae were validated also in this case with Monte
Carlo simulation where for each shell the noise realizations were
generated starting from the shot-noise power spectrum of the single
shell.\\*
In Figure~\ref{fig:shot_noise_a} we show the results obtained from
both these methods, for the B-mode power spectrum, which is the most sensitive to
the details of the lensing potential being entirely lens-induced in
our case (no primordial tensor modes). The TT and EE spectra are 
conversely quite insensitive to the
the shot-noise which impacts the results at the sub $0.1\%$ level
(see Figure~\ref{fig:tt_ee_lensed}).
Both the analytical and Monte Carlo estimates of the shot-noise
contribution in the Born approximation agree extremely well
at all angular scales. The shot-noise contribution in the ML approach
is comparable to the effective case at $\ell\lesssim 2000$ though the difference is less than $0.5\%$. % and slightly more important when going towards smaller angular scales.
%However, a full validation of the result with Monte Carlo
%simulations is too expensive from the numerical level as it
%would require to rerun multiple N-body simulation with random
%realization of the initial condition starting from the same
%cosmology. 

\begin{figure}[htbp]
\centering
%\subfloat[][\emph{B-modes shot-noise power spectrum}.]
{\includegraphics[width=.55\textwidth,angle=90]{Shot_noise_comp.ps}} %\\
%\subfloat[][\emph{Analytical Born and ML comparison}.]
%{\includegraphics[width=.55\textwidth,angle=90]{Immagini/Definitive_Demnuni/Shot_noise_comp_ratio.ps}} 
\caption{Angular power spectrum for the lensed B-modes induced by the
  simulation shot-noise. The red, solid line is computed
  using 100 realizations of the algorithm in the Born approximation
  scenario. The green-dotted and the blue-dashed lines are evaluated
  with the analytical formula for - respectively - the Born and the
  multiple plane scenario. 
}
\label{fig:shot_noise_a} 
\end{figure}
%\clearpage
\subsection{Angular power spectrum}
\label{sec:angular_power_spectrum}
%In this section we analyze the angular power spectra of the lensed
%CMB fields as obtained with our algorithm. 
Similarly to the case of the lensing potential extraction,
we now take into consideration two different approaches also 
for the evaluation of the angular power spectrum of the lensed CMB. 
The first set of primary CMB maps are lensed in the Born
approximation,  while the second set by mean of the ML approach.
In Figure~\ref{fig:tt_ee_lensed} we show the comparison between the expected CMB lensed temperature and the E-modes of polarization power spectra, ($C^{TT}_{\ell}$ and
$C^{EE}_{\ell}$), estimated using semi-analytical halo
mass function implemented in \lstinline!CAMB! \cite{Smith03, Takahashi12}, and the spectra
extracted from our lensed CMB maps. For both these cases the simulated power spectra
follows precisely the \lstinline!CAMB! signal. In particular, the
shot-noise-induced contribution (evaluated following the recipe of the
Section~\ref{sec:shotnoise}) for these two observables is negligible
given that the effect of lensing per-se is already minor. Thus,
changes introduced at percent variation in the lensing potential are
further mitigated and hidden in the numerical noise. After having
subtracted the shot-noise induced lensing contribution, the fractional
difference between \lstinline!CAMB! and the N-Body lensed spectra
shows no significant bias up to $\ell \approx 3000$ where we start
seeing effects due to the choice of $\ell_{max}^{CMB}$. The latter is
not high enough to properly resolve power on those scales with
high-accuracy. The difference between the results obtained with the
Born and ML method is negligible and important only towards scales
$\ell\approx\ell_{max}$ (see Figure~\ref{fig:tt_ee_lensed_ratio}). The
abrupt decrease in power observed on those scales for the ML approach
with respect to the Born approximation is due to the effect of polynomial
interpolation. As the latter tends to smooth the underlying signal,
the consecutive application effectively removes more power on small
angular scale with respect to the Born approach, for which the
interpolation is performed only once.
 %we found a slightly excess of power on all scales, consistent with the higher power recovered in the lensing potential angular power spectrum (see Figure~\ref{fig:phi_camb_subfig}).

\begin{figure}
\centering
%\subfloat[][\emph{TT angular power spectrum}.]
%{\includegraphics[width=.55\textwidth,angle=90]{Immagini/Definitive_Demnuni/NoiseSub_cmb_TT_power_spectum.ps}} \\
{\includegraphics[width=.55\textwidth,angle=90]{Referee_CoDECS_NoiseSub_TT_power_spectum.ps}} \\
%\subfloat[][\emph{EE angular power spectrum}.]
%{\includegraphics[width=.55\textwidth,angle=90]{Immagini/Definitive_Demnuni/NoiseSub_cmb_EE_power_spectum.ps}} 
{\includegraphics[width=.55\textwidth,angle=90]{Referee_CoDECS_NoiseSub_EE_power_spectum.ps}} 
\caption{
Angular power spectrum of the (lensed) total intensity T (top panel) and 
the E-mode polarization of the CMB (bottom panel). Black dashed lines are \lstinline!CAMB! realization of a
lensed spectra. Red dotted
line uses the lensing field as in the Born
approximation, while for the blue solid line the CMB is lensed through
multiple planes. Green lines in both panels show the shot-noise
contribution to the lensed TT and EE spectra; green,
dot-dashed lines represent the absolute value of this
contribution. Note that the shot-noise
power spectra is multiplied by a factor of 50 for visualization
purposes. }
\label{fig:tt_ee_lensed} 
\end{figure}
\begin{figure}
\centering
%\includegraphics[scale=0.7,angle=90]{Immagini/Definitive_Demnuni/NoiseSub_cmb_TT_EE_ratio.ps}
\includegraphics[scale=0.6,angle=90]{Referee_CoDECS_NoiseSub_TT_EE_ratio.ps}
\caption{
Fractional difference for the angular power spectrum of the temperature
(top panel) and E-mode polarization (bottom panel) with respect to
\lstinline!CAMB!. Red dotted
lines are obtained in the Born
approximation, while blue solid in the multiple lens plane
approach. The shot-noise has been subtracted in both cases.
}
\label{fig:tt_ee_lensed_ratio}
\end{figure}

The situation however is different for the B-modes of polarization, as
shown in  Figure~\ref{fig:bb_lensed_noise}. This signal is 
entirely caused by lensing as we have set the primordial tensor
modes to zero. Its behaviour is therefore a clear imprint of how the LSS
process the primary CMB field and thus we expect this observable to reflect more directly the features observed in the lensing potential. As expected from the analysis of the lensing
potential in Section~\ref{sec:lensing_pot}, the BB spectrum shows
a lack of of power at percent level with respect to \lstinline!CAMB! spectra,
in agreement with the matter power spectrum of the N-Body
simulation, though this effect is partially compensated by the increase in power at small scales in the lensing potential. This feature is not observable in the lensed T or E field, where power coming from primordial
anisotropies is dominant over the lensing-induced one. Moreover, 
while negligible in the TT and EE cases, we found the shot-noise contribution to be
important at the percent level for the BB power spectrum at small scales. This is
expected given that B-modes are very sensitive to non-linear power,
which is affected by shot-noise for $\ell >> 2000$ (see bottom panel
of Figure~\ref{fig:phi_camb_subfig}). The lack of power due to the
choice of $\ell_{max}^{CMB}$ starts to be important on angular scales
larger than the ones affected in T and E-modes power spectra. This can
be explained considering that at those scales a non-negligible
fraction of the contribution to the BB power spectrum starts to come
from progressively higher multipoles of both E and $\psi$. At $\ell_B
\sim 4000$, for instance, a 25\% contribution to the power in the
B-modes comes from scales in the E and $\psi$ fields at $\ell >
\ell_{max}=4096$ \citep{Fabbian13}. Since our algorithm is
band-limited to this $\ell_{max}$,  cutting power for those high
$\ell$, produce a loss of about 25\% in the BB power spectrum at that
particular multipole (as shown in Figure~\ref{fig:bb_lensed_noise}). 
%Finally, the gray area gives an idea of the
%uncertainty in the spectrum, estimated with the usual formula:
%\be
%\Delta C^{BB}_{\ell} = \sqrt{\frac{2}{2\ell+1}}\left(C^{BB}_{\ell} +
%N^{BB}_{\ell} \right),
%\ee
%where the $N^{BB}_{\ell}$ is the shot-noise-induced power spectrum for
%the B-modes. 
\begin{figure}
\centering
\includegraphics[scale=0.7,angle=90]{Referee_BB_noisesub_ps.ps}
\caption{Angular power spectrum for the B-mode polarization of CMB. Black dashed lines are
  \lstinline!CAMB! realization of a lensed spectra. Red dotted lines
  are obtained in the Born
approximation, while blue solid lines in the multiple lens
planes approach. The green dot-dashed line is
a lensed spectrum produced from a shot-noise-only lensing map. Note that the shot-noise
power spectra is multiplied by a factor of 10 for visualization
purposes. 
In the bottom panel, it is shown the fractional difference
with the reference \lstinline!CAMB! spectrum. In this panel, the noise power spectrum
  has been subtracted from the original signal.
}
\label{fig:bb_lensed_noise}
\end{figure}

%\begin{figure}
%\centering
%\includegraphics[scale=0.7,angle=90]{Immagini/Definitive_Demnuni/BB_integral.ps}
%\caption{
%Angular power spectrum for the BB modes. Black lines are \lstinline!CAMB!
%realization of lensed spectra. Crosses show the power spectra as
%computed lensing the CMB through multiple lensing planes, for different
%integration ranges in redshift. 
%}
%\label{fig:bb_lensed}
%\end{figure}
%
%\begin{figure}
%\centering
%\subfloat[][\emph{BB modes ($\ell < 2000$)}.]
%{\includegraphics[width=.45\textwidth]{Immagini/Giulio/bmodes_analytical.ps}} \quad
%\subfloat[][\emph{BB modes ($2000 < \ell < 4000$)}.]
%{\includegraphics[width=.45\textwidth]{Immagini/Giulio/bmodes_zoom.ps}} \\
%\caption{
%Analytical BB modes. }
%\label{fig:BBmodes_analytical}
%\end{figure}
%
%\begin{figure}
%\centering
%\subfloat[][\emph{BB modes ($\ell < 2000$)}.]
%{\includegraphics[width=.48\textwidth]{Immagini/Giulio/bmodes_noLS.ps}} \quad
%\subfloat[][\emph{BB modes ($2000 < \ell < 4000$)}.]
%{\includegraphics[width=.48\textwidth]{Immagini/Giulio/bmodes_noLS_zoom.ps}} \\
%\caption{
%Analytical BB modes with large scale power (for $\ell < 70$) set
%artificially to 0.}
%\label{fig:BBmodes_noLS}
%\end{figure}
%
%\begin{figure}
%\centering
%\includegraphics[scale=0.7,angle=+90]{Immagini/Definitive_Demnuni/BB_resolution_camb.ps}
%\caption{
%}
%\label{fig:bb_lensed_zoom}
%\end{figure} 
As argued in Section~\ref{sect:numerical-pixel-based}, one of the major numerical problem affecting the simulation of CMB B-modes is the power aliasing due to bandwidth extension induced by lensing.  In Figure~\ref{fig:resolution} we show the impact of this effect as a function of the map
resolution on the B-modes power spectrum recovery. For this tests we
extracted the lensing potential maps for both Born and ML approach
using two different HEALPix grid at $\textrm{\lstinline!NSIDE!}=2048,
4096$ and refer to these two setup as the low and high resolution case
respectively. We then synthesized on the same grid the CMB source
plane assuming the same bandlimit $\ell^{CMB}_{max}= 4096$, as done
for the results discussed above, and propagate it through the lensing
plane(s). As shown in Figure~\ref{fig:resolution}, the Born
approximation method is quite insensitive to the choice of
\lstinline!NSIDE! because the polynomial interpolation is effective in
removing most of the aliasing. However, for the ML scenario the
situation is worse as the aliasing generated by multiple deflection
can add up, becoming progressively more important.
%Comparing for example the Born approximation lines in both cases, we see that at higher
%resolution the power spectrum agrees better with \lstinline!CAMB!;
%basically with $\ell_{max}=$ \lstinline!NSIDE! the algorithm can efficiently
%resolve the simulated sky power all the way to its actual
%bandwidth. 
%In the low resolution case, instead, the signal suffers
%from residual aliasing which results in an higher power spectrum at
%small scales. 
This can then lead to a misinterpretation of the result obtained using
the ML, which seems to be significantly different from the once
obtained in the Born approximation. The fact that this difference
vanishes in the high-resolution case is a demonstration of the high
level of control of numerical effects which needs to be achieved for
this kind of algorithms. Even though these effects were limited in the
setup considered here, we expect those to become more important when
targeting accurate lensing simulations on scales $\ell>>2000$.
% In a future work we will fix this numerical
%aspect in order to better
%control aliasing, such that the need of keeping 
%$\ell_{max}=$ \lstinline!NSIDE! fixed can be dropped.
\begin{figure}
\centering
\includegraphics[scale=0.6,angle=+90]{resolution_effects.ps}
\caption{
Fractional difference of
angular power spectrum of the BB power spectrum with respect to the reference \lstinline!CAMB! spectrum. Dotted lines refers to the map at low resolution (\lstinline!NSIDE!=2048), while solid lines
to the map with \lstinline!NSIDE!=4096. Red lines plot the effective,
Born approximation case, blue lines are connected to the multiple
plane approach. Note that for this
comparison, the noise power spectrum
has not been subtracted from the original signal.
}
\label{fig:resolution}
\end{figure}
%
\\

Finally, we compare the differences in the angular power spectra between the Born approximation and the
multiple planes approach. First we define the quantity $O^{X}_{\ell}$
as the difference between the angular power spectra extracted with the
multiple lens approach and the one computed in the Born approximation, 
\be
O^{X}_{\ell} = C^{X, {\rm ML}}_{\ell} - C^{X, {\rm Born}}_{\ell}, 
\ee
where X = TT, EE, BB. Its uncertainty is given by the cosmic variance affecting both
spectra, or
\be
\sigma_{O^X_{\ell}} = \sqrt{\frac{2}{2\ell+1}\left( | C^{X, {\rm
      ML}}_{\ell}|^2 + |C^{X, {\rm Born}}_{\ell}|^2 \right)}.
\ee
Starting from these quantities we can define a reduced
chi-square $\tilde{\chi}^2$ statistics
\be
\tilde{\chi}^2 = \frac{1}{\ell_{max}-1} \sum^{\ell_{max}}_{\ell = 2}
\frac{O^2_{\ell}}{\sigma^2_{O_{\ell}}}, 
\label{eq:chi2-error}
\ee
to assess whether the two methods are inconsistent.
Since we expect $O^X_{\ell}/ \sigma^{2}_{{O_{\ell}}}$ to be a Gaussian
random variable,
we can we also perform a Kolmogorov-Smirnov (KS) to test
whether this hypothesis is verified or systematic differences exists
between the two methods. In defining both theses tests and the sample
variance of Eq.~(\ref{eq:chi2-error}), we assumed that the covariance
of the lensed power spectra is Gaussian. This assumption neglects the
fact that lensing introduces non-Gaussian correlations between
different modes \cite{benoit-levy12, KSmith04}, but this effect is
mainly important for B-modes, for which the gaussianity assumption
underestimates the sample variance.
 
In Table~\ref{table:tab1} we report the results of both those tests
expressed as the significance level probability.
In both cases we find that the power spectra obtained in the Born
approximation and with the ML method are statistically consistent. A
further possible test to compare the two methods would be to
reconstruct the effective integrated matter density from the simulated
lensed CMB maps as done in \cite{antolini}, but 
we leave this option to future work. 
\begin{table}
\centering
\begin{tabular}{lcc}
\toprule
$C^{X}_{\ell}$ & Significance $P_{KS}$ & Significance $P_{\tilde{\chi}^2}$\\
\midrule
TT & 0.47 & 0.70\\
EE & 0.19 & 0.51\\
BB & 0.21 & 0.19\\
%TE & 0.004 & 0.98\\
\bottomrule
\end{tabular}
\caption{\label{table:tab1} Results of statistical tests on difference
  between lensed CMB angular power spectra in the Born approximation
  and multiple lens planes approach. The significance level
  probability for the null hypothesis using a Kolmogorov-Smirnov test
  ($P_{KS}$) and a reduced chi-square $\tilde{\chi}^2$ statistics
  ($P_{\tilde{\chi}^2})$ show no difference between the power
  spectra computed with the two methods on a statistical level.
}
\end{table}


% {\bf GF: to remove and use a chi2 comparison with uncertainty given by the cosmic variance?}
% Finally, we compute the difference of the BB power spectrum to compare the results
% between the Born approximation and the multiple planes lensing:
% \be
% \Delta C^{BB}_{\ell} = C_{\ell}^{BB,ML} - C_{\ell}^{BB,Born}.
% \ee
% Figure~\ref{fig:tt_born_noborn} displays the difference relative to
% the Born spectrum, or $\Delta C^{BB}_{\ell} / C_{\ell}^{Born}$. In order to
% facilitate the comparison, we bin the raw power spectrum in $\ell$,
% and for each bin we compute the mean and the sample variance of the
% observations falling inside that bin (which will give, respectively,
% the value and the error bars shown in the figure). As we can see
% clearly from the Figure, we do not see any significant differences
% between the two approaches in the considered range of multipoles.
% \begin{figure}
% \centering
% \includegraphics[scale=0.6,angle=90]{lmax4096_CMB_BB_born_mlplane.ps}
% \caption{
% Relative fractional difference between the angular power spectrum of
% lensed B-modes computed in Born approximation and the CMB lensed
% through multiple planes. Similar results are obtained for the TT and
% EE angular power spectra. Error bars are given by the cosmic
% variance in each bin.
% }
% \label{fig:tt_born_noborn}
% \end{figure}



tests.tex
%!TEX root = jcap_article_MLcode.tex
\section{Test and Convergence}
\label{sec:test}

\subsection{Sanity Checks}
 
In this section we assess the reliability of our code by performing
 sanity checks similarly to \cite{DasBode} to ensure
that all the steps of the algorithm give stable and robust results. 
For the first test, we show that the total mass selected in each 3D
matter slice is equal to the theoretical
mass expected from the assumed cosmological model in the simulation, given by
\begin{equation}
M_{slice}^{theory} = 4 \pi \Omega_{m,0}  \rho_c \bar{\eta}^2 \Delta
\eta, 
\end{equation}
where $\Delta \eta$ is the comoving thickness of the slice at a comoving
(average) distance $\eta$, $\rho_c$ is the critical density and
$\Omega_{m,0}$ corresponds to the present
value of the matter density parameter. We compare this quantity with
the total mass obtained from the surface density maps
($\Sigma^{\theta}$) drawn with our procedure,
\begin{equation}
M_{slice} = \sum^{N_{pix}}_{p=1} \Sigma^{\theta}_p \: \Omega_{pix}, 
\end{equation}
by summing on all pixels of the spherical map.
\begin{figure}
\centering
\includegraphics[scale=0.5,angle=-90]{mass_diff.ps}
\caption{
The total mass for each shell compared with the one expected from
theory (fractional difference) as function of redshift.
 }\label{fig::mass_diff}
\end{figure}
Figure~\ref{fig::mass_diff} shows the fractional difference between the
two masses for the different redshifts at which each
spherical map is located. The agreement is very good within a few percents. As similarly found
by \cite{DasBode}, fluctuations respect to the theory appear at low $z$,
due to the tension between the small comoving volume as seen by the
observer, and a highly-clustered matter distribution at late times. Including or excluding large
dark matter halos in the selection process therefore leads to
differences between the mass extracted from the maps and the theoretical one. \\*
%
As a second test, we make sure that the projection from the simulation box onto the map
has been properly performed. Figure~\ref{fig:PDF} displays the probability density function
(PDF) as recovered from the surface mass density maps, compared
with analytic PDFs, drawn from the data, such as the Gaussian and
log-normal ones (as in \cite{Kayo01,Taruya02}).
%% \begin{align}
%% PDF(x) &= \frac{1}{\sigma\sqrt{2 \pi}} \exp \left( -\frac{(x - \mu)^2}{2
%%     \sigma^2} \right) \: \: \text{Gaussian}, \\
%% PDF(x) &= \frac{1}{\sigma x \sqrt{2 \pi}} \exp \left(  -\frac{\left(\ln{x} - \mu\right)^2}{2
%%     \sigma^2} \right) \: \: \text{Log-Normal}.
%% \end{align}
%PDF(x) &= \frac{N}{x} \exp \left( -\frac{(\ln{x} + \omega^2/2)^2
%    (1+A/x)}{2 \omega^2} \right) \: \: \text{Das \& Ostriker}
The extracted PDFs are quite similar to the ones found by \cite{DasBode},
even if - as already observed by the same authors - the analytical PDFs could not fit well
the data especially at high surface mass density where the
non-linearities becomes important and  where accurate models are yet
unknown.

%%  Our results fits the model proposed model proposed by
%% \cite{DasOstriker}, and in a good agreement
%% (regarding both the shape of the distribution and values of surface mass density)
%% with previous authors \cite{DasBode}.
\begin{figure}[!htb]
\centering
\includegraphics[scale=0.6,angle=-90]{SurfaceMass_hist.ps}
\caption{
The probability density function (PDF) of the surface
mass density in the lensing-planes (crosses) compared with the
Gaussian (red, dotted line) and the log-normal
(blue, dashed line).
 }\label{fig:PDF}
\end{figure}

\subsection{Lensing potential maps}
\label{sec:lensing_pot}
Once surface density maps have been validated, we can move one step
forward and verify lensing quantities.
As described in Section~\ref{sec:discr}, the effective convergence plane
is computed by weighting the surface
mass density planes with appropriate
geometrical factors. We validated such convergence maps by
comparing the extracted power
spectra to the theoretical expectations based on semi-analytical
computations of the matter perturbation evolution as implemented
in the publicly available Boltzmann code \lstinline!CAMB!\footnote{\url{http://camb.info}}.
%and by making a consistency
%test for the angular power spectrum.
Adopting the Born approximation, we drew an ``effective'' convergence map, as described in
Eq.~(\ref{eq:born_kappa}), using the matter shells at different
redshifts. We then compute, in Limber approximation of
Eq.~(\ref{eq:clkappa}), the theoretical convergence angular power
spectrum, exploiting directly the 3D matter power spectrum computed
with \lstinline!CAMB!. The comparison between the simulated and the
analytical power spectra is shown in Figure \ref{fig:born_kappa}.
In both cases we perform the integration up to a specific
redshift $z^{*}$ to assess the validity of the maps at different
times. We observe that the measured spectra agree at high accuracy with the
theoretical predictions on a large interval of angular scales,
indicating the validity of our map-making procedure. As expected, the
lack of power at small multipoles $\ell\lesssim 50$ is due to the finite box size of the
N-Body simulation. 
A source of possible contamination of
the signal is the so-called {\it shot-noise}, due to the finite
particle density in the N-Body simulation. The shot-noise power spectrum can be computed analytically substituting the shot noise power spectrum $P^{Shot}(k, z)
= 1/\bar{n}_k$ in Eq.~(\ref{eq:clkappa_dis}),  where $\bar{n}_k = N_k
/\left(4\pi \eta^2_k \Delta \eta \right)$ and $N_k$ is the total
number of particles in the $k$-th shell, we obtain the shot-noise
contribution to the convergence:
\be
\label{eq:clkappa_noise}
C_{\ell}^{\kappa\kappa,{\rm Shot}} = \frac{9 H^4_0 \Omega^2_{m,0}}{4 c^4} \sum_k
  \Delta \eta \:
  \frac{1}{\bar{n}_k} \left( \frac{D_A(\eta_s - \eta_k)}{D_A(\eta_s)
    a_k} \right)^2.
\ee
For the N-body simulation used in this code the shot-noise is small at all redshifts given the high spatial resolution and high number of particles employed (see Figure~\ref{fig:born_kappa}).
%In Figure~\ref{fig:born_kappa}, straight dashed lines show the
%shot-noise contribution to the convergence field at different times. 
\begin{figure}[!htb]
\centering
\includegraphics[scale=0.6,angle=90]{born_kappa.ps}
\caption{
Angular power spectrum from the ``effective'' convergence map
(solid lines) 
compared with theory predictions (dotted line) at three
different redshifts, $z^{*} \simeq 0.6$, 1.3, 4. Dashed lines quantify
the shot-noise contribution for each maps. Note that the shot-noise
power spectra is multiplied by a factor of 10 for visualization
purposes. 
}
\label{fig:born_kappa}
\end{figure}
\noindent
Figure~\ref{fig:phi_camb} shows a comparison of the partial
contributions to lensing potential
angular power spectrum computed at different redshifts
 with the corresponding analytical signal given by
Eqs.~(\ref{eq:clkappa_dis}) and (\ref{eq:kappa2lens})
in which we insert the 3D matter power spectrum extracted from
\lstinline!CAMB!. In this case, the
label $z$ refers to the redshift of the matter spherical map which
contributes to the 
lensing potential at that time. Each power spectrum
represents the ``real'' map given as input to \lstinline!LensPix! in
order to obtain the final CMB lensed maps in the multiple plane approach.
Here the Limber approximation is necessary 
to solve the Poisson equation using the
 transverse part of the Laplacian only, thus neglecting line-of-sight contributions as
previously discussed in Sec.~\ref{sec:discr}.
%% The grey area around each 
%% theoretical line
%% represents the uncertainty in the power spectrum given by the cosmic
%% variance, i.e.
%% \be
%% \Delta C^{\psi\psi}_{\ell} = \sqrt{\frac{2}{2\ell +1}} C^{\psi\psi,th}_{\ell},
%% \ee
%% where $th$ stands for the theoretical power spectrum.
%\footnote{http://camb.info/}.
The agreement between simulated and
analytical $C^{\psi\psi}_{\ell}$ as a function of the redshift is
clearly observable from Figure~\ref{fig:phi_camb}.
The recovered signal is stable and coherent on
a whole range of multipoles. As discussed in the following, we assume
a very conservative choice for the map resolution and power cut-off
$\ell_{max}$ ($\textrm{\lstinline!NSIDE!}=\ell_{max} =
4096$). Therefore, we do
not expect this result to be affected by power aliasing
given that an \lstinline!HEALPix! grid with resolution set by \lstinline!NSIDE!
parameter should be able to properly sample modes up to $\ell\approx
2\times\textrm{\lstinline!NSIDE!}$.  
%In general, our spectra follow the halofit prescription, even
%if we see an excess of power of about 2\%
%with respect to the \lstinline!CAMB! ones. The origin of this is in the
%% simulation itself: the code was run without the energy contribution of
%% the radiation component ($\Omega_{rad} = 0$), which results in an
%% excess at all scales in the matter power spectrum $P(k,z)$. In
%% Figure~\ref{fig:pk} we compare the $P(k)$ as computed by the
%% \lstinline!GADGET-2! code itself ($P(k)^{Gadget}$) and the analogue
%% \lstinline!CAMB! realization ($P(k)^{Camb}$). For different scales, the matter power
%% spectrum in the simulation is slightly higher than the semi-analytical one.
\begin{figure}[!htb]
\centering
\includegraphics[scale=0.6,angle=90]{lensingphi_various_z.ps}
\caption{
Comparison between the angular power spectrum of the lensing potential computed with
our algorithm (solid lines) and the
analytical results obtained using the Limber equation (dotted lines) at different redshift. The
spectra have been multiplied by a constant factor for displaying
purposes. The grey area displays the cosmic variance $1\sigma$
uncertainty for the theoretical spectra.}
\label{fig:phi_camb}
\end{figure}
\noindent
%% \begin{figure}
%% \centering
%% \includegraphics[scale=0.6,angle=90]{Analytical/NoiseImm/gadget_vs_camb_Feb90_plane.ps}
%% \caption{
%% Fractional difference $\Delta P(k,z) / P(k,z)$ between matter power
%% spectrum computed by \lstinline!GADGET! during the run of the
%% simulation ($P(k,z)^{Gadget}$) with the corresponding \lstinline!CAMB!
%% realization, at different redshifts (as done in Fig.~\ref{fig:phi_camb}). The top x-axis show the
%% corresponding range in multipole $\ell$ ($\ell = k \eta$ in Limber
%% approximation) for $\eta = \eta_{\star}$
%% (black) or $\eta = \eta(z=0.53)$ (red). }
%% \label{fig:pk}
%% \end{figure}
%s
An interesting and comprehensive way to see how the lensing process behaves at different scales is to look
at the integrated potential, as computed in Born approximation using
Eq.~(\ref{eq:born_phi_dis}). In Figure~\ref{fig:phi_camb_subfig} we show the angular power
spectrum for the effective lensing potential and its shot-noise
contribution, compared to the semi-analytical realizations by
\lstinline!CAMB!, where we fix the maximum
redshift of the integration, $z_{max}$, to be the same as the maximum
redshift used in our map-making procedure. Also in this case, we
find a very good agreement between the two methods, within the
1$\sigma$ uncertainty for the semi-analytical spectrum.
 As in Figures~\ref{fig:born_kappa} and ~\ref{fig:phi_camb} the spectrum
shows a lack of power due to the finite size of the simulation box for
$\ell < 50$. Note that the shot-noise contribution is negligible, as
we have multiplied it by a factor of 10 such that it could be compared with the
lensing potential signal.
In general, at intermediate scales our spectra show a small deficit of
power, within 3\% with respect to the HALOFIT prescription 
\cite{Takahashi12}, while at small scales, even after the shot-noise subtraction
(blue lines), the signal seems to increase towards $\ell\approx 3500$,
likely due to the underlying non-linear clustering underestimated by
the analytical models. Since the simulated spectra agree within percent level with the semi-analytical
realization of \lstinline!CAMB!, this means that our map-making
procedure traces with good accuracy the evolving matter distribution.
\begin{figure}
\centering
%\subfloat[][\emph{Lensing Potential with \lstinline!CAMB! }.]
\includegraphics[width=.58\textwidth,angle=90]{int_lens_crop.ps}\\
%\includegraphics[width=.38\textwidth,angle=90]{int_ratio_all_crop.ps}\\
%\subfloat[][\emph{Fractional difference between \lstinline!CAMB! and our lensing
%    potential}.]
%{\includegraphics[width=.4\textwidth,angle=90]{Immagini/Definitive/int_ratio_all_crop.ps}} \\
%\subfloat[][\emph{Fractional difference between \lstinline!CAMB! and our lensing
%    potential for $1000 < \ell < 2000 $}.]
\includegraphics[width=.38\textwidth,angle=90]{int_ratio_up1000_crop.ps}
\includegraphics[width=.38\textwidth,angle=90]{int_ratio_from2000_crop.ps}
%\subfloat[][\emph{Fractional difference between \lstinline!CAMB! and our lensing
%    potential for $2000 < \ell < 4000 $}.]
\caption{Top panel: Integrated lensing potential $\psi^{eff}$ obtained
  in the Born
  approximation from the simulation (green) and the \lstinline!CAMB! result based on
  semi-analytical approximation of the non-linear evolution (red). The
  fractional difference between the two is displayed for some intermediate and high multipoles
  scales in the bottom panels. The results obtained after the shot-noise
  subtraction are displayed in blue solid lines. Note that the
  shot-noise spectrum in the top panel is multiplied by an arbitrary factor of 10
  to be seen clearly in the Figure.} 
\label{fig:phi_camb_subfig}
\end{figure}


theory.tex
%!TEX root = /Users/fabbian/Dropbox/ray tracing for cmb lensing/Articolo/giulio/jcap_article_MLcode.tex
% !TEX spellcheck = en-EN

\section{Weak lensing in Cosmology}\label{sec:theory}
\subsection{Gravitational light deflection}
In this Section, we briefly recall the definition of the relevant
quantities concerning the weak lensing effect that will be used in the
rest of the paper. We refer the reader to more specialised reviews for a general
treatment of the weak lensing effect \cite{Bartelmann10,
  LewisChallinor06}. Following the approach in \cite{DasBode}, we assume a coordinate system
%$(t, \boldsymbol\theta, \eta)$ 
based on physical time $t$, two
angular coordinates $\boldsymbol\theta=(\theta_1,\theta_2)$, and the line-of-sight, radial comoving
distance $\eta$ relative to the observer. In a standard Universe with
a weakly perturbed Friedmann-Lemaitre-Robertson-Walker (FLRW) metric,
a light ray approaching a matter density distribution is deviated by an angle 
\be
d \boldsymbol\alpha = -2\nabla_{\perp} \Psi d \eta,
\label{eq:diff_defl}
\ee
where $d \boldsymbol\alpha$ is called deflection angle, $\Psi$ is the Newtonian
potential and $\nabla_{\perp}$ denotes the spatial gradient on a plane
perpendicular to the light propagation direction. The gradient in Eq.~(\ref{eq:diff_defl}) is defined in
the small-angle limit as $\nabla_{\perp} = (\partial / \partial \theta_1, \partial / \partial \theta_2)$
where $(\theta_1, \theta_2)$ describe a coordinate
system orthogonal to the light ray trajectory. 
The transverse shift of the
light ray position at $\eta$ due to a deflection at $\eta'$ can be thus written as
\be
d{\bf x}(\eta) = D_A(\eta - \eta')d\boldsymbol\alpha(\eta'),
\ee
where $D_A(\eta)$ is the comoving angular diameter distance.
%, which is
%defined as
%\be
%D_A(\eta) = 
%\begin{cases}
%1/\sqrt{K} \sin \left( \sqrt{K} \eta \right), & \mbox{if } K > 0 \\
%\eta, & \mbox{if } K = 0 \\
%1/\sqrt{-K} \sinh \left( \sqrt{-K} \eta \right), & \mbox{if } K < 0 \\
%\end{cases} 
%\ee
%with $K$ being the curvature of space.
%
In weak lensing calculations, an ``effective''
approach is commonly used, where the effect of
deflectors along the entire line of sight is approximated by a projected
potential computed along a
fiducial un-deflected ray (Born approximation).
%Such approximation is often referred to as the Born approximation. 
The final angular position $ \boldsymbol\theta(\eta) = {\bf
  x}(\eta)/D_A(\eta)$
is therefore given by 
\begin{align}
\label{eq:bornlens}
\boldsymbol\theta(\eta) & = \boldsymbol\theta(0) - \frac{2}{D_A(\eta)}\int_0^{\eta}
{\rm d} \eta' D_A(\eta - \eta') \nabla_{\perp}\Psi =\\
\nonumber &= \boldsymbol\theta(0) + \tilde{\boldsymbol\alpha} (\eta),
\end{align}
where $\tilde{\boldsymbol\alpha}$ is the total effective deflection
from the initial angular position of the light ray at the observer position $\boldsymbol\theta(\eta= 0)$. Note that the integral in Eq.~(\ref{eq:bornlens}) is evaluated along the
light ray trajectory and is thus an implicit equation for $\boldsymbol\theta(\eta)$. 
It is convenient to define $\tilde{\boldsymbol\alpha}$ as the gradient of an effective
potential including the contributions to the final deviation of all the structures present between the observer and the background source plane located at a comoving distance $\eta_s$, i.e. ${\bf \tilde{\boldsymbol\alpha}} = -  \nabla_{\perp} \psi^{eff}$, where
\be
\label{eq:born_phi}
\psi^{\rm eff}(\boldsymbol\theta) = \frac{2}{D_A(\eta_s)}
  \int^{\eta_s}_{0} {\rm d} \eta
  \frac{D_A(\eta_s - \eta)}{D_A(\eta)} \Psi(\boldsymbol\theta ,\eta).
\ee
The latter quantity is commonly known as the lensing potential.
If we consider in particular the case of the weak lensing of CMB anisotropies, $\eta_s$  is the comoving distance to the last scattering surface.
%In terms of the projected potential, the effective deflection
%is given by the angular gradient, ${\bf \tilde{\boldsymbol\alpha}} = -
 %\nabla_{\perp} \psi^{eff}$ . 
 An effective convergence field can be also defined in a similar manner:
\begin{equation}
\kappa(\boldsymbol\theta) \equiv \frac{1}{2} \nabla^2_{\perp}
  \psi^{\rm eff}(\boldsymbol\theta) %\\
   = \frac{2}{D_A(\eta_s)} \int^{\eta_s}_{0} {\rm d} \eta
  \frac{D_A(\eta_s - \eta)}{D_A(\eta)} \nabla^2_{\perp}\Psi(\boldsymbol\theta ,\eta).
\label{eq:born_kappa}
\end{equation}

%Equations describing the evolution of the lensing Jacobian $A_{ij}$
%(also called the magnification matrix) along the light rays
%trajectory can be derived by simply taking another set of derivatives
%of Equation (\ref{eq:bornlens}) with respect to the location of the ray at the observer
%$\theta(\eta = 0)$ (Jain {\it et al.}, 2000 \cite{Jain00}; Hilbert
%{\it et al.}, 2009 \cite{Hilbert09}):
%\be
%A_{ij} = \delta_{ij} - \frac{2}{d_A(\eta)} \int^{\eta}_0 d\tilde{\eta}
%  d_A(\eta - \tilde{\eta}) \frac{\partial^2 \Psi(\tilde{\boldsymbol\theta},\tilde{\eta})}{\partial \theta_i \partial \theta_k}
%  A_{kj}(\tilde{\eta})
%\label{eq::Aijcon}
%\ee
%The
%  distortion matrix $A$, also called shear-magnification matrix or
%  Jacobian matrix, is typically decomposed into four parts describing how the
%  light rays are transformed by the passage through the matter
%  distribution
%\be
%A_{ij} = 
%\begin{pmatrix}
%  1 - \kappa -\gamma_1 & -\gamma_2 + \omega \\
%  -\gamma_2 - \omega & 1-\kappa+\gamma_1
% \end{pmatrix}
%\ee
%where we assumed that the rotation angle $\omega$ which defines a
%rotation of the lensed image, is small. The
%parameter $\kappa$ is the convergence while $\gamma = \gamma_1 + i
%\gamma_2$ defines the complex shear, describing the shearing of the
%image.

\subsection{Discretised formalism}
\label{sec:discr}
The previous set of equations can be discretised by dividing
the interval between the observer and the source
into $N$ concentric shells, each of comoving thickness $\Delta \eta$.
The matter in the $k$-th shell is projected onto a spherical,
two-dimensional sheet which is positioned in the middle of 
the two edges of the matter shell\footnote{The shell index
  $k$ increases as moving away from the source plane.}, at comoving
distance $\eta_k$. To simplify the following calculations, it is common practice to use angular differential
operators defined on the sphere instead of spatial ones, since we will
be working with 2D spherical projections of the matter distribution. 
We thus can rewrite Eq.~(\ref{eq:diff_defl}) in terms
of the angular gradient $\nabla_{\hat{{\bf n}}}$ as 
\be
d \boldsymbol\alpha = -\frac{2}{D_A(\eta)}\nabla_{\hat{n}} \Psi d \eta.
\ee
Note that the versor ${\bf \hat{n}}$ refers to angular coordinates on the
sphere, or $\boldsymbol\theta = \eta {\bf \hat{n}}$. A photon in the $k$-th
shell at $\eta_k$ is deflected - due to the presence of matter - by an
angle $\boldsymbol\alpha^{(k)}$
which can be approximated by
\begin{align}
\label{eq:gradphi}
\boldsymbol\alpha^{(k)}(\hat{{\bf n}}) &= - \frac{2}{D_A(\eta_k)}
\frac{D_A(\eta_s - \eta_k)}{D_A(\eta_s)} \int^{\eta_k+\Delta \eta /2}
_{\eta_k-\Delta \eta /2} \nabla_{\hat{{\bf n}}}
\Psi(\tilde{\eta}\hat{{\bf n}};\tilde{\eta}) {\rm d}
\tilde{\eta}, \\
\nonumber & = - \nabla_{\hat{{\bf n}}} \psi^{(k)}(\hat{{\bf n}}),
\end{align}
where the 2-D gravitational potential on the sphere is defined as
\be
\Phi^{(k)}(\hat{{\bf n}}) = \frac{2}{D_A(\eta_k)} \int^{\eta_k+\Delta \eta /2}
_{\eta_k-\Delta \eta /2} \Psi(\tilde{\eta} \hat{ { \bf n
}};\tilde{\eta}){\rm d}\tilde{\eta},
\ee
and the corresponding contribution to the lensing potential is given by
\be
\psi^{(k)}(\hat{{\bf n}}) = \frac{D_A(\eta_s - \eta_k)}{D_A(\eta_s)}\Phi^{(k)}(\hat{{\bf n}}).
\ee
In the previous equations, the notation $( \eta \hat{{\bf n}}; \eta )$
indicates that the potential is 
evaluated when the
photon is at the position $\eta \hat{{\bf n}}$ on the sky, at the comoving distance $\eta$ from the observer.
We can
relate the gravitational potential to the mass overdensity in the
shell using the Poisson equation
\be
\nabla^2 \Psi= \frac{4 \pi G}{c^2} \frac{\rho - \bar{\rho}}{(1+z)^2},
\ee
where $\bar{\rho}$ is the mean matter density of the Universe at redshift
$z$. As in \cite{ValeWhite03}, we can integrate the above equation along the line
of sight to obtain the two dimensional version of
the Poisson equation for the $k$-th mass shell:
\be
\nabla^2_{\hat{{\bf n}}}\Phi^{(k)}(\hat{{\bf n}}) = \frac{8 \pi G}{c^2}
\frac{D_A(\eta_k)}{(1+z_k)^2} \Delta^{(k)}_\Sigma (\hat{{\bf n}}), 
\label{eq:poiss}
\ee
where the surface mass density is defined as
\be
\Delta^{(k)}_\Sigma (\hat{{\bf n}}) = \int^{\eta_k+\Delta \eta /2}
_{\eta_k-\Delta \eta /2} (\rho - \bar{\rho}){\rm d} \tilde{\eta}.
\label{eq::deltasigma}
\ee
In Eq~(\ref{eq:poiss}) we have dropped the term containing the
derivatives in the radial direction, ignoring thus 
long wavelength fluctuations along the line-of-sight via
the Limber approximation \cite{Jain00}.
%, therefore going from the three dimensional to the
%two dimensional version of the Laplacian.
However, as argued by \cite{DasBode,Becker12},
this is, at best, an approximation. In particular \cite{Li2010} showed that this assumption neglects
extra terms that become non-zero in presence of realistic
finite width lens plane. 
The problem arises because 
the matter distribution and, thus, the potential itself may 
become discontinuous at the boundaries if periodic conditions are not
enforced. 
Nevertheless, these corrections to the lens-plane approach adopting the  Born approximation which has been used for this work (see Sect.~\ref{sec:algorithm}) confine this problem to the single shells. In fact, partial
derivatives in the transverse plane
commute with the integral evaluated along the whole line of sight, resulting in the cancellation of
line-of-sight modes as required in the Limber approximation of the integral.
 \footnote{Note for example that assuming a flat-sky approximation, unlike what has been done in this work, is a stronger assumption with respect to the Limber approximation and results in a well-known excess of power on large scales as seen, e.g., by \cite{Hu2000}.}

\noindent
We use the following definition for the convergence field $K^{(k)}$ at
the $k$-th shell,
\be
K^{(k)}(\hat{{\bf n}}) = \frac{4 \pi
  G}{c^2}\frac{D_A(\eta_k)}{(1+z_k)^2}\Delta^{(k)}_{\Sigma}(\hat{{\bf
    n}}), 
\label{eq:kappa1}
\ee
to rewrite Eq.~(\ref{eq:poiss}) as
\be
\nabla^2_{\hat{{\bf n}}}\Phi^{(k)}(\hat{{\bf n}}) = 2 K^{(k)}(\hat{{\bf n}}).
\label{eq:poiss_little}
\ee
The lensing potential on the sphere is related to $K$ via Eq.~(\ref{eq:poiss_little}), and
it can be easily computed by expanding both sides of the Poisson
equation in spherical harmonics, obtaining
the following algebraic relation between the harmonic coefficients of
the two fields:
\be
\Phi_{\ell m} = \frac{2}{\ell(\ell+1)}K_{\ell m}.
\label{eq:kappa_to_phi}
\ee
The monopole term in the lensing potential does not contribute to the
deflection field: therefore to avoid any divergence in the above
equation we can safely set to zero $\Phi_{\ell m}$ for $\ell = 0$ in all calculations.
The quantity $K$ is directy computed when the matter distribution in
the shell is radially projected onto the spherical map; 
as discussed in Sect.~\ref{sec:maps2shell}, it is useful to define an angular surface mass density
$\Delta^{\theta}_{\Sigma}(\hat{{\bf n}})$ as the mass per steradians,
\be
\Delta^{\theta (k)}_{\Sigma}(\hat{{\bf n}}) = \int^{\eta_k+\Delta \eta /2}
_{\eta_k-\Delta \eta /2} (\rho -
\bar{\rho})\frac{D_A(\tilde{\eta})^2}{(1+\tilde{z})^3} {\rm d} \tilde{\eta}.
\ee
%The surface mass density defined in Equation (\ref{eq::deltasigma}) is related
%to this through the relation
%\be
%\Delta_{\Sigma} = \Delta^{\theta}_{\Sigma}\frac{(1+z)^3}{d_A(\eta)^2}.
%\ee
such that Eq.~(\ref{eq:kappa1}) can be rewritten as:
\be
K^{(k)}(\hat{{\bf n}}) = \frac{4 \pi
  G}{c^2}\frac{1+z_k}{D_A(\eta_k)}\Delta^{\theta (k)}_{\Sigma }(\hat{{\bf n}}).
\label{eq:key}
\ee
Finally, the vector field $\boldsymbol\alpha(\hat{{\bf n}})$ will be synthesised, 
as described in \cite{Zaldarriaga96,Kamionkowski96}, from the spherical harmonic components of the potential in
terms of spin-1 spherical harmonics.
The multiple-plane lens formalism can be also applied to exploit the
effective or single-plane approximation to lens the CMB.
Eqs.~(\ref{eq:born_phi}) and (\ref{eq:born_kappa}) can be discretised into the following sums,
\be
\label{eq:born_phi_dis}
\psi^{\rm eff}(\hat{{\bf n}}) = \sum_j \frac{D_A(\eta_s -
  \eta_j)}{D_A(\eta_s)} \Phi^{(j)},
\ee
\be
\kappa(\hat{{\bf n}}) = \sum_j \frac{D_A(\eta_s -
  \eta_j)}{D_A(\eta_s)} K^{(j)},
\ee
where we used the previous definitions of quantities on the $j$-th lens.
In the same framework, the convergence $\kappa$ can be seen as
just a weighted projected surface density \cite{Fosalba08,Hilbert09}:
\be
\label{eq:kappadelta}
\kappa(\boldsymbol\theta) = \frac{3 H^2_0 \Omega_{m,0}}{2 c^2} \int {\rm
  d}\eta \:
\delta(\eta, \boldsymbol\theta) \frac{D_A(\eta_s - \eta)
  D_A(\eta)}{D_A(\eta_s) a(\eta)},
\ee
where $\delta$ is the 3D matter density at radial distance $\eta$ and
angular position $\boldsymbol\theta$, $D_A(\eta_s)$ is the position
of the lensing source at the last scattering surface and $a(\eta)$ is
the scale factor at $\eta$.
Based on the definition in Eq.~(\ref{eq:kappadelta}), the
angular power spectrum of the convergence becomes
\be
\label{eq:clkappa}
C^{\kappa\kappa}_{\ell} = \frac{9 H^4_0 \Omega^2_{m,0}}{4 c^4}
  \int^{\eta_s}_0 {\rm d}\eta
  P(k,z)\left(\frac{D_A(\eta_s - \eta)}{D_A(\eta_s) a}\right)^2,
\ee
where $P(k,z)$ is 3D matter power spectrum, computed via the Limber
approximation at $k= \ell / D_A(\eta(z))$, valid for $\ell > 10$ within a few
percent accuracy \cite{ValeWhite03}. The discretised equation reads:
\be
\label{eq:clkappa_dis}
C_{\ell}^{\kappa\kappa} = \frac{9 H^4_0 \Omega^2_{m,0}}{4 c^4} \sum_k
  \Delta \eta \:
  P(\ell/D_A(\eta_k),z_k) \left( \frac{D_A(\eta_s - \eta_k)}{D_A(\eta_s)
    a_k} \right)^2,
\ee
summing all over the $k$ lens plane. Note that the convergence field can be converted into
lensing potential using the Poisson equation, or in terms of the angular
power spectrum:
\be
\label{eq:kappa2lens}
C^{\psi\psi}_{\ell} = \frac{4}{\ell^2(\ell+1)^2}C^{\kappa\kappa}_{\ell}.
\ee

jcap_article_mlcode.bbl
\providecommand{\href}[2]{#2}\begingroup\raggedright\begin{thebibliography}{10}

\bibitem{PlanckXVI}
{Planck Collaboration}, P.~A.~R. {Ade}, N.~{Aghanim}, C.~{Armitage-Caplan},
  M.~{Arnaud}, M.~{Ashdown}, F.~{Atrio-Barandela}, J.~{Aumont},
  C.~{Baccigalupi}, A.~J. {Banday}, and et~al., {\it {Planck 2013 results. XVI.
  Cosmological parameters}},  {\em ArXiv e-prints} (Mar., 2013)
  [\href{http://arxiv.org/abs/1303.5076}{{\tt arXiv:1303.5076}}].

\bibitem{EBEX}
B.~{Reichborn-Kjennerud}, A.~M. {Aboobaker}, P.~{Ade}, F.~{Aubin},
  C.~{Baccigalupi}, C.~{Bao}, J.~{Borrill}, C.~{Cantalupo}, D.~{Chapman},
  J.~{Didier}, M.~{Dobbs}, J.~{Grain}, W.~{Grainger}, S.~{Hanany},
  S.~{Hillbrand}, J.~{Hubmayr}, A.~{Jaffe}, B.~{Johnson}, T.~{Jones},
  T.~{Kisner}, J.~{Klein}, A.~{Korotkov}, S.~{Leach}, A.~{Lee}, L.~{Levinson},
  M.~{Limon}, K.~{MacDermid}, T.~{Matsumura}, X.~{Meng}, A.~{Miller},
  M.~{Milligan}, E.~{Pascale}, D.~{Polsgrove}, N.~{Ponthieu}, K.~{Raach},
  I.~{Sagiv}, G.~{Smecher}, F.~{Stivoli}, R.~{Stompor}, H.~{Tran},
  M.~{Tristram}, G.~S. {Tucker}, Y.~{Vinokurov}, A.~{Yadav}, M.~{Zaldarriaga},
  and K.~{Zilic}, {\it {EBEX: a balloon-borne CMB polarization experiment}},
  in {\em Society of Photo-Optical Instrumentation Engineers (SPIE) Conference
  Series}, vol.~7741 of {\em Society of Photo-Optical Instrumentation Engineers
  (SPIE) Conference Series}, July, 2010.
\newblock \href{http://arxiv.org/abs/1007.3672}{{\tt arXiv:1007.3672}}.

\bibitem{Polarbear}
Z.~D. {Kermish}, P.~{Ade}, A.~{Anthony}, K.~{Arnold}, D.~{Barron},
  D.~{Boettger}, J.~{Borrill}, S.~{Chapman}, Y.~{Chinone}, M.~A. {Dobbs},
  J.~{Errard}, G.~{Fabbian}, D.~{Flanigan}, G.~{Fuller}, A.~{Ghribi},
  W.~{Grainger}, N.~{Halverson}, M.~{Hasegawa}, K.~{Hattori}, M.~{Hazumi},
  W.~L. {Holzapfel}, J.~{Howard}, P.~{Hyland}, A.~{Jaffe}, B.~{Keating},
  T.~{Kisner}, A.~T. {Lee}, M.~{Le Jeune}, E.~{Linder}, M.~{Lungu},
  F.~{Matsuda}, T.~{Matsumura}, X.~{Meng}, N.~J. {Miller}, H.~{Morii},
  S.~{Moyerman}, M.~J. {Myers}, H.~{Nishino}, H.~{Paar}, E.~{Quealy}, C.~L.
  {Reichardt}, P.~L. {Richards}, C.~{Ross}, A.~{Shimizu}, M.~{Shimon},
  C.~{Shimmin}, M.~{Sholl}, P.~{Siritanasak}, H.~{Spieler}, N.~{Stebor},
  B.~{Steinbach}, R.~{Stompor}, A.~{Suzuki}, T.~{Tomaru}, C.~{Tucker}, and
  O.~{Zahn}, {\it {The POLARBEAR experiment}},  in {\em Society of
  Photo-Optical Instrumentation Engineers (SPIE) Conference Series}, vol.~8452
  of {\em Society of Photo-Optical Instrumentation Engineers (SPIE) Conference
  Series}, Sept., 2012.
\newblock \href{http://arxiv.org/abs/1210.7768}{{\tt arXiv:1210.7768}}.

\bibitem{SPTpol}
J.~E. {Austermann}, K.~A. {Aird}, J.~A. {Beall}, D.~{Becker}, A.~{Bender},
  B.~A. {Benson}, L.~E. {Bleem}, J.~{Britton}, J.~E. {Carlstrom}, C.~L.
  {Chang}, H.~C. {Chiang}, H.-M. {Cho}, T.~M. {Crawford}, A.~T. {Crites},
  A.~{Datesman}, T.~{de Haan}, M.~A. {Dobbs}, E.~M. {George}, N.~W.
  {Halverson}, N.~{Harrington}, J.~W. {Henning}, G.~C. {Hilton}, G.~P.
  {Holder}, W.~L. {Holzapfel}, S.~{Hoover}, N.~{Huang}, J.~{Hubmayr}, K.~D.
  {Irwin}, R.~{Keisler}, J.~{Kennedy}, L.~{Knox}, A.~T. {Lee}, E.~{Leitch},
  D.~{Li}, M.~{Lueker}, D.~P. {Marrone}, J.~J. {McMahon}, J.~{Mehl}, S.~S.
  {Meyer}, T.~E. {Montroy}, T.~{Natoli}, J.~P. {Nibarger}, M.~D. {Niemack},
  V.~{Novosad}, S.~{Padin}, C.~{Pryke}, C.~L. {Reichardt}, J.~E. {Ruhl}, B.~R.
  {Saliwanchik}, J.~T. {Sayre}, K.~K. {Schaffer}, E.~{Shirokoff}, A.~A.
  {Stark}, K.~{Story}, K.~{Vanderlinde}, J.~D. {Vieira}, G.~{Wang},
  R.~{Williamson}, V.~{Yefremenko}, K.~W. {Yoon}, and O.~{Zahn}, {\it {SPTpol:
  an instrument for CMB polarization measurements with the South Pole
  Telescope}},  in {\em Society of Photo-Optical Instrumentation Engineers
  (SPIE) Conference Series}, vol.~8452 of {\em Society of Photo-Optical
  Instrumentation Engineers (SPIE) Conference Series}, Sept., 2012.
\newblock \href{http://arxiv.org/abs/1210.4970}{{\tt arXiv:1210.4970}}.

\bibitem{ACTpol}
M.~D. {Niemack}, P.~A.~R. {Ade}, J.~{Aguirre}, F.~{Barrientos}, J.~A. {Beall},
  J.~R. {Bond}, J.~{Britton}, H.~M. {Cho}, S.~{Das}, M.~J. {Devlin},
  S.~{Dicker}, J.~{Dunkley}, R.~{D{\"u}nner}, J.~W. {Fowler}, A.~{Hajian},
  M.~{Halpern}, M.~{Hasselfield}, G.~C. {Hilton}, M.~{Hilton}, J.~{Hubmayr},
  J.~P. {Hughes}, L.~{Infante}, K.~D. {Irwin}, N.~{Jarosik}, J.~{Klein},
  A.~{Kosowsky}, T.~A. {Marriage}, J.~{McMahon}, F.~{Menanteau}, K.~{Moodley},
  J.~P. {Nibarger}, M.~R. {Nolta}, L.~A. {Page}, B.~{Partridge}, E.~D. {Reese},
  J.~{Sievers}, D.~N. {Spergel}, S.~T. {Staggs}, R.~{Thornton}, C.~{Tucker},
  E.~{Wollack}, and K.~W. {Yoon}, {\it {ACTPol: a polarization-sensitive
  receiver for the Atacama Cosmology Telescope}},  in {\em Society of
  Photo-Optical Instrumentation Engineers (SPIE) Conference Series}, vol.~7741
  of {\em Society of Photo-Optical Instrumentation Engineers (SPIE) Conference
  Series}, July, 2010.
\newblock \href{http://arxiv.org/abs/1006.5049}{{\tt arXiv:1006.5049}}.

\bibitem{Spider}
A.~A. {Fraisse}, P.~A.~R. {Ade}, M.~{Amiri}, S.~J. {Benton}, J.~J. {Bock},
  J.~R. {Bond}, J.~A. {Bonetti}, S.~{Bryan}, B.~{Burger}, H.~C. {Chiang}, C.~N.
  {Clark}, C.~R. {Contaldi}, B.~P. {Crill}, G.~{Davis}, O.~{Dor{\'e}},
  M.~{Farhang}, J.~P. {Filippini}, L.~M. {Fissel}, N.~N. {Gandilo},
  S.~{Golwala}, J.~E. {Gudmundsson}, M.~{Hasselfield}, G.~{Hilton},
  W.~{Holmes}, V.~V. {Hristov}, K.~{Irwin}, W.~C. {Jones}, C.~L. {Kuo}, C.~J.
  {MacTavish}, P.~V. {Mason}, T.~E. {Montroy}, T.~A. {Morford}, C.~B.
  {Netterfield}, D.~T. {O'Dea}, A.~S. {Rahlin}, C.~{Reintsema}, J.~E. {Ruhl},
  M.~C. {Runyan}, M.~A. {Schenker}, J.~A. {Shariff}, J.~D. {Soler},
  A.~{Trangsrud}, C.~{Tucker}, R.~S. {Tucker}, A.~D. {Turner}, and D.~{Wiebe},
  {\it {SPIDER: probing the early Universe with a suborbital polarimeter}},
  {\em \jcap} {\bf 4} (Apr., 2013) 47,
  [\href{http://arxiv.org/abs/1106.3087}{{\tt arXiv:1106.3087}}].

\bibitem{keck}
C.~D. {Sheehy}, P.~A.~R. {Ade}, R.~W. {Aikin}, M.~{Amiri}, S.~{Benton},
  C.~{Bischoff}, J.~J. {Bock}, J.~A. {Bonetti}, J.~A. {Brevik}, B.~{Burger},
  C.~D. {Dowell}, L.~{Duband}, J.~P. {Filippini}, S.~R. {Golwala},
  M.~{Halpern}, M.~{Hasselfield}, G.~{Hilton}, V.~V. {Hristov}, K.~{Irwin},
  J.~P. {Kaufman}, B.~G. {Keating}, J.~M. {Kovac}, C.~L. {Kuo}, A.~E. {Lange},
  E.~M. {Leitch}, M.~{Lueker}, C.~B. {Netterfield}, H.~T. {Nguyen}, R.~W.
  {Ogburn}, IV, A.~{Orlando}, C.~L. {Pryke}, C.~{Reintsema}, S.~{Richter},
  J.~E. {Ruhl}, M.~C. {Runyan}, Z.~{Staniszewski}, S.~{Stokes}, R.~{Sudiwala},
  G.~{Teply}, K.~L. {Thompson}, J.~E. {Tolan}, A.~D. {Turner}, P.~{Wilson}, and
  C.~L. {Wong}, {\it {The Keck Array: a pulse tube cooled CMB polarimeter}},
  {\em ArXiv e-prints} (Apr., 2011) [\href{http://arxiv.org/abs/1104.5516}{{\tt
  arXiv:1104.5516}}].

\bibitem{Zaldarriaga96}
M.~{Zaldarriaga} and U.~{Seljak}, {\it {All-sky analysis of polarization in the
  microwave background}},  {\em \prd} {\bf 55} (Feb., 1997) 1830--1840,
  [\href{http://arxiv.org/abs/astro-ph/9609170}{{\tt astro-ph/9609170}}].

\bibitem{Kamionkowski96}
M.~{Kamionkowski}, A.~{Kosowsky}, and A.~{Stebbins}, {\it {Statistics of cosmic
  microwave background polarization}},  {\em \prd} {\bf 55} (June, 1997)
  7368--7388, [\href{http://arxiv.org/abs/astro-ph/9611125}{{\tt
  astro-ph/9611125}}].

\bibitem{Seljak97}
U.~{Seljak} and M.~{Zaldarriaga}, {\it {Signature of Gravity Waves in the
  Polarization of the Microwave Background}},  {\em Physical Review Letters}
  {\bf 78} (Mar., 1997) 2054--2057,
  [\href{http://arxiv.org/abs/astro-ph/9609169}{{\tt astro-ph/9609169}}].

\bibitem{bicep2}
{BICEP2 Collaboration}, P.~A.~R. {Ade}, R.~W. {Aikin}, M.~{Amiri},
  D.~{Barkats}, S.~J. {Benton}, C.~A. {Bischoff}, J.~J. {Bock}, J.~A. {Brevik},
  I.~{Buder}, E.~{Bullock}, G.~{Davis}, C.~D. {Dowell}, L.~{Duband}, J.~P.
  {Filippini}, S.~{Fliescher}, S.~R. {Golwala}, M.~{Halpern}, M.~{Hasselfield},
  S.~R. {Hildebrandt}, G.~C. {Hilton}, V.~V. {Hristov}, K.~D. {Irwin}, K.~S.
  {Karkare}, J.~P. {Kaufman}, B.~G. {Keating}, S.~A. {Kernasovskiy}, J.~M.
  {Kovac}, C.~L. {Kuo}, E.~M. {Leitch}, N.~{Llombart}, M.~{Lueker}, C.~B.
  {Netterfield}, H.~T. {Nguyen}, R.~{O'Brient}, R.~W. {Ogburn}, IV,
  A.~{Orlando}, C.~{Pryke}, C.~D. {Reintsema}, S.~{Richter}, R.~{Schwarz},
  C.~D. {Sheehy}, Z.~K. {Staniszewski}, K.~T. {Story}, R.~V. {Sudiwala}, G.~P.
  {Teply}, J.~E. {Tolan}, A.~D. {Turner}, A.~G. {Vieregg}, P.~{Wilson}, C.~L.
  {Wong}, and K.~W. {Yoon}, {\it {BICEP2 II: Experiment and Three-Year Data
  Set}},  {\em ArXiv e-prints} (Mar., 2014)
  [\href{http://arxiv.org/abs/1403.4302}{{\tt arXiv:1403.4302}}].

\bibitem{flauger}
R.~{Flauger}, J.~C. {Hill}, and D.~N. {Spergel}, {\it {Toward an understanding
  of foreground emission in the BICEP2 region}},  {\em \jcap} {\bf 8} (Aug.,
  2014) 39, [\href{http://arxiv.org/abs/1405.7351}{{\tt arXiv:1405.7351}}].

\bibitem{Aghanim}
N.~{Aghanim}, S.~{Majumdar}, and J.~{Silk}, {\it {Secondary anisotropies of the
  CMB}},  {\em Reports on Progress in Physics} {\bf 71} (June, 2008) 066902,
  [\href{http://arxiv.org/abs/0711.0518}{{\tt arXiv:0711.0518}}].

\bibitem{LewisChallinor06}
A.~{Lewis} and A.~{Challinor}, {\it {Weak gravitational lensing of the CMB}},
  {\em \physrep} {\bf 429} (June, 2006) 1--65,
  [\href{http://arxiv.org/abs/astro-ph/0601594}{{\tt astro-ph/0601594}}].

\bibitem{ACT}
S.~{Das}, B.~D. {Sherwin}, P.~{Aguirre}, J.~W. {Appel}, J.~R. {Bond}, C.~S.
  {Carvalho}, M.~J. {Devlin}, J.~{Dunkley}, R.~{D{\"u}nner},
  T.~{Essinger-Hileman}, J.~W. {Fowler}, A.~{Hajian}, M.~{Halpern},
  M.~{Hasselfield}, A.~D. {Hincks}, R.~{Hlozek}, K.~M. {Huffenberger}, J.~P.
  {Hughes}, K.~D. {Irwin}, J.~{Klein}, A.~{Kosowsky}, R.~H. {Lupton}, T.~A.
  {Marriage}, D.~{Marsden}, F.~{Menanteau}, K.~{Moodley}, M.~D. {Niemack},
  M.~R. {Nolta}, L.~A. {Page}, L.~{Parker}, E.~D. {Reese}, B.~L. {Schmitt},
  N.~{Sehgal}, J.~{Sievers}, D.~N. {Spergel}, S.~T. {Staggs}, D.~S. {Swetz},
  E.~R. {Switzer}, R.~{Thornton}, K.~{Visnjic}, and E.~{Wollack}, {\it
  {Detection of the Power Spectrum of Cosmic Microwave Background Lensing by
  the Atacama Cosmology Telescope}},  {\em Physical Review Letters} {\bf 107}
  (July, 2011) 021301, [\href{http://arxiv.org/abs/1103.2124}{{\tt
  arXiv:1103.2124}}].

\bibitem{SPT}
R.~{Keisler}, C.~L. {Reichardt}, K.~A. {Aird}, B.~A. {Benson}, L.~E. {Bleem},
  J.~E. {Carlstrom}, C.~L. {Chang}, H.~M. {Cho}, T.~M. {Crawford}, A.~T.
  {Crites}, T.~{de Haan}, M.~A. {Dobbs}, J.~{Dudley}, E.~M. {George}, N.~W.
  {Halverson}, G.~P. {Holder}, W.~L. {Holzapfel}, S.~{Hoover}, Z.~{Hou}, J.~D.
  {Hrubes}, M.~{Joy}, L.~{Knox}, A.~T. {Lee}, E.~M. {Leitch}, M.~{Lueker},
  D.~{Luong-Van}, J.~J. {McMahon}, J.~{Mehl}, S.~S. {Meyer}, M.~{Millea}, J.~J.
  {Mohr}, T.~E. {Montroy}, T.~{Natoli}, S.~{Padin}, T.~{Plagge}, C.~{Pryke},
  J.~E. {Ruhl}, K.~K. {Schaffer}, L.~{Shaw}, E.~{Shirokoff}, H.~G. {Spieler},
  Z.~{Staniszewski}, A.~A. {Stark}, K.~{Story}, A.~{van Engelen},
  K.~{Vanderlinde}, J.~D. {Vieira}, R.~{Williamson}, and O.~{Zahn}, {\it {A
  Measurement of the Damping Tail of the Cosmic Microwave Background Power
  Spectrum with the South Pole Telescope}},  {\em \apj} {\bf 743} (Dec., 2011)
  28, [\href{http://arxiv.org/abs/1105.3182}{{\tt arXiv:1105.3182}}].

\bibitem{vanEngelen}
A.~{van Engelen}, R.~{Keisler}, O.~{Zahn}, K.~A. {Aird}, B.~A. {Benson}, L.~E.
  {Bleem}, J.~E. {Carlstrom}, C.~L. {Chang}, H.~M. {Cho}, T.~M. {Crawford},
  A.~T. {Crites}, T.~{de Haan}, M.~A. {Dobbs}, J.~{Dudley}, E.~M. {George},
  N.~W. {Halverson}, G.~P. {Holder}, W.~L. {Holzapfel}, S.~{Hoover}, Z.~{Hou},
  J.~D. {Hrubes}, M.~{Joy}, L.~{Knox}, A.~T. {Lee}, E.~M. {Leitch},
  M.~{Lueker}, D.~{Luong-Van}, J.~J. {McMahon}, J.~{Mehl}, S.~S. {Meyer},
  M.~{Millea}, J.~J. {Mohr}, T.~E. {Montroy}, T.~{Natoli}, S.~{Padin},
  T.~{Plagge}, C.~{Pryke}, C.~L. {Reichardt}, J.~E. {Ruhl}, J.~T. {Sayre},
  K.~K. {Schaffer}, L.~{Shaw}, E.~{Shirokoff}, H.~G. {Spieler},
  Z.~{Staniszewski}, A.~A. {Stark}, K.~{Story}, K.~{Vanderlinde}, J.~D.
  {Vieira}, and R.~{Williamson}, {\it {A Measurement of Gravitational Lensing
  of the Microwave Background Using South Pole Telescope Data}},  {\em \apj}
  {\bf 756} (Sept., 2012) 142, [\href{http://arxiv.org/abs/1202.0546}{{\tt
  arXiv:1202.0546}}].

\bibitem{pblens}
{POLARBEAR Collaboration}, P.~A.~R. {Ade}, Y.~{Akiba}, A.~E. {Anthony},
  K.~{Arnold}, M.~{Atlas}, D.~{Barron}, D.~{Boettger}, J.~{Borrill},
  S.~{Chapman}, Y.~{Chinone}, M.~{Dobbs}, T.~{Elleflot}, J.~{Errard},
  G.~{Fabbian}, C.~{Feng}, D.~{Flanigan}, A.~{Gilbert}, W.~{Grainger}, N.~W.
  {Halverson}, M.~{Hasegawa}, K.~{Hattori}, M.~{Hazumi}, W.~L. {Holzapfel},
  Y.~{Hori}, J.~{Howard}, P.~{Hyland}, Y.~{Inoue}, G.~C. {Jaehnig}, A.~{Jaffe},
  B.~{Keating}, Z.~{Kermish}, R.~{Keskitalo}, T.~{Kisner}, M.~{Le Jeune}, A.~T.
  {Lee}, E.~{Linder}, E.~M. {Leitch}, M.~{Lungu}, F.~{Matsuda}, T.~{Matsumura},
  X.~{Meng}, N.~J. {Miller}, H.~{Morii}, S.~{Moyerman}, M.~J. {Myers},
  M.~{Navaroli}, H.~{Nishino}, H.~{Paar}, J.~{Peloton}, E.~{Quealy},
  G.~{Rebeiz}, C.~L. {Reichardt}, P.~L. {Richards}, C.~{Ross}, I.~{Schanning},
  D.~E. {Schenck}, B.~{Sherwin}, A.~{Shimizu}, C.~{Shimmin}, M.~{Shimon},
  P.~{Siritanasak}, G.~{Smecher}, H.~{Spieler}, N.~{Stebor}, B.~{Steinbach},
  R.~{Stompor}, A.~{Suzuki}, S.~{Takakura}, T.~{Tomaru}, B.~{Wilson},
  A.~{Yadav}, and O.~{Zahn}, {\it {Measurement of the Cosmic Microwave
  Background Polarization Lensing Power Spectrum with the POLARBEAR
  experiment}},  {\em ArXiv e-prints} (Dec., 2013)
  [\href{http://arxiv.org/abs/1312.6646}{{\tt arXiv:1312.6646}}].

\bibitem{pbbb}
{The POLARBEAR Collaboration}, P.~A.~R. {Ade}, Y.~{Akiba}, A.~E. {Anthony},
  K.~{Arnold}, M.~{Atlas}, D.~{Barron}, D.~{Boettger}, J.~{Borrill},
  S.~{Chapman}, Y.~{Chinone}, M.~{Dobbs}, T.~{Elleflot}, J.~{Errard},
  G.~{Fabbian}, C.~{Feng}, D.~{Flanigan}, A.~{Gilbert}, W.~{Grainger}, N.~W.
  {Halverson}, M.~{Hasegawa}, K.~{Hattori}, M.~{Hazumi}, W.~L. {Holzapfel},
  Y.~{Hori}, J.~{Howard}, P.~{Hyland}, Y.~{Inoue}, G.~C. {Jaehnig}, A.~H.
  {Jaffe}, B.~{Keating}, Z.~{Kermish}, R.~{Keskitalo}, T.~{Kisner}, M.~{Le
  Jeune}, A.~T. {Lee}, E.~M. {Leitch}, E.~{Linder}, M.~{Lungu}, F.~{Matsuda},
  T.~{Matsumura}, X.~{Meng}, N.~J. {Miller}, H.~{Morii}, S.~{Moyerman}, M.~J.
  {Myers}, M.~{Navaroli}, H.~{Nishino}, H.~{Paar}, J.~{Peloton}, D.~{Poletti},
  E.~{Quealy}, G.~{Rebeiz}, C.~L. {Reichardt}, P.~L. {Richards}, C.~{Ross},
  I.~{Schanning}, D.~E. {Schenck}, B.~D. {Sherwin}, A.~{Shimizu}, C.~{Shimmin},
  M.~{Shimon}, P.~{Siritanasak}, G.~{Smecher}, H.~{Spieler}, N.~{Stebor},
  B.~{Steinbach}, R.~{Stompor}, A.~{Suzuki}, S.~{Takakura}, T.~{Tomaru},
  B.~{Wilson}, A.~{Yadav}, and O.~{Zahn}, {\it {A Measurement of the Cosmic
  Microwave Background B-Mode Polarization Power Spectrum at Sub-Degree Scales
  with POLARBEAR}},  {\em ArXiv e-prints} (Mar., 2014)
  [\href{http://arxiv.org/abs/1403.2369}{{\tt arXiv:1403.2369}}].

\bibitem{hanson13}
D.~{Hanson}, S.~{Hoover}, A.~{Crites}, P.~A.~R. {Ade}, K.~A. {Aird}, J.~E.
  {Austermann}, J.~A. {Beall}, A.~N. {Bender}, B.~A. {Benson}, L.~E. {Bleem},
  J.~J. {Bock}, J.~E. {Carlstrom}, C.~L. {Chang}, H.~C. {Chiang}, H.-M. {Cho},
  A.~{Conley}, T.~M. {Crawford}, T.~{de Haan}, M.~A. {Dobbs}, W.~{Everett},
  J.~{Gallicchio}, J.~{Gao}, E.~M. {George}, N.~W. {Halverson},
  N.~{Harrington}, J.~W. {Henning}, G.~C. {Hilton}, G.~P. {Holder}, W.~L.
  {Holzapfel}, J.~D. {Hrubes}, N.~{Huang}, J.~{Hubmayr}, K.~D. {Irwin},
  R.~{Keisler}, L.~{Knox}, A.~T. {Lee}, E.~{Leitch}, D.~{Li}, C.~{Liang},
  D.~{Luong-Van}, G.~{Marsden}, J.~J. {McMahon}, J.~{Mehl}, S.~S. {Meyer},
  L.~{Mocanu}, T.~E. {Montroy}, T.~{Natoli}, J.~P. {Nibarger}, V.~{Novosad},
  S.~{Padin}, C.~{Pryke}, C.~L. {Reichardt}, J.~E. {Ruhl}, B.~R. {Saliwanchik},
  J.~T. {Sayre}, K.~K. {Schaffer}, B.~{Schulz}, G.~{Smecher}, A.~A. {Stark},
  K.~T. {Story}, C.~{Tucker}, K.~{Vanderlinde}, J.~D. {Vieira}, M.~P. {Viero},
  G.~{Wang}, V.~{Yefremenko}, O.~{Zahn}, and M.~{Zemcov}, {\it {Detection of
  B-Mode Polarization in the Cosmic Microwave Background with Data from the
  South Pole Telescope}},  {\em Physical Review Letters} {\bf 111} (Oct., 2013)
  141301, [\href{http://arxiv.org/abs/1307.5830}{{\tt arXiv:1307.5830}}].

\bibitem{AcquavivaBacci06}
V.~{Acquaviva} and C.~{Baccigalupi}, {\it {Dark energy records in lensed cosmic
  microwave background}},  {\em \prd} {\bf 74} (Nov., 2006) 103510,
  [\href{http://arxiv.org/abs/astro-ph/0507644}{{\tt astro-ph/0507644}}].

\bibitem{Smith06}
K.~M. {Smith}, W.~{Hu}, and M.~{Kaplinghat}, {\it {Cosmological information
  from lensed CMB power spectra}},  {\em \prd} {\bf 74} (Dec., 2006) 123002,
  [\href{http://arxiv.org/abs/astro-ph/0607315}{{\tt astro-ph/0607315}}].

\bibitem{hu-dark}
W.~{Hu}, {\it {Dark energy and matter evolution from lensing tomography}},
  {\em \prd} {\bf 66} (Oct., 2002) 083515,
  [\href{http://arxiv.org/abs/astro-ph/0208093}{{\tt astro-ph/0208093}}].

\bibitem{Bleem12}
L.~E. {Bleem}, A.~{van Engelen}, G.~P. {Holder}, K.~A. {Aird}, R.~{Armstrong},
  M.~L.~N. {Ashby}, M.~R. {Becker}, B.~A. {Benson}, T.~{Biesiadzinski},
  M.~{Brodwin}, M.~T. {Busha}, J.~E. {Carlstrom}, C.~L. {Chang}, H.~M. {Cho},
  T.~M. {Crawford}, A.~T. {Crites}, T.~{de Haan}, S.~{Desai}, M.~A. {Dobbs},
  O.~{Dor{\'e}}, J.~{Dudley}, J.~E. {Geach}, E.~M. {George}, M.~D. {Gladders},
  A.~H. {Gonzalez}, N.~W. {Halverson}, N.~{Harrington}, F.~W. {High}, B.~P.
  {Holden}, W.~L. {Holzapfel}, S.~{Hoover}, J.~D. {Hrubes}, M.~{Joy},
  R.~{Keisler}, L.~{Knox}, A.~T. {Lee}, E.~M. {Leitch}, M.~{Lueker},
  D.~{Luong-Van}, D.~P. {Marrone}, J.~{Martinez-Manso}, J.~J. {McMahon},
  J.~{Mehl}, S.~S. {Meyer}, J.~J. {Mohr}, T.~E. {Montroy}, T.~{Natoli},
  S.~{Padin}, T.~{Plagge}, C.~{Pryke}, C.~L. {Reichardt}, A.~{Rest}, J.~E.
  {Ruhl}, B.~R. {Saliwanchik}, J.~T. {Sayre}, K.~K. {Schaffer}, L.~{Shaw},
  E.~{Shirokoff}, H.~G. {Spieler}, B.~{Stalder}, S.~A. {Stanford},
  Z.~{Staniszewski}, A.~A. {Stark}, D.~{Stern}, K.~{Story}, A.~{Vallinotto},
  K.~{Vanderlinde}, J.~D. {Vieira}, R.~H. {Wechsler}, R.~{Williamson}, and
  O.~{Zahn}, {\it {A Measurement of the Correlation of Galaxy Surveys with CMB
  Lensing Convergence Maps from the South Pole Telescope}},  {\em \apjl} {\bf
  753} (July, 2012) L9, [\href{http://arxiv.org/abs/1203.4808}{{\tt
  arXiv:1203.4808}}].

\bibitem{Sherwin12}
B.~D. {Sherwin}, S.~{Das}, A.~{Hajian}, G.~{Addison}, J.~R. {Bond},
  D.~{Crichton}, M.~J. {Devlin}, J.~{Dunkley}, M.~B. {Gralla}, M.~{Halpern},
  J.~C. {Hill}, A.~D. {Hincks}, J.~P. {Hughes}, K.~{Huffenberger}, R.~{Hlozek},
  A.~{Kosowsky}, T.~{Louis}, T.~A. {Marriage}, D.~{Marsden}, F.~{Menanteau},
  K.~{Moodley}, M.~D. {Niemack}, L.~A. {Page}, E.~D. {Reese}, N.~{Sehgal},
  J.~{Sievers}, C.~{Sif{\'o}n}, D.~N. {Spergel}, S.~T. {Staggs}, E.~R.
  {Switzer}, and E.~{Wollack}, {\it {The Atacama Cosmology Telescope:
  Cross-correlation of cosmic microwave background lensing and quasars}},  {\em
  \prd} {\bf 86} (Oct., 2012) 083006,
  [\href{http://arxiv.org/abs/1207.4543}{{\tt arXiv:1207.4543}}].

\bibitem{act-cfht}
N.~{Hand}, A.~{Leauthaud}, S.~{Das}, B.~D. {Sherwin}, G.~E. {Addison}, J.~R.
  {Bond}, E.~{Calabrese}, A.~{Charbonnier}, M.~J. {Devlin}, J.~{Dunkley},
  T.~{Erben}, A.~{Hajian}, M.~{Halpern}, J.~{Harnois-D{\'e}raps}, C.~{Heymans},
  H.~{Hildebrandt}, A.~D. {Hincks}, J.-P. {Kneib}, A.~{Kosowsky}, M.~{Makler},
  L.~{Miller}, K.~{Moodley}, B.~{Moraes}, M.~D. {Niemack}, L.~A. {Page},
  B.~{Partridge}, N.~{Sehgal}, H.~{Shan}, J.~L. {Sievers}, D.~N. {Spergel},
  S.~T. {Staggs}, E.~R. {Switzer}, J.~E. {Taylor}, L.~{Van Waerbeke}, and E.~J.
  {Wollack}, {\it {First Measurement of the Cross-Correlation of CMB Lensing
  and Galaxy Lensing}},  {\em ArXiv e-prints} (Nov., 2013)
  [\href{http://arxiv.org/abs/1311.6200}{{\tt arXiv:1311.6200}}].

\bibitem{pb-cib}
P.~A.~R. {Ade}, Y.~{Akiba}, A.~E. {Anthony}, K.~{Arnold}, M.~{Atlas},
  D.~{Barron}, D.~{Boettger}, J.~{Borrill}, S.~{Chapman}, Y.~{Chinone},
  M.~{Dobbs}, T.~{Elleflot}, J.~{Errard}, G.~{Fabbian}, C.~{Feng},
  D.~{Flanigan}, A.~{Gilbert}, W.~{Grainger}, N.~W. {Halverson}, M.~{Hasegawa},
  K.~{Hattori}, M.~{Hazumi}, W.~L. {Holzapfel}, Y.~{Hori}, J.~{Howard},
  P.~{Hyland}, Y.~{Inoue}, G.~C. {Jaehnig}, A.~{Jaffe}, B.~{Keating},
  Z.~{Kermish}, R.~{Keskitalo}, T.~{Kisner}, M.~{Le Jeune}, A.~T. {Lee},
  E.~{Linder}, E.~M. {Leitch}, M.~{Lungu}, F.~{Matsuda}, T.~{Matsumura},
  X.~{Meng}, N.~J. {Miller}, H.~{Morii}, S.~{Moyerman}, M.~J. {Myers},
  M.~{Navaroli}, H.~{Nishino}, H.~{Paar}, J.~{Peloton}, E.~{Quealy},
  G.~{Rebeiz}, C.~L. {Reichardt}, P.~L. {Richards}, C.~{Ross}, I.~{Schanning},
  D.~E. {Schenck}, B.~{Sherwin}, A.~{Shimizu}, C.~{Shimmin}, M.~{Shimon},
  P.~{Siritanasak}, G.~{Smecher}, H.~{Spieler}, N.~{Stebor}, B.~{Steinbach},
  R.~{Stompor}, A.~{Suzuki}, S.~{Takakura}, T.~{Tomaru}, B.~{Wilson},
  A.~{Yadav}, O.~{Zahn}, and {Polarbear Collaboration}, {\it {Measurement of
  the Cosmic Microwave Background Polarization Lensing Power Spectrum with the
  POLARBEAR Experiment}},  {\em Physical Review Letters} {\bf 113} (July, 2014)
  021301.

\bibitem{Laureijs11}
R.~{Laureijs}, J.~{Amiaux}, S.~{Arduini}, J.~. {Augu{\`e}res}, J.~{Brinchmann},
  R.~{Cole}, M.~{Cropper}, C.~{Dabin}, L.~{Duvet}, A.~{Ealet}, and et~al., {\it
  {Euclid Definition Study Report}},  {\em ArXiv e-prints} (Oct., 2011)
  [\href{http://arxiv.org/abs/1110.3193}{{\tt arXiv:1110.3193}}].

\bibitem{Amendola12}
L.~{Amendola}, S.~{Appleby}, D.~{Bacon}, T.~{Baker}, M.~{Baldi}, N.~{Bartolo},
  A.~{Blanchard}, C.~{Bonvin}, S.~{Borgani}, E.~{Branchini}, C.~{Burrage},
  S.~{Camera}, C.~{Carbone}, L.~{Casarini}, M.~{Cropper}, C.~{de Rham}, C.~{Di
  Porto}, A.~{Ealet}, P.~G. {Ferreira}, F.~{Finelli},
  J.~{Garc{\'{\i}}a-Bellido}, T.~{Giannantonio}, L.~{Guzzo}, A.~{Heavens},
  L.~{Heisenberg}, C.~{Heymans}, H.~{Hoekstra}, L.~{Hollenstein}, R.~{Holmes},
  O.~{Horst}, K.~{Jahnke}, T.~D. {Kitching}, T.~{Koivisto}, M.~{Kunz}, G.~{La
  Vacca}, M.~{March}, E.~{Majerotto}, K.~{Markovic}, D.~{Marsh}, F.~{Marulli},
  R.~{Massey}, Y.~{Mellier}, D.~F. {Mota}, N.~{Nunes}, W.~{Percival},
  V.~{Pettorino}, C.~{Porciani}, C.~{Quercellini}, J.~{Read}, M.~{Rinaldi},
  D.~{Sapone}, R.~{Scaramella}, C.~{Skordis}, F.~{Simpson}, A.~{Taylor},
  S.~{Thomas}, R.~{Trotta}, L.~{Verde}, F.~{Vernizzi}, A.~{Vollmer}, Y.~{Wang},
  J.~{Weller}, and T.~{Zlosnik}, {\it {Cosmology and Fundamental Physics with
  the Euclid Satellite}},  {\em Living Reviews in Relativity} {\bf 16} (Sept.,
  2013) 6, [\href{http://arxiv.org/abs/1206.1225}{{\tt arXiv:1206.1225}}].

\bibitem{pearson-zahn}
R.~{Pearson} and O.~{Zahn}, {\it {Cosmology from cross correlation of CMB
  lensing and galaxy surveys}},  {\em \prd} {\bf 89} (Feb., 2014) 043516,
  [\href{http://arxiv.org/abs/1311.0905}{{\tt arXiv:1311.0905}}].

\bibitem{Killedar11}
M.~{Killedar}, P.~D. {Lasky}, G.~F. {Lewis}, and C.~J. {Fluke}, {\it
  {Gravitational lensing with three-dimensional ray tracing}},  {\em \mnras}
  {\bf 420} (Feb., 2012) 155--169, [\href{http://arxiv.org/abs/1110.4894}{{\tt
  arXiv:1110.4894}}].

\bibitem{Hilbert07a}
S.~{Hilbert}, R.~B. {Metcalf}, and S.~D.~M. {White}, {\it {Imaging the cosmic
  matter distribution using gravitational lensing of pre-galactic HI}},  {\em
  \mnras} {\bf 382} (Dec., 2007) 1494--1502,
  [\href{http://arxiv.org/abs/0706.0849}{{\tt arXiv:0706.0849}}].

\bibitem{Couchman99}
H.~M.~P. {Couchman}, A.~J. {Barber}, and P.~A. {Thomas}, {\it {Measuring the
  three-dimensional shear from simulation data, with applications to weak
  gravitational lensing}},  {\em \mnras} {\bf 308} (Sept., 1999) 180--200,
  [\href{http://arxiv.org/abs/astro-ph/9810063}{{\tt astro-ph/9810063}}].

\bibitem{Carbone09}
C.~{Carbone}, C.~{Baccigalupi}, M.~{Bartelmann}, S.~{Matarrese}, and
  V.~{Springel}, {\it {Lensed CMB temperature and polarization maps from the
  Millennium Simulation}},  {\em \mnras} {\bf 396} (June, 2009) 668--679,
  [\href{http://arxiv.org/abs/0810.4145}{{\tt arXiv:0810.4145}}].

\bibitem{Carbone13}
C.~{Carbone}, M.~{Baldi}, V.~{Pettorino}, and C.~{Baccigalupi}, {\it {Maps of
  CMB lensing deflection from N-body simulations in Coupled Dark Energy
  Cosmologies}},  {\em \jcap} {\bf 9} (Sept., 2013) 4,
  [\href{http://arxiv.org/abs/1305.0829}{{\tt arXiv:1305.0829}}].

\bibitem{antolini}
C.~{Antolini}, Y.~{Fantaye}, M.~{Martinelli}, C.~{Carbone}, and
  C.~{Baccigalupi}, {\it {N-body lensed CMB maps: lensing extraction and
  characterization}},  {\em \jcap} {\bf 2} (Feb., 2014) 39,
  [\href{http://arxiv.org/abs/1311.7112}{{\tt arXiv:1311.7112}}].

\bibitem{hu-okamoto}
W.~{Hu} and T.~{Okamoto}, {\it {Mass Reconstruction with Cosmic Microwave
  Background Polarization}},  {\em \apj} {\bf 574} (Aug., 2002) 566--574,
  [\href{http://arxiv.org/abs/astro-ph/0111606}{{\tt astro-ph/0111606}}].

\bibitem{BlandfordNarayan}
R.~{Blandford} and R.~{Narayan}, {\it {Fermat's principle, caustics, and the
  classification of gravitational lens images}},  {\em \apj} {\bf 310} (Nov.,
  1986) 568--582.

\bibitem{Jain00}
B.~{Jain}, U.~{Seljak}, and S.~{White}, {\it {Ray-tracing Simulations of Weak
  Lensing by Large-Scale Structure}},  {\em \apj} {\bf 530} (Feb., 2000)
  547--577, [\href{http://arxiv.org/abs/astro-ph/9901191}{{\tt
  astro-ph/9901191}}].

\bibitem{Pace07}
F.~{Pace}, M.~{Maturi}, M.~{Meneghetti}, M.~{Bartelmann}, L.~{Moscardini}, and
  K.~{Dolag}, {\it {Testing the reliability of weak lensing cluster
  detections}},  {\em \aap} {\bf 471} (Sept., 2007) 731--742,
  [\href{http://arxiv.org/abs/astro-ph/0702031}{{\tt astro-ph/0702031}}].

\bibitem{Hilbert09}
S.~{Hilbert}, J.~{Hartlap}, S.~D.~M. {White}, and P.~{Schneider}, {\it
  {Ray-tracing through the Millennium Simulation: Born corrections and
  lens-lens coupling in cosmic shear and galaxy-galaxy lensing}},  {\em \aap}
  {\bf 499} (May, 2009) 31--43, [\href{http://arxiv.org/abs/0809.5035}{{\tt
  arXiv:0809.5035}}].

\bibitem{Becker12}
M.~R. {Becker}, {\it {CALCLENS: weak lensing simulations for large-area sky
  surveys and second-order effects in cosmic shear power spectra}},  {\em
  \mnras} {\bf 435} (Oct., 2013) 115--132.

\bibitem{DasBode}
S.~{Das} and P.~{Bode}, {\it {A Large Sky Simulation of the Gravitational
  Lensing of the Cosmic Microwave Background}},  {\em \apj} {\bf 682} (July,
  2008) 1--13, [\href{http://arxiv.org/abs/0711.3793}{{\tt arXiv:0711.3793}}].

\bibitem{Fosalba08}
P.~{Fosalba}, E.~{Gazta{\~n}aga}, F.~J. {Castander}, and M.~{Manera}, {\it {The
  onion universe: all sky lightcone simulations in spherical shells}},  {\em
  \mnras} {\bf 391} (Nov., 2008) 435--446,
  [\href{http://arxiv.org/abs/0711.1540}{{\tt arXiv:0711.1540}}].

\bibitem{Bartelmann10}
M.~{Bartelmann}, {\it {TOPICAL REVIEW Gravitational lensing}},  {\em Classical
  and Quantum Gravity} {\bf 27} (Dec., 2010) 233001,
  [\href{http://arxiv.org/abs/1010.3829}{{\tt arXiv:1010.3829}}].

\bibitem{ValeWhite03}
C.~{Vale} and M.~{White}, {\it {Simulating Weak Lensing by Large-Scale
  Structure}},  {\em \apj} {\bf 592} (Aug., 2003) 699--709,
  [\href{http://arxiv.org/abs/astro-ph/0303555}{{\tt astro-ph/0303555}}].

\bibitem{Li2010}
B.~{Li}, L.~J. {King}, G.-B. {Zhao}, and H.~{Zhao}, {\it {An analytic
  ray-tracing algorithm for weak lensing}},  {\em \mnras} {\bf 415} (July,
  2011) 881--892, [\href{http://arxiv.org/abs/1012.1625}{{\tt
  arXiv:1012.1625}}].

\bibitem{Hu2000}
W.~{Hu}, {\it {Weak lensing of the CMB: A harmonic approach}},  {\em \prd} {\bf
  62} (Aug., 2000) 043007, [\href{http://arxiv.org/abs/astro-ph/0001303}{{\tt
  astro-ph/0001303}}].

\bibitem{Baldi12}
M.~{Baldi}, {\it {The CoDECS project: a publicly available suite of
  cosmological N-body simulations for interacting dark energy models}},  {\em
  \mnras} {\bf 422} (May, 2012) 1028--1044,
  [\href{http://arxiv.org/abs/1109.5695}{{\tt arXiv:1109.5695}}].

\bibitem{Carbone08}
C.~{Carbone}, V.~{Springel}, C.~{Baccigalupi}, M.~{Bartelmann}, and
  S.~{Matarrese}, {\it {Full-sky maps for gravitational lensing of the cosmic
  microwave background}},  {\em \mnras} {\bf 388} (Aug., 2008) 1618--1626,
  [\href{http://arxiv.org/abs/0711.2655}{{\tt arXiv:0711.2655}}].

\bibitem{Springel01}
V.~{Springel}, M.~{White}, and L.~{Hernquist}, {\it {Hydrodynamic Simulations
  of the Sunyaev-Zeldovich Effect(s)}},  {\em \apj} {\bf 549} (Mar., 2001)
  681--687, [\href{http://arxiv.org/abs/astro-ph/0008133}{{\tt
  astro-ph/0008133}}].

\bibitem{Watson13}
W.~A. {Watson}, J.~M. {Diego}, S.~{Gottl{\"o}ber}, I.~T. {Iliev}, A.~{Knebe},
  E.~{Mart{\'{\i}}nez-Gonz{\'a}lez}, G.~{Yepes}, R.~B. {Barreiro},
  J.~{Gonz{\'a}lez-Nuevo}, S.~{Hotchkiss}, A.~{Marcos-Caballero},
  S.~{Nadathur}, and P.~{Vielva}, {\it {The Jubilee ISW project - I. Simulated
  ISW and weak lensing maps and initial power spectra results}},  {\em \mnras}
  {\bf 438} (Feb., 2014) 412--425, [\href{http://arxiv.org/abs/1307.1712}{{\tt
  arXiv:1307.1712}}].

\bibitem{Lewis05}
A.~{Lewis}, {\it {Lensed CMB simulation and parameter estimation}},  {\em \prd}
  {\bf 71} (Apr., 2005) 083008,
  [\href{http://arxiv.org/abs/astro-ph/0502469}{{\tt astro-ph/0502469}}].

\bibitem{Fabbian13}
G.~{Fabbian} and R.~{Stompor}, {\it {High-precision simulations of the weak
  lensing effect on cosmic microwave background polarization}},  {\em \aap}
  {\bf 556} (Aug., 2013) A109, [\href{http://arxiv.org/abs/1303.6550}{{\tt
  arXiv:1303.6550}}].

\bibitem{Kayo01}
I.~{Kayo}, A.~{Taruya}, and Y.~{Suto}, {\it {Probability Distribution Function
  of Cosmological Density Fluctuations from a Gaussian Initial Condition:
  Comparison of One-Point and Two-Point Lognormal Model Predictions with N-Body
  Simulations}},  {\em \apj} {\bf 561} (Nov., 2001) 22--34,
  [\href{http://arxiv.org/abs/astro-ph/0105218}{{\tt astro-ph/0105218}}].

\bibitem{Taruya02}
A.~{Taruya}, M.~{Takada}, T.~{Hamana}, I.~{Kayo}, and T.~{Futamase}, {\it
  {Lognormal Property of Weak-Lensing Fields}},  {\em \apj} {\bf 571} (June,
  2002) 638--653, [\href{http://arxiv.org/abs/astro-ph/0202090}{{\tt
  astro-ph/0202090}}].

\bibitem{Takahashi12}
R.~{Takahashi}, M.~{Sato}, T.~{Nishimichi}, A.~{Taruya}, and M.~{Oguri}, {\it
  {Revising the Halofit Model for the Nonlinear Matter Power Spectrum}},  {\em
  \apj} {\bf 761} (Dec., 2012) 152, [\href{http://arxiv.org/abs/1208.2701}{{\tt
  arXiv:1208.2701}}].

\bibitem{Smith03}
R.~E. {Smith}, J.~A. {Peacock}, A.~{Jenkins}, S.~D.~M. {White}, C.~S. {Frenk},
  F.~R. {Pearce}, P.~A. {Thomas}, G.~{Efstathiou}, and H.~M.~P. {Couchman},
  {\it {Stable clustering, the halo model and non-linear cosmological power
  spectra}},  {\em \mnras} {\bf 341} (June, 2003) 1311--1332,
  [\href{http://arxiv.org/abs/astro-ph/0207664}{{\tt astro-ph/0207664}}].

\bibitem{benoit-levy12}
A.~{Benoit-L{\'e}vy}, K.~M. {Smith}, and W.~{Hu}, {\it {Non-Gaussian structure
  of the lensed CMB power spectra covariance matrix}},  {\em \prd} {\bf 86}
  (Dec., 2012) 123008, [\href{http://arxiv.org/abs/1205.0474}{{\tt
  arXiv:1205.0474}}].

\bibitem{KSmith04}
K.~M. {Smith}, W.~{Hu}, and M.~{Kaplinghat}, {\it {Weak lensing of the CMB:
  Sampling errors on B modes}},  {\em \prd} {\bf 70} (Aug., 2004) 043002,
  [\href{http://arxiv.org/abs/astro-ph/0402442}{{\tt astro-ph/0402442}}].

\bibitem{Hagstotz14}
S.~{Hagstotz}, B.~M. {Sch{\"a}fer}, and P.~M. {Merkel}, {\it {Born-corrections
  to weak lensing of the cosmic microwave background temperature and
  polarisation anisotropies}},  {\em ArXiv e-prints} (Oct., 2014)
  [\href{http://arxiv.org/abs/1410.8452}{{\tt arXiv:1410.8452}}].

\bibitem{CoorayHu02}
A.~{Cooray} and W.~{Hu}, {\it {Second-Order Corrections to Weak Lensing by
  Large-Scale Structure}},  {\em \apj} {\bf 574} (July, 2002) 19--23,
  [\href{http://arxiv.org/abs/astro-ph/0202411}{{\tt astro-ph/0202411}}].

\bibitem{ShapiroCooray06}
C.~{Shapiro} and A.~{Cooray}, {\it {The Born and lens lens corrections to weak
  gravitational lensing angular power spectra}},  {\em \jcap} {\bf 3} (Mar.,
  2006) 7, [\href{http://arxiv.org/abs/astro-ph/0601226}{{\tt
  astro-ph/0601226}}].

\bibitem{PeacockDodds96}
J.~A. {Peacock} and S.~J. {Dodds}, {\it {Non-linear evolution of cosmological
  power spectra}},  {\em \mnras} {\bf 280} (June, 1996) L19--L26,
  [\href{http://arxiv.org/abs/astro-ph/9603031}{{\tt astro-ph/9603031}}].

\end{thebibliography}\endgroup
