Lensing Reconstruction in Post-Born Cosmic Microwave Background Weak Lensing
LLensing reconstruction in post-Born cosmic microwave background weak lensing
Dominic Beck, ∗ Giulio Fabbian, † and Josquin Errard AstroParticule et Cosmologie, Universit´e Paris Diderot,CNRS/IN2P3,CEA/Irfu, Obs de Paris, Sorbonne Paris Cit´e, France Institut d’Astrophysique Spatiale, CNRS (UMR 8617), Univ. Paris-Sud,Universit´e Paris-Saclay, bˆatiment 121, 91405 Orsay, France (Dated: September 5, 2018)The study of the Cosmic Microwave Background (CMB) lensing potential has established itselfby now as a robust way of probing the physics of large-scale structure growth. The most commonestimators of the lensing potential are derived under the assumption of Gaussianity of the matterdistribution and in the Born approximation of the photon diffusion. In this paper we study theperformance of quadratic estimators when applied to realistic sky maps extracted from multiple-lensray tracing techniques in cosmological N -body simulations. These are expected to model accuratelythe effects due to the non-Gaussianity of the matter distribution induced by its nonlinearity andthe deviation from the Born approximation. We show that both these effects on their own leadto reconstruction biases, but these tend to partially cancel each other when both these effects areconsidered together. We forecast the impact of these biases on the estimation of cosmologicalparameters for future high-sensitivity CMB experiments like CMB-S4. We find that the cold darkmatter density, Ω cdm , the optical depth to reionization τ , the amplitude of primordial inflationaryperturbations, A s and the sum of neutrino masses, M ν , could be biased at the 1-2 σ level, if noexternal data set is used. We also observe a reduction of the bias if external data like baryonacoustic oscillations (BAO) is included. I. INTRODUCTION
The cosmic microwave background (CMB) photonsdetected today have interacted with the matter distri-bution in the universe throughout their journey fromthe last scattering surface towards us. Such interactionsresult in the generation of the so-called secondaryanisotropies, i.e. fluctuations generated after the epochof matter-radiation decoupling (see e.g [1] for a review).These can be either due to scattering between CMBphotons and free electrons, such as inverse Compton orvelocity-induced scatterings (the thermal and kineticSunyaev-Zel’dovich effect and the Ostriker-Vishniaceffect) or to interactions of the photons with gravita-tional potential wells (e.g. the integrated Sachs-Wolfe[2] and Rees-Sciama [3] effects). Within this last classof secondary anisotropies the weak gravitational lensingof CMB anisotropies in temperature and polarizationis one of the key signals exploited by current andfuture experiments to obtain constraints on cosmologicalmodels.CMB lensing is sourced by the growth of all matterlocated between z = 0 and the last-scattering surface( z ≈ M ν ) and the properties of the darkenergy (see [4] for a review). ∗ Corresponding author. [email protected] † Corresponding author. [email protected]
The effect of lensing on the CMB manifests itself in ascale-dependent smoothing of the acoustic oscillations inthe angular power spectrum of temperature and E -modepolarization as well as in an increase of power in thedamping tails. The first evidence of this effect wasreported by the ACBAR experiment [5] and measuredwith high significance by SPT [6]. In addition, lensinginduces correlations in the harmonic coefficients of CMBanisotropies that can be used to reconstruct the distribu-tion of the line-of-sight integrated gravitational potentialthat lensed the CMB, i.e. the so-called lensing potential.The first attempts to measure the latter in CMB datafrom WMAP using cross-correlation techniques with ex-ternal LSS tracers were performed by [7, 8] and the firstsignificant direct detections were reported by the ACTand SPT Collaborations [9, 10]. Currently, the most pre-cise measurement of the CMB lensing potential has beenachieved by the Planck Collaboration, who measuredthis signal with a significance higher than 40 σ on nearlythe full sky [11]. The effect of gravitational lensing onthe CMB polarization anisotropies has been recentlyisolated by the current generation of ground-based CMBpolarization experiments POLARBEAR [12], SPTpol[13], and ACTpol [14] using CMB data alone and incross-correlation with LSS tracers [14–16]. Additionally,limits on the CMB B -mode power on subdegree scaleshave been obtained [17–19]. The B -mode signal of CMBpolarization on these scales is largely sourced by thelensing distortion of the primordial E -mode polarization.Hence, achieving high-sensitivity measurements of thelensed CMB polarization is a crucial step to increase theprecision of the CMB lensing potential reconstruction.With decreasing noise levels, higher angular resolutionsand larger areas observed by future experiments (e.g. a r X i v : . [ a s t r o - ph . C O ] S e p CMB-S4 [20]), the accuracy of reconstruction techniquesand theoretical modeling of the measurements has toimprove alike.To date, CMB lensing potential reconstruction anal-yses commonly rely on the assumption of Gaussianityof the unlensed CMB temperature field and of thelensing potential itself. The lensing potential, however,becomes non-Gaussian due to nonlinear structureformation mainly at late times. Although the level ofnon-Gaussianity is expected to be small due to the largenumber of potential wells that deflect CMB photons,the impact of this effect has to be quantified in light offuture high-precision measurements [21–23].Moreover, the Born approximation (i.e. the evaluationof the deflections of the photons with respect to theoriginal unperturbed line of sight), usually employedfor modeling CMB lensing, does not account accuratelyfor all features of the actual deflection process (e.g.the correlation between subsequent lensing events)neglecting therefore some of the sources of non-Gaussianstatistics in the lensing potential. Earlier attempts tomodel CMB lensing including the effects of nonlinearstructure formation were presented in [21, 24, 25].Recent works investigated the effect of the relaxation ofthe Born approximation on lensed CMB power spectraand CMB lensing power spectra, from both an analyticaland a numerical point of view [23, 26–28]. Similaranalytical studies were previously performed also inthe context of the weak lensing shear power spectrum[29–32]. While the most recent studies showed that themain post-Born effects are observed on the higher-orderstatistics of the CMB lensing potential rather than onits power spectrum, the impact of such effects on lensingreconstruction has not yet been evaluated. Recenttheoretical works further suggested that the presenceof non-Gaussianities in the CMB lensing potentialcould lead to percent level biases in the reconstructedCMB lensing potential power spectrum if they are leftunaccounted for [33]. This could in turn lead to a biasedestimation of cosmological parameters. In this paper we evaluate the impact of the non-Gaussian statistics of the CMB lensing potential on thecommonly employed quadratic estimator techniques forthe CMB lensing reconstruction. As these effects are of-ten too complex to model analytically, we use the simula-tions of [28] that include both the nonlinear evolution ofLSS and post-Born effects to model and investigate thisproblem numerically. The paper is organized in the fol-lowing way. In Sec. II we review the theoretical aspectsof weak lensing in the Born and post-Born regimes, andin Sec. III we review the properties of the statistical es-timators to extract this effect in the CMB. In Sec. IV wereview the details of the modeling implemented in thesimulations used in this work. In Sec. V we show theresults of our numerical experiments as their impact onthe lensing potential power spectrum and in Sec. VI wedescribe the impact of our findings on the estimation ofseveral cosmological parameters with a particular focuson the total mass of neutrinos, M ν , which is one of themain science targets of future CMB experiments. Finally,conclusions are made in Sec. VII. II. GRAVITATIONAL LENSING FORMALISM
In the weak lensing formalism the effect of deflectionsof light rays coming from a source plane is described bythe lens equation. This maps the final position ( β , χ ) ofa ray to the angular position of its source θ , i.e., β i ( θ , χ ) = θ i − c (cid:90) χ D A ( χ − χ (cid:48) ) D A ( χ ) D A ( χ (cid:48) ) Ψ , β i ( β ( θ , χ (cid:48) ) , χ (cid:48) ) d χ (cid:48) , (1)where χ is the conformal time, Ψ( β , χ ) is a gravitationalpotential located on the photon path, Ψ( β , χ ) , β i theirangular derivatives and D A ( χ ) is the comoving angulardiameter distance. The linearized mapping between animage at the source plane and the lensed image at agiven lens plane is described by the lensing magnificationmatrix (or lensing Jacobian). This can be computed asthe simple derivative of the equation above. A ij ( θ , χ ) ≡ ∂β i ( θ , χ ) ∂θ j = δ Kij − c (cid:90) χ D A ( χ − χ (cid:48) ) D A ( χ ) D A ( χ (cid:48) ) Ψ , β i β k ( β ( θ , χ (cid:48) ) , χ (cid:48) ) A kj ( θ , χ (cid:48) )d χ (cid:48) , (2) We warn the reader that the findings of [26] disagrees with laterstudies of [23, 28], probably due to numerical error in the evalu-ation of their analytical expressions. We refer the reader to thediscussion in [23] for more details. The derivatives in the small angle limit should be computed usinga coordinate system orthogonal to the current light ray’s direc-tion of travel. Numerical tests have shown that using angular derivatives causes a negligible error (see e.g. [34] and referencestherein). We note that the following formula can be extended to the full-sky case by promoting the partial derivatives to covariant deriva-tives.
In the weak lensing regime the magnification matrix isusually decomposed into A ij ≈ (cid:18) − κ − γ − γ + ω − γ − ω − κ + γ (cid:19) , (3)where κ is the lensing convergence, γ , are the compo-nents of the lensing shear, and ω is the lensing rotationangle [35]. The components of the magnification matrixare not independent and are connected through a seriesof consistency relations [36, 37].In the leading-order computations of the lensing effect,the photon path is approximated by the unperturbedphoton geodesic x ( χ ) ≈ θ χ , such that the line integral ofthe Newtonian potential Ψ simplifies to β i ( θ , χ ) = θ i − c (cid:90) χ D A ( χ − χ (cid:48) ) D A ( χ ) D A ( χ (cid:48) ) Ψ , β i ( θ , χ (cid:48) ) d χ (cid:48) . (4)At linear order in Ψ, the overall deflection of a photon α is then given by α ( θ ) = 2 D A ( χ ∗ ) (cid:90) χ ∗ dχ D A ( χ ∗ − χ ) D A ( χ ) ∇ Ψ( θ , χ ) , (5)where χ ∗ is the distance to the source plane. In the caseof CMB lensing, for instance, it is the distance to the lastscattering surface. The lens equation is usually rewrittenin terms of the lensing potential φ , which is connected tothe total photon deflection, as β ( θ , χ ∗ ) = θ − α ( θ ) ≡ θ − ∇ φ ( θ ) . (6)We note that the lensing potential and the lensing con-vergence can be connected in the weak lensing regimethrough the relations κ = − ∇ φ, (7) C κκL = [ L ( L + 1)]4 C φφL . (8)If we want to evaluate the lens equation at higher order,i.e. beyond the Born approximation (post-Born), wehave to account for the fact that photons do not travelalong the unperturbed background geodesics. Higher-order corrections are typically introduced perturbativelyin Eq. (1) by Taylor expanding the potential Ψ aroundthe unperturbed geodesic position.The distinct additional couplings that arise reflectthe change in the shape of a light ray bundle by one Despite being derived in the Born approximation, these rela-tions hold in the post-Born regime at sub-percent accuracy asdiscussed in [28]. lensing event affecting the amount of lensing generatedby a later lensing event (lens-lens coupling) as well asby changing gravitational potentials in the directionin which the ray path is bent. We refer the reader to[23, 29, 31, 38] for further details. Post-Born correctionsaffect the angular power spectrum of CMB lensingobservables in a minor way. In particular, the amplitudeof C κκL is suppressed on scales L (cid:46) L (cid:38) κ field, such asthe bispectrum, are, however, more affected and we willdiscuss these effects in the coming sections.A characteristic signature of post-Born corrections isthe appearance of curl-like modes in the overall lensingdeflection angle [23, 39], such that β ( θ ) = θ − ∇ φ ( θ ) − ∇ × Ω( θ ) . (9)Here we define ( ∇ × Ω) i ≡ (cid:15) ij ∂ j Ω,where (cid:15) ij is the Levi-Civita symbol in two dimensions and Ω is a pseudo-scalarfield. In analogy to the case of κ and ψ , Ω is related tothe lensing rotation ω as ω = − ∇ Ω , (10) C ωωL = [ L ( L + 1)] C ΩΩ L . (11) III. CMB LENSING RECONSTRUCTIONWITH QUADRATIC ESTIMATORSA. Formalism
Weak lensing by the large-scale structure of the Uni-verse remaps the primary CMB anisotropies according toEq. (6) such that its observed lensed Stokes parameter X along the θ direction is given by X ( θ ) = ˜ X ( θ − ∇ φ ) X ∈ ( T, Q, U ) (12)where ˜ X is the primordial unlensed CMB and φ the lensing potential. In the harmonic domain theremapping operation acts as a convolution that mixespower in different multipoles and therefore correlatesmodes across a band determined by the power in thelensing potential [40]. The lensing potential itself can beextracted statistically using the observed CMB, assum-ing that the underlying, unlensed CMB is on averagehomogeneous and isotropic. This operation commonlyinvolves the use of the so-called quadratic estimator[41, 42], which relies on the lensing information in thetwo-point correlation of the CMB fields. Althoughhigher-order correlations will become more important inreconstructing the CMB lensing potential to exploit thefull power of future data sets to a more optimal precision[43–45], the quadratic estimator is proven to be a veryrobust tool thanks to its well understood biases andcapability of quick forward modeling of instrumentaland systematic effects. Furthermore we note that, todate, none of the proposed alternative estimators candispense easily with the assumption of Gaussianity ofthe CMB lensing potential. In the following, we willuse an implementation of the quadratic estimator, thepublicly available code lensquest . In the quadraticestimator context we assume the primordial unlensedCMB to be a Gaussian field such that the harmoniccoefficients for its temperature, E -mode and B -modeanisotropies, a X(cid:96)m , X ∈ T, E, B , have a variance givenby the four nonzero power spectra C TT (cid:96) , C TE (cid:96) , C EE (cid:96) and C BB (cid:96) . Likewise, for the harmonic coefficients afterlensing, ˜ a X(cid:96)m , X ∈ T, E, B , we can write the variance ofthe harmonic coefficients, computed taking the ensembleaverage over primordial CMB and matter realizations,as (cid:68) ˜ a X † (cid:96)m ˜ a Y(cid:96) (cid:48) m (cid:48) (cid:69) = δ (cid:96)(cid:96) (cid:48) δ mm (cid:48) ˜ C XY(cid:96) .For a given realization of the lensing potential, thisvariance acquires non-diagonal terms due to the charac-teristic introduction of correlations in the harmonic coef-ficients due to gravitational lensing. This can be used toconstruct an estimator for the lensing potential [41, 42].On the full sky, this estimator takes the formˆ φ XYLM = (cid:88) (cid:96) m (cid:96) m K XYLM(cid:96) m (cid:96) m ˆ a X(cid:96) m ˆ a Y(cid:96) m , (13)where the convolution kernel is given by K XYLM(cid:96) m (cid:96) m = A XYL L ( L + 1) ( − M (cid:18) L (cid:96) (cid:96) − M m m (cid:19) g XYL(cid:96) (cid:96) . (14)This kernel has cosmology and experiment-dependentweights, which read g XYL(cid:96) (cid:96) = f XY ∗ L(cid:96) (cid:96) C XXn(cid:96) ˜ C Y Y n(cid:96) or g XYL(cid:96) (cid:96) = f XY ∗ L(cid:96) (cid:96) ˜ C XXn(cid:96) ˜ C Y Y n(cid:96) , (15)if X = Y or X (cid:54) = Y , respectively. These weights arechosen such that the variance of ˆ φ XYLM is minimal and we adopt the measured power spectra includingthe instrumental noise power spectrum N XY(cid:96) , i.e.˜ C XY n(cid:96) = ˜ C XY(cid:96) + N XY(cid:96) .The response functions f XYL(cid:96) (cid:96) used in this workare those of [42], with the distinction of using thelensed CMB power spectra ˜ C (cid:96) to mitigate the biases http://github.com/doicbek/lensquest Correlation between T and E is neglected in the estimatorweights, causing the estimator for XY = T E to be slightly sub-optimal in favor of computational time. of O (cid:16) C φφL (cid:17) that arise in the lensing potential powerspectrum calculation using this estimator [46]. We notethat this choice of weights might still be suboptimal andlead to biased results from very small-scale CMB temper-ature signal [47]. This bias could be mitigated replacingthe temperature autopower spectrum with the lensedtemperature-gradient power spectrum C ˜ T ∇ ˜ T(cid:96) , appearingin the nonperturbative response function calculation [48].Because in the following we will compare lensed CMBrealizations among each other and do not compare toa specific model, this has a negligible effect on our results.The normalization vector A XYL in Eq. (13) is given by A XYL = L ( L + 1)(2 L + 1) (cid:32)(cid:88) (cid:96) (cid:96) g XY ∗ L(cid:96) (cid:96) f XYLl l (cid:33) − (16)and ensures that the quadratic estimator is unbiased.The three CMB anisotropy fields allow for six separateestimators of φ . The estimator for XY = BB hasa vanishing signal-to-noise in cosmological scenarioswhere gravitational wave perturbations are negligiblecompared to scalar perturbations, effectively reducingthe number of estimators to 5. We will thus ignore it inthe following without introducing an appreciable loss inthe overall sensitivity.The power spectrum of the quadratic estimate of thelensing potential is then a contraction of the CMB four-point function, which includes three terms up to firstorder in φ (cid:42) L + 1 (cid:88) M (cid:16) ˆ φ ABLM (cid:17) † ˆ φ CDLM (cid:43) ≈≈ N (0) ,ABCDL + C φφL + N (1) ,ABCDL . (17)The biases N (0) , ABCDL and N (1) , ABCDL arise fromdisconnected Gaussian two-point contractions of theCMB fields and — in the case of the latter — of thelensing potential up to first order in C φφL [49]. Ananalytic expression for the zero-order bias, N (0) ABCDL ,can be found in [42]. In practice, the computation ofthe realization-dependent zero-order bias [50, 51] withthe help of Monte Carlo simulations is preferred to theevaluation of the analytic formula, since it accountsfor small mismatches in the two-point statistic betweensimulation and data. Analytic expressions for N (1) L canbe found in [49, 52] and an analog method to computeit using Monte Carlo simulations in [13].The different estimators for φ can be combined into anoptimal minimum-variance estimator asˆ φ mv LM ≡ (cid:88) XY w AB ˆ φ ABLM , (18)with weights w ABL = N mvL (cid:88) CD (cid:0) N − L (cid:1) ABCD (19)and minimum variance lensing noise N mv L = (cid:32) (cid:88) ABCD (cid:0) N − L (cid:1) ABCD (cid:33) − . (20) B. Effect of non-Gaussianities on quadraticestimators
The formalism derived in the previous section assumesthat all the non-Gaussianity in the CMB is entirely dueto the lensing effect and that the lensing potential is aGaussian field. However, this is just an approximationand if the lensing potential (or equivalently the lensingconvergence) has nonzero higher-order correlations,there are additional terms involving four-point functionsof lensed CMB fields that create distinct biases. Thisproblem was first studied in [33] in the context ofassessing the impact of the nonlinear evolution of thematter distribution in the lensing reconstruction. Inthis work the authors derived expressions for the biasinduced by a nonzero bispectrum in the lensing potentialcaused by the nonlinear gravitational collapse that is oforder O (cid:18)(cid:16) C φφL (cid:17) / (cid:19) and is referred to as N (3 / L . TheTT reconstruction channel was found to be the mostsensitive on angular scales (cid:96) (cid:46) .
5% for low noiseand large sky coverage experiments. This level of bias issignificant in light of the expected future experimentalsensitivity. Understanding the amplitude and natureof higher-order biases and their effect on our ability ofconstraining the cosmology is therefore crucial.Modeling these effects analytically becomes cumber-some very quickly. Therefore, we decide to adopt a nu-merical approach and assess the impact of these biasesthrough accurate and realistic numerical simulations. Inorder to tackle the problem in its full complexity we de-cide to use simulations that include not only the non-linear evolution of matter studied in [33] but also non-Gaussianity induced by post-Born effects. Analyticalpredictions of the shape and amplitude of these non-Gaussian correlations have recently been computed in[23, 27, 53].
IV. MODELING CMB LENSING AT HIGHERORDER
To test the bias in the lensing reconstruction in-duced by non-Gaussian evolution of the matter and post-Born effect we need to simulate the lensing of CMBanisotropies including both these effects. For this pur-pose we use the simulation method and results of [28](hereafter FCC18). This work produced a collection oflensing observables κ, ω, φ,
Ω, derived from a ΛCDM N -body simulation of the DEMNUni suite [54, 55] in theBorn approximation and using multiple-lens ray trac-ing techniques. The N -body simulation employed inFCC18 used 2048 dark matter particles and a box sizeof 2 Gpc/ h from z = 99 to z = 0. This redshiftrange cover allows to reproduce the CMB lensing ker-nel D A ( χ ∗ − χ ) /D A ( χ ) D A ( χ ∗ ) with sub-percent preci-sion. The mass resolution of the simulation at z = 0 is M CDM = 8 . × M (cid:12) /h and the gravitational soft-ening length is set to (cid:15) s = 20 kpc /h corresponding to0.04 times the mean linear inter-particle separation. Be-low, we briefly summarize the specificities of the light-cone construction and ray tracing algorithm adopted inthese simulations as well as further tests complementaryto the one presented in FCC18 and specifically performedfor this work. We refer the reader to FCC18 for a moredetailed discussion. A. Ray tracing algorithm for CMB lensing
Starting from a series of snapshots in time of an N -body simulation, the algorithm adopted in FCC18 re-constructs the full-sky past light cone of the observerfrom redshift z = 0 to the maximum redshift covered bythe simulation z max (in our case z max = 99). Becausethe universe volume simulated in the N -body is finite,we replicate the box volume in space to fill the wholeobservable volume between 0 ≤ z ≤ z max . To avoid re-peating the same structures along the line of sight andto recover (at least partially) structures on scales com-parable to the box size, the algorithm employs a specificrandomization procedure for the particle positions as in[24, 56]. The light cone is then sliced into spherical shellsof thickness ∆ χ = 150 Mpc /h . The particles inside eachof these volumes are then projected onto spherical planesof surface mass density, as described in [57]. The algo-rithm then converts the surface mass density planes intoconvergence fields. With this discrete version of the lightcone at hand, it is convenient to discretize the geodesicand lens equation of Eqs. (1) and (2) [58–60] β ( θ , χ ) = θ − N − (cid:88) k =0 D A ( χ − χ k ) D A ( χ ) α ( k ) ( β ( k ) ) . (21)Here k is the shell index and we define the gradient of thetwo-dimensional (2D) projected gravitational potentialas α ( k ) ( β ( k ) ) = 2 c (cid:90) χ k +∆ χχ k − ∆ χ dχ ∇ Ψ( β ( k ) , χ ) D A ( χ k ) . (22) α ( k ) is easily computed starting from the convergencefield of each shell ( k ) using a spin-1 spherical harmonictransform [34, 61] in the E and B decomposition α ( k ) ,E(cid:96)m = 2 κ ( k ) (cid:96)m (cid:112) (cid:96) ( (cid:96) + 1) α ( k ) ,B(cid:96)m = 0 . (23)The latest operation requires the computation of thespherical harmonic coefficients κ ( k ) (cid:96)m using a fast sphericalharmonic transform up to a given cut-off in power (cid:96) max .The choice of (cid:96) max for each different shell is optimized toensure the total deflection is computed with sub-percentprecision for scales (cid:96) (cid:46) A Nij ( θ , χ N ) = δ Kij − N − (cid:88) k =0 D k,N D N U ( k ) ip ( β ( k ) , χ k ) A ( k ) pj ( θ , χ k ) , (24)where N is the number of planes necessary to reachthe source at comoving distance χ N and U ij is the ma-trix of the second derivatives of the gravitational poten-tial, ∂ Ψ /∂β i ∂β j . U ij can be computed easily as deriva-tives of the component of the spin-1 field α ( k ) (see Ap-pendix A of FCC18). In Eq. (24) we use the notation D k,N ≡ D A ( χ N − χ k ) and D k ≡ D A ( χ k ) for simplicity.Implementing Eq. (24) in numerical simulations becomesquickly prohibitive for a large number of lens planes andlarge sky fraction. FCC18 adopted the multiple lens ap-proach of [62], who showed that the equation can berewritten in a more efficient form that requires one tostore in a memory for a given k th iteration just the posi-tion of the light rays at the two previous positions β ( k − and β ( k − , β ( k ) = (cid:18) − D k − D k D k − ,k D k − ,k − (cid:19) β ( k − (25)+ D k − D k D k − ,k D k − ,k − β ( k − − D k − ,k D k α ( k − ( β ( k − ) . By differentiating with respect to θ as in Eq. (2), we ob-tain the recurrence relation for the magnification matrix A ( k ) ij = (cid:18) − D k − D k D k − ,k D k − ,k − (cid:19) A ( k − ij (26)+ D k − D k D k − ,k D k − ,k − A ( k − ij − D k − ,k D k U ( k − ip A ( k − pj . This algorithm was originally developed in the contextof galaxy lensing, but adapted to spherical geometry in[61] and developed first in [28, 56] for CMB lensing. Thisapproach is also convenient to derive the magnificationmatrix and lensing observables in the Born approxima-tion that we will use later to isolate the contribution com-ing from post-Born effects. Assuming the background distortion, the first-order magnification matrix is A ( N ) , ij ( θ , χ s ) = δ Kij − N − (cid:88) k =0 D k,N D N U ( k ) ij ( θ , χ k ) . (27)We note that the U ij matrix is symmetric because mixedderivatives commute and thus the rotation, ω , is identi-cally zero. B. Impact of the LSS bispectrum
FCC18 carried out an accurate characterization of thepost-Born corrections on κ , ω and lensed CMB powerspectra and compared extensively with their analyticalpredictions derived in [23, 27, 38, 53, 63]. However,the analysis did not investigate in detail the impact ofthe nonlinear evolution of large-scale structures and howsimulation properties match with analytical predictionsof the higher-order statistics of the κ field. Below wepresent additional validation tests performed to assessthe reliability of these simulations in modeling higher-order statistics of post-Born corrections and nonlinearLSS evolution. We limit our analysis to the statistics ofthe κ field and its cross-correlation with ω . Higher-orderstatistics of the curl mode of the deflection field beyondthe mixed κκω bispectrum [23], which appear at higherorder in the perturbative expansion, are lacking theoret-ical predictions. The measurement of the κκω and κωω bispectrum in the simulations used in this work throughits effect on lensed B -modes power spectrum was pre-sented in FCC18, together with the measurement of thepost-Born induced curl mode on lensed CMB power spec-tra. We refer the reader to that work for a more in-depthdiscussion and comparison with theoretical predictions.
1. Higher-order statistics of the CMB convergence
To verify the accuracy of the simulations in reproduc-ing the expected level of non-Gaussianity in κ , we com-pare its skewness as measured in the simulations withthe values obtained by contracting the predicted theoreti-cal bispectrum including LSS nonlinearity and post-Borncorrections. The definition of skewness given a pixelizedmap of a scalar field, X , is S [ X ] = (cid:104) XXX (cid:105) = 1 N pix N pix (cid:88) p X p , (28)where p is the pixel index and N pix the total number ofpixels in the map. Following [64, 65], we compute theskewness in terms of the reduced bispectrum b L L L as S [ b L L L ] = L max (cid:88) L L L (2 L + 1)(2 L + 1)(2 L + 1)(4 π ) (cid:18) L L L (cid:19) b L L L , (29)with corresponding variance dominated by the disconnected six-point function σ S (cid:39) π L max (cid:88) L L L (2 L + 1)(2 L + 1)(2 L + 1)(4 π ) (cid:18) L L L (cid:19) C L C L C L . (30)In particular, the skewness of the Born-approximatedconvergence, κ F , obtained from the first-order mag-nification matrix, provides a measurement of theLSS-induced bispectrum. The bispectrum of theconvergence computed using the multiple lens raytracing algorithm, κ R , receives a contribution from theLSS-induced bispectrum as well as from the post-Borncorrections induced bispectrum. The difference ofthe skewness of κ R and κ F gives thus a direct mea-surement of the collapsed post-Born-induced bispectrum.We use the formulas presented in [66] and [23] tocompute the bispectrum of κ due to LSS nonlinearity(at tree level in density perturbations or adopting thenonlinear fitting formula from [67]) and post-Borneffects, respectively. In Fig. 1 we show a comparisonbetween the skewness measured in the low-pass-filteredsimulations and their expected theoretical value as afunction of the maximum multipole cutoff used in thecalculations. We find a good agreement between sim-ulation and theoretical expectations for the post-Bornbispectrum part, confirming the findings of FCC18 onthe level of the lensed CMB B -mode power spectrum.For this observable, the post-Born κκκ bispectrum isthe dominant correction while the contribution of thecurl mode in terms of the κκω bispectrum is negligiblysmall (see also [27]). The LSS skewness agrees well withtheoretical expectation on scales 75 (cid:46) L max (cid:46) N -body simulation used forthis work and found an excess of power at low valuesof k (cid:46) . − h for both squeezed and equilateralconfigurations. These scales contribute significantly tothe signal on angular scales (cid:96) (cid:46)
100 (see e.g. [4]) andcould be responsible of the excess of skewness observedwhen only such scales are included. Although in FCC18the replication procedure was shown to produce accurateresults on the large scales of C κκL and no significantspurious excess of power was observed, we tested the FIG. 1. Comparison of the skewness for different cutoff valuesof the convergence multipoles. The theory curves are com-puted using the tree-level expression of the LSS convergencebispectrum including the Scoccimarro & Couchman fit of [67],as well as post-Born corrections of the b κκκ and b κκω (+) bis-pectra of [23]. Only the absolute values are shown; negativevalues are marked by a dashed line or triangular marker. stability of our results on lensing reconstruction withrespect to the minimum multipole employed in theanalysis. We found negligible differences when excludingCMB angular scales (cid:96) ≤ L max (cid:38) L ≈ k (cid:38) − h [4, 28]and on these angular scales uncertainties on the matterpower spectrum are already of the order of 15% [69]. Theuse of nonlinear fitting formulas improves the agreementwith simulation results with respect to the tree-level bis-pectrum. We note that we do not investigate possible FIG. 2. Impact of CMB convergence bispectrum on lensed temperature (left) and E -mode (right) power spectra. The toppanel shows the total correction accounting for the LSS and post-Born induced bispectrum, while the bottom panel shows thecorrection due to only nonlinear LSS evolution. The theoretical predictions of [27] are shown in black and simulation results inred. The green curves show the values of the nonperturbative corrections computed in [27] for the temperature and E -modepower spectra. Binned theoretical predictions are shown with empty markers. The error bars include only the uncertainty onthe average over the Gaussian MC realizations and do not include the sample variance of the convergence bispectrum. improvement using alternative nonlinear bispectrum fit-ting formulas, as, for example, the one introduced in [70].The validity of these equations at high redshifts was notvalidated and the differences with respect to the Scocci- FIG. 3. Same as Fig. 2 for the BB power spectrum. marro & Couchman formulas formulas [67] were shownto be marginal and relevant only for a subset of the bis-pectrum configurations (see discussion in [23]).
2. LSS bispectrum effect on lensed CMB
Non-Gaussianity in the lensing potential can affect theshape of the lensed CMB power spectra. The authors of[27] (hereafter LP16) computed the effect on the CMBpower spectrum induced by the bispectrum of the CMBconvergence due to the nonlinear evolution of matter(hereafter LSS term) and the one due to post-Borncorrections (hereafter PB term). FCC18 showed thatthe corrections computed by LP16 for the PB termmatch very well the results extracted from ray tracingsimulations. As a validation test for this work, wefocused on measuring the corrections to lensed CMBpower spectra generated by the LSS term alone, as wellas those due to the combination of LSS and PB terms.We then compared the results of the simulations withthe theoretical prediction of LP16. To isolate the LSSterm, we lens 100 Gaussian realizations of unlensedCMB maps with a deflection field extracted from the κ F map as performed in FCC18. From the averageof the power spectra of these maps we subtract theaverage power spectrum of the 100 CMB realizationsthat were lensed with a deflection field computed froma Gaussian realization of the lensing convergence κ G with power spectrum C κ F κ F L . Similarly, to measure thetotal correction, we repeat the same procedure with κ R and C κ R κ R L to produce the Gaussian realizations of thedeflection field. In Figs. 2 and 3 we show the resultsof this analysis together with a comparison with theprediction of LP16.The theoretical predictions for both the total and LSSbispectrum (which is the dominant term) agree quite wellwith the simulation results on the relevant angular scales,especially the ones implementing the nonperturbativeformalism for the TT and EE power spectrum as dis-cussed in PL16 and FCC18. This approach accounts forthe fact that even in the Gaussian approximation lensingis a O (1) effect at small scales and therefore treatingthe corrections due to non-Gaussianity as perturba-tions around an unlensed field leads to inaccurate results.Despite the overall good agreement, however, some dif-ferences can be observed. This is expected because, un-like the analytical approximations, simulations includethe effects of non-Gaussianity nonperturbatively and theexact shape of the correction depends on the detailedshape of the bispectrum. In particular, simulation resultsshow an excess of power on the B -mode power spectrumcompared to analytical predictions. This is consistentsince B -modes are more sensitive to small scale lenses andthus non-Gaussianities due to strongly nonlinear densityfields are expected to give larger corrections where theperturbative expansion becomes less accurate. The dis-crepancies at scales (cid:96) (cid:46)
100 could conversely arise dueto the excess of skewness discussed in the previous sec-tion, although we stress that a larger skewness does notseem to affect significantly the temperature and E -modepower spectrum, where the corrections are dominated bystructures at (cid:96) (cid:46) V. RESULTSA. Numerical setup
To measure the N (3 / L bias, we produce several setsof lensed CMB maps using the Lenspix code . Theseare later combined in different ways to isolate differentcontributions to this bias and to perform consistency androbustness tests. A subset of these simulations are brieflydescribed in Sec. IV B, and here we review the procedure We found consistent results when analyzing maps simulated withthe LenS HAT code [40] which implements a different interpo-lation scheme to resample the unlensed CMB realization at thedisplaced ray position. in more detail. First, we simulate 100 Gaussian realiza-tions of the primordial CMB. Each of these simulationsis then lensed using seven different simulated deflectionfields α eff = ∇ φ + ∇× Ω and adopting the effective remap-ping for the CMB photons as in Eq. (9). The φ and Ωpotentials are obtained from the κ and ω field of FCC18using the consistency relations in Eqs. (7) and (10) inthe harmonic domain. For this operation as well as inthe synthesis of the unlensed CMB, we adopted a ban-dlimit parameter (cid:96) max = 6200. According to the findingsof [40], this setup allows us to recover lensed CMB witha precision of O (10 − ) on scales (cid:96) (cid:46) O (10 − )at (cid:96) ≈ • κ G . A Gaussian realization of convergence withpower spectrum C κ F κ F L . • ± κ F . These simulations measure the bias includingonly the effects of the nonlinear LSS evolution. • ± κ R alone. These simulations measure the bias dueto LSS nonlinearity and PB effects in the conver-gence field. • ± κ R and ± ω R ( ± κ Rω hereafter). They include thefull set of nonlinearity of LSS and PB corrections,including the so-called mixed bispectrum correla-tions κκω and κωω (we refer the reader to [23, 28]for further discussion).We denote the resulting lensed CMB simulations witha given deflection field by a superscript G , ± F , ± R or ± Rω , respectively. For the results described in this pa-per we use maps having an angular resolution of 52 arcsecin Healpix pixelization, corresponding to N side = 4096.On each of these sets we run the lensing reconstructionusing a quadratic estimator and compare them to ex-tract different sources of biases. Each simulation set isdesigned to contain a lensing potential with the samemean power spectrum C φφL . Remaining relative devia-tions from the fiducial C φφL due to post-Born correctionsare below 0 .
2% on the relevant scales considered in thispaper. Hence, in the following, we assume N (0) L and N (1) L to be equal for all simulations. Under this assumptionwe can writeˆ C φφL [ κ ] = 12 L + 1 (cid:88) M ˆ φ † LM ˆ φ LM (31) ≈ C φφL + N (0) L + N (1) L + N (3 / L [ κ ] + O (cid:0) φ , Ω (cid:1) , C φφL extracted from a N -body simulation has a potential bias atsmall angular scales due to the presence of shot noise due to thefinite number of particles in the N -body simulation. Accordingto the estimates of FCC18, the shot noise accounts for roughly15% of the amplitude of the power spectrum on the maximummultipole relevant for this analysis. Because in the followingtext we compare simulated quantities, all including the shot-noise term, the impact of the shot-noise term on the results isexpected to be highly reduced. N (3 / L bias depends on the specificstatistic of the κ field used to lens a specific simulation.We will test the validity of this assumption in Sec. V C.In order to evaluate the bias in a specific experi-mental configuration we add Gaussian noise realizationswith corresponding power spectrum N (cid:96) = σ n B (cid:96) , withwhite noise level, σ n , and a circular Gaussian beam withFHWM size, θ [71], B (cid:96) = exp (cid:18) (cid:96) ( (cid:96) + 1) θ
16 log 2 (cid:19) . (32) B. Measurements of N (3 / L bias To measure the N (3 / L biases from the simulations anddistinguish the contributions to the biases originatingfrom all the different contributions of the κ bispectrumand correlations involving curl modes ( κκω + κωω , PB ω hereafter), we combine the reconstructed CMB lensingpotential power spectrum on each set of lensed CMB re-alizations as follows:LSS: N (3 / L = (cid:68) ˆ C φφL [ κ F ] − ˆ C φφL [ κ G ] (cid:69) Lensed CMB
PB: N (3 / L = (cid:68) ˆ C φφL [ κ R ] − ˆ C φφL [ κ F ] (cid:69) Lensed CMB PB ω : N (3 / L = (cid:68) ˆ C φφL [ κ Rω ] − ˆ C φφL [ κ R ] (cid:69) Lensed CMB
Total: N (3 / L = (cid:68) ˆ C φφL [ κ Rω ] − ˆ C φφL [ κ G ] (cid:69) Lensed CMB ,where we denote in squared brackets the correspondingset of CMB realizations used in the lensing reconstruc-tion. The total bias is equal to the sum of the formerthree, well within the uncertainties shown later in thetext.We report the measurement of the N (3 / L as the av-erage over the 100 lensed CMB simulations at our dis-posal for each deflection field configuration. The errorbars shown in the following figures are computed fromthe dispersion of the lensed CMB simulations and repre-sent the uncertainty on the mean of the simulations. Dueto the fact that the realizations of primordial CMB arethe same for all sets of simulations, we avoid realization-dependent biases (up to bispectrum terms) and cosmicvariance noise. In the following we discuss the impact of N (3 / L bias in terms of the ratio between the bias and thelensing potential power spectrum measured in the FCC18simulations. The reported signal-to-noise ratio (SNR) iscomputed as the ratio between N (3 / L and the error barexpected for a specific experimental configuration (cid:115) L + 1 1 f sky ∆ L (cid:16) N (0) L + N (1) L (cid:17) , (33) where we assume the observed sky fraction to be f sky = 40%, to match the expected sky coverage ofCMB-S4, and the bin size ∆ L ≈ (cid:96) min = 2.In Fig. 4 we show the total N (3 / L bias for theminimum-variance quadratic estimator due to non-Gaussianity in the lensing deflection field, along withthe breakdown of the contribution of each source ofnon-Gaussianity (LSS, PB, and PB ω ). These resultsare derived performing lensing reconstruction using asharp cutoff in harmonic space that removed all theCMB harmonic coefficients having (cid:96) ≥ (cid:96) max = 4000 andassuming an experiment with 1 . µK -arcmin white noisein polarization (1 µK -arcmin in temperature) and a1 arcmin beam size to match CMB-S4 configuration. Wefind that post-Born effects produce a positive bias in thelensing reconstruction, while LSS effects suppress powerin the reconstructed potential. This leads to an impor-tant cancellation of the two effects and, in fact, the total N (3 / L bias becomes a subpercent effect. The amplitudeof N (3 / L , however, changes quite significantly dependingon which combination of the quadratic estimator is usedfor the lensing reconstruction. At low multipoles theindividual relative contributions to the biases induced byLSS and PB can reach up to 7% in the autopower spec-trum of the T T estimator. Generally, the bias amplitudegrows with the number of contributing temperaturefields used in the estimator. For polarization-basedestimators the overall bias can reach 2% for both LSSand PB terms when considered separately. In ourexperimental setup, the polarization-based estimatorsprovide the most important contribution to the mini-mum variance combination below L ≈ L ≈ L + L , while for equilateral configurations, i.e. L ≈ L ≈ L , the bispectrum gets enhanced. Simple1 FIG. 4. Relative biases in the estimated lensing potential power spectrum induced by non-Gaussian statistics of the underlyinglensing potential (black curves) as measured in the FCC18 simulations. This case included lensed CMB modes up to (cid:96) max = 4000and CMB-S4-like experimental configuration. We differentiate the effects caused by nonlinearities of large-scale structures (LSS,purple curve), post-Born lensing effects (PB, orange curve) as well as post-Born mixed bispectrum terms (PB ω , yellow curve)accounting for higher-order correlation between the lensing gradient and curl potential. The shaded areas show the uncertaintyon the bias computed from the dispersion of 100 lensed CMB simulations. arguments can be made to understand why there is asign difference in the bispectrum when all the conver-gence modes are aligned, i.e. in the flattened limit [23].In this case lens-lens deflection, i.e. the deflection of alight ray bundle off two consecutive lenses, dominates.In this case, the first lens induces a contraction of thelight bundle area. This in turn causes the second lens tohave a smaller effect than it would have without the firstlens. This results in an anti-correlation between largeand small scale convergence modes, leading to a negativesign of the bispectrum in the flattened limit. The pos-itive contributions conversely represent a change in thedeflection field along the direction in which the ray is de-flected. A ray passing the edge of an overdensity couldbe deflected towards the center, where the potential gra-dients are larger. This generates more lensing than if thetwo contributions had been added independently and apositive correlation between angular scales. The fact thatthe post-Born and LSS contributions roughly match inamplitude is coincidental and not anymore the case whenthe source plane is at low redshifts [23]. We note, however, that due to the complex convolu-tion of the bispectrum configurations in the quadraticestimator, the details of the cancellations happening onthe N (3 / L bias are nontrivial and their analytical model-ing for the different combinations of quadratic estimatorsis challenging. A more detailed discussion can be foundin [33, 72]. The important cancellation effects betweenthe LSS and PB terms observed for CMB lensing mightnot be as effective in the case of lensing of other diffusebackground emissions that have a redshift kernel peakingat lower redshift, such as the cosmic infrared backgroundor line intensity mapping data [73]. With a shorter lineof sight integration, the relative importance of the post-Born effect is in fact decreased and the LSS term for the N (3 / L bias will become the leading one, thus increasingthe impact of N (3 / L on the reconstructed power spec-trum.The shape of the N (3 / L biases depends not only onthe type of reconstruction channel used, but also on therange of multipoles included in the reconstruction. We2 FIG. 5. The relative suppression or enhancement of the con-vergence bispectrum from large-scale structure nonlinearity(derived using the fitting formula of [67]) due to post-Borneffects, for L = 200 and L = 2000. perform the lensing reconstruction using different cutoffvalues (cid:96) max for the harmonic coefficient used in thelensing reconstruction and show the value of N (3 / L forthe minimum-variance estimator for a cosmic-variancelimited experiment in Fig. 6. Because of the differenceswith analytical predictions discussed in Sec. IV B, wetest the stability of our results with respect to thechoice of (cid:96) min and verified that increasing the cutoff to (cid:96) min = 200 did not affect our results. As expected, wecan observe that the non-Gaussian effects become moreprominent when we include progressively smaller angularscales in the lensing reconstruction. For (cid:96) max = 2000the bias is not detectable and its signal-to-noise ratio issmaller than one. In the case of (cid:96) max = 3000, at smallscales the LSS bias becomes positive, such that the totalbias includes positive contributions from LSS and thepost-Born gradient and curl fields, which causes the pre-viously detected cancellation to fail. The total bias cantherefore reach levels up to 4%, although at multipoles FIG. 6. Dependence of the N (3 / L bias for the minimum-variance lensing estimator on the maximum lensed CMB mul-tipole used in the reconstruction algorithm in the limit of noinstrumental noise. The shaded areas show the uncertaintyon the bias, computed from the dispersion of the 100 simula-tions. with poor SNR. Including progressively smaller scalescauses the LSS terms to increase in amplitude fasterthan the PB term, and as a result, the cancellationbecome less effective, causing the N (3 / L bias to grow.In this scenario the bias becomes very significant and itsSNR could be larger than 10. We warn the reader thatsuch an extreme case serves an illustrative purpose andshould be taken with a grain of salt. In fact, the matterdistribution on scales k ≥ − h affects significantlythe CMB lensing signal at (cid:96) (cid:39) N -bodysimulations used to model the deflection field and theabsence of baryonic effects. These might become moreimportant when analyzing non-Gaussian effects (see,e.g., [74, 75]). Furthermore, one can observe that thecross-bispectrum contribution from the curl potentialdominates at scales (cid:96) (cid:46) N (3 / L gets suppressedcompared to the cosmic-variance limit case. Reducingthe cutoff in power for the reconstruction to (cid:96) max = 3000has a net effect of making the bias practically disappear-ing, despite that the individual LSS and PB effects canbe of order of the error bar.3 FIG. 7. The bias in the reconstructed minimum-variance,TTTT and EBEB lensing power spectrum for an instrumentwith 1 . µK -arcmin white noise in polarization, 1 arcminbeam, and including only lensed CMB multipoles up to (cid:96) max = 3000 in the reconstruction. For comparison the signal-to-noise ratios of the biases in this case ( (cid:96) max = 3000) and thecase with (cid:96) max = 4000 (cf. Fig. 4) are shown in the bottomtwo rows. In Fig. 7 we show a comparison of the SNR obtainedusing these two cutoffs in CMB multipoles. As can beseen in this figure, we observe a rapid increase in thebias amplitudes between the two cases, in particularin the temperature reconstruction channels. Usingpolarization-only lensing reconstruction and comparingthe results with temperature-only reconstruction can beappropriate tools to identify and potentially mitigatethe N (3 / L biases. Since the TTTT reconstruction is themost sensitive for the CMB-S4 experimental configura-tions for L (cid:38) N (3 / L biases while mitigating theloss of sensitivity. The contamination by unresolvedextragalactic foreground residual might in any caseprevent the use of multipoles (cid:96) (cid:29) N (3 / L in the cross-correlation power spectrum between the reconstructedlensing potential and an external large-scale structuretracer. The bias of the cross-spectrum, induced by anonzero CMB lensing potential bispectrum, is mainlycaused by the correlation of the external large-scale struc-ture tracer with the second-order response of the recon-structed lensing potential to the true lensing potential.For the sake of simplicity we limit our analysis to the caseof the cross-correlation with a perfect tracer of the CMB FIG. 8. Summary of cumulative signal-to-noise ratio of thebias for an instrument with 1 . µK -arcmin white noise inpolarization, 1 arcmin beam, and different CMB multipolecutoffs (cid:96) max , comparing temperature (T), polarization (P),and minimum-variance (T+P) estimators. lensing potential, i.e., the lensing potential directly ex-tracted from the FCC18 simulations. Since in the cross-correlation case the tracer is almost uncorrelated withthe CMB, there are fewer contractions of the matter fieldthat contribute to the N (3 / L bias, and thus we shouldsee a reduction in the amplitude of N (3 / L by a factor ofroughly 2 with respect to the bias on the autospectrum,in particular for the TTTT estimator [33]. We verify thatthis prediction holds, as we show in Fig. 9. A similar levelof suppression is observed also for other estimators and,in particular, for EBEB we saw a reduction of more thana factor of 4 for L (cid:38) N (3 / L bias in the au-tospectrum analysis. However, we stress that due to thedistinctive impact of the post-Born term with respect tothe LSS one in the case of CMB lensing, the overall vari-ation in amplitude of the bias in cross-correlation mightchange significantly if a tracer of structures at lower red-shift is considered. Nevertheless, these techniques mightbe affected by other type of biases, such as those dueto the galaxy intrinsic alignements in the case of galaxyweak lensing [76–78]. In addition, the tracers at lowerredshift are in fact more sensitive to the non-Gaussianitydue to matter nonlinearity and less sensitive to post-Borneffects. Therefore we expect to observe an increase in the N / bias as the cancellation between LSS and post-Bornbecomes less effective in this case. We leave the investi-gation of this topic to future work.4 FIG. 9. Top: N (3 / L bias for the reconstructed CMB lensingpotential autospectrum (dashed lines) and cross-correlationwith the input CMB lensing potential of FCC18 simulations(solid lines) for a CMB-S4 experiment and a cutoff in power (cid:96) max = 4000 for the lensing reconstruction with the temper-ature estimator. Bottom: The ratio of the N (3 / L biases forthe cross-correlation and autopower spectrum compared withthe leading order predictions of [33] (dashed black lines). C. Consistency checks
To ensure that the reported biases were not causedby a mismatch in the CMB and lensing potential powerspectra and therefore are not residual N (0) L and N (1) L bi-ases, we check the consistency of our measurements withan alternative method to extract the N (3 / L bias. In par-ticular, we compare the spectra∆ C φφ, L [ κ X ] = (cid:68) ˆ C φφL [ κ X ] − ˆ C φφL [ κ G ] (cid:69)
100 sims (34)∆ C φφ, L [ κ X ] = 12 (cid:68) ˆ C φφL [ κ X ] − ˆ C φφL [ − κ X ] (cid:69)
100 sims , (35)where X ∈ { F, R } . The averaging in Eq. (34) is per-formed over the 100 realizations of lensed CMB derivedwith the set of simulations including a Gaussian conver-gence, and the averaging in Eq. (35) is computed overthe 100 realizations of lensed CMB lensed with the non-Gaussian convergence κ X . We have that (cid:68) ˆ C φφL [ κ X ] (cid:69) ≈ N (0) L (cid:2) C CMB(cid:96) (cid:3) + C φφL (36)+ N (1) L (cid:104) C CMB(cid:96) , C φφL (cid:105) +sgn( κ X ) N (3 / L (cid:104) C CMB(cid:96) , C φφL , b φφφL L L (cid:105) , where we denote in squared brackets the functional de-pendencies of the biases for clarity. Hence both tech-niques in Eqs. (34) and (35) isolate in principle the N (3 / L TABLE I. Global p values of null spectra of noiseless config-uration and (cid:96) max = 3000. p Value [%] LSS TotalMVMV 13.8 47.2TTTT 38.9 52.7TTTE 16.0 17.8TTEE 87.8 96.0TTTB 76.3 87.5TTEB 67.6 99.6TETE 43.6 47.6TEEE 5.3 9.2TETB 97.8 86.0TEEB 83.6 98.9EEEE 30.5 27.9EETB 21.7 20.7EEEB 46.1 49.0TBTB 45.2 20.4TBEB 60.1 73.4EBEB 17.3 51.5 bias. However, a mismatch of N (0) L and N (1) L betweensimulations lensed with κ F , κ R and κ G or correlations atorder higher than the bispectrum should manifest them-selves in a discrepancy between the two spectra. We con-structed null spectra and computed Welch’s t -test statis-tics for both the κ F and the κ R sets of simulations to testseparately LSS effects alone and LSS and PB together.In both cases we use the spectra from the three most rel-evant reconstruction channels (TTTT, EBEB, and theautopower spectrum of the minimum-variance estimator(MVMV)) binned in 21 bins within L ∈ [30 , C φφ, L − ∆ C φφ, L in Fig. 10. The deviations fromzero in the high signal-to-noise regions are subdominant,while small deviations at mostly large multipoles are wellwithin the 1 σ error bar. We furthermore obtained global p values by averaging over the bins for each estimator andfind no PTE lower than 5%, as summarized in Table I.These results made us conclude that the simulation andreconstruction pipeline up to the lensing power spectrumstep are internally consistent, increasing our confidencein the results shown in Sec. V B. VI. N / L IMPACT ON COSMOLOGICALPARAMETER ESTIMATION
Future sensitive measurements of the CMB lensing po-tential will provide important constraints on cosmologi-cal parameters. Therefore a biased reconstruction of thelensing potential power spectrum could affect their es-timation. For example, we find that at the high sensi-tivities envisioned for CMB-S4 measurements the total5
FIG. 10. The null spectra obtained taking the difference be-tween ∆ C φφ, L and ∆ C φφ, L as defined in Eqs. (34) and (35)for the minimum-variance, TTTT, EEEE, and EBEB lensingreconstruction in the limit of no instrumental noise. The re-constructions on κ F -lensed CMB fields are shown in purple(LSS only contribution), the same with κ R are shown in or-ange [LSS and PB (total) contributions]. The error bars showthe uncertainties as measured from the scatter in the simu-lations while the shaded area show the expected statisticaluncertainty in the respective bin. N (3 / L bias could produce deviations of more than 3 σ from the fiducial value of 1 when fitting the lensing am-plitude parameter A lens . In Table II we show the fit-ted A lens parameter for different CMB multipole cutoffsobtained by maximizing the simple one-parameter likeli-hood defined by − L = (cid:88) L (2 L +1) f sky (cid:18) ln (cid:18) C L D L (cid:19) + D L C L − (cid:19) , (37)where C L = A lens × C fid. L + N φφ, tot. (cid:96) , D L = C fid. L + N φφ, tot. (cid:96) + N (3 / L , and N φφ, tot. L = N (0) L + N (1) L .Because of the nontrivial scale dependence of the N / L bias, we expand our cosmological parameter estimationstudy to the exploration of a broader parameter spaceusing Markov chain Monte Carlo (MCMC) techniques.The goal is to quantify the significance of possible biasesin parameters like the total neutrino mass, M ν , or the TABLE II. Fitted A lens parameter of the biased reconstructedlensing power spectrum with a fiducial value of A lens = 1 fortemperature-only (T), polarization-only (P), and minimumvariance (T+P) lensing estimators and no noise in the CMB.Cases with significant bias are marked in bold. Total bias (cid:96) max = 3000 (cid:96) max = 4000 (cid:96) max = 5000T 0 . ± . . ± .
003 0 . ± . P . ± .
002 1 . ± .
001 1 . ± . T+P . ± .
002 1 . ± .
001 0 . ± . amplitude of primordial inflationary perturbations, A s ,if N / L is unaccounted for in the power-spectra modelingand cosmological parameters sampling. For this purposewe use the publicly available package MontePython [79, 80] based on the Metropolis-Hastings sampling algo-rithm. In this analysis we consider the CMB and lensinglikelihood for a set of parameters θ given the measuredpower spectra of CMB temperature, E -modes and lens-ing potential as Gaussian in the respective fields. Un-der these assumptions the likelihood function is given by(e.g., [81]) − L ( θ | ˆ C ) = (cid:88) (cid:96) (2 (cid:96) +1) f sky (cid:32) ln | C (cid:96) || ˆ C (cid:96) | + C − (cid:96) ˆ C (cid:96) − (cid:33) , (38)where the covariance matrix for the fiducial model ˆ C (cid:96) and the theoretical signal C (cid:96) are constructed as C (cid:96) = C T T(cid:96) + N T T(cid:96) C T E(cid:96) C T φ(cid:96) C T E(cid:96) C EE(cid:96) + N EE(cid:96) C T φ(cid:96) C φφ(cid:96) + N φφ, tot. (cid:96) , where N TT (cid:96) and N EE (cid:96) are the white noise power spec-tra for the temperature and the E -modes and N φφ, tot. (cid:96) = N (0) (cid:96) + N (1) (cid:96) . All these quantities are computed assum-ing the fiducial cosmology with CMB-S4 sensitivities andconsidered to be independent of the cosmological param-eters in order to simplify and speed up the sampling. Inthe evaluation of the fiducial ˆ C (cid:96) we use the biased lens-ing potential power spectrum, ˜ C φφL , which includes the N (3 / L bias measured in the simulations and depends onthe cosmological parameters of the fiducial model as˜ C φφL ≡ C φφL [ θ fid. ] + N (0) L + N (1) L + C φφL [ θ fid. ] C φφ sims L N (3 / L . (39)This definition allows one to mitigate the impact ofthe shot-noise term and the difference in the modeling ofthe nonlinear evolution between the simulation resultsand the Boltzmann solvers which typically employ theHalofit fitting formulas [69]. This enables us to have aconsistent modeling of nonlinearity between the fiducialand the fitted model, reducing the chance to obtainspurious results in the fitting that are driven by thedifferences in the CMB lensing potential power spectrummodeling. We note, however, that the uncertaintiesin the modeling of nonlinearity on the CMB lensingpower spectrum reach the 10 – 15% level on the scalesconsidered in this work [23, 69] and might becomenon-negligible. http://baudren.github.io/montepython.html TABLE III. The cosmological parameters from Planck 2015[85, 86] together with their 1 σ proposal scale or parameterbounds used in the cosmological parameter inference.Ω b h . ± . c h . ± . τ . ± . A s . ± . n s . ± . θ s . ± . M ν [meV] [0 , In the construction of the covariance we neglect the φE correlation because it is confined at very large angularscales and carries little information on the parametersof interest in our analysis. For the sake of simplicitywe do not include the B -mode power spectrum in ˆ C (cid:96) and C (cid:96) to avoid the need to model the non-Gaussiancovariance between C BB(cid:96) and C φφL [82]. We note thatmore optimal formalisms to deal with the non-Gaussiancorrelations between CMB and lensing power spectrahave been discussed in the literature [47, 83, 84]. Asthe present analysis is intended to quantify biases oncosmological parameters estimation due to mismodelingof the lensing potential bias rather than to provide anaccurate forecast of future CMB experiment constraints,the approximations adopted here are not expected toaffect our conclusions at the level of accuracy consideredin this work.In the likelihood construction we assume a fiducialΛCDM cosmology taken from Planck 2015 results [85, 86]devoid of massive neutrinos, while we allow for a singleneutrino to be massive in the parameter fit. We includeangular scales 30 ≤ (cid:96) ≤ f sky = 40% to mimic a CMB-S4-like surveywith 1 . µK -arcmin white noise in polarization and a1 arcmin beam size in the likelihood. We summarize thevalues of our fiducial cosmology as well as the detailsof the priors adopted for the cosmological parameterssampled in our analysis in Table III.We neglect the effects of the LSS non-Gaussianity andpost-Born corrections on the lensed CMB TT and EEpower spectra since the cumulative signal-to-noise ratiofor these corrections is below the detection thresholdseven for CMB-S4 sensitivity. In Fig. 11 we show the 2Dposteriors obtained for the parameter combinations ofΩ c h , ln (cid:0) A s (cid:1) , and M ν for the minimum-variancelensing estimator and CMB-S4 sensitivity. The figureshows an example of the main biases in the parameter es-timation induced by different sources (LSS, PB, total) ofunaccounted N (3 / L bias. Similar to what was observedin Sec.V B, the compensating effect between the LSSand PB biases observed at the level of the lensing powerspectrum is also visible in the cosmological parameterestimation, where we find a cancellation of the parameter biases when both these terms are included. Each sourceof N (3 / L bias might considerably affect the estimationof the cosmological parameters when considered alone atthe level of CMB-S4 sensitivity. Assuming we can modelthese biases analytically we need to include both theterms in the modeling as the inclusion of only one of theLSS or PB term would lead to an overcorrection of theeffect. This is clearly visible in the case the LSS-induced N (3 / L for A s and M ν , where the large negative biasover a large range of scales in the power spectrum causesa significant false detection of a 169 +50 − meV neutrinomass. The cancellation due to post-Born correctionsmitigates this bias, reducing it to 83 +40 − , and hence stillcompatible with zero neutrino mass only at the 2 σ level. FIG. 11. The 2D posteriors for the cold dark matter den-sity, Ω c h , the amplitude of primordial inflationary perturba-tions, A s and the neutrino mass, M ν , including biases fromLSS nonlinearities and post-Born effect in ˆ C φφL , reconstructedusing the minimum variance estimator, CMB modes up to (cid:96) max = 4000 and CMB-S4 experimental specifications. The same analysis carried out adding only the N (3 / L biases of polarization-based estimators indicates thatusing these reconstruction channels leads to more robustconstraints on cosmological parameters, even whenincluding the smaller angular scales in the lensing recon-struction. In Table IV and Fig. 12 we show the best-fitvalues and marginalized posteriors obtained includingthe total N (3 / L computed varying the CMB multipolecutoff used in the reconstruction for the two differentcases including and excluding temperature data whenforming the minimum-variance estimator. Includingmultipoles up to (cid:96) = 5000 in the reconstruction leadsto a neutrino mass bias larger than 1 σ , even after7 TABLE IV. This table shows the deviation of the best-fit from the fiducial values (bias) and 68% confidence level (1 σ )uncertainties for the cold dark matter density, Ω c h , the optical depth to reionization, τ , the amplitude of primordial inflationaryperturbations, A s and the neutrino mass M ν . A configuration with 1 . µK -arcmin white noise and 1 arcmin beam with differentCMB multipole cutoff and estimator combinations was used. We show biases using minimum-variance lensing reconstructionincluding CMB temperature (T+P) and using polarization only (P). Upper limits are given in terms of 95% confidence level. (cid:96) max = 3000 (cid:96) max = 4000 (cid:96) max = 5000 (cid:96) max = 5000+DESI Bias 1 σ (stat.) Bias 1 σ (stat.) Bias 1 σ (stat.) Bias 1 σ (stat.)T+P Ω c h ×
25 85 14 88 −
45 85 −
66 55 τ × (cid:0) A s (cid:1) ×
11 15 18 18 27 16 16 14 M ν [meV] 0 79 90 60 110 50 0 55 P Ω c h ×
16 84 26 82 25 80 −
37 56 τ × (cid:0) A s (cid:1) ×
12 16 13 16 14 15 15 16 M ν [meV] 0 75 0 84 65 60 0 44FIG. 12. The one-dimensional posteriors for the total neutrino mass M ν for different CMB multipole cutoffs used in the lensingreconstruction. (cid:96) max = 3000 case is shown as a solid line, while (cid:96) max = 4000 and (cid:96) max = 5000 are shown as dashed and dottedlines, respectively. The left figure shows the results obtained including all reconstruction estimators including temperature(T+P), while the right figure uses only polarization-based estimators (P). Each figure also includes the posterior after includinga prior using DESI BAO data [87] in the sampling in green, for the most extreme case of (cid:96) max = 5000. excluding temperature data. Nevertheless, on the levelof the parameter estimation we can observe that thepolarization lensing estimator is more robust to thesekinds of biases, which can be attributed in part tothe slightly worse reconstruction lensing noise whenexcluding small-scale temperature data and in part tothe smaller amplitude of the N (3 / L bias for polarizationestimators.We note, however, that the error on the total neu-trino mass does not decrease significantly with decreas-ing noise in the CMB lensing potential power spectrum.This is due to the degeneracy of the total neutrino masswith the A s parameter and the sensitivity of the con-straint on the latter (or, more precisely, on the combi-nation A s e − τ ). Since we are assuming future data fromground-based CMB-S4 instruments, which are limited tomultipoles (cid:96) ≥
30, we are not able to push the uncertainty on τ to the cosmic-variance limit. However, accessingthe reionization bump at (cid:96) ≤
20 down to cosmic-varianceprecision could be achieved by proposed all-sky polarizedCMB surveys like CLASS [88], CORE [81], LiteBIRD [89]or PIXIE [90]. This would provide a tighter constraint on τ [91], and would lead to the expected decrease in statis-tical uncertainty with increasing multipole cut-off in thelensing reconstruction. Furthermore, a N (3 / -bias in thelensing potential estimation would bias A s and τ high.We observe that being able to include the constrainingpower of the reionization bump at large-scales would re-duce the bias on τ , and consequently significantly reducethe bias on the total neutrino mass. This would occur atthe expense of a ≤ − σ total bias on cold-dark matterdensity Ω m and a negligibly larger χ goodness of fit.8 VII. CONCLUSIONS
In this work we investigate the properties of higher-order correlations of the CMB lensing deflection fieldarising from nonlinear evolution of the matter as wellas post-Born corrections, modeled through numericalsimulations, and their impact on the CMB lensingpotential reconstruction using quadratic estimators( N (3 / L bias). We validate the numerical simulationsused to model these effects comparing the expectedcorrections on the lensed CMB power spectrum dueto both LSS nonlinearity and post-Born correctionsmodeled analytically, finding a good agreement. We findthat both the matter nonlinearity and post-Born non-Gaussianity cause significant biases of the reconstructedCMB lensing potential power spectrum. However,when these effects are analyzed jointly, the amplitudeof the total N (3 / L bias is greatly reduced both on theCMB lensing autospectrum and in the cross-correlation.This is directly related to the different shape and signproperties of the post-Born bispectrum and the matterbispectrum. The cancellation is more effective in thepresence of experimental noise. Despite this fact, wefind that the estimation of the A lens parameter from theCMB lensing potential could be biased by more than 3 σ for future high-sensitivity experiments like CMB-S4.We further perform a MCMC analysis to evaluate theimpact of the residual N (3 / L bias on the estimation ofother cosmological parameters at the CMB-S4 sensi-tivity. We find that the best-fit value of cosmologicalparameters such as M ν and A s could be biased dueto the N (3 / L bias by up to 2 σ , but the significance ofthese biases greatly depends on the type of quadraticestimator and the maximum multipole used for thelensing reconstruction. Using multipoles (cid:96) ≤ N (3 / L effects. The latter, in particular, is causedby a consistently more effective cancellation of LSS andpost-Born effects. As an illustrative case, we perform thecosmological parameter analysis including multipoles up to (cid:96) = 5000. In this case we find a shift of the likelihoodpeak causing a detection of a nonzero neutrino mass atthe 2 σ level when including all the lensing reconstructionchannels. The inclusion of external data sets such asDESI BAO seems to help remove the biases, though 1 σ tensions might still remain. Nevertheless, based on theresults above, we could expect inconsistencies betweenthe inferred neutrino mass estimates from differentdata sets, if the N (3 / L bias is not accounted for in theparameter estimation for future, high-sensitivity/high-resolution CMB experiments. Finally, we find that the N (3 / L bias in the cross-correlation with a perfect tracerof the CMB lensing potential is reduced by a factor ofroughly 2 with respect to the bias on the autospectrum,in agreement with the prediction of [33]. The biasobserved in cross-correlation with lower-redshift tracersmight, however, be different due to the different weightthat the post-Born term has for lower redshift probes,but we leave the investigation of this aspect to futurework.During the final stage of this work we compared ourresults on the N (3 / L -bias with those of [72], who alsoperformed a similar analysis using a CMB lensing fieldextracted from different N -body simulations. Despitetheir N -body simulations differing in resolution and boxsize, and the simulated sky area used for the lensingreconstruction being smaller than the full-sky resultsof our work, we find similar conclusions. This suggeststhat despite some quantitative conclusion of this workmight still be simulation dependent and more complexphysical effects are excluded from our modeling, thehigher-order effects in CMB lensing should be treatedcarefully in future analyses in order to exploit the fullscientific capacity of a CMB-S4-like observation. ACKNOWLEDGMENTS
We thank Vanessa B¨ohm and Blake Sherwin for use-ful discussions and the comparison of numerical results.We thank Antony Lewis for providing some of the ana-lytical results used in the comparison with the resultsof this work. We thank Radek Stompor, Silvia Galliand Luca Pagano for useful discussions and commentson the manuscript. G.F. acknowledges the support of theCNES postdoctoral program. We acknowledge the use of camb [95, 96], class [97], MontePython [98, 99], GetDist , Lenspix and Healpix [100] software http://camb.info http://github.com/lesgourg/class_public http://github.com/brinckmann/montepython_public http://github.com/cmbant/getdist http://cosmologist.info/lenspix http://healpix.sourceforge.net [1] N. Aghanim, S. Majumdar, and J. Silk, Rep. Prog.Phys. , 066902 (2008), arXiv:0711.0518.[2] R. K. Sachs and A. M. Wolfe, Astrophys. J. , 73(1967).[3] M. J. Rees and D. W. Sciama, Nature (London) ,511 (1968).[4] A. Lewis and A. Challinor, Physics Reports , 1(2006), arXiv:astro-ph/0601594.[5] C. L. Reichardt et al. , Astrophys. J. , 1200 (2009),arXiv:0801.1491.[6] R. Keisler et al. , Astrophys. J. , 28 (2011),arXiv:1105.3182.[7] K. M. Smith, O. Zahn, and O. Dor´e, Phys. Rev. D ,043510 (2007), arXiv:0705.3980.[8] C. M. Hirata, S. Ho, N. Padmanabhan, U. Seljak,and N. A. Bahcall, Phys. Rev. D , 043520 (2008),arXiv:0801.0644.[9] S. Das et al. , Phys. Rev. Lett. , 021301 (2011),arXiv:1103.2124.[10] A. van Engelen et al. , Astrophys. J. , 142 (2012),arXiv:1202.0546.[11] P. A. R. Ade et al. (Planck Collaboration), A&A ,A15 (2016), arXiv:1502.01591.[12] P. A. R. Ade et al. (POLARBEAR Collaboration),Phys. Rev. Lett. , 021301 (2014), arXiv:1312.6646.[13] K. T. Story et al. , Astrophys. J. , 50 (2015),arXiv:1412.4760.[14] B. D. Sherwin et al. , Phys. Rev. D , 123529 (2017).[15] P. A. R. Ade et al. (POLARBEAR Collabora-tion), Physical Review Letters , 131302 (2014),arXiv:1312.6645.[16] D. Hanson et al. , Phys. Rev. Lett. , 141301 (2013),arXiv:1307.5830.[17] P. A. R. Ade et al. (POLARBEAR Collaboration), As-trophys. J. , 171 (2014), arXiv:1403.2369.[18] P. A. R. Ade et al. (POLARBEAR Collaboration), As-trophys. J. , 121 (2017), arXiv:1705.02907.[19] R. Keisler et al. , Astrophys. J. , 151 (2015),arXiv:1503.02315.[20] K. N. Abazajian et al. , (2016), arXiv:1610.02743.[21] P. M. Merkel and B. M. Sch¨afer, Mon. Not. R. Astron.Soc. , 1067 (2011), arXiv:1007.1408.[22] T. Namikawa, Phys. Rev. D , 121301 (2016),arXiv:1604.08578.[23] G. Pratten and A. Lewis, J. Cosmol. Astropart. Phys. , 047 (2016), arXiv:1605.05662.[24] C. Carbone, V. Springel, C. Baccigalupi, M. Bartel-mann, and S. Matarrese, Mon. Not. R. Astron. Soc. , 1618 (2008), arXiv:0711.2655.[25] C. Carbone, C. Baccigalupi, M. Bartelmann, S. Matar-rese, and V. Springel, Mon. Not. R. Astron. Soc. ,668 (2009), arXiv:0810.4145.[26] S. Hagstotz, B. M. Sch¨afer, and P. M. Merkel, Mon.Not. R. Astron. Soc. , 831 (2015), arXiv:1410.8452.[27] A. Lewis and G. Pratten, J. Cosmol. Astropart. Phys. , 003 (2016), arXiv:1608.01263. [28] G. Fabbian, M. Calabrese, and C. Carbone, J. Cosmol.Astropart. Phys. , 050 (2018), arXiv:1702.03317.[29] A. Cooray and W. Hu, Astrophys. J. , 19 (2002),arXiv:astro-ph/0202411.[30] C. Shapiro and A. Cooray, J. Cosmol. Astropart. Phys. , 007 (2006), arXiv:astro-ph/0601226.[31] E. Krause and C. M. Hirata, A&A , A28 (2010),arXiv:0910.3786.[32] A. Petri, Z. Haiman, and M. May, Phys. Rev. D ,123503 (2017), arXiv:1612.00852.[33] V. B¨ohm, M. Schmittfull, and B. D. Sherwin, Phys.Rev. D , 043519 (2016), arXiv:1605.01392.[34] M. R. Becker, Mon. Not. R. Astron. Soc. , 115(2013), arXiv:arXiv:1210.3069.[35] M. Bartelmann and P. Schneider, Physics Reports ,291 (2001), arXiv:astro-ph/9912508.[36] A. Stebbins, (1996), arXiv:astro-ph/9609149.[37] W. Hu, Phys. Rev. D , 043007 (2000), arXiv:astro-ph/0001303.[38] G. Marozzi, G. Fanizza, E. Di Dio, and R. Dur-rer, J. Cosmol. Astropart. Phys. , 028 (2016),arXiv:1605.08761.[39] C. M. Hirata and U. Seljak, Phys. Rev. D , 043001(2003), arXiv:astro-ph/0209489.[40] G. Fabbian and R. Stompor, A&A , A109 (2013),arXiv:1303.6550.[41] W. Hu and T. Okamoto, Astrophys. J. , 566 (2002),arXiv:astro-ph/0111606.[42] T. Okamoto and W. Hu, Phys. Rev. D , 083002(2003), arXiv:astro-ph/0301031.[43] C. M. Hirata and U. Seljak, Phys. Rev. D , 083002(2003), arXiv:astro-ph/0306354.[44] J. Carron and A. Lewis, Phys. Rev. D , 063510 (2017),arXiv:1704.08230.[45] M. Millea, E. Anderes, and B. D. Wandelt, (2017),arXiv:1708.06753.[46] D. Hanson, A. Challinor, G. Efstathiou, andP. Bielewicz, Phys. Rev. D , 043005 (2011),arXiv:1008.4403.[47] J. Peloton, M. Schmittfull, A. Lewis, J. Carron,and O. Zahn, Phys. Rev. D , 043508 (2017),arXiv:1611.01446.[48] A. Lewis, A. Challinor, and D. Hanson, J. Cosmol.Astropart. Phys. , 018 (2011), arXiv:1101.2234.[49] M. Kesden, A. Cooray, and M. Kamionkowski, Phys.Rev. D , 123507 (2003), arXiv:astro-ph/0302536.[50] T. Namikawa, D. Hanson, and R. Takahashi, Mon. Not.R. Astron. Soc. , 609 (2013), arXiv:1209.0091.[51] T. Namikawa and R. Takahashi, Mon. Not. R. Astron.Soc. , 1507 (2014), arXiv:1310.2372.[52] E. Anderes, Phys. Rev. D , 083517 (2013),arXiv:1301.2576.[53] G. Marozzi, G. Fanizza, E. Di Dio, and R. Durrer, Phys.Rev. D , 023535 (2018), arXiv:arXiv:1612.07263.[54] C. Carbone, M. Petkova, and K. Dolag, J. Cosmol.Astropart. Phys. , 034 (2016), arXiv:1605.02024. [55] E. Castorina, C. Carbone, J. Bel, E. Sefusatti, andK. Dolag, J. Cosmol. Astropart. Phys. , 043 (2015),arXiv:1505.07148.[56] M. Calabrese, C. Carbone, G. Fabbian, M. Baldi, andC. Baccigalupi, J. Cosmol. Astropart. Phys. , 049(2015), arXiv:1409.7680.[57] P. Fosalba, E. Gazta˜naga, F. J. Castander, andM. Manera, Mon. Not. R. Astron. Soc. , 435 (2008),arXiv:0711.1540.[58] B. Jain, U. Seljak, and S. White, Astrophys. J. ,547 (2000), arXiv:astro-ph/9901191.[59] T. Hamana and Y. Mellier, Mon. Not. R. Astron. Soc. , 169 (2001), arXiv:astro-ph/0101333.[60] P. Schneider, J. Ehlers, and E. E. Falco, Astronomyand Astrophysics Library (Springer-Verlag, New York,1992) p. 112.[61] S. Das and P. Bode, Astrophys. J. , 1 (2008),arXiv:0711.3793.[62] S. Hilbert, J. Hartlap, S. D. M. White, and P. Schnei-der, A&A , 31 (2009), arXiv:0809.5035.[63] A. Lewis, A. Hall, and A. Challinor, J. Cosmol. As-tropart. Phys. , 023 (2017), arXiv:1706.02673.[64] M. Srednicki, ApJ Letters , L1 (1993), arXiv:astro-ph/9306012.[65] E. Komatsu and D. N. Spergel, Phys. Rev. D , 063002(2001), arXiv:astro-ph/0005036.[66] F. Bernardeau, S. Colombi, E. Gazta˜naga, and R. Scoc-cimarro, Physics Reports , 1 (2002), arXiv:astro-ph/0112551.[67] R. Scoccimarro and H. M. P. Couchman, Mon. Not. R.Astron. Soc. , 1312 (2001), arXiv:astro-ph/0009427.[68] R. Ruggeri, E. Castorina, C. Carbone, and E. Se-fusatti, J. Cosmol. Astropart. Phys. , 003 (2018),arXiv:1712.02334.[69] R. Takahashi, M. Sato, T. Nishimichi, A. Taruya,and M. Oguri, Astrophys. J. , 152 (2012),arXiv:1208.2701.[70] H. Gil-Mar´ın, C. Wagner, F. Fragkoudi, R. Jimenez,and L. Verde, J. Cosmol. Astropart. Phys. , 047 (2012),arXiv:1111.4477.[71] L. Knox, Phys. Rev. D , 4307 (1995), arXiv:astro-ph/9504054.[72] V. B¨ohm, B. D. Sherwin, J. Liu, J. C. Hill, M. Schmit-tfull, and T. Namikawa, (2018), arXiv:1806.01157.[73] E. Schaan, S. Ferraro, and D. N. Spergel, Phys. Rev.D , 123539 (2018), arXiv:1802.05706.[74] T. Castro, M. Quartin, C. Giocoli, S. Borgani, andK. Dolag, Mon. Not. Roy. Astron. Soc. , 1305 (2018),arXiv:1711.10017.[75] A. Natarajan, A. R. Zentner, N. Battaglia, and H. Trac,Phys. Rev. D , 063516 (2014), arXiv:1405.6205.[76] P. M. Merkel and B. M. Sch¨afer, Mon. Not. R. Astron.Soc. , 2431 (2017), arXiv:1709.04444.[77] P. Larsen and A. Challinor, Mon. Not. R. Astron. Soc. , 4343 (2016), arXiv:1510.02617.[78] M. A. Troxel and M. Ishak, Phys. Rev. D , 063528 (2014), arXiv:1401.7051.[79] T. Brinckmann and J. Lesgourgues, (2018),arXiv:1804.07261.[80] B. Audren, J. Lesgourgues, K. Benabed, andS. Prunet, J. Cosmol. Astropart. Phys. , 001(2013), arXiv:1210.7183.[81] E. Di Valentino, T. Brinckmann, M. Gerbino, V. Poulin,F. R. Bouchet, J. Lesgourgues, A. Melchiorri, J. Chluba,S. Clesse, J. Delabrouille, et al. , J. Cosmol. Astropart.Phys. , 017 (2018), arXiv:1612.00021.[82] K. M. Smith, W. Hu, and M. Kaplinghat, Phys. Rev.D , 043002 (2004), astro-ph/0402442.[83] M. M. Schmittfull, A. Challinor, D. Hanson,and A. Lewis, Phys. Rev. D , 063012 (2013),arXiv:1308.0286.[84] A. Benoit-L´evy, K. M. Smith, and W. Hu, Phys. Rev.D , 123008 (2012), arXiv:1205.0474.[85] P. A. R. Ade et al. (Planck Collaboration), Astron. As-trophys. , A13 (2016), arXiv:1502.01589.[86] R. Adam et al. (Planck), Astron. Astrophys. , A108(2016), arXiv:1605.03507.[87] M. Archidiacono, T. Brinckmann, J. Lesgourgues, andV. Poulin, J. Cosmol. Astropart. Phys. , 052(2017), arXiv:1610.09852.[88] D. J. Watts et al. , Astrophys. J. , 121 (2018),arXiv:1801.01481 [astro-ph.CO].[89] A. Suzuki et al. (2018) arXiv:1801.06987.[90] E. Calabrese, D. Alonso, and J. Dunkley, Phys. Rev. D , 063504 (2017), arXiv:1611.10269.[91] R. Allison, P. Caucal, E. Calabrese, J. Dunkley,and T. Louis, Phys. Rev. D , 123535 (2015),arXiv:1509.07471.[92] A. van Engelen, S. Bhattacharya, N. Sehgal, G. P.Holder, O. Zahn, and D. Nagai, Astrophys. J. ,13 (2014), arXiv:1310.7023.[93] S. J. Osborne, D. Hanson, and O. Dor´e, J. Cosmol.Astropart. Phys. , 024 (2014), arXiv:1310.7547.[94] S. Ferraro and J. C. Hill, Phys. Rev. D , 023512(2018), arXiv:1705.06751.[95] A. Lewis, A. Challinor, and A. Lasenby, Astrophys. J. , 473 (2000), arXiv:astro-ph/9911177.[96] C. Howlett, A. Lewis, A. Hall, and A. Challi-nor, J. Cosmol. Astropart. Phys. , 027 (2012),arXiv:1201.3654.[97] D. Blas, J. Lesgourgues, and T. Tram, J. Cosmol. As-tropart. Phys. , 034 (2011), arXiv:1104.2933.[98] T. Brinckmann and J. Lesgourgues, (2018),arXiv:1804.07261.[99] B. Audren, J. Lesgourgues, K. Benabed, andS. Prunet, J. Cosmol. Astropart. Phys. , 001 (2013),arXiv:1210.7183.[100] K. M. G´orski, E. Hivon, A. J. Banday, B. D. Wandelt,F. K. Hansen, M. Reinecke, and M. Bartelmann, As-trophys. J.622