CMB lensing bi-spectrum: assessing analytical predictions against full-sky lensing simulations

TL;DR

This paper evaluates the accuracy of analytical models predicting the CMB lensing bi-spectrum by comparing them with full-sky simulations, highlighting their strengths and limitations across different scales and configurations.

Contribution

It provides a detailed comparison of analytical predictions and simulations for the CMB lensing bi-spectrum, including the effects of perturbation theory and fitting formulas, and discusses implications for modified gravity theories.

Findings

Tree-level perturbation theory matches simulations up to 200.

One-loop order captures results up to 600.

Fitting formula predictions agree within 10% for equilateral configurations up to 2000.

Abstract

Cosmic microwave background (CMB) lensing is an integrated effect whose kernel is greater than half the peak value in the range $1<z<5$. Measuring this effect offers a powerful tool to probe the large-scale structure of the Universe at high redshifts. With the increasing precision of ongoing CMB surveys, other statistics than the lensing power spectrum, in particular the lensing bi-spectrum, will be measured at high statistical significance. This will provide ways to improve the constraints on cosmological models and lift degeneracies. Following on an earlier paper, we test analytical predictions of the CMB lensing bi-spectrum against full-sky lensing simulations, and discuss their validity and limitation in detail. The tree-level prediction of perturbation theory agrees with the simulation only up to $\ell\sim 200$, but the one-loop order allows capturing the simulation results up to $\ell\sim 600$. We also show that analytical predictions based on fitting formulas for the matter bi-spectrum agree reasonably well with simulation results, although the precision of the agreement depends on the configurations and scales considered. For instance, the agreement is at the $10\%$-level for the equilateral configuration at multipoles up to $\ell\sim2000$, but the difference in the squeezed limit raises to more than a factor of two at $\ell\sim2000$. This discrepancy appears to come from limitations in the fitting formula of the matter bi-spectrum. We also find that the analytical prediction for the post-Born correction to the bi-spectrum is in good agreement with the simulation. We conclude by discussing the bi-spectrum prediction in some theories of modified gravity.

Figures

Equations

$$B_{{\ell }_1{\ell }_2{\ell }_3}\equiv \sum _{m_1{m}_2{m}_3}\left(\begin{array}{ccc}{\ell }_1& {\ell }_2& {\ell }_3\\ m_1& m_2& m_3\end{array}\right)⟨a_{{\ell }_1{m}_1}{a}_{{\ell }_2{m}_2}{a}_{{\ell }_3{m}_3}⟩.$$

$$b_{{\ell }_1{\ell }_2{\ell }_3}\equiv h_{{\ell }_1{\ell }_2{\ell }_3}^{-1}{B}_{{\ell }_1{\ell }_2{\ell }_3}({\ell }_1+{\ell }_2+{\ell }_3\text{ is even}),$$

$$h_{{\ell }_1{\ell }_2{\ell }_3}\equiv \sqrt{\frac{(2{\ell }_1+1)(2{\ell }_2+1)(2{\ell }_3+1)}{4\pi }}\left(\begin{array}{ccc}{\ell }_1& {\ell }_2& {\ell }_3\\ 0& 0& 0\end{array}\right).$$

$$⟨a_{{\mathbf{\ell }}_1}{a}_{{\mathbf{\ell }}_2}{a}_{{\mathbf{\ell }}_3}⟩=(2\pi {)}^2\delta ({\mathbf{\ell }}_1+{\mathbf{\ell }}_2+{\mathbf{\ell }}_3)B_{{\mathbf{\ell }}_1{\mathbf{\ell }}_2{\mathbf{\ell }}_3},$$

$$B_{{\ell }_1{\ell }_2{\ell }_3}$$

$$\varphi (\hat{\mathbf{n}})=-2{\int }_0^{{\chi }_{\ast }}\mathrm{d}\chi W(\chi ,{\chi }_{\ast })\Psi (\chi ,\hat{\mathbf{n}}).$$

$$W(\chi ,{\chi }_{\ast })=\frac{{\chi }_{\ast }-\chi }{\chi {\chi }_{\ast }}\Theta ({\chi }_{\ast }-\chi ),$$

$$B_{{\mathbf{\ell }}_1{\mathbf{\ell }}_2{\mathbf{\ell }}_3}^{\mathrm{LSS}}={\int }_0^{{\chi }_{\ast }}\mathrm{d}\chi$$

$$\times$$

$$B_{\delta }({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3,\chi )$$

$$+2perms.,$$

$$F_2({\mathbf{k}}_1,{\mathbf{k}}_2,z)$$

$$+\frac{1}{2}\frac{k_1^2+k_2^2}{k_1{k}_2}b(k_1,z)b(k_2,z)\cos \theta$$

$$+\frac{2}{7}c(k_1,z)c(k_2,z){\cos }^2\theta ,$$

$$B_{{\mathbf{\ell }}_1{\mathbf{\ell }}_2{\mathbf{\ell }}_3}^{\mathrm{pb}}$$

$$+cyc.perm.,$$

$$M_{\ell {\ell }^{\prime }}$$

$$\times {\int }_0^{\chi }\mathrm{d}{\chi }^{\prime }\frac{W({\chi }^{\prime },\chi )W({\chi }^{\prime },{\chi }_{\ast })}{({\chi }^{\prime }{)}^2}{P}_{\Psi }\left(\frac{{\ell }^{\prime }}{{\chi }^{\prime }},{\chi }^{\prime }\right),$$

$$b_{b_1{b}_2{b}_3}\equiv \frac{{\sum }_{{\ell }_i}^{b_i}{h}_{{\ell }_1{\ell }_2{\ell }_3}^2{b}_{{\ell }_1{\ell }_2{\ell }_3}}{{\sum }_{{\ell }_i}^{b_i}{h}_{{\ell }_1{\ell }_2{\ell }_3}^2},$$

$$b_{{\ell }_1{\ell }_2{\ell }_3}$$

$$=h_{{\ell }_1{\ell }_2{\ell }_3}^{-1}\sum _{m_1{m}_2{m}_3}{h}_{{\ell }_1{\ell }_2{\ell }_3}^{-1}\int {\mathrm{d}}^2\hat{\mathbf{n}}{Y}_{{\ell }_1{m}_1}(\hat{\mathbf{n}})Y_{{\ell }_2{m}_2}(\hat{\mathbf{n}})Y_{{\ell }_3{m}_3}(\hat{\mathbf{n}})⟨{\kappa }_{{\ell }_1{m}_1}{\kappa }_{{\ell }_2{m}_2}{\kappa }_{{\ell }_3{m}_3}⟩$$

$$=h_{{\ell }_1{\ell }_2{\ell }_3}^{-2}⟨\int {\mathrm{d}}^2\hat{\mathbf{n}}{\kappa }_1(\hat{\mathbf{n}}){\kappa }_2(\hat{\mathbf{n}}){\kappa }_3(\hat{\mathbf{n}})⟩,$$

$$h_{{\ell }_1{\ell }_2{\ell }_3}^2$$

$$=\sqrt{\frac{(2{\ell }_1+1)(2{\ell }_2+1)(2{\ell }_3+1)}{4\pi }}\int {\mathrm{d}}^2\hat{\mathbf{n}}{Y}_{{\ell }_{1}0}(\hat{\mathbf{n}})Y_{{\ell }_{2}0}(\hat{\mathbf{n}})Y_{{\ell }_{3}0}(\hat{\mathbf{n}})$$

$$=\frac{1}{\sqrt{4\pi }}\int {\mathrm{d}}^2\hat{\mathbf{n}}{h}_1(\hat{\mathbf{n}})h_2(\hat{\mathbf{n}})h_3(\hat{\mathbf{n}}),$$

$$b_{b_1{b}_2{b}_3}$$

$$a(\hat{\mathbf{n}})\equiv g(\hat{\mathbf{n}})+[g^2(\hat{\mathbf{n}})-⟨g^2⟩],$$

$$a_{\ell m}$$

$$=g_{\ell m}+\int {\mathrm{d}}^2{\hat{\mathbf{n}}}^{\prime }\sum _{{\ell }^{\prime }{\ell }^{\prime \prime }{m}^{\prime }{m}^{\prime \prime }}{Y}_{\ell m}^{\ast }({\hat{\mathbf{n}}}^{\prime })Y_{{\ell }^{\prime }{m}^{\prime }}({\hat{\mathbf{n}}}^{\prime })Y_{{\ell }^{\prime \prime }{m}^{\prime \prime }}({\hat{\mathbf{n}}}^{\prime })g_{{\ell }^{\prime }{m}^{\prime }}{g}_{{\ell }^{\prime \prime }{m}^{\prime \prime }}$$

$$=g_{\ell m}+(-1{)}^m\sum _{{\ell }^{\prime }{\ell }^{\prime \prime }{m}^{\prime }{m}^{\prime \prime }}{G}_{-m,m^{\prime },m^{\prime \prime }}^{\ell {\ell }^{\prime }{\ell }^{\prime \prime }}{g}_{{\ell }^{\prime }{m}^{\prime }}{g}_{{\ell }^{\prime \prime }{m}^{\prime \prime }},$$

$$G_{m_1{m}_2{m}_3}^{{\ell }_1{\ell }_2{\ell }_3}$$

$$=h_{{\ell }_1{\ell }_2{\ell }_3}\left(\begin{array}{ccc}{\ell }_1& {\ell }_2& {\ell }_3\\ m_1& m_2& m_3\end{array}\right),$$

$$G_{m_1{m}_2{m}_3}^{{\ell }_1{\ell }_2{\ell }_3}=G_{m_2{m}_3{m}_1}^{{\ell }_2{\ell }_3{\ell }_1}=G_{m_3{m}_1{m}_2}^{{\ell }_3{\ell }_1{\ell }_2}.$$

$$\sum _{m_1{m}_2{m}_3}[G_{m_1{m}_2{m}_3}^{{\ell }_1{\ell }_2{\ell }_3}{]}^2=h_{{\ell }_1{\ell }_2{\ell }_3}^2.$$

$$b_{{\ell }_1{\ell }_2{\ell }_3}^{gg{g}^2}$$

$$=h_{{\ell }_1{\ell }_2{\ell }_3}^{-2}\sum _{m_i}{G}_{m_1{m}_2{m}_3}^{{\ell }_1{\ell }_2{\ell }_3}(-1{)}^{m_3}\sum _{{\ell }^{\prime }{\ell }^{\prime \prime }{m}^{\prime }{m}^{\prime \prime }}{G}_{-m_3,m^{\prime },m^{\prime \prime }}^{{\ell }_3{\ell }^{\prime }{\ell }^{\prime \prime }}⟨g_{{\ell }_1{m}_1}{g}_{{\ell }_2{m}_2}{g}_{{\ell }^{\prime }{m}^{\prime }}{g}_{{\ell }^{\prime \prime }{m}^{\prime \prime }}⟩+(\text{2 perms.}).$$

$$b_{{\ell }_1{\ell }_2{\ell }_3}^{gg{g}^2}$$

Full text

cosmology, cosmic microwave background, weak gravitational lensing, large-scale structure

††preprint: YITP-18-127

I Introduction

The path of cosmic microwave background (CMB) photons between the last scattering surface and an observer is deflected by the gravitational potential of the large-scale structure (LSS) of the Universe. This leads to distortions of the observed CMB anisotropies, which are correlated across scales. By extracting the lensing effect in these anisotropies, we can directly probe the underlying gravitational potential of the LSS, and constrain cosmological properties which affect it, such as characteristics of dark matter, dark energy, and massive neutrinos. By its nature, the lensing of CMB is sensitive to higher redshifts than most other cosmological observations of LSS, and it offers one of the most powerful probes of fundamental issues in cosmology and physics in the near future.

Recent CMB experiments have already detected the angular power spectrum of the lensing potential very precisely ACT16:phi ; BKVIII ; P18:phi ; PB14:phi ; SPT:phi . For instance, the Planck detection (P18:phi, ) has a 40 $\sigma$ significance. The detection and precise determination of the CMB lensing bi-spectrum is therefore an obvious and important next step in CMB scientific analyses. Recent studies indeed showed that the lensing potential bi-spectrum and other higher-order statistics due to non-linearity are detectable in near future CMB experiments Namikawa:2016b ; Pratten:2016 ; Liu:2016nfs ; Namikawa:2018b .

In previous works Namikawa:2016b ; Pratten:2016 , the signal-to-noise ratio of the bi-spectrum was obtained by using the fitting formulae of Refs. Scoccimarro:2001 ; Gil-Marin:2012 . However, these fitting formulae were tested in other contexts and it is unclear how they perform in this case. This paper therefore assesses the validity of these fitting formulae to calculate the CMB lensing bi-spectrum by comparing with a full-sky CMB lensing simulation. In addition to allowing refined expectations, accurate analytic predictions for the bi-spectrum are also necessary to reduce the computational cost when analyzing real data.

This paper is organized as follows. Section II reviews theoretical predictions of the bi-spectrum of CMB lensing and Section III describes the simulations used in our analysis. Section IV presents our main results. Section V discusses theoretical predictions in modified gravity theories. Section VI summarizes our work. Appendix A verifies our bi-spectrum measurements. Appendix B briefly summarizes the fitting formulas of the matter bi-spectrum. Appendix C explains the bi-spectrum in modified gravity theories. Appendix D further discusses the discrepancy between fitting formulae and simulations appearing in the squeezed lensing bi-spectrum.

Throughout this paper we assume the fiducial flat $\Lambda$CDM cosmology already used in Takahashi:2017 . The cosmological parameters are $\Omega_{\rm cdm}=0.233$, $\Omega_{\rm b}=0.046$, $h=0.7$, $\sigma_{8}=0.82$ and $n_{\rm s}=0.97$. We also use natural units where $c=G=1$.

II CMB lensing bi-spectrum

In this section, we review briefly the general formalism of the bi-spectrum and the theoretical prediction of the bi-spectrum of CMB lensing as shown in, e.g., Refs. Namikawa:2016b ; Pratten:2016 .

II.1 Bi-spectrum

Denoting the harmonic coefficients of the fluctuations (e.g., the CMB temperature anisotropies) as $a_{\ell m}$, the (angular) bi-spectrum is in general given by (e.g., Fergusson:2012 )

[TABLE]

It is convenient to use the reduced bi-spectrum,

[TABLE]

which removes the geometrical factor

[TABLE]

In the flat-sky approximation, the bi-spectrum is given by

[TABLE]

where $a_{\bm{\ell}_{i}}$ is a Fourier mode of fluctuations. The relationship between the full- and flat-sky bi-spectrum is given by Hu:2000ee ,

[TABLE]

This equation shows that the full-sky reduced bi-spectrum is equivalent to the flat-sky bi-spectrum. In the following, we use the flat-sky bi-spectrum to compute the reduced bi-spectrum and compare it with the results from the full sky simulation.

II.2 Bi-spectrum from CMB lensing

In a CMB lensing analysis, we measure the lensing potential which is defined as (see, e.g., Lewis:2006fu ; Hanson:2009kr )

[TABLE]

Here $\chi_{*}$ is the comoving distance to the CMB last-scattering surface and $\Psi$ is the Weyl potential. The lensing kernel, $W(\chi,\chi_{*})$, is defined (for a spatially flat cosmology) as

[TABLE]

where $\Theta(x)$ denotes the Heaviside step function. The bi-spectrum of the lensing potential is then given by Eq. (1). Hereafter, instead of using the lensing potential, we frequently use the lensing convergence, $\kappa_{\ell m}=\ell(\ell+1)\phi_{\ell m}/2$.

In general, a non-zero bi-spectrum is induced by non-Gaussian fluctuations, and in CMB lensing, it is sourced by two effects. One source is the non-Gaussian matter distribution induced by the nonlinear gravitational evolution Namikawa:2016b , and the other is the next-to-leading order correction to the Born approximation for computing the deflection of the path, which is referred to as the post-Born correction Pratten:2016 . We denote these contributions respectively by $B^{\rm LSS}$ and $B^{\rm pb}$, and present their explicit expressions below.

The bi-spectrum of the lensing convergence which arises from the nonlinear growth of the density perturbations is given in the flat-sky limit by Namikawa:2016b

[TABLE]

where $\bm{k}_{i}=\bm{\ell}_{i}/\chi$. Here, $B_{\delta}$ is the matter bi-spectrum which results from the nonlinear growth of structure. In the weakly nonlinear regime, it can be obtained by using perturbation theory. The result at the tree-level order is of the general form

[TABLE]

where $P_{\delta}(k,z)$ is the matter power spectrum at redshift $z(\chi)$, and the function $F_{2}$ is the second-order perturbation theory kernel (e.g., Bernardeau:2001qr ). Writing $\bm{k}_{1}\cdot\bm{k}_{2}=k_{1}k_{2}\cos\theta$, it is given by

[TABLE]

where $a(k,z)$, $b(k,z)$ and $c(k,z)$ are unity at the tree-level in perturbation theory. In order to capture the deviation from the tree-level prediction, one may calculate the correction at next-to-leading order (one-loop) to the prediction. We will briefly discuss its validity in subsequent analysis (see Sec. IV.3). In this paper, our primary focus is to compare the lensing simulations with analytical predictions based on fitting formulae of the matter bi-spectrum. These fits provide an analytical functional form for $a(k,z)$, $b(k,z)$ and $c(k,z)$. We will assess the most recent results for these fitting formulae, as given by Gil-Marin et al. in Ref. Gil-Marin:2012 , hereafter denoted by “GM”, and the earlier results of Scoccimaro and Couchman Scoccimarro:2001 , denoted below by “SC”.

Furthermore, the post-Born correction to the CMB lensing bi-spectrum, $B^{\rm pb}$ is expressed, as in Pratten:2016 , by

[TABLE]

where

[TABLE]

$P_{\Psi}(k,\chi)$ being the power spectrum of the Weyl potential at a comoving distance $\chi$.

Fig. 1 shows the theoretical predictions of the CMB lensing bi-spectrum given by Eqs. 8 and 11, assuming either the GM or SC fitting formulae, for the following four specific configurations:

1. C1.

equilateral, with $\ell_{1}=\ell_{2}=\ell_{3}=\ell$, 2. C2.

folded, with $\ell_{1}=2\ell_{2}=2\ell_{3}=\ell$, 3. C3.

squeezed, with $\ell_{1}=50$, $\ell_{2}=\ell_{3}=\ell$, and 4. C4.

isosceles, with $\ell_{1}=\ell$, $\ell_{2}=\ell_{3}=1000$.

We also show the contribution from the matter bi-spectrum alone. This figure clearly shows that the post-Born correction is important and changes the shape of the bi-spectrum in all of the above configurations.

The figure also shows the effect of multipole binning, which is discussed further in Sec. III.2. Binning is not necessary to constrain cosmological parameters, but in practice, it is rather essential, e.g., to reduce the computational cost. Here, we adopt binning in order to reduce the dispersion of the simulation measurements, as follows:

[TABLE]

where $\sum_{\ell_{i}}^{b_{i}}$ ($i=1,2,3$) means the summation over $\ell_{i}$ within the multipole bin, $b_{i}$. As shown in Fig. 1, the binned bi-spectrum (dot-dashed) can be very different from the unbinned bi-spectrum (solid). For example, the equilateral binned bi-spectrum, $b\equiv b_{1}=b_{2}=b_{3}$, contains not only the unbinned equilateral bi-spectrum but also all of the other bi-spectra which satisfies $\ell_{1},\ell_{2},\ell_{3}\in b$. The figures show that the binning effect is not significant in the equilateral case but crucial in other configurations.

III CMB lensing bi-spectrum in simulations

III.1 Full-sky ray-tracing simulations

In this subsection, we briefly describe the full-sky ray-tracing simulations which we used to construct the lensing maps. Details may be found in Ref. Takahashi:2017 . First, we performed cosmological $N$-body simulations of dark matter to produce an inhomogeneous mass distribution in the Universe, which properly describe the non-Gaussian character of the matter distribution. The simulations followed the gravitational evolution of $2048^{3}$ particles in cubic simulation boxes by using the public code Gadget2 Springel:2001 ; Springel:2005 . The cubic boxes have side lengths of $450,\,900,\,1350,\dots,6300\,h^{-1}{\rm Mpc}$ from low to high redshift $(z\leq 7.1)$. We chose an arbitrary point in the box as the observer’s position and constructed spherical lens shells with a thickness of $150h^{-1}\,{\rm Mpc}$ around the observer (who is located at the center of the shells). We construct three lens shells (having a $450h^{-1}\,{\rm Mpc}$ total thickness) in each box. The inner shell is taken from the smaller box at lower redshift. Then, we projected the $N$-body particles onto each lens shell, and calculated the surface mass density and deflection angle. For higher redshifts ($z=7.1-1100$), we prepared the lens shells based on the Gaussian fluctuations (instead of the $N$-body method), which is a good approximation in the linear regime. Light rays, emitted from the observer, are deflected at each lens shell, and these ray paths are numerically evaluated up to the last scattering surface. Therefore, the simulation includes the post-Born correction. We used the public code GRayTrix111http://th.nao.ac.jp/MEMBER/hamanatk/GRayTrix/ Hamana:2015 ; Shirasaki:2015 for this ray-tracing computation. The code relies on the Healpix scheme Gorski:2004by . We adopt the angular resolution of $N_{\rm side}=4096$ in our main results (but we also increase the resolution to $N_{\rm side}=8192$ in order to see the dependence of our results on the map resolution). We numerically obtain the CMB deflection angle $\bm{d}(\hat{\bm{n}})$ on the full sky maps, and then compute the harmonic coefficients of the lensing potential, $\phi_{\ell m}$, using the spin-1 harmonic transform. The harmonic coefficients of the lensing convergence are finally obtained by $\kappa_{\ell m}=\ell(\ell+1)\phi_{\ell m}/2$. We prepared $108$ such maps in total. We checked that the angular power spectrum of the lensing potential agrees with the theoretical prediction calculated by CAMB Lewis:1999bs within $5\%$ up to $\ell=3000$ (see Section 3.5 of Ref. Takahashi:2017 ).

In addition to the full-sky CMB lensing maps which include both the nonlinear growth and the post-Born effect, we also created CMB lensing maps only including the latter. Such maps are useful to study each contribution separately. To create the maps, we construct lens shells based on the Gaussian fluctuations without modifying the density power spectrum. Then, we performed the ray-tracing simulation, and repeated the same simulation with a different realization to obtain 10 maps.

III.2 Measurement of the lensing bi-spectrum

Next we describe how we measure the binned lensing reduced bi-spectrum in practice. The full-sky reduced (unbinned) bi-spectrum, $b_{\ell_{1}\ell_{2}\ell_{3}}$, is rewritten with

[TABLE]

where $\kappa_{i}(\hat{\bm{n}})$ is obtained by the inverse harmonic transform of $\delta_{\ell\ell_{i}}\kappa_{\ell m}$. From the first to the second line, we use the formula of the Gaunt integral shown in Eq. (19). The square of the geometrical factor is rewritten as

[TABLE]

where $h_{i}(\hat{\bm{n}})$ is given by the inverse harmonic transform of $\sqrt{2\ell_{i}+1}\delta_{\ell\ell_{i}}\delta_{m0}$.

To reduce the simulation error from the cosmic variance, we measure the binned reduced bi-spectrum of Eq. (13) from the simulation data. Substituting Eqs. 14 and 15 into Eq. (13), we find

[TABLE]

where $p$ is the pixel index, and $\sum_{p}$ is the sum over all pixels. The discretization of the sphere applied to the above equation is tested in Appendix A.

IV Results

In this section, we compare analytic predictions of the lensing bi-spectrum with measured results from lensing simulations, focusing mainly on the predictions based on fitting formula. We first present the results of the total lensing bi-spectrum in Sec. IV.1. Then, the contribution from post-Born correction is measured separately and compared with analytic predictions in Sec. IV.2. Based on this, possible reasons for the discrepancy seen in Sec. IV.1 are discussed. In Sec. IV.3, we examine the perturbation theory calculation, and the applicable range of higher-order (one-loop) predictions is discussed.

IV.1 Total contributions

Fig. 2 compares theory and simulations for the same four configurations presented in Fig 1 (equilateral, folded, squeezed and isosceles). The analytic prediction is computed based on Eqs. 8 and 11. Note that, in the simulation, we compute the bi-spectrum of the lensing convergence instead of the lensing potential. The top panel in each figure shows the bi-spectrum, the middle panel shows the fractional difference, $\Delta b/b\equiv(b^{\rm sim}-b^{\rm theory})/b^{\rm theory}$, and the bottom panel shows the difference against the simulation dispersion, $\Delta b/\sigma\equiv(b^{\rm sim}-b^{\rm theory})/\sigma$.

In the equilateral case, the analytic prediction with the SC fitting formula agrees with the simulation result, at least, up to $\ell=2048$ within the simulation error. The analytic result with the GM fitting formula slightly overestimates the simulation results. In the folded case, the analytic results with both the SC and GM fitting formulas overestimate the simulation at almost all scales. The most significant discrepancy is found in the squeezed configuration. The analytic predictions start to overestimate the simulation results as we go to small scale and the discrepancy is much larger than the simulation error. The fractional difference is roughly $30-40$% at small scales ($\ell\geq 1000$). In the isosceles case, the analytic prediction becomes consistent with the simulation if $\ell_{1}$ becomes close to $\ell_{2}=\ell_{3}=1000$.

In order to see the dependence on the size of the multipole bins, Fig. 3 shows the bi-spectrum with an increased number of multipole bins. If the number of multipole bins increases, the binned bi-spectrum is less contaminated by different configurations of the bi-spectrum while the simulation error increases. The consistency between the analytic and simulated bi-spectrum is improved compared to the case with $20$ multipole bins, mostly due to the increase of the error bars. However, a huge discrepancy in the squeezed configuration still remains at small scales, and its size is almost the same as that shown in Fig. 2 (see also Fig. 12 in Appendix D).

IV.2 Post-Born correction

As shown in Refs. Pratten:2016 ; Lewis:Pratten:2016 , the post-Born correction changes the shape of the bi-spectrum. Nevertheless, their discussions are based on the leading-order contributions to the bi-spectrum, and the validity and accuracy of their predictions have not yet been tested. Here, we compare the leading-order analytic prediction with the post-Born simulation. Fig. 4 shows the bi-spectrum measurements from the post-Born simulation to check the consistency between the analytic predictions and simulation results. The analytic prediction of the post-Born bi-spectrum is given by Eq. (11), with $P(k)$ computed using the HALOFIT prescription Smith:2002dz ; Takahashi:2012em . The leading-order predictions are indeed in good agreement with the simulation results, and no notable discrepancy is found over all scales, suggesting that the leading-order contribution is sufficient to describe post-Born correction to the lensing bi-spectrum.

As we discuss in Sec. IV.1, although the prediction by the analytic fitting formula is consistent with the simulation results in most configurations and scales, we find a significant discrepancy in the squeezed configuration (a similar result was found in Ref. Coulton:2018ebd using a flat-sky simulation). The good match of the post-Born prediction with simulation measurements seen in Fig. 4 suggests that the discrepancy mostly comes from the LSS contribution, as theoretically predicted by Eq. (8). Further, the close match of the post-Born prediction and simulation measurements also suggests the discrepancy observed is not a failing of the Limber approximation. Indeed, we find a discrepancy between the squeezed matter bi-spectra obtained from simulation and fitting formula especially at low $z$, i.e., the fitting formula is not accurate in the squeezed limit. This may be because of the configurations and scales the formula was fit to. This may also be partly because the fitting formula is derived using the unbinned calculation while the binning effect significantly changes the squeezed bi-spectrum prediction. Ref. Coulton:2018ebd found that the poor resolution of the simulation map leads to a discrepancy in the small-scale squeezed bi-spectrum, but our results are not significantly affected by the map resolution. While a thorough exploration of the discrepancy requires the validation of the fitting formula at all redshifts and scales and is beyond the scope of this paper, we further discuss possible sources of the discrepancy in Appendix D.

IV.3 Comparison with perturbation theory

The analytic predictions of the lensing bi-spectrum have up to this point been computed with fitting formula for the matter bi-spectrum. Here, we examine the range of validity of perturbation theory calculations, specifically we compare the prediction based on one-loop order calculations with our simulation measurements. Fig. 5 shows the results for the equilateral and squeezed configurations. Note that only the LSS contribution as given by Eq. (8) is shown, with the black points being computed by subtracting the post-Born correction from the simulation measurement. The analytic predictions for the equilateral configuration (left panel) are computed without the binning effect since the impact of the multipole binning is negligible in this case as shown in Fig. 1. We find the one-loop prediction agrees with the simulation up to $\ell\sim 600$ and does substantially better than linear theory although valuable non-linear information is still missing from this prediction. In the right panel, the predictions for the squeezed limit configuration are plotted, with the binned predictions depicted as dot-dashed lines. We find that all predictions fail to trace the simulation measurements for $300\leq\ell$ (see Appendix D for more details). While the binned predictions from the fitting formula and one-loop calculation are indistinguishable in the binned case at $\ell\leq 1000$, the unbinned results show a significant discrepancy at higher $\ell$ when unbinned.

Albeit its limited applicability range, one advantage of the perturbation theory approach is its flexibility. Within the framework of perturbative calculation, one can construct predictions for general theories of gravity and dark energy with relatively few assumptions thrown into the modelling. In the next section, based on the range of applicability deduced here, we use one-loop perturbation theory to investigate the signal of modified gravity on the CMB lensing bi-spectrum.

V Signal of modified gravity

In this section, we investigate the impact of modified gravity on the CMB lensing bi-spectrum. Beyond GR, currently available predictions of lensing bispectrum, especially for LSS contribution, are from the perturbation theory calculation through the prediction of matter bi-spectrum that has been recently extended by Ref. Bose:2018zpk , based on the general framework of perturbation theory treatment (Koyama_etal2009 ; Taruya2016 ; Bose:2016qun ). No fitting formula for matter bi-spectrum has been calibrated from modified gravity simulations. Nevertheless, a simple treatment proposed by Ref. Namikawa:2018b is to stick to the the fitting formula of matter bi-spectrum in Eq. (9) with the kernel function given by Eq. (10), but introducing some freedoms in the kernel function. With this modification, the fitting formula allows us to predict the lensing bi-spectrum in a rather wide class of modified theories of gravity, including Horndeski theoryHorndeski1974 and beyond-Horndeski class theory Gleyzes:2015a ; Gleyzes:2015b ; Langlois:2016a ; Langlois:2016b . Indeed, Ref. Bose:2018zpk has checked with modified gravity simulations that this treatment provide a reasonable description of matter bi-spectrum at large scales in several modified gravity models.

Here, we restrict our analysis to two representative models of modified gravity, both with distinct screening mechanisms which are essential theoretical implementations that recover GR at small scales. This allows them to be consistent with solar system tests of gravity. The representative models we choose are the Vainshtein-screened DGP model of gravity Dvali:2000hr and the chameleon-screened Hu-Sawicki $f(R)$ model of gravity Hu:2007nk . The DGP model also serves as a representative of the Horndeski class of models and so we may look at two distinct predictions, that given by perturbation theory and that given by the fitting formula described in Ref. Namikawa:2018b . On the other hand, the $f(R)$ case has no non-linear matter bi-spectrum fitting formula available in the literature, and so only a perturbative treatment will be considered. For the perturbative treatment we consider both the tree level and one-loop prediction (see Appendix C). The theory parameters we choose are $\Omega_{\rm rc}=0.438$222A common parametrization of the cross-over scale, see Bose:2016qun for example. and $|\bar{f}_{R0}|=2.5\times 10^{-6}$ which is a recent constraint coming from the LSS observations Smith:2009fn .

In Fig. 6, we show the ratio of the LSS contribution (to the CMB lensing bi-spectrum) in DGP (left) and $f(R)$ gravity (right) to the prediction of GR. We consider the equilateral configuration as it has been shown to offer the largest signal of modified gravity in the representative models chosen Bose:2018zpk as well as the theoretical predictions giving good consistency with simulation measurements (see Fig. 2). In both models, the linear power spectrum is identical to GR at large scales.333Note that in the DGP case we also have a modification to the linear growth factor which we normalize to that of GR. This factor is completely degenerate with $\sigma_{8}$ and so such scale independent modifications to the linear power spectrum are not very useful to test modified gravity unless one can break such a degeneracy using other probes. For the DGP case, the deviation shown here comes only from non-linear screening effects while for $f(R)$ gravity we have effects both from screening and scale-dependent modifications to the linear growth of structure. In both panels, we show the tree-level prediction (orange) and one-loop prediction (green).444The noise seen in the plots is numerical. For the DGP case we also show the prediction coming from the non-linear fitting formula of Ref. Namikawa:2018b (red).

Note that the effect of screening mechanism is perturbatively incorporated into the perturbation theory predictions through the modified higher-order kernels. Because of this, the screening effect is ineffective at tree-level order, and the predictions do not approach unity (i.e., GR) at small scales. On the other hand, the one-loop predictions depicted as green curves include the screening effect up to the fourth-order kernels in perturbation theory. These higher-order corrections give distinct scale-dependent shapes and dominate the signal at $\ell>100$. They also work to reduce the linear modification. This is not evident in the DGP case (left panel) as we have normalized the linear growth of DGP to that of GR. Rather, this is clearer in the $f(R)$ case (right). Additional orders in perturbation theory would work to complete the screening at small scales. On this point, the fitting formula shown for the DGP case deviates from the one-loop calculation noticeably at smaller scales. This can be expected as the effect of screening mechanism is not fully incorporated into the fitting formula.

Fig. 6 indicates that the impact of the modification to gravity we considered is of the order of a few percent at most. The CMB Stage-IV experiments would not be able to probe such small signals of modified gravity at the multipoles considered Namikawa:2018b . However, this is only the case if we restrict the scales to the applicable range of one-loop perturbation theory. At smaller scales, a rather notable deviation from GR is still allowed in beyond-Horndeski class theory Hirano:2018 . In any case, for a decisive detection of modified gravity signal, small-scale information will be key. In doing so, additional complication may arise from the post-Born correction (Eq.11), for which we must also rely on a non-linear prescription for the power spectrum in order to get an accurate prediction. For DGP model, a good approximation is to use the fitting formula in GR (i.e., HALOFIT) and simply replace the linear power spectrum in GR with that in DGP. For $f(R)$ gravity, no such a simple treatment can work well, but there exists an extension of HALOFIT formula Zhao:2013dza . Since such extensions are not readily available in general models, one may alternatively consider the prediction based on perturbation theory Taruya2016 ; Bose:2016qun , although it would not capture accurately the small-scale behaviors.

Another interesting point is that in the DGP case (left plot of Fig. 6) we find the fitting formula gives a significantly different prediction from the one-loop. This may hint at a lack of screening effect included in the fitting formula as only the second-order kernel is modified in this case. To investigate this further, one would require modified gravity simulation measurements and is beyond the scope of this work.

VI Summary

The CMB lensing bi-spectrum will be soon detected from upcoming CMB experiments, and will become an important observable, complementary to the lensing power spectrum. In this paper, making use of full-sky lensing simulations, we have tested the analytic predictions of the CMB lensing bi-spectrum based on existing fitting formulae of the dark matter bi-spectrum and on perturbation theory. We found that the agreement between the simulation and fitting formulae depends on the configuration and scales of the bi-spectrum. In the equilateral case, the bi-spectrum obtained from the simulation is in good agreement with that derived from the fitting formulae. On the other hand, the squeezed bi-spectrum from the simulation deviates significantly from that predicted by the fitting formulae.

We also found that a significant discrepancy also appears in the squeezed matter bi-spectrum at low $z$, indicating that the discrepancy in the squeezed lensing bi-spectrum probably comes mostly from the inconsistency between the analytic prediction and simulation results of the matter bi-spectrum.

We have also compared the measured equilateral and squeezed configurations with one-loop perturbation theory. For the equilateral case we find that one-loop perturbation theory gives a much better prediction than tree-level calculation, but the scales of applicability is still restricted to rather low multipoles ($\ell\leq 600$). A notable point may be that the one-loop prediction becomes comparable to the fitting formula for the squeezed configuration up to $\ell\leq 1000$, though both fail to model the measured bi-spectrum as previously mentioned.

Using one-loop perturbation theory, we have also discussed the effects of modified gravity on the lensing bi-spectrum. Considering DGP and $f(R)$ gravity models, we find that the signal of possible modification in the equilateral configuration is small at multipoles $\ell\lesssim 800$, with a maximum of $\sim 1\%$ deviation from GR in DGP model and $\sim 4\%$ in $f(R)$ gravity. Such a signal cannot be distinguished even with stage-IV surveys, however, a large deviation from GR is still allowed in beyond-Horndeski class theory. Accurately modeling bi-spectrum in general theories of gravity is thus crucial especially at small scales.

Further, we discussed the accuracy of the lensing bi-spectrum predictions while leaving aside possible practical issues in measuring the lensing bi-spectrum with real data. For example, with the quadratic estimator to reconstruct the lensing convergence, the measured bi-spectrum is a six-point correlation in CMB anisotropies. Similar to the lensing power spectrum measurement, an accurate subtraction of the disconnected six-point correlation is required. More practical issues on the lensing bi-spectrum measurement will be investigated elsewhere.

Another interesting direction would be to revise the fitting formula of the matter bi-spectrum as we see the discrepancy between the simulated and analytic matter bi-spectrum in the squeezed configuration. The current fitting formula of the matter bi-spectrum is derived within a limited range of scale, configuration, and cosmological parameters. Since the bi-spectrum contains additional information on cosmology, an accurate fitting formula of the matter bi-spectrum will be necessary for a fast computation of the bi-spectrum prediction in future cosmological implications. The updated fitting formula will be also useful for studies on not only the CMB lensing but also, e.g., the galaxy weak lensing and galaxy clustering. We leave these for our future work.

Acknowledgements.

TN acknowledges the support from the Ministry of Science and Technology (MOST), Taiwan, R.O.C. through the MOST research project grants (no. 107-2112-M-002-002-MY3). This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. BB acknowledges support from JSPS International Research Fellowship PE17043 and from the Swiss National Science Foundation (SNSF) Professorship grant No.170547. This work was supported in part by MEXT/JSPS KAKENHI Grant Number JP15H05899 (AT), JP16H03977 (AT), 15H05893 (RT) and 17H01131 (RT). Numerical computations were in part carried out on Cray XC50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan.

Appendix A Measurement of the bi-spectrum from full-sky simulation

Here we show the verification of our pipeline for measuring the bi-spectrum from the simulated CMB lensing convergence map. To test the algorithm of Eq. (16), we consider the following non-Gaussian fluctuations:

[TABLE]

where $g(\hat{\bm{n}})$ is a random Gaussian map generated by a power spectrum, $C_{\ell}$. We use the lensing convergence power spectrum to generate $g(\hat{\bm{n}})$. The harmonic coefficients at $\ell>0$ are obtained as

[TABLE]

with $g_{\ell m}$ being the harmonic coefficients of $g(\hat{\bm{n}})$. The Gaunt integral is defined as

[TABLE]

where $h_{\ell_{1}\ell_{2}\ell_{3}}$ is defined in Eq. (3). The Gaunt integral vanishes if $\ell_{1}+\ell_{2}+\ell_{3}$ is an odd integer. Taking the complex conjugate of the above equation, we obtain $G^{\ell_{1}\ell_{2}\ell_{3}}_{m_{1}m_{2}m_{3}}=(-1)^{m_{1}+m_{2}+m_{3}}G^{\ell_{1}\ell_{2}\ell_{3}}_{-m_{1},-m_{2},-m_{3}}$. By definition, the Gaunt integral satisfies the following symmetric property:

[TABLE]

Using the properties of the Wigner 3j symbol,

[TABLE]

The expected reduced bi-spectrum is decomposed into the contribution from the four and six point correlations. The contribution from the four point correlation is given by

[TABLE]

Since $g_{\ell m}$ is a random Gaussian fields, the four point correlation, $\langle g_{\ell_{1}m_{1}}g_{\ell_{2}m_{2}}g_{\ell^{\prime}m^{\prime}}g_{\ell^{\prime\prime}m^{\prime\prime}}\rangle$, is decomposed into the three terms by the Wick theorem. However, the term $\langle g_{\ell_{1}m_{1}}g_{\ell_{2}m_{2}}\rangle\langle g_{\ell^{\prime}m^{\prime}}g_{\ell^{\prime\prime}m^{\prime\prime}}\rangle$ provides the disconnected three-point correlation ($\langle a_{\ell_{1}m_{1}}a_{\ell_{2}m_{2}}\rangle\langle a_{\ell_{3}m_{3}}\rangle$), and should be ignored in the above equation. Then we obtain

[TABLE]

From the second to the third equation, we use Eq. (21).

The contribution from the six point correlation contains the product of the three angular power spectrum. Thus, we choose the amplitude of the angular power spectrum so that the terms involving the six point correlation is negligible.

Fig. 7 shows the validity of our code for measuring the reduced bi-spectrum from simulations. The binned analytic bi-spectrum is obtained from Eq. (16) where the reduced bi-spectrum at each multipole is given by Eq. (23). In measuring the bi-spectrum from simulations, a resolution parameter, $N_{\rm side}$, in obtaining $\kappa_{i}(\hat{\bm{n}})$ and $h_{i}(\hat{\bm{n}})$ in Eqs. (14) and (15) is chosen for each $i$ (multipole bin). We show the cases with $N_{\rm side}=\ell^{i}_{\rm max}$ or $3\ell^{i}_{\rm max}$ where $\ell^{i}_{\rm max}$ is the maximum multipole at the $i$th bin.

Appendix B Fitting formulas of the matter bi-spectrum

Here we briefly summarize the fitting formulas of the matter bi-spectrum used in this paper. Ref. Scoccimarro:2001 obtains the fitting formula of the matter bi-spectrum as

[TABLE]

with $Q(x)=(4-2^{x})/(1+2^{x+1})$. Here, $q=k/k_{\rm NL}$ with $k_{\rm NL}$ where $4\pi k^{3}_{\rm NL}P_{\rm m}^{\rm lin}(k_{\rm NL})=1$. The quantity $\sigma_{8}(z)$ is the variance of the matter density fluctuations smoothed with a top-hat sphere of radius $8h^{-1}$Mpc at redshift $z$. $n_{\rm eff}\equiv d\ln P_{\rm m}^{\rm lin}(k)/d\ln k$, is the effective spectral index of the linear power spectrum, $P_{\rm m}^{\rm lin}(k)$. We apply the smoothing to compute $n_{\rm eff}$ in order to remove un-physical oscillation as discussed in Ref. Gil-Marin:2012 . The parameters, $a_{i}$, are determined by fitting results of N-body simulations, which yields Scoccimarro:2001

[TABLE]

Later on, Ref. Gil-Marin:2012 proposed an improved fit given by

[TABLE]

Appendix C One-loop matter bi-spectrum

The one-loop matter bi-spectrum is defined as

[TABLE]

where \theta=\cos^{-1}{(\hat{\mbox{\boldmathk}}_{1}\cdot\hat{\mbox{\boldmathk}}_{2})}. The tree level ($B^{112}$) and one-loop terms (terms in square brackets) are defined in the usual way

[TABLE]

where we must add all permutations of the over-density perturbations, $\delta_{n}$, on the LHS, $n$ denoting the order of the perturbation. The angled brackets denote an ensemble average and under the assumptions of perturbation theory these averages can be decomposed into a product of linear power spectra convolved with perturbative kernels (see Bernardeau:2001qr for a review). These kernels are determined by solving energy and momentum conservation equations order by order. We direct the reader to Bose:2018zpk for a full description of this procedure for general models of gravity, including the nDGP and $f(R)$ models considered in this paper.

Appendix D Discrepancy in the squeezed lensing bi-spectrum

Here we discuss possible sources for the huge discrepancy which we found in the squeezed configuration between the theory based on fitting formulae and the direct simulation measurements

D.1 Role of non-linearity

The first check we perform is to see which configurations encode the most non-linear information. Those configurations with the most non-linearity will naturally be modelled the worst by the theoretical predictions. Fig.8 shows a contour of the ratio of the non-linear prediction to the tree level prediction for the LSS contribution to the CMB lensing bi-spectrum for various configurations. Specifically we show the GM (left) and 1-loop (right). We fix $\ell_{2}=500$ and vary $\ell_{1}$ and the angle between $\ell_{1}$ and $\ell_{2}$. The equilateral configuration has a good deal of non-linearity ($\ell_{1}=500,\cos(\theta)=-0.5$) which is expected as all modes are at their maximum multipole. We find in this case we have a $45\%$ enhancement and $33\%$ enhancement over the tree level prediction for GM and 1-loop respectively. In the squeezed limit we also find a lot of non-linear information ($\ell_{1}=50,\cos(\theta)=0.05$) with a $33\%$ and $30\%$ enhancement for GM and one-loop respectively. Further, in the very squeezed limit ($\ell_{1}\rightarrow 10$) we get enhancements up to $60\%$. When binning, as the simulation measurement is, these highly non-linear configurations will be included and so failure in the theoretical predictions is likely to become an issue. Finally, we also note that the GM formula systematically predicts significantly more non-linearity than the one-loop for these limits which further suggests inaccuracy of the fit when considering the bottom left plot of Fig.2.

Of course to perform a more accurate investigation one would ideally plot the same contour but with simulation measurements in the numerator of the ratio, but this would be very computationally expensive. Despite this, this does serve as a good indication as to which configurations are the most non-linear.

D.2 Matter bi-spectrum

Given the good match of the post-Born prediction with simulation measurements seen in Fig.4), it is strongly suggested that the discrepancy comes from the LSS contribution, theoretically predicted by Eq. (8). To investigate this issue we first have plotted the kernel of Eq. (8) as a function of redshift in Fig.9 for the equilateral (left) and squeezed (right) configurations. Each curve denotes a different multipole with the dashed lines being the tree level prediction and the solid the GM prediction. We find that in the squeezed limit the kernel with the GM fitting formula peaks in the range $z\sim 0.5$, and for multipoles $\ell\geq 600$ the peak is an order of magnitude larger than the equilateral shape. Given this we check the consistency of the squeezed matter bi-spectrum obtained from simulation and fitting formula at $z\sim 0.5$. Figs. 10 and 11 plot the matter bi-spectra at $z=0.541$ and $0.990$ measured from the N-body data used to construct the CMB maps. The box-sizes are $L=1800\,h^{-1}$Mpc (for $z=0.541$) and $2700\,h^{-1}$Mpc (for $z=0.990$) with $2048^{3}$ particles. The matter bi-spectrum at $z\sim 0.5$ gives the most contribution to the CMB bi-spectrum from Fig. 9. The black circles with bars are mean with errors from $6$ realizations. The GM fitting formula gives a good agreement for equilateral case but slightly over predicts for the squeezed case at small scales. This discrepancy is seen even for half $k$-bin width. This trend seems consistent with the CMB lensing bi-spectrum. Such a discrepancy can lead to the significant over-prediction noted at small scales in Fig. 2 once we integrate over all redshifts. The discrepancy noted at large scales in Fig. 2 could be a combination of effects including post-Born effects.

As mentioned in the main text this discrepancy may be ascribed to an inaccuracy in the GM fitting formula. Indeed, the formula is not fit to very squeezed configurations and was only fit in the range of $0.03\,h$ Mpc*-1*$\leq k\leq 0.4\,h$ Mpc*-1* Gil-Marin:2012 , corresponding to $46\leq\ell\leq 570$ at $z=0.541$. Further, the squeezed limit considered here for $100\leq\ell$ has $0.9\leq\theta/\pi$, $\theta$ being the angle between the ‘non-squeezed’ wave vectors, while the formula is only fit for $\theta/\pi\leq 0.9$. Despite this possibility, we perform a number of further tests to clarify the issue.

D.3 Resolution of simulation

As discussed in Ref. Coulton:2018ebd , the resolution of the simulation map also changes the small scale behavior of the squeezed lensing bi-spectrum. In the main text, we show the measurements of the bi-spectrum from the simulation of $N_{\rm side}=4096$. Here we check how the results depend on the simulation map resolution, $N_{\rm side}$.

Fig. 12 shows the fractional difference of the bi-spectrum at each configuration, $\Delta b/b$, where $\Delta b/b\equiv b^{N_{\rm side}=8192}/b^{N_{\rm side}=4096}-1$ is the fractional difference between the bi-spectra measured from the simulation of $N_{\rm side}=8192$ and $N_{\rm side}=4096$. In the squeezed case, the increase of the resolution enhances the amplitude of the bi-spectrum by $5-10\%$ at smaller scales. This increase is, however, too small compared to the discrepancy we found.

D.4 Other possibilities

Fig. 13 shows the dependence of the discrepancy in the squeezed bi-spectrum on the multipole bin for $\ell_{1}$. As the multipole range of $\ell_{1}$ goes to smaller scales, the discrepancy between the analytic and simulated bi-spectrum is reduced. However, the analytic prediction always overestimates the simulation. Indeed, the discrepancy is still large when we choose a larger multipole bin to $\ell_{1}$ which is much less effected by sample variance or possible large scale systematics.

Fig. 14 plots the fractional difference of the bi-spectrum with varying the number of multipole bins ($20,40$ and $60$). The discrepancy is almost independent from the number of multipole bins.