CMB weak-lensing beyond the Born approximation: a numerical approach

TL;DR

This paper investigates the effects of beyond-Born gravitational lensing on CMB anisotropies using advanced numerical simulations, revealing small but significant corrections to lensing observables and polarization rotation, with implications for future experiments.

Contribution

It introduces a comprehensive numerical approach to quantify beyond-Born lensing effects on the CMB, extending previous analytical predictions to smaller scales and including polarization rotation.

Findings

Good agreement with analytical predictions for lensing corrections.

Beyond-Born effects are detectable in temperature spectra with foreground control.

Polarization corrections are challenging to measure with upcoming experiments.

Abstract

We perform a complete study of the gravitational lensing effect beyond the Born approximation on the Cosmic Microwave Background (CMB) anisotropies using a multiple-lens raytracing technique through cosmological N-body simulations of the DEMNUni suite. The impact of second order effects accounting for the non-linear evolution of large-scale structures is evaluated propagating for the first time the full CMB lensing jacobian together with the light rays trajectory. We carefully investigate the robustness of our approach against several numerical effects in the raytracing procedure and in the N-body simulation itself, and find no evidence of large contaminations. We discuss the impact of beyond-Born corrections on lensed CMB observables, and compare our results with recent analytical predictions that appeared in the literature, finding a good agreement, and extend these results to smaller angular scales. We measure the gravitationally-induced CMB polarization rotation that appears at second order, and compare this result with the latest analytical predictions. We then present the detection prospect of beyond-Born effects with the future CMB-S4 experiment. We show that corrections to the temperature power spectrum can be measured only if a good control of the extragalactic foregrounds is achieved. Conversely, the beyond-Born corrections on E and B-modes power spectra will be much more difficult to detect.

Tables

Table 1. Table 1 : Moments of the first-order and beyond-Born convergence and rotation maps (together with moments of the corresponding shot-noise maps) obtained filtering modes ℓ > 8000 ℓ 8000 \ell>8000 .

Map	Variance (${10}^{-3}$)	Skewness	Kurtosis
${\kappa }^{1st}$	8.465	0.370	0.554
${\kappa }^{\mathrm{ML}}$	8.472	0.266	0.450
$\omega$	$6.87\cdot {10}^{-3}$	$3.184\cdot {10}^{-3}$	3.185
${\kappa }^{\mathrm{ML}}$ shot-noise	1.070	$1.289\cdot {10}^{-3}$	$-3.865\cdot {10}^{-4}$
$\omega$ shot-noise	$6.357\cdot {10}^{-5}$	$6.570\cdot {10}^{-3}$	$1.337$

Equations

$${\beta }_i(\mathbf{\theta },\chi )={\theta }_i-\frac{2}{c^2}{\int }_0^{\chi }\frac{f_K(\chi -{\chi }^{\prime })}{f_K(\chi )f_K({\chi }^{\prime })}\Psi {,}_{{\beta }_i}\left(\mathbf{\beta }(\mathbf{\theta },{\chi }^{\prime }),{\chi }^{\prime }\right)\mathrm{d}{\chi }^{\prime },$$

$$\begin{array}{rl}{A}_{ij}(\mathbf{\theta },\chi )& \equiv \frac{\partial {\beta }_i(\mathbf{\theta },\chi )}{\partial {\theta }_j}=\\ & ={\delta }_{ij}^K-\frac{2}{c^2}{\int }_0^{\chi }\frac{f_K(\chi -{\chi }^{\prime })}{f_K(\chi )f_K({\chi }^{\prime })}\Psi {,}_{{\beta }_i{\beta }_k}\left(\mathbf{\beta }(\mathbf{\theta },{\chi }^{\prime }),{\chi }^{\prime }\right)A_{kj}(\mathbf{\theta },{\chi }^{\prime })\mathrm{d}{\chi }^{\prime },\end{array}$$

$$A_{ij}\equiv \left(\begin{array}{cc}\cos \omega & \sin \omega \\ -\sin \omega & \cos \omega \end{array}\right)\left(\begin{array}{cc}1-\kappa -{\gamma }_1& -{\gamma }_2\\ -{\gamma }_2& 1-\kappa +{\gamma }_1\end{array}\right)\approx \left(\begin{array}{cc}1-\kappa -{\gamma }_1& -{\gamma }_2+\omega \\ -{\gamma }_2-\omega & 1-\kappa +{\gamma }_1\end{array}\right),$$

$$\mathbf{\beta }(\mathbf{\theta },{\chi }_s)=\mathbf{\theta }-\nabla \psi -\nabla \times \Omega ,$$

$${\mathbf{\alpha }}^{(k)}=\frac{2}{c^2{f}_K({\chi }_k)}{\int }_{{\chi }_k-\Delta \chi /2}^{{\chi }_k+\Delta \chi /2}{\nabla }_{\beta }\Psi (\mathbf{\beta }(\mathbf{\theta },\tilde{\chi }),\tilde{\chi })d\tilde{\chi }={\nabla }_{\beta }{\Phi }^{(k)}(\mathbf{\beta }(\mathbf{\theta },{\chi }_k),{\chi }_k),$$

$${\Phi }^{(k)}(\mathbf{\beta }(\mathbf{\theta },{\chi }_k),{\chi }_k)=\frac{2}{c^2{f}_K({\chi }_k)}{\int }_{{\chi }_k-\Delta \chi /2}^{{\chi }_k+\Delta \chi /2}\Psi (\mathbf{\beta }(\mathbf{\theta },\tilde{\chi }),\tilde{\chi })d\tilde{\chi }.$$

$$U_{ij}^{(k)}=\frac{{\partial }^2{\Psi }^{(k)}({\mathbf{\beta }}^{(k)})}{\partial {\beta }_i^{(k)}\partial {\beta }_j^{(k)}}=\frac{\partial {\alpha }_i^{(k)}({\mathbf{\beta }}^{(k)})}{\partial {\beta }_j^{(k)}},$$

$${\nabla }_{\hat{\mathbf{n}}}^2{\Phi }^{(k)}({\mathbf{\beta }}^{(k)})=2K^{(k)}({\mathbf{\beta }}^{(k)}),$$

$$K^{(k)}({\mathbf{\beta }}^{(k)})=\frac{4\pi G}{c^2}\frac{D_A({\chi }_k)}{(1+z_k{)}^2}{\Delta }_{\Sigma }^{(k)}({\mathbf{\beta }}^{(k)}).$$

$$\mathbf{\beta }(\mathbf{\theta },{\chi }_s)=\mathbf{\theta }-\sum _{i=0}^{N-1}\frac{f_K({\chi }_s-{\chi }_k)}{f_K({\chi }_s)}{\mathbf{\alpha }}^{(k)}({\mathbf{\beta }}^{(k)}).$$

$$A_{ij}^N(\mathbf{\theta },{\chi }_N)={\delta }_{ij}^K-\sum _{k=0}^{N-1}\frac{D_{k,N}}{D_N}{U}_{ip}^{(k)}({\mathbf{\beta }}^{(k)},{\chi }_k)A_{pj}^{(k)}(\mathbf{\theta },{\chi }_k),$$

$${\mathbf{\beta }}^{(k)}=\left(1-\frac{D_{k-1}}{D_k}\frac{D_{k-2,k}}{D_{k-2,k-1}}\right){\mathbf{\beta }}^{(k-2)}+\frac{D_{k-1}}{D_k}\frac{D_{k-2,k}}{D_{k-2,k-1}}{\mathbf{\beta }}^{(k-1)}-\frac{D_{k-1,k}}{D_k}{\mathbf{\alpha }}^{(k-1)}({\mathbf{\beta }}^{(k-1)}),$$

$$A_{ij}^{(k)}=\left(1-\frac{D_{k-1}}{D_k}\frac{D_{k-2,k}}{D_{k-2,k-1}}\right)A_{ij}^{(k-2)}+\frac{D_{k-1}}{D_k}\frac{D_{k-2,k}}{D_{k-2,k-1}}{A}_{ij}^{(k-1)}-\frac{D_{k-1,k}}{D_k}{U}_{ip}^{(k-1)}{A}_{pj}^{(k-1)}.$$

$$A_{ij}^{(N),1st}(\mathbf{\theta },{\chi }_s)={\delta }_{ij}^K-\sum _{k=0}^{N-1}\frac{D_{k,N}}{D_N}{U}_{ij}^{(k)}(\mathbf{\theta },{\chi }_k).$$

$$\{{\Omega }_{dm},{\Omega }_b,{\Omega }_{\Lambda },n_s,{\sigma }_8,H_0\}=\{0.27,0.05,0.68,0.96,0.83,67Km/s/Mpc\}.$$

$${}_1{\alpha }_{\ell m}^E{}^{(k)}=\sqrt{\ell (\ell +1)}{\psi }_{\ell m}^{(k)}{}_1{\alpha }_{\ell m}^B{}^{(k)}=0.$$

$${\mathbf{\beta }}^{(k)}={\mathbf{\beta }}^{(k-1)}+{\mathbf{d}}^{(k)}$$

$${\mathbf{\beta }}^N=\mathbf{\theta }-{\mathbf{d}}^{\mathrm{eff}}$$

$${\mathbf{d}}_{\ell m}^{\mathrm{eff},\mathrm{E}}=\sqrt{\ell (\ell +1)}{\psi }_{\ell m}^{\mathrm{eff}}{\mathbf{d}}_{\ell m}^{\mathrm{eff},\mathrm{B}}=\sqrt{\ell (\ell +1)}{\Omega }_{\ell m}^{\mathrm{eff}},$$

$$C_{\ell }^{\kappa ϵ}$$

$$C_{\ell }^{\omega \beta }$$

$$C_{\ell }^{\kappa ϵ}$$

$$C_{\ell }^{\omega \beta }$$

$$C_{\ell }^{\kappa \kappa ,\text{Shot-noise}(k)}=\frac{9H_0^4{\Omega }_{m,0}^2}{4c^4}\delta \chi \frac{1}{{\bar{n}}_k}{\left(\frac{{\chi }_s-{\chi }_k}{{\chi }_s{\chi }_k}\right)}^2$$

$$\Delta {\sigma }_{\kappa }^{2,\mathrm{Theory}}\equiv {\left({\sigma }_{{\kappa }^{\mathrm{ML}}}^2-{\sigma }_{{\kappa }^{1st}}^2\right)}^{\mathrm{Theory}}=\sum _{\ell =0}^{{\ell }_{\max }}\frac{(2\ell +1)}{4\pi }\Delta C_{\ell }^{\kappa \kappa }.$$

$${}_1{d}_{\ell m}^{\mathrm{E},\mathrm{eff}}=-2\frac{{\kappa }_{\ell m}^{\mathrm{ML}}}{\sqrt{\ell (\ell +1)}}{}_1{d}_{\ell m}^{\mathrm{B},\mathrm{eff}}=-2\frac{{\omega }_{\ell m}}{\sqrt{\ell (\ell +1)}}.$$

$$\tilde{P}(\mathbf{\theta })=e^{-2i{\beta }^{\mathrm{rotation}}}\left[e^{-2i{\delta }^{\mathrm{transport}}}P(\mathbf{\theta }-{\mathbf{d}}^{\mathrm{eff}})\right].$$

$$\tilde{B}(\mathbf{\ell }{)}^{d^{\mathrm{eff}}}$$

$${\Delta }^{\psi \Omega }(\mathbf{\ell })$$

$$C_{\ell }^{BB,\mathrm{corr}}$$

$$C_{\ell }^{BB,\mathrm{mixed}}\equiv C_{\ell }^{BB,\psi \psi \mathrm{\Omega }}$$

$$C_{\ell }^{BB,\psi \Omega \Omega }$$

$$\begin{array}{rl}{A}_{ij}& =(1-\kappa ){\delta }_{ij}^K-{\gamma }_{ij}+\omega {ϵ}_{ij}\\ & ={\delta }_{ij}^K-\sum _{\ell m}\left({\psi }_{\ell m}{\nabla }_i{\nabla }_j{Y}_{\ell m}+{\Omega }_{\ell m}{ϵ}_j^k{\nabla }_i{\nabla }_k{Y}_{\ell m}\right),\end{array}$$

Full text

1 Introduction

The Cosmic Microwave Background (CMB) is now firmly established as one of the most important cosmological probe, and has been spectacularly exploited to very high sensitivity and accuracy by the recent Planck satellite mission [1, 2] and current generation of suborbital experiments. In parallel to the studies on the primordial signal for cosmological applications, the attention of the community has swiftly moved towards the effects imprinted at later stages of the Universe evolution on the CMB, the so-called secondary anisotropies. The weak gravitational lensing of CMB anisotropies stands as the primary example of such group and consists in a modification of the geodesic path of the CMB light coming from the last scattering surface induced by growing matter inhomogeneities.

The first robust detections of such effect have been achieved using CMB temperature data only by ACT [3] and SPT [4, 5], and later confirmed by Planck with a significance greater than 25$\sigma$ [6]. Only recently, however, the evidence of lensing was detected for the first time in polarization data by POLARBEAR [7, 8], SPTpol [9], ACTpol [10] and recently BICEP2/Keck Array [11]. Being sensitive to the whole matter distribution along the line of sight, CMB lensing can be used for cosmological analysis to infer information about the Large Scale Structure (LSS) distribution and thus on the parameters that govern the physics of structure formation at intermediate and late time like (e.g. the dark energy (DE) and massive neutrinos properties). Direct measurements of CMB lensing can be improved and complemented by the cross correlation analysis with observations of the actual lenses in LSS surveys as independent tracers of the matter distribution. This approach has already been exploited to obtained astrophysical and cosmological information [12, 13, 14, 15, 9, 16, 17, 18]. However, a major improvement is expected for this field in the forthcoming years when the next generation of high sensitivity CMB polarization experiments (Simons Array [19], AdvACTpol [20], SPT-3G [21] and ultimately CMB-S4 [22]) together with the next generation of galaxy surveys DESI [23], LSST [24] and the ESA Euclid satellite [25] will start observing the sky.

The exponentially growing quality of CMB and galaxy surveys data requires a great effort to simulate and make theoretical predictions for the cosmological observables with the highest possible accuracy. In particular it will be crucial to take into account second-order effects and non-linear evolution of the large-scale structures at the same time. Several recent theoretical works [26, 27, 28, 29, 30] studied the impact of the relaxation of the Born approximation to CMB lensing and lensed CMB anisotropies adopting the perturbative approach to the lens equation presented in [31, 32, 33]. Similar studies have also been performed in the context of weak lensing (see e.g. [34, 35]). Analytical and numerical studies on the impact of the non-linear evolution of the matter distribution on several aspect of CMB lensing have also made significant progresses [29, 36, 37, 38]. At the numerical level, ray-tracing through large, high-resolution N-Body numerical simulations is still the best tool we have to analyze the signal up to the full non-linear scales. Despite the problem can be solved exactly [39], even though on limited sky fractions, a cheaper yet accurate and popular approach consists in extracting the lensing observables by photon ray-tracing along “unperturbed”, i.e. undeflected light paths in the so-called Born approximation (e.g. [40, 41, 42]). Approximated methods based on halo model formalism have also been recently proposed in this context [43]. In particular, [44, 45] applied techniques based on the Born approximation to study a set of N-Body simulation with different cosmologies and dark energy and/or massive neutrinos scenarios and investigate the variation of the lensing pattern with respect to the standard $\Lambda$CDM model on the full sky.

However, when facing the challenge of producing accurate and realistic simulations of the lensing effect, we must take into account that each light ray undergoes several deflections due to matter inhomogeneities. In this context, we should therefore replace the modeling through a single effective deflection adopted in the Born approximation by a multiple-lens (ML) approach. In the latter case, large volumes of matter are projected onto a series of lens planes [46, 47, 48, 49, 50, 51, 52, 53] so that the continuous deflection experienced by a light ray is approximated by finite deflections at each of the mass planes. Faster approaches based on A ML algorithm for CMB lensing application were first sketched in [54]. In [55] we presented the first generalized implementation of this algorithm through large N-body simulations and applied it to study the lensing effect on both CMB temperature and polarization anisotropies.

Motivated by the recent theoretical results and to address concerns on the capability of our numerical setup to resolve beyond Born corrections in CMB lensing expressed in [26], in this work we expanded and improved the algorithm of [55]. In particular we modified the ray-tracing code to explore and characterize the full set of lensing observables derived with the ML method up to arcminute scales, where corrections and impact of non-linear effects are most noticeable and beyond the limit of [55].

The paper is organized as follows. In Sec. 2 we introduce the theoretical background and notation used for our lensing algorithm while in Sec. 3 we discusse our ray-tracing technique emphasizing the improvements made from the previous version of the algorithm. In Sec. 4 we discusse the properties of the lensing observables extracted from our simulation with particular emphasis on stability, accuracy and reliability of the signal for different simulation setups. In Sec. 5 and 6 we compare the results of our simulations with analytical predictions for the lensing observables beyond Born approximation as well as their impact on lensed CMB observables and discuss perspective for the measurements of these corrections in future data sets. The last section draws the conclusions.

2 Theory

2.1 Gravitational light deflection

In this Section we briefly introduce the relevant quantities of weak lensing used in the rest of this paper. A detailed description of the weak lensing formalism is summarized in several reviews [56, 57] and recent articles [47, 49, 54, 50, 55] therefore we refer the reader to these papers for further details.

In weak lensing formalism, the effect of deflections of light rays along the entire line of sight is described by the lens equation, which maps the final position $(t,\boldsymbol{\beta},\chi)$ of the ray to the position of its source $\boldsymbol{\theta}$, i.e.

[TABLE]

where $\Psi(t,\boldsymbol{\beta},\chi)$ is a gravitational potential located on the photon path and $,_{\beta_{i}}$ its spatial derivative with respect to the photon’s position $\boldsymbol{\beta}$. Here $f_{K}(\chi)$ is the standard angular diameter distance for a universe with curvature $K$. The relative position of nearby light rays is quantified by the derivative of the equation above

[TABLE]

where $\delta^{K}_{ij}$ is the Kronecker delta. The image distortions of light sources are described by the magnification matrix111We will refer to the A matrix as magnification matrix, lensing jacobian or distortion tensor as synonyms in the following. ${\bf A}(\boldsymbol{\theta},\chi)\equiv\left\{A_{ij}(\boldsymbol{\theta},\chi)\right\}$, which holds the information of the mapping induced by lensing between the original image and the one at the current position on the lens plane. Note that the gravitational potential $\Psi$ is evaluated at the ray angular position $\boldsymbol{\beta}(\boldsymbol{\theta},\chi)$, while the distortion itself - which is present at the r.h.s of Eq. (2.2) and describes the lens-lens coupling - is computed at the “background” position $\boldsymbol{\theta}$. The magnification matrix ${\bf A}$ is typically decomposed into four fields describing how the light rays coming from a source at $\chi\equiv\chi_{s}$ are transformed by the passage through the matter distribution,

[TABLE]

where we assumed that the image rotation angle $\omega(\boldsymbol{\theta},\chi_{s})$ which defines the rotation of the lensed image, is small, and we can work in the weak lensing regime. The field $\kappa(\boldsymbol{\theta},\chi_{s})$ is referred to as the convergence while $\gamma(\boldsymbol{\theta},\chi_{s})=\gamma_{1}(\boldsymbol{\theta},\chi_{s})+i\gamma_{2}(\boldsymbol{\theta},\chi_{s})$ defines the complex shear, describing the shearing of the image along the two orthogonal directions of the basis, and can be decomposed into a curl-free part, the shear E-modes $\gamma_{\epsilon}(\boldsymbol{\theta},\chi_{s})$ and a divergence-free one, the shear B-modes $\gamma_{\beta}(\boldsymbol{\theta},\chi_{s})$.

In the so-called Born approximation [58, 49], the lens equation can be integrated over the unperturbed, unlensed photon paths $\left(\boldsymbol{\theta},\chi\right)$, therefore dropping in Eqs. 2.1, 2.2 the dependence over the lensed position $\boldsymbol{\beta}$. While absent in first order weak lensing approximation, shear B-modes are generated if we take into account the full non-linear equation at higher order [59, 26]. Modern surveys [60] can detect the presence of B-modes in the shear field and although these are mainly used as a monitor for systematic effects in weak-lensing data, this signal could also be used for cosmological applications such as e.g. testing anisotropic cosmological models [61].

An additional assumption usually made is to consider negligible the coupling between lenses at two different redshifts. Unlike the correction to the Born approximation, the lens-lens coupling results in both B-mode generation and a net rotation of galaxy images, which has also been measured in numerical simulations [47].

The presence of these effects authorizes the introduction of an auxiliary curl potential $\Omega(\boldsymbol{\theta})^{\rm eff}$, such that the deflection angle may be expressed as a combination of the standard gradient contribution from a scalar field $\psi(\boldsymbol{\theta})^{\rm eff}$, and a curl contribution [59]:

[TABLE]

where we defined the two-dimensional curl $(\nabla\times\Omega)_{i}=\epsilon_{ij}\partial\Omega$ and dropped the spatial dependence of the potentials for sake of clarity. This potential $\Omega^{\rm eff}$ has to be intended as an “effective” Born-like potential, integrated along the line of sight, that encodes all the information about the rotation of the image as curl-like patterns are originated in the signal through multiple deflections. All the quantities above can be treated on the sphere using the spin-$s$ spherical harmonic decomposition of the full-sky [62, 50]. In addition, combining results from [62] and [63], we can derive consistency relations between the components of the magnification matrix and the effective potentials (see Appendix A for further details).

2.2 Multiple-lens-plane approach

Following [54, 55] we note that the Eqs. (2.1),(2.2) can be discretized by dividing the interval between the observer and the source into concentric $N$ spherical shells, each of comoving thickness $\Delta\chi$, denoted by $\chi_{k}$ as the comoving distance to the middle of the $k$-th shell and its related redshift $z_{k}$. A photon incoming on the $k$-th shell at $\chi_{k}$ is deflected due to the presence of matter by an angle $\boldsymbol{\alpha}^{(k)}$, which can be approximated by

[TABLE]

where we have defined the 2-D gravitational potential on the sphere as

[TABLE]

Here, the notation $(\boldsymbol{\beta}(\boldsymbol{\theta},\chi_{k}),\chi_{k})$ means that the potential is evaluated at the conformal look-back time $\chi_{k}$, when the photon, coming from a source at distance $\chi_{s}$ from the observer at position $\boldsymbol{\theta}$ on the celestial sphere, was at the position $\boldsymbol{\beta}^{(k)}$. The second derivatives can be combined into the shear matrix ${\bf U}$:

[TABLE]

being $\boldsymbol{\alpha}^{(k)}$ the lensing angle for the $k$-th shell. In case of the full-sky analysis the partial derivative operators have to be promoted to covariant derivatives [50]. The lensing potential for each matter shell $k$ is the solution of the Poisson equation, i.e.

[TABLE]

where the convergence field $K^{(k)}$ at the $k$-th shell is

[TABLE]

and $\Delta^{(k)}_{\Sigma}(\hat{{\bf n}})$ is the (projected) surface mass overdensity, as in [54, 55]. In Eq. (2.8) we dropped the term containing the derivatives in the radial direction, ignoring thus long wavelength fluctuations along the line-of-sight via the Limber approximation [47, 50]. As noted in [55], the results of this approximation are particularly evident if we look at the angular power spectrum for the lensing potential of each single matter shell. However, if we look at the overall effect after a sufficiently long photon path, the 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 [54, 64, 50, 55]. The lensing potential on the sphere is related to $K^{(k)}$ via Eq. (2.8), and it can be easily computed by expanding both sides of the Poisson equation in spherical harmonics. The quantity $K^{(k)}$ is directly computed when the matter distribution in the shell is radially projected onto the spherical map, as explained in Sec. 3.2.

Given the deflection angle at each lens-plane, we can trace back the light coming from a source at position $(\boldsymbol{\theta},\chi_{s})$ to the observer after $N$ deflections:

[TABLE]

We can easily discretize the Eq. (2.2) as

[TABLE]

where we defined for simplicity $D_{k,N}=f_{K}(\chi_{N}-\chi_{k})$, $D_{k}=f_{K}(\chi_{k})$ while $N$ is the number of planes necessary to reach the source at comoving distance $\chi_{N}$. The method in Eq. (2.11) becomes computationally unfeasible very quickly, especially when we have a large number of lens planes covering a wide sky fraction. It, in fact, requires that for each $k\textrm{-th}$ iteration all the information of the $k-1$ deflections is kept. This becomes particularly problematic in the case of CMB where the source plane is located at very high redshift and at least 50 (or more) iterations are required to model the path of CMB photons. [65, 49] proposed a more efficient method that requires the combination of only two previous lens-planes instead of the whole set as Eq. (2.11). The method has been validated on the full-sky in [50]. The angular position $\boldsymbol{\beta}^{(k)}$ at the $k$-th shell is a function of its two previous positions $\boldsymbol{\beta}^{(k-2)}$ and $\boldsymbol{\beta}^{(k-1)}$ as

[TABLE]

and, by differentiating with respect to $\boldsymbol{\theta}$ as in Eq. (2.2), we obtain the recurrence relation for the magnification matrix [49, 50] as well:

[TABLE]

These relations require fewer arithmetic operations and memory usage than the standard discretization of Eq. (2.11), therefore allowing us to compute iteratively the magnification matrix for each light-rays from the observer to any source. In the following we will also make use of the so-called magnification matrix in the first order approximation [49] to assess the impact of the second-order effects and to distinguish those from numerical effects. In the multiple-lens formalism this takes the form

[TABLE]

3 The Algorithm

A detailed outline of the algorithm, including the construction of the past light-cone and the map-making procedure, has been given in [55]. In the following, we summarized the map-making procedure to produce lensing planes and the numerical and improvements specific to this work.

3.1 N-body simulation

The results reported in this were derived using the reference $\Lambda$CDM simulation belonging to the “Dark Energy and Massive Neutrino Universe” (DEMNUni) simulation project. We refer the reader to [66, 45] for a more detailed description of the simulations and for an extended discussion on the physical results issued by the project. The DEMNUni suite consists in a baseline $\Lambda$CDM model, and three modified $\Lambda$CDM cosmologies characterized by three active neutrinos222The simulations do not account for an effective neutrino number $N_{e}ff>3$. For example, neutrino isocurvature perturbations could produce larger $N_{eff}$ and thus affect galaxy and CMB power spectra [67]. However, these are excluded by present data (see e.g. [68]). with different values of their total mass $\sum m_{\nu}$. All these simulations share the same total matter density $\Omega_{m}$, as well as the same amplitude of primordial curvature perturbations and are based on a Planck-2013 cosmology [1]:

[TABLE]

The employed $\Lambda$CDM simulation follows the evolution of $2048^{3}$ CDM particles in a cubic comoving volume $($2 $h^{-1}\text{Gpc})^{3}$ from redshift $z=99$ to the present epoch (see e.g. Fig. 3). **Being $z=99$ the highest redshift probed by the simulation, the theoretical curves used for comparisons with numerical results in the following include contributions from matter perturbations $z\leq 99$. ** The mass resolution of the simulation at $z=0$ is $M_{\rm CDM}=8.27\times 10^{10}M_{\odot}/h$ and the gravitational softening length is set to $\epsilon_{s}=20$ kpc$/h$ corresponding to 0.04 times the mean linear inter-particle separation. The simulation was carried out using a modified TreePM version of GADGET-3333http://www.mpa-garching.mpg.de/gadget [69], specifically developed to include all the additional physical effects that characterize different neutrino models [70]. This version of GADGET-3 follows the evolution of CDM and neutrino particles, treating them as two distinct sets of collisionless particles. Moreover, to test the effect the N-Body resolution has on our results (see Sec. 4.4), we considered a Millennium-like numerical simulation, with the same background cosmology and amplitude of primordial curvature perturbations as described before. In this case the $2048^{3}$ CDM particles populate a box with comoving side of 500 Mpc$/h$ and a mass resolution, for a DM particle at $z=0$, of $M_{\rm CDM}=8.27\times 10^{8}M_{\odot}/h$.

3.2 Map making

Starting from different snapshots in time of the N-body simulation we reconstructed the full-sky past lightcone of the observer back to the maximum redshift available in the simulation (in our case $z_{\rm max}=99$). Because the size of the simulation box is limited, we need to replicate the box volume several times in space to fill the entire observable volume between the observer and $z_{\rm max}$. In order to avoid the repetition of the same structures along the line of sight and to mitigate the deficit of lensing power on angular scales comparable to the boxsize, we employed a specific randomization framework. We refer the reader to [71, 55] and references therein for a further discussion on this topic. Following [58, 42, 55], we divided the volume up to $z_{\rm max}$ into large spherical shells, each of thickness comparable to the simulation boxsize. All the simulation boxes falling into the same larger 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 3D matter shell. The coordinate transformation however changes from (large) shell to shell. We enforced the periodic boundary conditions when performing box randomization to guarantee a smooth matter distribution at the box edges and minimize numerical artifacts due to potential matter discontinuities at the boundaries. The particles that fall within the radius of the $k$-th spherical shell of width $\delta\chi$, ($\chi_{k}-\delta\chi/2,\chi_{k}+\delta\chi/2$) are considered part of the shell, where $\delta\chi$ is the chosen average comoving thickness of the spherical shells. For our baseline setup $\delta\chi=150\textrm{Mpc}/h$.

Following the scheme proposed in [54, 55], we then converted the position of a particle of mass $m$ belonging to a given matter shell into its angular position on the 2D sphere assuming a HEALPix 444http://healpix.sourceforge.net pixelization scheme. These maps are then projected into surface matter density maps $\Sigma^{\theta}$ assuming that in each pixel $p$ $\Sigma^{\theta}_{p}=\sum_{0}^{n}m_{i}/\Delta\Omega_{\rm pix}$, where $n$ is the number of particles per pixel and $\Delta\Omega_{\rm pix}$ its area in steradians. For the results described in this paper we used maps having an angular resolution of 52 arcsec corresponding to a NSIDE parameter $4096$.

3.2.1 Systematics tests

Similarly to what was done in [55], we performed several systematic test to assess the impact of possible systematic effect in the map-making procedure. The results of these tests are shown in Fig. 1. We checked that the amount of mass present in each mass sheet was consistent with the amount of total mass expected from the cosmology and found differences of sub percent level on almost all shells and no greater than 2% at the lowest redshifts bin. We also checked the distribution of the surface mass density and found consistency with previous numerical work. The mass distribution is in fact well described by Gaussian or log-normal probability density function (PDF) at high redshift with a tail in the ultra-nonlinear regime well fitted by the phenomenological Das and Ostriker model [72] with smoothing scale set to the variance of each map.

Once the surface mass density had been computed, we extracted the $k$-th shell surface mass density contrast $\Delta^{(k)}_{\Sigma^{\theta}}({\bf\hat{n}})=\Sigma^{\theta}({\bf\hat{n}})-\langle\Sigma^{\theta}\rangle$ and the convergence maps $K^{(k)}({\bf\hat{n}})$ (see Eq. (2.9)) that will be used in the following.

3.3 Raytracing and remapping

In order to propagate the CMB photons through the different shells we used in [55] a pixel-based approach exploiting the publicly available code LensPix [73]. For this work, we changed this part of the algorithm and adopted a modified version of the lenS2HAT code [74]. In particular we implemented the raytracing algorithm of Eq. (2.12) and Eq. (2.13) while carefully preserving the optimal data and workload distribution of the original algorithm.

The algorithm operates as follow. Starting from the convergence harmonic coefficients $K^{(k)}_{\ell m}$ it extracts the lensing potential of each shell solving the Poisson equation in the harmonic domain as $\psi^{(k)}_{\ell m}=-2K^{(k)}_{\ell m}/[\ell(\ell+1)]$. It then computes the deflection field for each lensing plane $\boldsymbol{\alpha}^{(k)}$ assuming the deflection field as a purely gradient field (i.e. a spin-1 curl-free vector field) having E and B decomposition

[TABLE]

Note that the lensing potential $\psi$ is simply the 2-D gravitational potential of Eq (2.6) multiplied by a geometrical lensing weight $f_{K}(\chi_{s}-\chi_{k})/f_{K}(\chi_{s})$. Once the deflection field is obtained, it remaps the CMB field as

[TABLE]

where $\boldsymbol{\beta}^{(k-1)}$ represents the position of an incoming CMB photons on the $k$-th lensing plane and ${\bf d}^{(k)}$ the total deflection field at each $k$-th lensing plane. This takes into account the deflection field $\boldsymbol{\alpha}^{(k)}$ and the incoming direction of the photon with respect to the normal to the surface $\boldsymbol{\xi}^{(k)}$ so that ${\bf d}^{(k)}=\boldsymbol{\alpha}^{(k)}+\boldsymbol{\xi}^{(k)}$. In the code, we adopted the algorithm of [54] to compute $\boldsymbol{\xi}^{(k)}$, which can be shown to be equivalent to the recurrence equation of Eq. (2.12)555We renamed $\boldsymbol{\xi}$ the $\boldsymbol{\beta}$ of [54] to avoid confusion with the notation adopted in the text so far.. In the case $k=1$, $\boldsymbol{\beta}^{(0)}=\boldsymbol{\theta}$ i.e. the unlensed CMB photon position. Following [50] we accounted for the change of the local tensor basis when parallel transporting the lensing jacobian along the geodesic connecting the old and displaced position at each deflection plane adopting the formalism of [75]. For sake of clarity we show in Fig. 2 the definition of the angles adopted in the algorithm.

Although algorithmically similar to LensPix, lenS2HAT offers several advantages when performing full-sky raytracing simulation for application that require the use of the highest possible angular resolution. The advantages lie mainly in the nearly perfect scalability and the very low memory overhead inherited from the S2HAT library666http://www.apc.univ-paris7.fr/APC_CS/Recherche/Adamis/MIDAS09/software/s2hat/s2hat.htm [76, 77, 78]. These allow to manipulate maps and compute spherical harmonics at very high resolution required by this specific application. The code is also very flexible as it works for multiple pixelization schemes beyond the commonly used HEALPix. In particular, for this work, we resampled the displacement field and performed the raytracing part on grids based on the Equidistant Cylindrical Projection (ECP) pixelization [79]. Thanks to its symmetry properties, this pixelization supports fast spherical harmonic transforms and nearly perfect quadrature for band-limited functions, minimizing thus aliasing effects that may become important in codes of this type.

lenS2HAT implements the raytracing using a simple yet powerful pixel-based lensing method adopting an efficient Nearest Grid Points (NGP) assignment scheme to evaluate the source plane (or, in this case the lensing plane) along the displaced direction. We thus assign to every deflected direction a value of the sky signal computed at the nearest grid points defined by the centers of a pixel of the assumed pixelization scheme. The NGP assignment is quick and allows to control the accuracy of the method and the signal smoothing scale through the parameters defining the grid resolution alone. Interpolation methods on the sphere, though helpful in cutting numerical costs, might modify the underlying signal and alter its statistical properties especially at the smallest scales, where we expect the effects of our interest to lie. We refer the reader to [74] for a further discussion on this topic together with an accurate benchmark of lenS2HAT.

3.4 Derivatives computation on the sphere

In order to propagate the lensing jacobian together with the displaced photon direction it is necessary to compute second derivatives together with first derivatives of the projected gravitational potential (see Eq. 2.11, 2.10). It is possible to compute efficiently the derivatives of this field on the sphere using the routine alm2mapspin_derv included in the S2HAT library. This routine can compute the first partial derivatives of an arbitrary spin-$s$ field on the sphere, $\partial_{\theta},\partial_{\phi}$, with respect to the spherical coordinates ($\theta,\phi$). The routine implements an efficient recurrence relation for the Legendre polynomials part of the harmonics and for their derivatives starting from their expression in terms of Wigner D-matrices [80]. This recurrence is similar to the one implemented in the HEALPix library, where however only the spin-2 field derivatives are available. Numerical tests conducted so far have shown that derivatives of the same field on a sphere, computed with S2HAT or HEALPix, are equal within the numerical error. Thanks to the capability of the routine to compute derivative for arbitrary spin-$s$ field it is possible to compute the second derivatives of the matter distribution at the same time when computing the deflection field (see Eq. (3.1)). As shown in Appendix B in fact, the second derivatives of the lensing potential can be naturally related to the first derivatives of the spin-1 deflection field. Thus, the deflection field and the second derivative of the lensing potential can be obtained with a single call to the spin-$1$ Legendre polynomial recurrence. This minimizes the extra computational overhead required to compute the second derivatives given that the Legendre polynomial recurrence is the heaviest part of each spherical harmonic transform algorithm (see, e.g., [77, 78]). Once partial derivatives are available it is then straightforward to extract the second covariant derivatives, which are the key quantity for a full-sky formalism adopted in this work.

4 Simulation results

4.1 Lensing observables angular power spectra

The first problem that we addressed using our raytracing algorithm was the evaluation of the full set of components of the CMB lensing jacobian. We thus set the position of our source at the CMB last scattering surface ($z^{*}\approx 1100$) and propagated the light ray position to redshift $z=0$. For these runs we assumed that no lensing took place between the maximum redshift covered by the DEMNUni simulation and the last scattering surface. One of the free parameters of the algorithm, together with the spatial resolution of the simulation, is the band-limit value $\ell_{\rm max}$ adopted to solve the Poisson equation and extract the lensing potential $\psi^{(k)}_{\ell m}$ from each convergence map $K^{(k)}$ (see Sec. 3.3). The choice of this band-limit defines effectively the Fourier modes $\Psi(k)$ of the matter distribution that are included in the simulation and thus can affect the precision of the results itself.

We decided to adopt a pragmatic approach to identify a suitable choice of band-limit for this application and proceeded as follows. We first computed the expected theoretical CMB lensing power spectrum for the Planck 2013 cosmology with CAMB and then evaluated the error introduced with respect to this reference quantity when the same calculation was performed including only the modes $k<k_{\rm ref}$ of the matter power spectrum. We report the results of this calculations for reference in Fig. 3. We set our goal to be the recovery of the effective CMB lensing potential with roughly 1% precision up to $\ell\approx 8000$, thus we had to include modes of the matter distribution up to $k_{\rm max}\approx 5h/Mpc$. Following this reasoning, we then translated this constraint in the value of the band-limit parameter to be used at a given i-th shell using the Limber approximation, i.e. $\ell^{(i)}_{\rm max}=\ell_{\rm max}^{(i),\rm theory}\equiv k_{\rm max}\chi_{i}$. However, it turned out that for several shells $\ell_{\rm max}^{(i)}$ became higher than the band-limit supported by the HEALPix grid where we sampled the $K^{(k)}$ field ($\ell_{\rm max}^{grid}=3\cdot$NSIDE). For these cases we set $\ell_{\rm max}^{(i)}=\ell_{\rm max}^{grid}$ while we chose to fix $\ell_{\rm max}^{(i)}$ to a conservative reference value of 8192 if $\ell_{\rm max}^{(i),\rm theory}<8192$. In Fig. 3 we show the modes of the matter distribution effectively included in the raytracing simulation as a function of redshift given the choice of $\ell_{\rm max}^{(i)}$ just presented. In Fig. 3 we also report the corresponding precision on the lensing potential obtained with the same choice of band-limit. As it can be seen from these figures, the optimized procedure improves the achievable precision at $\ell=8000$ from the 6% obtained assuming the same band-limit $\ell_{\rm max}=8192$ for all shells to the targeted 1% level. We, in fact, cannot resolve modes $k<5$ at all redshift because we are limited by the band-limit imposed by the resolution of our convergence maps and the optimized band-limit mitigates the impact of lack of power on those scales. This problem could be bypassed either fixing a band-limit to higher value (increasing the computational cost) or producing convergence maps at higher HEALPix resolution. The latter option would in fact allow to solve the Poisson equation at higher multipoles. As it will be shown in the next sections, however, we found no indication that this baseline setup suffers from a lack of power at small angular scales, and therefore we decided not to pursue this option due to its higher computational cost.

Once $\psi^{(k)}_{\ell m}$ had been obtained, we converted them into Fourier modes of the deflection field and resampled the deflection field and all its second derivatives required by the algorithm on a very high resolution ECP grid on which we later performed the raytracing step. We note that this part of the algorithm is the most computationally expensive step and therefore a more careful optimization of the choice of the band-limit parameter according to the users needs is crucial to reduce the computational cost of these kind of simulations.

In the following sections we will show results derived employing ECP grids characterized by the same number of isolatitude rings on the sphere $N_{\theta}$ and number of pixel per isolatitude ring $N_{\phi}$. The map resolution is therefore non uniform across the sky (with smaller pixels close to the poles and bigger close to the equator) but the area of the coarser pixels is well approximated by $\theta_{\rm res}=\sqrt{\frac{4\pi}{N_{\rm pix}}}$. We thus express the map resolution in terms of the parameter $\sqrt{N_{\rm pix}}=\sqrt{N_{\theta}\cdot N_{\phi}}$. In Fig. 4 we show the angular power spectrum of the different lensing observables extracted from the CMB lensing jacobian at $z=0$ – convergence, $C^{\kappa\kappa}_{\ell}$, rotation $C^{\omega\omega}_{\ell}$ and E- and B-modes of lensing shear, $C^{\epsilon\epsilon}_{\ell}$, $C^{\beta\beta}_{\ell}$. For these simulations we employed an ECP grid with $\sqrt{N_{\rm pix}}=262144$ which corresponds to $\theta_{\rm res}\approx 3^{\prime\prime}$ to derive these results. The previous set of parameters represents the highest setup in terms of resolution and band-limit shown in this paper, although we tested that setups with $\sqrt{N_{\rm pix}}=524288$ are achievable with the current numerical implementation of the code.

The simulated convergence and E-modes power spectra are in good agreement with the HALOFIT based results derived with CAMB. For the purpose of comparing the two quantities we restricted the integration of the matter evolution in the redshift range covered by the N-body simulation. The accuracy of the results is compatible with the uncertainties proper of the HALOFIT fitting formulae [81]. Matter perturbation on scales larger than the box-size of the simulation used to construct the lightcone should only be partially recovered by the stacking procedure employed in the map-making step [55]. However simulations having a box-size of 2Gpc like the ones of the DEMNUni suite seem to reproduce well the convergence power spectrum at large angular scales for all practical purposes.

In Fig. 4 we can also observe a non-zero image rotation $C^{\omega\omega}_{\ell}$ and shear B-modes $C^{\beta\beta}_{\ell}$ power spectrum which are due to multiple deflections of photons and are a pure imprint of beyond-Born corrections properly resolved by our simulations. We cross-checked this fact performing the raytracing with the same numerical setup but propagating the CMB lensing jacobian at first odder, i.e. in the Born approximation (Eq. (2.14)). The shear B-modes extracted from the first-order simulations are shown in Fig. 4 together with the beyond-Born result and they are consistent with numerical noise. We note that the second mixed covariant derivatives of a scalar field on the sphere commute, therefore we did not show the image rotation at first-order because it is identically zero by construction. All these results answer and confute the statement made in [26] according to which the multiple lens approach is not able to properly resolve beyond-Born corrections.

4.2 Consistency tests

We assessed the consistency of our results and the potential presence of numerical artifacts in two ways. First, we tested if the consistency relations between the components of the lensing jacobian were satisfied (see Appendix A). In addition, we tested if these also satisfied the consistency relation that can be derived between the E and B-modes of the effective displacement field of the rays. Following the notation in Eq. (3.2) we can define the effective displacement field ${\bf d}^{\rm eff}$ after $N$ deflections as

[TABLE]

where $\boldsymbol{\theta}$ is the photon un-displaced direction. Because we know the initial and final position of the simulated light rays, it is possible to invert the equation above using spherical triangles identities [73] and derive ${\bf d}^{\rm eff}$. We can extract the effective potentials of 2.4 from the E/B decomposition of this vector (i.e. spin-1) field as

[TABLE]

In Fig. 5 we show the maps of the effective potentials extracted from our simulations777From now on we will drop the eff superscript for the effective scalar potentials for sake of brevity..

The second test we performed was computing cross power spectra of lensing observables that are supposed to be null in absence of systematics (see Appendix A). In particular we looked at combinations of $\kappa\times\beta$, $\epsilon\times\beta$, $\omega\times\epsilon$, $\kappa\times\omega$ for the magnification matrix, and $\psi\times\Omega$ for the lensing potentials fields themselves. In Fig. 6 we show a selection of these cross-spectra omitting the ones that can be derived using the lensing consistency relations. All are clearly consistent with a null spectrum with no significant trend at all scales.

In order to test the consistency relations between lensing observables we adopted a cross-spectra technique and focused on cross-spectra combinations $\kappa\times\epsilon$ and $\omega\times\beta$. As it can be derived using the formulae in Appendix A, these cross-spectra can be constructed from the angular power spectra of the lensing fields themselves, as

[TABLE]

where the angular power spectra on the r.h.s. of the equations is computed directly from simulated maps. A similar construction can be obtained instead using the effective lensing potentials, as

[TABLE]

In Fig. 7, we compare cross-spectra as extracted from simulations with the two analytical constructions of Eqs. (4.3),(4.5). We can see that the consistency relation involving the gradient-like components is recovered with exquisite accuracy (below 0.1%) for both $\kappa$ and $\epsilon$ and $\epsilon$ and $\psi$. Conversely, the consistency relation between the shear B-modes and the rotation field does not reach the same level of accuracy, nevertheless it is satisfied at the 1% level at small angular scale and is better than numerical results presented so far in the literature (see e.g. [50]). The consistency relation derived using the $\Omega$ effective potential is satisfied at the 10% level and we can see a clear excess of power at small-scales ($\ell>1000$). As discussed in the following, the stability of the results of the ${\bf d}^{\rm eff,B}$ and shear B-modes with respect to the choice of numerical setup gives an indication that the cause of this discrepancy is likely purely numerical. A residual of aliasing effect might in fact be present in the shear and ${\bf d}^{\rm eff,B}$ affecting thus the consistency relation involving these two quantities. As explained in [74], the B-modes of spin fields are much more sensitive to numerical effects than the scalar quantities, and controlling residual resolution effects may become crucial.

4.3 Testing resolution effects

We additionally tested the stability of the simulated $\kappa,\gamma$ and $\omega$ maps with respect to different numerical parameters used for the raytracing. In Fig. 8 we show the effect of the map resolution used in the raytracing on the reconstructed angular power spectra of the lensing observables. The panels show the fractional difference between results of raytracing performed at progressively higher spatial resolution but keeping the same band-limits choice discussed in Sec. 4.1. For the kind of resolution considered for this test we recovered a percent level precision on $\kappa$ and $\omega$ already using a resolution of $\sqrt{N_{\rm pix}}=65538$ at scales $\ell\approx 5000$. This level of precision reaches the sub-percent level on $\ell<8000$ using our default numerical setup. However, a compromise solution between accuracy and computing sustainability could be identified for the case ECP $\sqrt{N_{\rm pix}}=131072$, where the reconstruction is well within percent level with respect to the highest resolution case up to $\ell\lesssim 8000$.

As it can be expected from the results of the previous sections, similar conclusions can be drawn for the shear E and B-modes and ${\bf d}^{\rm E,eff}$, showing that the code is able to simulate the expected signal with a high precision. We note however that the B-modes of the effective displacement did not reach a numerical stability with the same level of precision. For this field the numerical convergence between the cases adopting a resolution of $\sqrt{N_{\rm pix}}=131072$ and $\sqrt{N_{\rm pix}}=262144$ is of the order of 10%. As confirmed by the consistency relation analysis of the previous section, this field seems to be the most sensitive to resolution effects. Even though this level of precision is likely enough for all practical applications, we conservatively decided not to use this quantity for the scientific analysis of the following sections and derive the effective displacement from the $\psi-\kappa$ and $\Omega-\omega$ consistency relations.

4.4 Testing the N-Body simulation

Since beyond-Born effects are expected to show up as small angular scales where non-linear scales of the matter distribution contribute the most, it is important to quantify the impact of the resolution of the N-Body itself on the final result. In analogy to what we did in our previous work [55], we first characterized the impact of the shot-noise generated from the discrete matter distribution inside the box after projection on the 2D sphere. We used for this purpose Monte Carlo Gaussian realization of the surface mass density field derived from the shot-noise power spectra in each shell [51]

[TABLE]

where $\bar{n}_{k}$ is the particle density in the k-th shell. We then used these maps as an input of the raytracing algorithm and propagated the shot-noise induced uncertainties on all the simulated lensing fields. In Fig. 9 we show the shot-noise contribution to $\kappa$ and $\omega$. From this analysis it is evident that the excess of power at small angular scales that is observed in the recovered convergence power spectrum when compared with the CAMB results based on HALOFIT, is compatible with a noise bias induced by the simulation shot noise. We note however that even if the shot-noise affects at 10% level the convergence and shear E-modes, it is negligible for the rotation and shear B-modes. This can be explained as the physical origin of those signals is the coupling of subsequent deflections along the line of sight. A random realization of the shot noise has in fact a very small probability to replicate such lens-lens configuration, which is indeed a peculiar and interesting feature of the matter distribution from which the rotation field arises.

In addition, we also investigated whether the simulation resolution for high value of $k$ significantly affected our results. In particular we used as a benchmark case a Millennium-like simulation (500 Mpc boxsize, $2048^{3}$ dark matter particles) that adopts the same background cosmology and amplitude of matter perturbations of the baseline N-Body simulation described so far but having different initial conditions. In Fig. 10, we show the convergence and rotation power spectrum extracted from the Millennium-like and the baseline simulation. Even though it was not possible to run the Millennium-like simulation with the same initial conditions of our baseline setup, the agreement between the results recovered with these two simulations is remarkable. In particular we can see that dark matter halos of lower mass resolved in the Millennium-like simulation impact the convergence and rotation power spectrum below the 1% level for $\ell>1000$. The only noticeable difference appears at large angular scales where the Millennium-like based results display a lack of power. This is expected since the box size is four times smaller than the standard case, resulting in a more pronounced lack of power on scales larger than the boxsize, as noted in [55, 42]

4.5 Algorithm robustness

A standard back-of-the-envelope calculation predicts that a CMB photon undergoes a total of about 50 independent deflections from large scale lenses (i.e. of $300h^{-1}$Mpc) during its journey from the last scattering surface towards us [57]. However, when considering lensing on smaller scales, there could be more deflections along the line of sight and thus it is important to test the stability of the simulation results against the number of employed lensing planes. We note that for an accurate modeling of post-Born correction, a relatively low number of planes could be sufficient. Post-Born effects require, in fact, a fairly large separation between the background and forward planes to have any impact on the lensing signal. There are no post-Born effects for lenses right next to each other, i.e the $f_{K}(\chi_{k+1}-\chi_{k})$ factor of (2.10) approaches zero.

In the previous sections we presented the results of raytracing simulations employing 62 lensing planes as a reference setup. This number corresponds to the number of N-body snapshot at our disposal and thus limits the number of randomization the we have to perform in order to fill the full past lightcone. To assess the robustness of the results with respect to this specific choice in the number of lensing planes, we constructed the same lightcone realization employed throughout this work but then adopted a thinner matter shells slicing. In the benchmark case, the 62 outputs of the simulation have been translated into 62 matter shells of (average) comoving width of about 150 Mpc$/h$. Instead, for this test, we sliced the lightcone into 141 matter shells of (average) comoving width of about 50 Mpc$/h$. This scenario traces exactly the same matter distribution and evolution as the benchmark case, but it models the lensing signal with an increased number of deflections. In Fig. 11 we compared the results obtained with these two approaches. We found that the convergence signal is stable at better than 0.5% level, however, the rotation signal changes by a factor of 1.5 % on scales close to $\ell=8000$. The latter is expected, because if we increase the number of deflections, we expect the signal that depends on the lens-lens configuration to increase as well. However such difference is barely above the amplitude of the cosmic variance at those scales and the number of lens planes increased at the same time by more than a factor of 2. The number of lens planes employed in our baseline simulation setup seems therefore sufficient to resolve beyond-Born corrections, induced by the lens-lens coupling, with reasonable accuracy.

5 Comparison with theoretical predictions

After having confirmed the robustness of our result from a numerical point of view, in this section we compared the recovered signal with theoretical predictions based on the perturbative approach to the lens equation and beyond-Born correction [59, 31, 33, 27, 29, 28, 30].

5.1 Convergence corrections

As noted in [55] and confirmed analytically by [27] (hereafter PL16), the impact of beyond-Born corrections on the CMB lensing convergence signal is small and the first-order result dominates. A direct and accurate comparison between analytical predictions and numerical results is not straightforward because differences might be easily hidden in the cosmic variance scatter. The latter cannot be reduced given the impossibility of performing a full Monte Carlo analysis with N-body simulation of the size employed in this work. In Fig. 12 we compared the fractional difference of the lensing convergence power spectrum beyond-Born and the first-order result with the theoretical prediction computed with the numerical implementation of the perturbative corrections presented in PL16888http://cmbant.github.io/notebooks. Despite the fact that it was not possible to investigate quantitatively the agreement $\ell$ by $\ell$ for the reasons outlined above, the binned version of these curves shows a good agreement between simulations and analytical predictions, with a clear trend of an increase of power at small angular scales and a deficit of power at multipoles of few hundreds.

For sake of completeness we investigated the impact of beyond-Born corrections also in the real domain. In particular, we looked at the moments of the convergence maps obtained in the first-order approximation $\kappa^{\rm 1st}$ and in the beyond-Born regime $\kappa^{\rm ML}$. The quantity $\Delta\sigma^{2}_{\kappa}=\sigma^{2}_{\kappa^{\rm ML}}-\sigma^{2}_{\kappa^{\rm 1st}}$ can in fact be directly connected to the beyond-Born correction on the convergence power spectrum $\Delta C^{\kappa\kappa}_{\ell}$ using the well-known relation between the variance of a scalar field and its power spectrum

[TABLE]

We estimated the same quantity from the simulated $\kappa^{\rm 1st}$, $\kappa^{\rm ML}$ maps, $\Delta\sigma_{\kappa}^{2,\rm Sim}$, and compared it to the predictions of (5.1) derived using the analytical expectations for $\Delta C^{\kappa\kappa}_{\ell}$. To minimize the impact of numerical effects towards the end of the band-limit of the signal, we filtered the simulated maps setting to zero all the $\ell>8000$ modes and fixed the upper limit of the sum in Eq. (5.1) accordingly. We found that $\Delta\sigma_{\kappa}^{2,\rm Sim}=7.25\cdot 10^{-6}$ and $\Delta\sigma_{\kappa}^{2,\rm Theory}=7.00\cdot 10^{-6}$ respectively, thus the analytical predictions of beyond-Born correction match to a precision of roughly $3\%$ the results derived in the simulations for this specific realization of the matter distribution.

The recent work of PL16 noted that beyond-Born corrections generate a non-negligible amount of non-Gaussianities, in particular in the $\kappa\kappa\kappa$ bispectrum, and that those might be detected in the near future unlike corrections on $C_{\ell}^{\kappa\kappa}$. For this reason, we measured the variation of skewness and kurtosis of the 1-point PDF of $\kappa^{\rm 1st}$ and $\kappa^{\rm ML}$. The 1-point PDF of $\kappa^{\rm ML}$ and $\kappa^{\rm 1st}$ are shown in Fig. 13 together with the Gaussian approximation of $\kappa^{\rm 1st}$ defined as a Gaussian PDF having mean and variance equal to the one measured on the simulated map itself. Both $\kappa^{\rm ML}$ and $\kappa^{\rm 1st}$ PDF follow a Das and Ostriker model as it could have been expected from the analysis of Sec. 3.2.1. As discussed above, the variances of the two distributions are very similar. However their skewness and kurtosis change by a factor of roughly 30% when beyond-Born corrections are considered (see table 1 for a summary of these values). In particular, as shown in Fig. 13, the beyond-Born correction tends to modify mainly the high and low ends of the distribution. In Fig. 14 we show the measurements of the higher order moments of $\kappa^{\rm ML}$ and $\kappa^{\rm 1st}$ and $\kappa^{\rm shot-noise}$ obtained filtering the maps in order to retain only the modes $\ell\leq\ell_{\rm cut-off}$. We performed this test to assess whether large-scales artifacts induced by the N-body boxsize replication (especially in the high redshift shells) or shot noise effects dominated the values of the skewness and kurtosis of the map. As it is visible from the figure, the higher order moments are dominated by non-Gaussianities located at progressively smaller angular scale, consistent with non-linearities arising in the gravitational evolution of the matter distribution. In particular comparing the values of skewness and kurtosis of the first-order and ML maps we can clearly see the onset of the change in the amount of non-Gaussianities due to beyond-Born corrections starting at scales $\ell\approx 250$ and a negligible shot-noise contribution to these values for scales $\ell>1000$. Although a quantitative assessment of the precision of the analytical predictions for the amount of non-Gaussianities produced by beyond-Born corrections is beyond the scope of this paper, we qualitatively confirm the findings of PL16: beyond-Born corrections affect more the non-Gaussian part of $\kappa$ than its power spectrum. PL16 also noted that the specific signature of beyond-Born corrections to the $\kappa\kappa\kappa$ bispectrum depends on the specific triangle configuration and redshift. The effect consists mainly in a combination of a reduced lensing efficiency generated by lens coupling (which suppresses the bispectrum of the matter distribution) and additional distortions due to deflections across the direction where the gravitational potential changes. These two contributions have an opposite sign and their relative amplitude changes with redshift. The results of our simulations suggest that the leading effect of beyond-Born corrections consists in an overall suppression of the amount of non-Gaussianities in the CMB lensing convergence. This result seems to agree with the naive expectations from the results of PL16 (see, e.g., their Fig. 4 where the beyond-Born $\kappa\kappa\kappa$ bispectrum is mostly negative). However, being the skewness the sum of all the configurations of the $\kappa\kappa\kappa$ bispectrum, we cannot make a quantitative statement on the validity of the analytical predictions of this signal. Thus, we postpone a more quantitative analysis of this aspect on our simulated map to a future work.

5.2 Lensing rotation corrections

One of the most characteristic signs of beyond-Born corrections is the introduction of an image rotation field. In Fig. 12 we show the comparison between the analytical expectation for the rotation power spectrum and the result extracted from the simulation. We note that, despite the procedure devised to define the band-limits for the Poisson equation resolution (see Sec. 4.1) was gauged on the estimation the $\kappa$ field, $\omega$ is recovered with a similar level of accuracy. The results of the simulation agree to better than 5% at all scales and to better than 1% on scales around the rotation power spectrum peak. The largest discrepancy with analytical results can be observed at large angular scales ($\ell\lesssim 50$). Few factors may be the reason of this discrepancy on top of the uncertainties in modeling of the non-linear matter power spectrum. Analytical results have in fact been derived in flat-sky and using the Limber approximation, which becomes less accurate in these regimes. At the same time additional correlation may be introduced by boxsize replication during the map-making step of the algorithm (see Sec. 3.2) and the excess might be related to this procedure.

Similarly to what we did in the previous section, we investigated whether the variance measured in the map matched the one computed from the theoretical $C_{\ell}^{\omega\omega}$. We found that this number agrees within 10%, with a contribution of the shot noise of 1% or less. On the contrary, if we look at the 1-point PDF of the $\omega$ map in Fig. 15 we can see that the field is significantly non-Gaussian and the observed distribution is peaked and with long tails (see table 1 for the moments of the distribution). Unlike the convergence case, where higher order moments never exceed a value of 1, the level of non-Gaussianities is significant. This is an additional source of possible discrepancies with the analytical calculations that so far have not addressed the impact of non-Gaussianities in the $\omega$ field alone. The skewness of the $\omega$ map is much smaller than its kurtosis. For this reason we suppose that the major contribution to the rotation field non-Gaussianities may come from the trispectrum. Since a large and positive trispectrum is usually sign of highly non-linear and local effects [82], this signature might be related to the presence of massive and clustered halos in the maps that cause deviations from the purely weak lensing regime. This hypothesis should however be investigated quantitatively with dedicated tests and we leave this analysis to future work.

We finally note that the typical level of kurtosis expected for a shot-noise realization following the model used in Sec. 4.4 can account for 30-40% of the amplitude of the kurtosis (see table 1). However, a simple analysis of the scale dependency of the kurtosis analogous to the one of the previous section showed that the shot-noise contribution is negligible (few percent) up to scales $\ell\approx 2000$ where the kurtosis is already significant and equal to 2 (see Fig. 14).

Despite the prospect of measuring $\omega$ from the reconstruction of curl-deflections of CMB lensing are not promising [27, 83] - and thus analysis of its non-Gaussianities might not seem useful for upcoming experiments - we note that a curl-like deflection is very efficient in generating B-modes [84] and might leave a detectable imprint on top of the standard lensing B-modes signal generated by the gradient mode. In particular, in the context of CMB lensing, a significant level of non-Gaussianities in the curl mode may affect the estimators used to reconstruct the CMB lensing and the curl deflection potential from CMB maps [85, 86]. These in fact look for non-Gaussian signatures in the CMB maps left by lensing assuming the potentials are Gaussian. A first analysis of this aspect for the CMB lensing potential reconstruction has been reported in [36]. A more detailed analysis might be important to understand the ultimate precision limit with which the lensing-induced contribution to the CMB B-modes signal can be subtracted from the measured signal to enhance the detection of their primordial counterpart. This is in fact ultimately related to the presence of curl deflections and to their peculiar statistical properties. [59, 87]. We will discuss these issues in more details in the following section.

We finally note that despite $\kappa$ and $\omega$ are uncorrelated at first-order and thus $C_{\ell}^{\kappa\omega}=0$, these fields are not independent. They in fact have a non-null $\kappa\kappa\omega$ mixed bispectrum, as explained in PL16. We discuss details of this correlation in Sec. 6.

6 Lensed CMB beyond-Born approximation

6.1 Methodology

In the previous sections we described the properties of the simulated CMB lensing convergence and rotation fields including beyond-Born corrections. In this section we use these quantities to evaluate the impact of beyond-Born corrections on the CMB observables commonly used to estimate cosmological parameters, i.e. the angular power spectra of its temperature and polarization anisotropies.

For this purpose, it is convenient to adopt an effective approach where CMB photons undergo a single deflection as done in the standard Born approach but using a deflection field $\mathbf{d}^{\rm eff}$ synthesized from the $\kappa^{\rm ML}$ and $\omega$. The spin-1 E and B harmonic coefficients of this field are

[TABLE]

Because the consistency relations between $\kappa^{\rm ML}$ and $\psi$, as well as $\omega$ and $\Omega$, are satisfied with good accuracy (see Sec. 4.2), this is equivalent to using the ${\bf d}^{\rm eff}$ field but minimizing the impact of residual aliasing. Moreover this procedure is analogous to propagating the CMB sky through the shells as done in [55]. However, it allows to minimize the numerical effects induced by the pixel remapping operation that otherwise would have to be performed each time the ray trajectory is deflected by a lens plane in the multiple lens approach (see discussion in [55, 74]). A snapshot of the $\mathbf{d}^{\rm eff}$ field is shown in Fig. 16.

For sake of completeness, we note that when we lens CMB polarization we have to take into account the parallel transport of the basis of the tangent plane between the unlensed and lensed direction on the 2-sphere each time a light ray is deflected. This is a direct consequence of the fact that CMB polarization is a spin-2 field defined on the sphere. The parallel transport operation can be recast as a rotation of the components of the polarization field after each deflection by a proper angle $\delta$, which is used to rotate the polarization components back into the initial reference frame. This angle can be expressed in terms of the deflection field components and the incoming direction of the light ray [73, 75]. Every time we computed a ray deflection at each lensing plane in the raytracing, we evaluated at the same time the corresponding rotation angle $\delta^{(i)}$, and then updated subsequently a global effective rotation angle $\delta^{\rm raytracing}$ that describes the total rotation of the polarization components after the i-th deflection. To compute $\delta^{(i)}$, we adopted the formulae presented in the appendix of [73] and used the deflection field ${\bf d}^{(i)}$ defined in Sec. 3.3. This takes into account the incidence angle and the actual deflection at the plane crossing and it is equal to the effective displacement of the angular position of the ray seen by an observer. In the following we will compare $\delta^{\rm raytracing}$ with the angle derived solely from the parallel transport of the ray from its unlensed to its final direction. We compute the latter using the effective displacement described above, and denote this quantity as $\delta^{\rm transport}$ in the following. We will comment more on this point in the following sections.

Finally, we extracted the deflection field $\mathbf{d}^{\rm eff-1st}$ using the effective lensing potential recovered from the first-order convergence. This deflection field is used to derive the lensed CMB power spectra that serve as a benchmark for the results based on the Born approximation.

6.2 Impact on lensed CMB power spectra

Once we extracted the effective displacement field, we lensed the incoming CMB photons using the pixel based method implemented in the lenS2HAT code [74] with an oversampling factor of 8. For these simulations we filtered all the power in the displacement field for scales $\ell>8000$ and set the band-limit for the unlensed CMB to $\ell^{\rm CMB}_{\rm max}={8000}$. This setup is sufficient to obtain 0.1% accuracy on the lensed temperature, E and B-modes power spectra up to $\ell\approx 6000$ [74].

In order to measure the impact of the beyond-Born corrections, we produced a set of simulated CMB maps where we lensed 100 independent unlensed Gaussian CMB realizations using four different deflection fields: $\mathbf{d}^{\rm eff}$, $\mathbf{d}^{\rm eff-grad}$, $\mathbf{d}^{\rm eff-curl}$ and $\mathbf{d}^{\rm eff-1st}$. $\mathbf{d}^{\rm eff-grad}$ and $\mathbf{d}^{\rm eff-curl}$, in particular, were computed setting either B or E modes to zero respectively in Eq. (6.1). The CMB maps lensed with $\mathbf{d}^{\rm eff-1st}$ are instead used as reference quantities for results based on the Born approximation. We then compared each of the lensed CMB power spectra computed with the four different displacements in order to isolate the impact of specific beyond-Born corrections.

Following [27, 29, 30] we can classify the corrections to the lensed CMB power spectra into two types. The first involves corrections to the deflection field, which can be further divided into the contributions coming from the beyond-Born curl and gradient components of the deflection, or contributions from the higher order correlations of the two effective lensing and curl potentials. Since the consistency relations are valid, and following the terminology of PL16, we refer to these higher order effects as $\kappa\kappa\kappa$ and $\kappa\kappa\omega$ (or mixed) bispectra. In our simulations, lensed CMB maps obtained with $\mathbf{d}^{\rm eff}$ hold all beyond-Born effects coming from the combined action of both gradient and curl deflections as well as their higher order correlations $\kappa\kappa\kappa$ and $\kappa\kappa\omega$. $\mathbf{d}^{\rm eff-grad}$ and $\mathbf{d}^{\rm eff-curl}$ contain only beyond-Born corrections on the gradient deflections (including the effect of $\kappa\kappa\kappa$ bispectrum) and curl deflections respectively.

The second type of beyond-Born corrections affects only the polarization field and it consists in the rotation of the polarization tensor about its propagation direction by an angle $\beta^{\rm rotation}$. This effect appears when considering gravitational perturbations induced by higher-order scalar perturbations or vector and tensor perturbations at all perturbative orders [88]. The total beyond-Born lensed polarization field $\tilde{P}\equiv\tilde{Q}+i\tilde{U}$, where Q and U are the CMB Stokes parameters, is then related to the unlensed field $P$ as

[TABLE]

In the following, we denote with the superscript rotation the quantities including this additional rotation of the lensed polarization field. We measured $\beta^{\rm rotation}$ from our simulation as the additional rotation of the coordinate basis on the sphere observed when performing raytracing with respect to the rotation angle corresponding to the simple parallel transport of the ray. Following the notation of the previous section, this reads $\beta^{\rm rotation}=\delta^{\rm raytracing}-\delta^{\rm transport}$.

Because each CMB realization differs only by the lensing deflection used to lens the CMB, we can isolate the effect of each type of correction computing simply fractional differences. For the sake of brevity, we defined $\Delta^{x,y}_{z}(C_{\ell}^{XY})$ as the difference between the lensed CMB power spectrum $C_{\ell}^{XY}$ obtained using the $x$ and $y$ deflection fields described above, divided by the same CMB power spectrum obtained with the $z$ deflection field. In particular, to display all the beyond-Born corrections we focused on the following quantities:

$\Delta^{\rm eff,eff-1st}_{\rm eff-1st}(C_{\ell}^{XY})$ shows the amplitude of all beyond-Born corrections on the deflection field. These are shown as gray dots in Fig. 17. 2. 2.

$\Delta^{\rm eff-grad,eff-1st}_{\rm eff-1st}(C_{\ell}^{XY})$ isolates the impact of the beyond-Born corrections on the deflection field coming from the convergence and $\kappa\kappa\kappa$ bispectrum. Red diamonds mark this contribution in Fig. 17. 3. 3.

$\Delta^{\rm eff,eff-grad}_{\rm eff-1st}(C_{\ell}^{XY})$ isolates the impact of the beyond-Born curl deflection corrections as well as the mixed $\kappa\kappa\omega$ bispectrum. They are shown as blue squares in Fig. 17. 4. 4.

$\Delta^{\rm eff+rotation,eff}_{\rm eff}(C_{\ell}^{XY})$ shows the impact of the polarization rotation only. Fig. 17 shows this correction as pink dots. Here the superscript eff+rotation means that polarization rotation effects have been applied on maps lensed with the $\mathbf{d}^{\rm eff}$ field. 5. 5.

$\Delta^{\rm eff+rotation,eff-1st}_{\rm eff-1st}(C_{\ell}^{XY})$ gives the impact of all the beyond-Born corrections with respect to the first-order calculation. These are shown as black dots in Fig. 17.

We note that the $\Delta^{x,y}_{z}(C_{\ell}^{XY})$ quantities do not account for the effects of the convergence bispectrum due to LSS non-linearities. Because all the deflection fields are extracted from the same fully non-linear matter distribution of the DEMNUni N-Body simulation, they all contain the same amount of non-linearities. They differ only by how we model the bending of the light-ray. Thus, the majority of this contribution of LSS non-Gaussianity is canceled when computing any $\Delta^{x,y}_{z}(C_{\ell}^{XY})$ quantity.

In Fig. 17 we show the average of the fractional differences of the lensed CMB power spectra for each component of the beyond-Born corrections computed over 100 different realizations of unlensed CMB. The beyond-Born deflection corrections on the power spectra reach the $0.2\%$ level at $\ell\approx 4000$ for TT and TE power spectra. Measuring in detail those corrections is challenging due to their oscillatory shape and our impossibility to run MC simulations on a large set of N-body realizations. However, this is not the case for the corrections on the B-modes power spectrum, because the amplitude of the corrections is larger and reaches the $0.6\%$ level at $\ell\approx 4000$. The amplitude of the beyond-Born convergence correction is consistent in amplitude and shape with the one computed in [29, 30] for the angular scales considered in these papers.

The contribution induced by curl-like deflection is negligible for the TT, EE and TE power spectra but is important for the B-modes and accounts for $30\%$ of the signal at small angular scales. Thus, it cannot be neglected for an accurate analysis of the signal at $\ell\gtrsim 1000$. Moreover, we confirm the findings of [29] concerning the impact of the $\kappa\kappa\omega$ correction to the B-modes power spectrum. In particular the difference between $\Delta^{\rm eff,eff-grad}_{\rm eff-1st}(C_{\ell}^{BB})$ and the quantity $C_{\ell}^{BB,\rm eff-curl}/C_{\ell}^{BB,\rm eff-1st}$999We denote here with $C_{\ell}^{BB,\rm eff-curl}$ and $C_{\ell}^{BB,\rm eff-1st}$ the B-modes power spectrum obtained with $\mathbf{d}^{\rm eff-curl}$ and $\mathbf{d}^{\rm eff-1st}$ respectively., which accounts for the curl deflection without the contribution of $\kappa\kappa\omega$, is negligible for all practical purposes as discussed in the next section. Furthermore, the amplitude of the $\beta^{\rm rotation}$ correction on the polarization power spectra looks negligible up to $\ell\approx 6000$, similarly to the recent analytical findings of [89]. We postpone a more detailed discussion on this subject to Sec. 6.4.

In Fig. 18 and Fig. 19 we show a direct comparison of the results of our simulations with the analytical predictions of [30] (hereafter MFDD). We found a good agreement between those and our results for the beyond-Born deflection corrections including both the curl and gradient components on all TT, EE and, in particular, for BB power spectra. In Fig. 18 we also show a comparison to the predictions adopting the non-perturbative method of [29] for the TT and EE power spectra. This approach leads to a suppression of roughly a factor of 2 of the corrections at small angular scales. However, on the TT power spectrum the difference between the two analytical results is practically unobservable on all scales. On the other hand, on the EE power spectrum, we can clearly observe a better agreement with the non-perturbative theoretical predictions of [29] at $\ell\gtrsim 2200$. We note that both the theoretical predictions shown in Fig. 19 include all the matter perturbations at $z\leq 1100$, while simulations include contributions from structures located at $z\leq 99$. The good agreement of the theoretical and numerical results shows that the contribution of structures located at $z\geq 99$ to the beyond-Born signal is negligible.

Because beyond-Born corrections become more important at progressively smaller angular scales, via numerical simulations we were able to investigate them in more detail at $\ell\gg 3000$. At these scales, the signal receives contributions from matter perturbations on scales where the fitting formulae for the matter power spectrum and bispectrum used in the analytic calculations are less accurate. We found that the beyond-Born deflection corrections reach the $1\%$ level at $\ell\approx 6000$ for all TT, EE and BB power spectra.

We finally note that the beyond-Born corrections do not modify the TB and EB power spectrum at a sensible level. This is expected because lensing preserves parity, and we verified that this is true even when beyond-Born corrections are included.

6.3 Zooming on lensed B-modes

The B-modes are the most sensitive signal to beyond-Born corrections and, despite some contribution might become negligible at the power spectrum level, we tried to isolate them on the map level. The most interesting one is the contribution of the mixed bispectrum $\kappa\kappa\omega$ (or, equivalently, $\psi\psi\Omega$). As noted in PL16, two specific contractions of the mixed bispectrum can contribute to the large scales B-modes power spectrum. Following PL16 and adopting the flat-sky approximation plus a second order Taylor expansion around the unlensed direction of Eq. (6.2), we can write the lensed B-modes in Fourier space as the superposition of different contributions

[TABLE]

where $\varphi_{\boldsymbol{\ell}\boldsymbol{\ell}^{\boldsymbol{\prime}}}$ is the angle between the two Fourier modes $\boldsymbol{\ell}$ and $\boldsymbol{\ell}^{\boldsymbol{\prime}}$. The last equation corrects a sign for the first-order $\Omega(\boldsymbol{\ell})$ term reported in the published version of PL16 that was later corrected in the latest version of the paper submitted to the arXiv101010version 3 or later at https://arxiv.org/abs/1605.05662. We identify with $\tilde{B}(\boldsymbol{\ell})^{d^{\rm eff-grad}}$, $\tilde{B}(\boldsymbol{\ell})^{d^{\rm eff-curl}}$, $\tilde{B}(\boldsymbol{\ell})^{d^{\rm eff-1st}}$ the B-modes obtained setting to zero the $\Omega$ or $\psi$ respectively in the deflection field or using the first-order lensing potential $\psi^{\rm 1st}$, as explained in the previous section. With these quantities in mind and assuming $\psi$ and $\Omega^{\rm eff}$ are nearly Gaussian, we can isolate the two contributions of the mixed bispectrum $\kappa\kappa\omega$ to the B-modes power spectrum computing the following power spectra

[TABLE]

where we dropped the specific $\boldsymbol{\ell}$-dependencies in the higher order terms to simplify the notation. Following the convention adopted in PL16, we refer to $C_{\ell}^{BB,\rm corr}$ as the correlation term, and to $C_{\ell}^{BB,\rm mixed}$ as the mixed term. In addition to these terms, we isolated an additional contribution to the B-modes, $C_{\ell}^{BB,\psi\Omega\Omega}$, which is related to the $\psi\Omega\Omega$ bispectrum. All these terms are shown in Fig. 20 together with the B-modes signal induced by the curl deflection alone. $C_{\ell}^{BB,\rm corr}$ and $C_{\ell}^{BB,\rm mixed}$ can be compared with the results shown in Fig. 11 of [27]. The mixed bispectrum contribution exceeds the signal induced by the curl deflections alone, and displays an amplitude and a characteristic sign inversion at around $\ell\approx 1500$ in agreement with the findings of PL16. The $C_{\ell}^{BB,\rm corr}$ term shows a good agreement with PL16 for $150\leq\ell\leq 1500$. At $\ell\leq 150$ we get a much higher contribution with respect to the analytical results but we cannot disentangle if the discrepancy is due to numerical errors in the analytical calculations or boxsize effects in our lightcone construction. The sign of $C_{\ell}^{BB,\rm corr}$ at $1500\lesssim\ell\lesssim 2000$ disagrees with analytical prediction but the discrepancy is driven by the higher order terms contributing to our measurements. The contribution coming from $\psi\Omega\Omega$ bispectrum is in fact non negligible and if we properly subtract it, we recover a qualitative agreement with PL16. The situation is different for $C_{\ell}^{BB,\rm mixed}$. In fact, if we consider $\psi$ and $\Omega$ as Gaussian fields and recall that $C_{\ell}^{\psi\Omega}=0$, the higher order term $\langle\psi^{\rm eff}\Omega^{\rm eff}\psi^{\rm eff}\psi^{\rm eff}$ is approximately 0, providing thus a cleaner measurement of the mixed bispectrum contribution. The same reasoning applies to $C_{\ell}^{BB,\psi\Omega\Omega}$. Given their amplitude, we do not expect any of these terms to affect significantly the inflationary B-modes estimation. In fact, they partly cancel at large scales, besides their contribution being further suppressed if a delensing procedure reducing the contribution of the $\psi$-related signal is applied. The dominant contribution that can affect the estimation of $r$ after delensing comes from the presence of curl mode in the deflection. This generates a signal as high as the one produced by primordial tensor perturbations with an amplitude $r\lesssim 10^{-6}$ at $\ell\approx 80$, where the primordial B-modes signal is expected to peak. This level of contamination is far from the sensitivity of any future experiment proposed so far, but might affect any attempt to extract the tilt of the primordial tensor perturbations power spectrum $n_{T}$ because those measurements require an accurate lensing B-modes subtraction at smaller angular scales.

6.4 The impact of polarization rotation

Recently, two independent groups investigated the effect of the polarization rotation corrections on lensed CMB power spectra. MFDD adopted a next-to-leading order formalism in the light-cone gauge, and found that $\beta^{\rm rotation}=\omega$ up to second order in scalar perturbations. This result seems to confirm the guess of [88] that derived the same result if deflections are induced by vector or tensor perturbations at linear order, but speculated on its validity at higher order in scalar perturbations. Conversely, [89] (hereafter LHC) adopted a different approach based on the geometry of parallel transport, and found that $\beta^{\rm rotation}$ is too small to have a significant impact on lensed CMB power spectra.

In Sec. 6.2 we showed that the polarization rotation correction computed using the $\beta^{\rm rotation}$ extracted from our simulations is essentially an unobservable effect. In Fig. 20 we show a comparison of the power spectrum of $\beta^{\rm rotation}$ with the power spectrum of the lensing rotation field $\omega$ extracted from our simulations. These should be equal in the MFDD approach. On the other hand, LHC predicted the power spectrum $\beta^{\rm rotation}$ to have a negative power-law shape and an amplitude between $10^{-16}$ and $10^{-22}$ going from large to small angular scales (see Fig. 6 in [89] or Fig. 21). The $C_{\ell}^{\beta^{\rm rotation}\beta^{\rm rotation}}$ of our simulations differs from both these predictions when computed on the full sky. It is in fact dominated by numerical noise at large angular scales but displays an excess of power at $\ell\gtrsim 1000$ approaching the result of MFDD while keeping a systematically lower amplitude. We verified that this result was stable with respect to the choice of the raytracing resolution at 15% level. However, an accurate inspection of our $\beta^{\rm rotation}$ maps showed us that the technique used to compute $\beta^{\rm rotation}$ produces values of the polarization rotation angle roughly 1000 times higher than the typical r.m.s of the map for pixels located close to the poles. This effect is clearly unphysical, because the $\beta^{\rm rotation}$ field should be isotropic. Such systematic effect at the poles generated ringing in harmonic domain, and was the responsible of the highly oscillatory behaviour of the full sky $C_{\ell}^{\beta^{\rm rotation}\beta^{\rm rotation}}$. In order to mitigate these artifacts, we decided to mask the pixels having colatitude $0<\theta<\theta_{\rm cut}$ and $\theta_{\rm cut}<\theta<\pi-\theta_{\rm cut}$ when computing $C_{\ell}^{\beta^{\rm rotation}\beta^{\rm rotation}}$. As it can be seen in Fig. 20, once we masked a modest sky fraction around the poles (e.g. $\theta_{\rm cut}=0.5{\rm deg}$ ), $C_{\ell}^{\beta^{\rm rotation}\beta^{\rm rotation}}$ changed significantly. This value remained stable at $\ell\lesssim 4000$ if we masked larger areas around the poles. This new result disagrees clearly with the results of MFDD. On the other hand, in Fig. 21 we show a comparison between this new trend of $C_{\ell}^{\beta^{\rm rotation}\beta^{\rm rotation}}$ and the predictions of LHC. Despite a difference in the amplitude of $\beta^{\rm rotation}$ of roughly one order of magnitude, our measurement of $C_{\ell}^{\beta^{\rm rotation}\beta^{\rm rotation}}$ shows a qualitative agreement with the shape of LHC results. We also recover with good agreement the prediction of the cross-correlation power spectrum between $\beta^{\rm rotation}$ and $\omega$ (see Fig. 21). Despite this good qualitative agreement, the fact that the results are sensitive to the details of the processing of the maps suggests that a more precise measurement of $\beta^{\rm rotation}$ might require a redefinition of the technique used to isolate this effect. We plan to investigate this problem in more detail in future work.

Keeping in mind the caveat just outlined and for the sake of completeness, we also checked explicitly the validity of the calculations of MFDD for the lensed CMB power spectra. To this purpose, we rotated the lensed CMB simulations described in the previous sections by the lensing rotation angle $\omega$ extracted from the simulations, instead of $\beta^{\rm rotation}$. In Fig. 20 we show that this might indeed become the dominant effect on the power spectrum level, and could increase the amplitude of the correction to 6% on the E-modes and to 11% on the B-modes power spectrum at $\ell\approx 6000$. In Fig. 20 we also show a comparison of these results with the theoretical prediction derived in MFDD for the polarization rotation correction. In this scenario, the residual contribution to the large scale B-modes increased with respect to the results shown in Fig. 19, and could reach the level of $r\approx 10^{-5}$. This value is still below any foreseeable experimental target. We stress that this test does not validate the $\beta^{\rm rotation}=\omega$ result of MFDD (which is disfavoured in light of the results discussed above), but only their calculations of the lensed CMB power spectra at the correct perturbative level if $\beta^{\rm rotation}=\omega$.

Finally we note that the polarization rotation correction does not contribute to any isotropic parity breaking that would create non-zero TB and EB correlations, because both $\omega$ and $\beta^{\rm rotation}$ have effectively zero mean111111An analogous mechanism can arise in the context of CMB polarization experiments if the orientation of polarized detectors on the sky is misestimated [90]..

6.5 Detection perspectives

The future generation of CMB experiments, such as CMB-S4, will produce high resolution maps of CMB polarization over a large fraction of the sky (nearly 50%). Assuming the level of signal extracted from our simulations, these corrections might not be negligible [91]. We expect that for an experiment similar to [22], with $1.4\mu\textrm{K-arcmin}$ polarization sensitivity, 1 arcmin FWHM angular resolution, and observing a range of angular scales $30\leq\ell\leq 5000$, beyond-Born correction on TT, EE and BB power spectra could be measured with a cumulative signal to noise ratio (S/N) of 4.8, 0.5, 0.7 respectively, if we consider both the $\beta^{\rm rotation}$ and deflection corrections extracted from the simulations (as well as the same input cosmology used for the N-body simulation121212We assumed a Gaussian covariance for the B-modes power spectrum degraded by a scale-independent factor 1.2 to take into account the non-Gaussian component of the covariance introduced by weak lensing [92]).. Even in the most optimistic case where $\beta^{\rm rotation}=\omega$ (and thus the impact of the total beyond-Born corrections on the polarization power spectra is higher), the S/N for the corrections on EE and BB power spectra would be 1.4 and 2.2, respectively. Thus, in practice, a robust detection of post-Born corrections on the lensed CMB power spectra will be challenging as will highly depend on the capability to control the impact of extragalactic foreground to the required level of precision and accuracy. This is true especially for the TT power spectrum, where both the signal and the extragalactic contaminations are more important. We note that a realistic forecast of the detectability of these corrections is also complicated by the uncertainties in the modeling of the non-linear matter power spectrum and by degeneracies with other cosmological parameters affecting the CMB small-scale power (see e.g. [91] for a discussion on post-Born correction and $N_{eff}$ measurement).

A direct detection of the corrections on the $C_{\ell}^{\kappa\kappa}$, or a direct reconstruction of the lensing curl potential $\Omega$, will be extremely challenging for upcoming experiments [27, 83]. For the curl potential extracted from the simulations, in particular, the cumulative signal to noise for its direct reconstruction will only be 1.2, assuming a simplified Gaussian reconstruction noise as in [93]. Given the amplitude of the effect measured in the simulation, the detection of any $\beta^{\rm rotation}$ contribution on lensed CMB power spectra is not feasible. We discuss this issue in more detail in appendix C.

7 Conclusions

In this paper we presented an improvement of the raytracing algorithm through N-body simulations for CMB lensing presented in [55]. The modified method propagates the full CMB lensing jacobian, deflecting light rays trajectories using a similar multiple-lens approach. The improved numerical efficiency of the algorithm made possible a detailed analysis of the second-order lensing effects due to the relaxation of the Born approximation on the convergence, shear and lensing rotation fields at arcsecond resolution on the full sky. We also improved the setup of the N-body simulation of [55] employing a simulation with larger boxsize and resolution included in the DEMNUni suite [66, 45].

We tested the robustness of our method against resolution effects in the raytracing and in the N-body simulation used in this work, as well as against the number of employed lens planes finding no evidence of numerical artifacts.

We showed that the effect of beyond-Born corrections on the statistical properties of both the $\kappa$ and $\omega$ field 1-point PDF and power spectra, as well as on non-Gaussian statistics. We found a good agreement with theoretical predictions of [27] for the beyond-Born $C_{\ell}^{\kappa\kappa}$ and $C_{\ell}^{\omega\omega}$ power spectra and a qualitative confirmation that the beyond-Born corrections tend to suppress the amount of non-Gaussianity in the $\kappa$ field. However, we decided to postpone a quantitative characterization of the $\kappa\kappa\kappa$ and $\kappa\kappa\omega$ bispectrum to future works.

We then used the lensing observables extracted from our raytracing simulation to evaluate the effect of beyond-Born corrections on the lensed CMB power spectra up to $\ell\approx 6000$. We compared these findings with recent analytical predictions including beyond-Born deflection corrections of [30, 29], finding a very good agreement in the range of angular scales considered in that work ($\ell\lesssim 3500$) in particular for the B-modes power spectrum. We also found that the non-perturbative approach of [29], which predicts a suppression of the amplitude of beyond-Born corrections on lensed temperature and E-modes power spectra at $\ell\gtrsim 2000$, describes the results of the simulation better than the calculations of [30]. We showed, for the first time in the literature, the impact of the $\kappa\omega\omega$ bispectrum on the lensed B-modes power spectrum, and found it to be comparable to the contribution of the curl-mode in the deflection field at $\ell\gtrsim 1000$.

In addition, we measured the second-order gravitational rotation of CMB polarization, $\beta^{\rm rotation}$ directly from the simulations. Our results indicate that the impact of $\beta^{\rm rotation}$ on lensed CMB power spectra is negligible as found recently by [89]. We compared $C_{\ell}^{\beta^{\rm rotation}\beta^{\rm rotation}}$ and the cross-correlation $C_{\ell}^{\beta^{\rm rotation}\omega}$ measured in our simulations with theoretical predictions of [30] and [89]. Despite some differences in amplitude for both these observables still remain, we found a good qualitative agreement with the prediction of [89]. We recovered the shape and sign inversion of $C_{\ell}^{\beta^{\rm rotation}\omega}$ and found a good qualitative agreement with the shape and amplitude of $C_{\ell}^{\beta^{\rm rotation}\beta^{\rm rotation}}$ once we masked the region around the poles of the $\beta^{\rm rotation}$ map in the power spectrum computation. These pixels present in fact an anomalously high contribution which we interpreted as a consequence of a limitation of our approach to extract $\beta^{\rm rotation}$. Because this effect is small and a direct reconstruction of $\beta^{\rm rotation}$ field is essentially unfeasible, we decided to postpone the design of a more accurate simulation method of $\beta^{\rm rotation}$ to future work.

Finally, we discussed the possibility to detect the beyond-Born signatures with future high sensitivity and high resolution CMB polarization experiments, such as CMB-S4. We showed that the corrections to the temperature power spectrum could be detected with good significance, although the extragalactic foregrounds will act as the main (probably definitive) obstacle. The situation for the detection of beyond-Born corrections to E and B-modes power spectra is conversely more challenging as the expected signal to noise is less than 2.

Acknowledgments

We thank Carlo Baccigalupi, Anthony Challinor, Enea Di Dio, Giuseppe Fanizza, Julien Grain, Stefan Hilbert, Antony Lewis and Giovanni Marozzi for useful discussion and comments. We thank Enea Di Dio, Giuseppe Fanizza , Giovanni Marozzi and Antony Lewis for providing the analytical results used for comparison with the simulations of this work. We thank Julien Carron for providing the fitting code for the Das & Ostriker PDF model.

GF warmly thanks the hospitality of the University of Cambridge DAMTP and Institute of Astronomy where part of this work has been carried out. GF acknowledges the support of the Dennis Sciama Legacy fellowship and the CNES postdoctoral program. This work was partially supported by the RADIOFOREGROUNDS grant of the European Union Horizon 2020 research and innovation program (COMPET-05-2015, grant agreement number 687312). MC and CC acknowledge financial support to the “INAF Fellowships Programme 2010” and to the European Research Council through the Darklight Advanced Research Grant (# 291521). CC acknowledges the support from the grant MIUR PRIN 2015 “Cosmology and Fundamental Physics: illuminating the Dark Universe with Euclid”. MC has carried out part of the work for this project while supported by a grant (CUP B36G15002310006) of the European Union-European Social Fund, the Autonomous Region of the Aosta Valley and the Italian Ministry of Labour and Social Policy.

This research used resources of the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. Additional computing resources for this project have been provided by several class-C calls on the Fermi and Marconi machines of the Centro Interuniversitario del Nord-Est per il Calcolo Elettronico (CINECA, Bologna, Italy). The DEMNUni simulations were carried out on the Fermi machine at CINECA via the five million cpu-hrs budget provided by the Italian SuperComputing Resource Allocation (ISCRA) to the class-A proposal entitled “The Dark Energy and Massive-Neutrino Universe”.

We acknowledge the use of the publicly available CAMB code and HEALPix package.

Appendix A Weak lensing consistency relations

The full-sky magnification matrix can be written as

[TABLE]

where $\gamma_{ij}$ is a symmetric traceless tensor,

[TABLE]

$\omega\epsilon_{ij}$ is the anti-symmetric part of $A_{ij}$ and $(1-\kappa)\delta^{K}_{ij}$ is the non-zero trace part. The last equality follows from the definition of Eq. (2.4). Comparing the identities in Eq. (A.1) and applying the relations between derivatives of scalar and spin-$s$ spherical harmonics of [62], we get

[TABLE]

[TABLE]

[TABLE]

Being $\boldsymbol{\gamma}$ a spin-2 field, we can decompose it into E and B modes and therefore the following relations between the harmonic coefficients hold:

[TABLE]

We can then easily obtain the relations between convergence ($\kappa$), rotation ($\omega$), shear E ($\epsilon$) and B-modes ($\beta$) and lensing and curl-potential angular power spectra [62]:

[TABLE]

The cross-spectra $C^{\epsilon\beta}_{\ell}$, $C^{\omega\kappa}_{\ell}$, $C^{\kappa\beta}_{\ell}$, $C^{\omega\epsilon}_{\ell}$ and $C_{\ell}^{\psi^{*}\Omega}$ are all zero at first order in perturbation assuming that the Universe is statistically parity invariant [59]. Other cross spectra combination can be computed from the relations of harmonic coefficients provided above. Finally, note that at small scales the factor

[TABLE]

tends to one, thus the $\kappa$ and $\epsilon$ and $\beta$ and $\omega$ describe the same quantity [63, 59, 33].

Appendix B Partial derivatives of an arbitrary

spin-$s$ field

Starting from spin-$s$ spherical harmonics we can define their even and odd combination as harmonics

[TABLE]

If we consider the case of a spin-1 field and dropping the $(\theta,\phi)$ dependency for clarity we can write an explicit expression for ${}_{1}F^{\pm}$ applying the spin raising and lowering operators defined in [94]

[TABLE]

where $P_{\ell m}(\theta)$ are the scalar associated Legendre polynomials. In the following we dropped the $s$ index to identify scalar quantities. We recall that the harmonics ${}_{s}F_{\ell m}^{\pm}$ can be recast in terms of spin-s associated Legendre polynomials ${}_{s}P_{\ell m}(\theta)$

[TABLE]

where ${}_{s}P_{\ell m}^{\pm}(\theta)$ are odd and even combinations of the polynomials in analogy with Eq. (B.1).

We recall that for a spin-s complex field ${}_{s}\eta(\theta,\phi)=R(\theta,\phi)+iI(\theta,\phi)$ we have

[TABLE]

where H defines the arbitrary right-handed and left handed polarization convention. We assume H=1 in the following. In case of spin-1 field and assuming harmonic coefficients of Eq. 3.1 we get

[TABLE]

thus the derivatives of this spin-1 field give

[TABLE]

These are the expression of the partial derivatives of a spin-1 field computed by the S2HAT routines used in this work. If the quantities in Eq. B.10 and B.12 are renormalized by an additional factor $1/\sin\theta$ we obtain the normalized version of the derivatives needed to compute the magnification matrix.

Appendix C Detection of polarization rotation in future experiments

The direction dependent polarization rotation $\beta^{\rm rotation}$ mimics the effect of a direction-dependent cosmological birefringence process [95, 96], which generates coupling of CMB harmonic coefficients $E_{\ell m}$ and $B_{\ell^{\prime}m^{\prime}}$ with $\ell\neq\ell^{\prime},m\neq m^{\prime}$. Thus, future experiments aiming at constraining parity violating processes - or primordial magnetic fields that can generate anisotropic polarization rotation in the CMB - might have to deal with the beyond-Born signal as a potential contaminant. The level of $\beta^{\rm rotation}$ expected from the raytracing does not seem to be of any concern. However, we checked if this was the case if $\beta^{\rm rotation}=\omega$. In order to compare the beyond-Born rotation signal with the expected sensitivity of future experiments to anisotropic birefringence processes, we converted the observed level of variance of our $\omega$ map into the power spectrum amplitude $K_{\alpha\alpha}$ of an equivalent scale invariant anisotropic birefringence process $\alpha$, i.e. $\ell(\ell+1)C_{\ell}^{\alpha\alpha}/2\pi=K_{\alpha\alpha}={\rm constant}$. This class of models is commonly used in the literature to provide upper limits on anisotropic birefringence and can be easily related to the amplitude of a scale invariant primordial magnetic field that could source these kind of effect through Faraday rotation as [97, 98]

[TABLE]

In particular we found

[TABLE]

where we included in the calculation only the scales measured by our CMB-S4-like setup. The best constraint to date on $K_{\alpha\alpha}$ has been set by the POLARBEAR experiment [99] who reported a null detection of $K_{\alpha\alpha}=0.33\textrm{ deg}^{2}$ and a 95% confidence level upper limit of $K_{\alpha\alpha}\leq 1.0\textrm{ deg}^{2}$. CMB-S4 is expected to improve the error bars on $C_{\ell}^{\alpha\alpha}$ by almost three orders of magnitude [22] and constrain $B_{\textrm{1Mpc}}\leq 0.6\textrm{nG}$ at 95% CL. Thus, this level of signal might affect birefringence measurements in the upcoming years.

We also considered the possibility of accurately measuring $\beta^{\rm rotation}$ in this case through dedicated EB and TB anisotropic birefringence estimators [96, 100] and found that the cumulative S/N of those will be around 2. However, recently proposed complementary techniques to measure rotation through radio galaxies polarization might also help in this process [101].