CMB lensing bi-spectrum: assessing analytical predictions against full-sky lensing simulations
Toshiya Namikawa Benjamin Bose François R. Bouchet Ryuichi Takahashi Atsushi Taruya








Leung Center for Cosmology and Particle Astrophysics, National Taiwan University, Taipei, 10617, Taiwan


Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan
Departement de Physique Theorique, Universite de Geneve, 24 quai Ernest Ansermet, 1211 Geneve 4, Switzerland


Institut d'Astrophysique de Paris, UMR 7095, Sorbonne Université & CNRS, 75014 Paris, France


Faculty of Science and Technology, Hirosaki University, 3 Bunkyo-cho, Hirosaki, Aomori 036-8561, Japan


Center for Gravitational Physics, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan
Kavli Institute for the Physics and Mathematics of the Universe, The University of Tokyo Institutes for Advanced Study (UTIAS), The University of Tokyo, Chiba 277-8583, Japan





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 ł∼ 200, but the one-loop order allows capturing the simulation results up to ł∼ 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 ł∼2000, but the difference in the squeezed limit raises to more than a factor of two at ł∼2000. 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. 



YITP-18-127

CMB lensing  bi-spectrum: assessing analytical predictions against full-sky lensing simulations
    Atsushi Taruya
    December 30, 2023
===============================================================================================








§ 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 <cit.>. For instance, the Planck detection <cit.> has a 40 σ 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 <cit.>. 


In previous works <cit.>, the signal-to-noise ratio of the bi-spectrum was obtained by using the fitting formulae of Refs. <cit.>. 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 <ref> reviews theoretical predictions of the bi-spectrum of CMB lensing and  
Section <ref> describes the simulations used in our analysis. 
Section <ref> presents our main results. 
Section <ref> discusses theoretical predictions in modified gravity theories. 
Section <ref> summarizes our work. 
Appendix <ref> verifies our bi-spectrum measurements. 
Appendix <ref> briefly summarizes the fitting formulas of the matter bi-spectrum.  
Appendix <ref> explains the bi-spectrum in modified gravity theories. 
Appendix <ref> further discusses the discrepancy between fitting formulae and simulations appearing in the squeezed lensing bi-spectrum. 


Throughout this paper we assume the fiducial flat ΛCDM cosmology already used in <cit.>. The cosmological parameters are Ω_ cdm = 0.233, Ω_ b = 0.046, h = 0.7, σ_8 = 0.82 and n_ s = 0.97. We also use natural units where c = G = 1.




§ 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. <cit.>. 




 §.§ Bi-spectrum
 


Denoting the harmonic coefficients of the fluctuations (e.g., the CMB temperature anisotropies) as a_ł m, the (angular) bi-spectrum is in general given by (e.g., <cit.>)


	B_ł_1ł_2ł_3 ≡∑_m_1m_2m_3ł_1ł_2ł_3m_1m_2m_3
    	a_ł_1m_1a_ł_2m_2a_ł_3m_3 
	 . 


It is convenient to use the reduced bi-spectrum,  


	b_ł_1ł_2ł_3 ≡h^-1_ł_1ł_2ł_3B_ł_1ł_2ł_3   (ł_1+ł_2+ł_3  is even)  ,


which removes the geometrical factor


	h_ł_1ł_2ł_3 ≡√((2ł_1+1)(2ł_2+1)(2ł_3+1)/4π)
    	ł_1ł_2ł_3000 . 
    


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


	a__1a__2a__3 
		= (2π)^2δ(_1+_2+_3) B__1_2_3
	 ,


where a__i is a Fourier mode of fluctuations. The relationship between the full- and flat-sky bi-spectrum is given by <cit.>,


	B_ł_1ł_2ł_3    ≃h_ł_1ł_2ł_3 B__1_2_3
	 .


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.





 §.§ Bi-spectrum from CMB lensing
 


In a CMB lensing analysis, we measure the lensing potential which is defined as (see, e.g., <cit.>)


	ϕ() = -2χ0χ_* W(χ,χ_*) Ψ(χ,)  .


Here χ_* is the comoving distance to the CMB last-scattering surface and Ψ is the Weyl potential. The lensing kernel, W(χ,χ_*), is defined (for a spatially flat cosmology) as


	W(χ,χ_*) = χ_*-χ/χχ_*Θ(χ_*-χ)  ,


where Θ(x) denotes the Heaviside step function. The bi-spectrum of the lensing potential is then given by Eq:bispec. Hereafter, instead of using the lensing potential, we frequently use the lensing convergence, κ_ł m=ł(ł+1)ϕ_ł 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 <cit.>, 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 <cit.>. We denote these contributions respectively by B^ LSS and B^ 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 <cit.>


	B^LSS__1_2_3 = χ0χ_*    
		[3Ω_m,0H_0^2/2a(χ)]^3
		

		×    χ^2 W^3(χ,χ_*)  B_δ(_1,_2,_3,χ)
	 , 


where _i=_i/χ. Here, B_δ 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


	B_δ(_1,_2,_3,χ) 
		   = 2F_2(_1,_2,z) P_δ(k_1,z)P_δ(k_2,z) 
	

		   + 2 perms.
	  ,
	


where P_δ(k,z) is the matter power spectrum at redshift z(χ), and the function F_2 is the second-order perturbation theory kernel (e.g., <cit.>). Writing _1·_2=k_1k_2cosθ, it is given by


	F_2(_1,_2,z)    = 5/7 a(k_1,z)a(k_2,z) 
		

		   + 1/2 k_1^2+k_2^2/k_1k_2 b(k_1,z)b(k_2,z)cosθ

		   + 2/7 c(k_1,z)c(k_2,z)cos^2θ , 


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. <ref>). 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. <cit.>, hereafter denoted by “GM", and the earlier results of Scoccimaro and Couchman <cit.>, denoted below by “SC". 

Furthermore, the post-Born correction to the CMB lensing bi-spectrum, B^ pb is expressed, as in <cit.>, by


	B^pb__1_2_3 
		   = 2_1·_2/ł^2_1ł^2_2[_1·_3  M_ł_1ł_2+_2·_3 M_ł_2ł_1] 
	

		   + cyc. perm. 
	 , 


where 


	M_łł'    = (łł')^4χ0χ_* [W(χ,χ_*)]^2/χ^2P_Ψ(ł/χ,χ) 
		

		   ×χ'0χW(χ',χ)W(χ',χ_*)/(χ')^2
		P_Ψ(ł'/χ',χ')
	 , 


P_Ψ(k,χ) being the power spectrum of the Weyl potential at a comoving distance χ. 





Fig. <ref> shows the theoretical predictions of the CMB lensing bi-spectrum given by Eq:bisp:lss,Eq:bisp:pb, assuming either the GM or SC fitting formulae, for the following four specific configurations: 
    
  C1. equilateral, with ł_1=ł_2=ł_3=ł, 
    
  C2. folded, with ł_1=2ł_2=2ł_3=ł, 
    
  C3. squeezed, with ł_1=50, ł_2=ł_3=ł, and
    
  C4. isosceles, with ł_1=ł, ł_2=ł_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. <ref>. 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: 


	b_b_1b_2b_3 ≡∑_ł_i^b_ih^2_ł_1ł_2ł_3b_ł_1ł_2ł_3/∑_ł_i^b_ih^2_ł_1ł_2ł_3 
     , 


where ∑_ł_i^b_i (i=1,2,3) means the summation over ł_i within the multipole bin, b_i. As shown in Fig. <ref>, the binned bi-spectrum (dot-dashed) can be very different from the unbinned bi-spectrum (solid). For example, the equilateral binned bi-spectrum, b≡ b_1=b_2=b_3, contains not only the unbinned equilateral bi-spectrum but also all of the other bi-spectra which satisfies ł_1,ł_2,ł_3∈ b. The figures show that the binning effect is not significant in the equilateral case but crucial in other configurations. 





§ CMB LENSING BI-SPECTRUM IN SIMULATIONS
 





 §.§ 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. <cit.>. 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 <cit.>. The cubic boxes have side lengths of 450,   900,   1350,…, 6300   h^-1 Mpc from low to high redshift (z ≤ 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  Mpc around the observer (who is located at the center of the shells). 
We construct three lens shells (having a 450h^-1  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 GRayTrix[<http://th.nao.ac.jp/MEMBER/hamanatk/GRayTrix/>] <cit.> for this ray-tracing computation. The code relies on the Healpix scheme <cit.>. We adopt the angular resolution of N_ side=4096 in our main results (but we also increase the resolution to N_ side=8192 in order to see the dependence of our results on the map resolution). We numerically obtain the CMB deflection angle d() on the full sky maps, and then compute the harmonic coefficients of the lensing potential, _ł m, using the spin-1 harmonic transform. The harmonic coefficients of the lensing convergence are finally obtained by κ_ł m=ł(ł+1)_ł 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 <cit.> within 5% up to ł=3000 (see Section 3.5 of Ref. <cit.>). 

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.




 §.§ 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_ℓ_1ℓ_2ℓ_3, is rewritten with 



	b_ł_1ł_2ł_3    = h^-1_ł_1ł_2ł_3 ∑_m_1m_2m_3 
		ł_1ł_2ł_3m_1m_2m_3 κ_ł_1m_1κ_ł_2m_2κ_ł_3m_3 
        

 	      = h^-1_ł_1ł_2ł_3 ∑_m_1m_2m_3
       	h^-1_ł_1ł_2ł_32 Y_ł_1m_1()Y_ł_2m_2()Y_ł_3m_3()
        κ_ł_1m_1κ_ł_2m_2κ_ł_3m_3
        

 	   
       
       
       
    
    
 	      = h^-2_ł_1ł_2ł_3 2 κ_1()κ_2()κ_3()
     , 



where κ_i() is obtained by the inverse harmonic transform of δ_łł_iκ_ł m. From the first to the second line, we use the formula of the Gaunt integral shown in Eq:gaunt. The square of the geometrical factor is rewritten as 



	h^2_ł_1ł_2ł_3    = (2ł_1+1)(2ł_2+1)(2ł_3+1)/4πł_1ł_2ł_3000^2
    	

       = √((2ł_1+1)(2ł_2+1)(2ł_3+1)/4π)
        
    	2 Y_ł_10()Y_ł_20()Y_ł_30()
    	

       = 1/√(4π) 2 h_1()h_2()h_3()
     , 



where h_i() is given by the inverse harmonic transform of √(2ł_i+1)δ_łł_iδ_m0. 

To reduce the simulation error from the cosmic variance, we measure the binned reduced bi-spectrum of Eq:bispec:binned:def from the simulation data. Substituting Eq:bispec:sim,Eq:bispec:norm into Eq:bispec:binned:def, we find 


	b_b_1b_2b_3    ≃√(4π)∑_p κ_1(_p)κ_2(_p)κ_3(_p)/∑_p h_1(_p)h_2(_p)h_3(_p)
     , 


where p is the pixel index, and ∑_p is the sum over all pixels. The discretization of the sphere applied to the above equation is tested in Appendix <ref>. 





§ 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. <ref>. Then, the contribution from post-Born correction is measured separately and compared with analytic predictions in Sec. <ref>. Based on this, possible reasons for the discrepancy seen in Sec. <ref> are discussed. In Sec. <ref>, we examine the perturbation theory calculation, and the applicable range of higher-order (one-loop) predictions is discussed. 




 §.§ Total contributions
 


Fig. <ref> compares theory and simulations for the same four configurations presented in Fig <ref> (equilateral, folded, squeezed and isosceles). The analytic prediction is computed based on Eq:bisp:lss,Eq:bisp:pb. 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, Δ b/b≡ (b^ sim-b^ theory)/b^ theory, and the bottom panel shows the difference against the simulation dispersion, Δ b/σ≡ (b^ sim-b^ theory)/σ. 

In the equilateral case, the analytic prediction with the SC fitting formula agrees with the simulation result, at least, up to ł=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 (ł≥ 1000). In the isosceles case, the analytic prediction becomes consistent with the simulation if ł_1 becomes close to ł_2=ł_3=1000. 





In order to see the dependence on the size of the multipole bins, Fig. <ref> 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. <ref> (see also Fig. <ref> in Appendix <ref>). 








 §.§ Post-Born correction



As shown in Refs. <cit.>, 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. <ref> 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:bisp:pb, with P(k) computed using the HALOFIT prescription <cit.>. 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. <ref>, 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. <cit.> using a flat-sky simulation). 


The good match of the post-Born prediction with simulation measurements seen in Fig. <ref> suggests that the discrepancy mostly comes from the LSS contribution, as theoretically predicted by Eq:bisp:lss. 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. <cit.> 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 <ref>. 








 §.§ 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. <ref> shows the results for the equilateral and squeezed configurations. Note that only the LSS contribution as given by Eq:bisp:lss 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. <ref>. We find the one-loop prediction agrees with the simulation up to ł∼ 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 ≤ℓ (see Appendix <ref> for more details). While the binned predictions from  the fitting formula and one-loop calculation are indistinguishable in the binned case at ℓ≤ 1000, the  unbinned results show a significant discrepancy at higher ℓ 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.









§ 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. <cit.>, based on the general framework of perturbation theory treatment (<cit.>). No fitting formula for matter bi-spectrum has been calibrated from modified gravity simulations. Nevertheless, a simple treatment proposed by Ref. <cit.> is to stick to the the fitting formula of matter bi-spectrum in Eq. (<ref>) with the kernel function given by Eq. (<ref>), 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 theory<cit.> and beyond-Horndeski class theory <cit.>. Indeed, Ref. <cit.> 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 <cit.> and the chameleon-screened Hu-Sawicki f(R) model of gravity <cit.>. 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. <cit.>. 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 <ref>). The theory parameters we choose are Ω_ rc=0.438[A common parametrization of the cross-over scale, see <cit.> for example.] and |f̅_R0|=2.5×10^-6 which is a recent constraint coming from the LSS observations <cit.>. 

In Fig. <ref>, 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 <cit.> as well as the theoretical predictions giving good consistency with simulation measurements (see Fig. <ref>). In both models, the linear power spectrum is identical to GR at large scales.[Note 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 σ_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).[The 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. <cit.> (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 ł>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. <ref> 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 <cit.>. 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 <cit.>. 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.<ref>), 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 <cit.>. Since such extensions are not readily available in general models, one may alternatively consider  the prediction based on perturbation theory  <cit.>, although it would not capture accurately the small-scale behaviors.

Another interesting point is that in the DGP case (left plot of Fig. <ref>) 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. 
 



§ 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 (ℓ≤ 600). A notable point may be that the one-loop prediction becomes comparable to the fitting formula for the squeezed configuration up to ℓ≤ 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 ℓ≲800, with a maximum of ∼ 1% deviation from GR in DGP model and ∼ 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. 





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.








§ 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:bispec:binned, we consider the following non-Gaussian fluctuations: 


	a() ≡g() + [g^2()-g^2]  ,


where g() is a random Gaussian map generated by a power spectrum, C_ł. 
We use the lensing convergence power spectrum to generate g(). 
The harmonic coefficients at ł>0 are obtained as 


	a_łm    = g_łm + 2' Y^*_łm(') g^2(') 
    	

    	   = g_łm + 2'∑_ł'ł”m'm”
    		Y^*_łm(')Y_ł'm'(')Y_ł”m”(')g_ł'm'g_ł”m”
        

    	   = g_łm + (-1)^m∑_ł'ł”m'm”G^łł'ł”_-m,m',m”g_ł'm'g_ł”m”
	 ,


with g_ł m being the harmonic coefficients of g(). The Gaunt integral is defined as


	G^ł_1ł_2ł_3_m_1m_2m_3    ≡2 Y_ł_1m_1()Y_ł_2m_2()Y_ł_3m_3()
		

		   = h_ł_1ł_2ł_3ł_1ł_2ł_3m_1m_2m_3 
     , 


where h_ł_1ł_2ł_3 is defined in Eq:hlll. The Gaunt integral vanishes if ł_1+ł_2+ł_3 is an odd integer. Taking the complex conjugate of the above equation, we obtain G^ł_1ł_2ł_3_m_1m_2m_3=(-1)^m_1+m_2+m_3G^ł_1ł_2ł_3_-m_1,-m_2,-m_3. By definition, the Gaunt integral satisfies the following symmetric property:


	G^ł_1ł_2ł_3_m_1m_2m_3 = G^ł_2ł_3ł_1_m_2m_3m_1 = G^ł_3ł_1ł_2_m_3m_1m_2
     . 


Using the properties of the Wigner 3j symbol, 


	∑_m_1m_2m_3[G^ł_1ł_2ł_3_m_1m_2m_3]^2 = h^2_ł_1ł_2ł_3  . 
    



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 


	b^ggg^2_ł_1ł_2ł_3
 	      = h^-2_ł_1ł_2ł_3 ∑_m_iG^ł_1ł_2ł_3_m_1m_2m_3
        a_ł_1m_1a_ł_2m_2a_ł_3m_3
        

 	      = h^-2_ł_1ł_2ł_3 ∑_m_iG^ł_1ł_2ł_3_m_1m_2m_3
        (-1)^m_3∑_ł'ł”m'm”G^ł_3ł'ł”_-m_3,m',m”
        g_ł_1m_1g_ł_2m_2g_ł'm'g_ł”m”
        + (2 perms.) 
     . 


Since g_ł m is a random Gaussian fields, the four point correlation, g_ł_1m_1g_ł_2m_2g_ł'm'g_ł”m”, is decomposed into the three terms by the Wick theorem. However, the term g_ł_1m_1g_ł_2m_2g_ł'm'g_ł”m” provides the disconnected three-point correlation (a_ł_1m_1a_ł_2m_2a_ł_3m_3), and should be ignored in the above equation. Then we obtain


	b^ggg^2_ł_1ł_2ł_3
 	      = h^-2_ł_1ł_2ł_3 ∑_m_iG^ł_1ł_2ł_3_m_1m_2m_3
        (-1)^m_3∑_ł'ł”m'm”G^ł_3ł'ł”_-m_3,m',m” 
       

            ×(-1)^m_1+m_2C_ł_1C_ł_2
       [δ_ł_1ł'δ_ł_2ł”δ_m_1,-m'δ_m_2,-m”
    	+ δ_ł_1ł'δ_ł_2ł”δ_m_1,-m'δ_m_2,-m”]
        + (2 perms.) 
       

 	      = 2C_ł_1C_ł_2h^-2_ł_1ł_2ł_3 ∑_m_iG^ł_1ł_2ł_3_m_1m_2m_3
        (-1)^m_1+m_2+m_3G^ł_3ł_1ł_2_-m_3,-m_1,-m_2
        + (2 perms.) 
       

 	      = 2[C_ł_1C_ł_2+(2 perms.)]h^-2_ł_1ł_2ł_3
     . 


From the second to the third equation, we use Eq:gaunt:sum. 

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. <ref> shows the validity of our code for measuring the reduced bi-spectrum from simulations. The binned analytic bi-spectrum is obtained from Eq:bispec:binned where the reduced bi-spectrum at each multipole is given by Eq:bispec:test. In measuring the bi-spectrum from simulations, a resolution parameter, N_ side, in obtaining κ_i() and h_i() in Eqs. (<ref>) and (<ref>) is chosen for each i (multipole bin). We show the cases with N_ side=ł^i_ max or 3ł^i_ max where ł^i_ max is the maximum multipole at the ith bin. 




§ FITTING FORMULAS OF THE MATTER BI-SPECTRUM
 


Here we briefly summarize the fitting formulas of the matter bi-spectrum used in this paper. Ref. <cit.> obtains the fitting formula of the matter bi-spectrum as 


	a(k,z)    = 1+{σ_8(z)}^a_6√(0.7Q(n_eff))(qa_1)^n_eff+a_2/1+(qa_1)^n_eff+a_2 

	b(k,z)    = 1+0.2a_3(n_eff+3)(qa_7)^n_eff+3+a_8/1+(qa_7)^n_eff+3.5+a_8 

	c(k,z)    = 1+[4.5a_4/(1.5+(n_eff+3)^4)](qa_5)^n_eff+3+a_9/1+(qa_5)^n_eff+3.5+a_9
	 , 


with Q(x)=(4-2^x)/(1+2^x+1). Here, q=k/k_ NL with k_ NL where 4π k^3_ NLP_ m^ lin(k_ NL)=1. The quantity σ_8(z) is the variance of the matter density fluctuations smoothed with a top-hat sphere of radius 8h^-1Mpc at redshift z. n_ eff≡ dln P_ m^ lin(k)/dln k, is the effective spectral index of the linear power spectrum, P_ m^ lin(k). We apply the smoothing to compute n_ eff in order to remove un-physical oscillation as discussed in Ref. <cit.>. The parameters, a_i, are determined by fitting results of N-body simulations, which yields <cit.>


	a_1    = 0.250     a_2    = 3.50     a_3    = 2.00     a_4    = 1.00     a_5    =2.00     a_6    = -0.200    
    a_7    = 1.00     a_8    = 0.00     a_9    = 0.00  .
    


Later on, Ref. <cit.> proposed an improved fit given by  


	a_1    = 0.484     a_2    = 3.74     a_3    = -0.849     a_4    = 0.392     a_5    = 1.01     a_6    = -0.575    
    a_7    = 0.128     a_8    = -0.722     a_9    = -0.926  .
    







§ ONE-LOOP MATTER BI-SPECTRUM
 


The one-loop matter bi-spectrum is defined as 

    B^1- loop(k_1,k_2,θ;a) =      B^112(k_1,k_2,θ;a) 
        + [B^222(k_1,k_2,θ;a)  +  B^321(k_1,k_2,θ;a) + B^114(k_1,k_2,θ;a)] ,

where θ = cos^-1(_1·_2). The tree level (B^112) and one-loop terms (terms in square brackets) are defined in the usual way 

    ⟨δ_n_1(_1) δ_n_2(_2) δ_n_3(_3) ⟩   =
    (2π)^3δ_ D(_1+_2 + _3 ) B^n_1n_2n_3(_1,_2,_3),

where we must add all permutations of the over-density perturbations, δ_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 <cit.> for a review). These kernels are determined by solving energy and momentum conservation equations order by order. We direct the reader to <cit.> for a full description of this procedure for general models of gravity, including the nDGP and f(R) models considered in this paper.




§ 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 





 §.§ 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.<ref> 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 ℓ_2 = 500 and vary ℓ_1 and the angle between ℓ_1 and ℓ_2. The equilateral configuration has a good deal of non-linearity (ℓ_1 = 500, cos(θ) = -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 (ℓ_1=50, cos(θ)=0.05) with a 33% and 30% enhancement for GM and 
one-loop respectively. Further, in the very squeezed limit (ℓ_1 → 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.<ref>. 

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. 








 §.§ Matter bi-spectrum















Given the good match of the post-Born prediction with simulation measurements seen in Fig.<ref>), it is strongly suggested that the discrepancy comes from the LSS contribution, theoretically predicted by Eq:bisp:lss. To investigate this issue we first have plotted the kernel of Eq:bisp:lss as a function of redshift in Fig.<ref> 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∼ 0.5, and for multipoles ł≥ 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∼ 0.5. Figs. <ref> and <ref>  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^-1Mpc (for z=0.541) and 2700 h^-1Mpc (for z=0.990) with 2048^3 particles. The matter bi-spectrum at z∼0.5 gives the most contribution to the CMB bi-spectrum from Fig. <ref>. 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. <ref> once we integrate over all redshifts. The discrepancy noted at large scales in Fig. <ref> 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≤ k ≤ 0.4 h Mpc^-1 <cit.>, corresponding to 46 ≤ℓ≤ 570 at z=0.541. Further, the squeezed limit considered here for 100 ≤ℓ has 0.9 ≤θ /π, θ being the angle between the `non-squeezed' wave vectors, while the formula is only fit for θ /π≤ 0.9. Despite this possibility, we perform a number of further tests to clarify the issue. 




 §.§ Resolution of simulation



As discussed in Ref. <cit.>, 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_ side=4096. Here we check how the results depend on the simulation map resolution, N_ side. 





Fig. <ref> shows the fractional difference of the bi-spectrum at each configuration, Δ b/b, where Δ b/b≡ b^N_ side=8192/b^N_ side=4096 - 1 is the fractional difference between the bi-spectra measured from the simulation of N_ side=8192 and N_ 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.  





 §.§ Other possibilities







Fig. <ref> shows the dependence of the discrepancy in the squeezed bi-spectrum on the multipole bin for ł_1. As the multipole range of ł_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 ł_1 which is much less effected by sample variance or possible large scale systematics. 





Fig. <ref> 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. 











mybst


