# The reconstructed CMB lensing bispectrum

Alba Kalaja  Giorgio Orlando  Aleksandr Bowkis  Anthony Challinor  P. Daniel
Meerburg  Toshiya Namikawa

###### Abstract

Weak gravitational lensing by the intervening large-scale structure (LSS) of
the Universe is the leading non-linear effect on the anisotropies of the
cosmic microwave background (CMB). The integrated line-of-sight mass that
causes the distortion – known as lensing convergence – can be reconstructed
from the lensed temperature and polarization anisotropies via estimators
quadratic in the CMB modes, and its power spectrum has been measured from
multiple CMB experiments. Sourced by the non-linear evolution of structure,
the bispectrum of the lensing convergence provides additional information on
late-time cosmological evolution complementary to the power spectrum. However,
when trying to estimate the summary statistics of the reconstructed lensing
convergence, a number of noise-biases are introduced, as previous studies have
shown for the power spectrum. Here, we explore for the first time the noise-
biases in measuring the bispectrum of the reconstructed lensing convergence.
We compute the leading noise-biases in the flat-sky limit and compare our
analytical results against simulations, finding excellent agreement. Our
results are critical for future attempts to reconstruct the lensing
convergence bispectrum with real CMB data.

##  1 Introduction

The photons of the cosmic microwave background (CMB) are deflected by the
large-scale structure (LSS) of the universe as they travel from the last
scattering surface ($z\sim 1100$) towards our detectors. This effect, called
weak gravitational lensing of the CMB [1], traces the distribution of matter
over a wide range of redshifts ($0.1<z<5$), making it a powerful probe of
late-time cosmological evolution.

Lensing alters the CMB power spectra and introduces non-Gaussianities in the
distribution of the anisotropies. It also generates a $B$-mode polarization
signal that represents the main source of confusion for the signal from
primordial gravitational waves generated during inflation (see, e.g., [2] and
references therein). From this perspective, lensing can be a nuisance that
needs to be removed (by delensing), e.g., to extract the primordial $B$-mode
signal [3, 4]. On the other hand, the induced statistics of lensed CMB modes
can be used to reconstruct the line-of-sight integrated gravitational
potential, also known as the lensing potential ($\phi$), whose gradient
determines the lensing deflections and whose second derivatives determine the
lensing shear and convergence (denoted with $\kappa$).

The first measurement of the CMB lensing signal was achieved in cross-
correlation with radio galaxy counts [5], whereas the detection with CMB data
only was done for the first time by ACT [6, 7], followed by Planck [8],
POLARBEAR [9], the South Pole Telescope (SPT) [10, 11] and BICEP/Keck [12].
The most recent measurement by Planck [13] detected the lensing power spectrum
at $40\,\sigma$ significance.

Given that CMB lensing is most sensitive to structure growth at redshift
$z\sim 2$, its power spectrum is particularly useful to constrain the content
of the Universe around that redshift, such as the dark energy equation-of-
state [14, 10, 8] and the sum of the neutrino masses [15, 16]. In addition,
the cross-correlation between CMB lensing and other tracers of large-scale
structure (LSS), e.g., galaxy surveys, can help circumvent cosmic variance
[17, 18] and break degeneracies between cosmological parameters, for example,
galaxy bias and growth in the LSS [19, 20, 21].

Upcoming experiments including Simons Observatory (SO) [22] and CMB-S4 [23]
will provide measurements of the CMB sky with unprecedented precision. In this
context, CMB lensing plays a major role in refining our knowledge of the late-
time Universe over a wide range of redshifts, both as a cosmological probe on
its own as well as in combination with correlated LSS tracers. Thus far,
cosmological applications of CMB lensing have relied on the assumption that
the lensing potential is Gaussian, namely that the underlying matter
distribution is well described by linear theory. However, at late times,
nonlinear clustering in the matter density fluctuations generates large non-
Gaussianities, which in turn induce a small, but non-zero, bispectrum in the
CMB lensing convergence field. While this effect is negligible for the
sensitivity reached by, e.g., Planck, future CMB experiments should be
sensitive to the non-Gaussian convergence field [24, 25, 26]. Interestingly,
it has been shown that ongoing Stage-III experiments, like the Atacama
Cosmology Telescope (ACT) [27] and Simons Array [28], may already be able to
measure the lensing convergence bispectrum arising from the non-linear growth
of the LSS – the next higher-order correlation function – and that in the
future Stage-IV experiments, like CMB-S4, have the potential to detect it at
$50\,\sigma$ statistical significance [24].

A simple and efficient (albeit somewhat suboptimal) method to reconstruct the
lensing potential employs estimators that are quadratic in the CMB
anisotropies. Here, we work with the widely used Hu-Okamoto (HO) quadratic
estimator [29, 30]. Because of the nature of the quadratic estimator,
additional terms called noise-biases arise when estimating the $n$-point
functions of the lensing potential, and are comparable to or even larger than
the signal in the recent power spectrum measurements. These noise-biases have
been worked out in detail for the power spectrum in Refs. [31, 32] and methods
developed for their practical subtraction in survey data [33, 11]. However,
when including non-Gaussian effects in the convergence field, new
reconstruction bias terms contribute to the measured lensing power spectrum
and must be taken into account in future analysis (see, e.g., Ref. [34, 35]).
This non-Gaussian bias is also important for the cross-correlation between CMB
lensing and other LSS tracers [36].

Besides these additional biases in the measurement of the lensing power
spectrum, non-Gaussianity of the lensing convergence field encodes
cosmological information that could be used to tighten the bounds on neutrino
masses and dark energy models. Fisher forecasts show that the expected
constraints on the dark energy equation-of-state and sum of neutrino masses
from the lensing bispectrum111Note that from here on we will use the
terminology lensing bispectrum, for which we intend the bispectrum of the
lensing convergence field, not to be confused with the CMB bispectrum induced
by lensing, such as ISW-lensing [37, 38]. alone are comparable to those
obtained from the power spectrum [24]. Additionally, the combination of the
power spectrum and the bispectrum suggests the potential to improve
constraints by around $30\,\%$, making the lensing bispectrum an appealing
cosmological tool comparable and complementary to the lensing power spectrum.

In this work, we aim to characterize for the first time the leading-order
noise-biases222The largest noise-biases for the lensing bispectrum are
Gaussian, therefore they do not introduce a bias in the bispectrum signal.
However, we keep this nomenclature to highlight the analogy with the power
spectrum reconstruction. in the reconstructed lensing convergence bispectrum.
We limit our analysis to temperature modes only, but the extension to
polarization is straightforward. The procedure to compute these biases follows
the same perturbative methods that were applied to obtain the noise-biases for
the power spectrum. However, the bispectrum case can be significantly more
challenging compared to the power spectrum, as the bispectrum of the estimated
convergence field depends on the 6-pt correlation function of the lensed CMB
modes. While in principle the calculations can still be carried out using
Wick’s theorem for the products of (Gaussian) unlensed CMB fields and low-
order cumulant expansions for the (non-Gaussian) lensing convergence fields
that appear in the perturbative expansion of the lensed CMB $n$-pt functions,
keeping track of all terms in the expansion can become increasingly difficult.
Therefore we will additionally use and adapt a more convenient way to organize
these calculations based on Feynman diagrams [39, 40], which allows one to
write down the noise-biases without resorting to lengthy calculations. The
noise-bias term can be identified by combining the proper propagators,
similarly to how these diagrams are used in quantum field theory. In this way,
we will be able to characterize the noise-bias terms $N_{B}^{(0)}$,
$N_{B}^{(1)}$, $N_{B}^{(3/2)}$ in the flat-sky limit (the superscripts on the
bispectrum noise-biases denote the leading-order dependence on the lensing
potential power spectrum). Finally, in order to validate our analytical and
numerical computations of the noise-biases, we will compare the results of the
largest noise-biases, $N_{B}^{(0)}$ and $N_{B}^{(1)}$, with simulations of the
lensed temperature CMB maps.

This paper is organized as follows. In Sec. 2, we introduce the HO quadratic
estimator for the lensing convergence field and establish our notation. In
Sec. 3, we discuss the origin of noise-biases that arise in the bispectrum of
the reconstructed convergence, starting from the cumulant expansion of the
lensed CMB 6-pt function. We also present the direct (“Wick expansion”)
calculation of these biases along with the diagrammatic approach. In Sec. 4,
we show our numerical results and in Sec. 5 we validate $N_{B}^{(0)}$ and
$N_{B}^{(1)}$ using simulations of the lensed temperature CMB maps. We draw
our conclusions and outline perspectives for future work in Sec. 6. Finally,
in App. A we give the full technical details of the computation of the
$N_{B}^{(3/2)}$ bias, and in App. B we describe an alternative set of
simulations to test our theoretical findings.

##### Notation and conventions

Throughout the paper, we adopt the flat-sky approximation. We denote 2D
Fourier wavenumbers by $\bm{\ell}$ for CMB temperature fields and $\bm{L}$ for
the lensing potential/convergence. In general, the modulus of a vector is
defined as, e.g., $\ell\equiv|\bm{\ell}|$.

The power spectrum and the bispectrum of the relevant observables are defined,
respectively, as

$\langle
X(\bm{\ell})X(\bm{\ell}^{\prime})\rangle=(2\pi)^{2}\,\delta^{(2)}(\bm{\ell}+\bm{\ell}^{\prime})C_{\ell}^{XX},$
(1.1)

and

$\langle
X(\bm{\ell}_{1})\,X(\bm{\ell}_{2})\,X(\bm{\ell}_{3})\rangle=(2\pi)^{2}\delta^{(2)}(\bm{\ell}_{1}+\bm{\ell}_{2}+\bm{\ell}_{3})\,B_{X}(\ell_{1},\ell_{2},\ell_{3})\,,$
(1.2)

where $X$ can be either $T$, $\phi$ or $\kappa$. We further use the notation
$n_{C}$ to denote connected $n$-pt functions. We use the Fourier convention
(in 2D)

$\displaystyle\mathcal{F}[Y(\hat{\bm{n}})](\bm{\ell})$
$\displaystyle=\int\mathrm{d}^{2}\hat{\bm{n}}\,Y(\hat{\bm{n}})\,e^{-i\bm{\ell}\cdot\hat{\bm{n}}},$
(1.3) $\displaystyle\mathcal{F}^{-1}[Y(\bm{\ell})](\hat{\bm{n}})$
$\displaystyle=\int\frac{\mathrm{d^{2}}\bm{\ell}}{(2\pi)^{2}}\,Y(\bm{\ell})\,e^{i\bm{\ell}\cdot\hat{\bm{n}}}.$
(1.4)

For brevity, we introduce the compact notations $\bm{\ell}_{ij\dots
n}=\bm{\ell}_{i}+\bm{\ell}_{j}+\dots+\bm{\ell}_{n}$,
$\bar{\bm{L}}_{ij}=\bm{L}_{i}-\bm{\ell}_{j}$ (where either $i=j$ or $i\neq j$)
and

$\int_{\bm{\ell}_{i},\dots,\,\bm{\ell}_{n}}=\int\frac{\mathrm{d^{2}}\bm{\ell}_{i}}{(2\pi)^{2}}\ldots\int\frac{\mathrm{d^{2}}\bm{\ell}_{n}}{(2\pi)^{2}}\,.$

Additionally, the lensed temperature fields
$\widetilde{T}(\bm{L}_{i}-\bm{\ell}_{j})$ and $\widetilde{T}(\bm{\ell}_{i})$
are sometimes represented as $\widetilde{T}_{ij}$ and $\widetilde{T}_{i}$,
respectively. In our numerical computations, we consider a flat $\Lambda$CDM
cosmology, with cosmological parameters in accordance with the latest Planck
results [41], summarized in Tab. 1.

|  $\bm{\Lambda}$CDM parameters |
---|---|---
$H_{0}=67.32\,\text{km\,s}^{-1}\,\text{Mpc}^{-1}$ | $\Omega_{\mathrm{b}}h^{2}=0.02238$ | $\sum m_{\nu}=0.06\,\text{eV}$
$\Omega_{K}=0$ | $\Omega_{\mathrm{c}}h^{2}=0.12010$ | $\tau=0.0543$
$n_{s}=0.9660$ | $A_{\mathrm{s}}=2.1005\times 10^{-9}$ | $r=0$
Table 1: Best-fit Planck parameters (specifically, Tab. 1 of Ref. [41] with
$TT,TE,EE\text{+low$E$+lensing+BAO}$) used in our numerical computations.

##  2 The Hu & Okamoto quadratic estimator

In this section, we briefly review the effect of weak gravitational lensing on
the CMB temperature anisotropies, and introduce the HO quadratic estimator
[29] in the flat-sky approximation.

The effect of lensing is usually expressed as a remapping of the primary CMB
anisotropies

$\widetilde{T}(\hat{\bm{n}})=T(\hat{\bm{n}}^{\prime})=T(\hat{\bm{n}}+\nabla\phi(\hat{\bm{n}}))\,,$
(2.1)

where $\hat{\bm{n}}$ is the line-of-sight direction on the sky, $T$ refers to
the unlensed temperature field and $\widetilde{T}$ indicates the lensed
counterpart. We shall adopt this notation from here on.

The lens remapping in Eq. (2.1) is given by the deflection angle
$\bm{d}(\hat{\bm{n}})=\nabla\phi(\hat{\bm{n}})$, where the lensing potential
$\phi(\hat{\bm{n}})$ is determined by the projection of the gravitational
potential along the line-of-sight:

$\phi(\hat{\bm{n}})=-2\,\int_{0}^{\chi_{*}}\mathrm{d}\chi\,\frac{\chi_{*}-\chi}{\chi_{*}\chi}\,\Psi(\chi\hat{\bm{n}};\eta_{0}-\chi)\,.$
(2.2)

Here, $\chi$ denotes the comoving distance, $\eta$ the conformal time (with
present-day value $\eta_{0}$) and $\Psi$ is the (Weyl) gravitational
potential. The definition in Eq. (2.2) assumes the Born approximation, whereby
the integration is carried out along the undeflected photon path, which is
valid when the deflection angles are small compared to the scales on which the
lenses vary.333Post-Born corrections to lensing are, however, an important
source of the convergence bispectrum along with the contribution from non-
Gaussian matter clustering [25]. Post-Born effects also introduce a curl-mode
to the lensing deflections, so they are no longer describable by the lensing
potential alone. We ignore the small effect of curl-modes in the deflections
throughout this paper and leave their analysis to future work.

In the following we will also use the lensing convergence
$\kappa(\hat{\bm{n}})$, which captures the isotropic distortion due to weak
gravitational lensing and is related to the lensing potential via

$\kappa(\hat{\bm{n}})=-\frac{1}{2}\nabla^{2}\phi(\hat{\bm{n}})\,.$ (2.3)

We can expand the lensed temperature modes Eq. (2.1) (up to second order in
$\phi$) as

$\widetilde{T}(\hat{\bm{n}})=T(\hat{\bm{n}})+\nabla
T(\hat{\bm{n}})\cdot\nabla\phi(\hat{\bm{n}})+\frac{1}{2}\nabla_{i}\nabla_{j}T(\hat{\bm{n}})\,\,\nabla^{i}\phi(\hat{\bm{n}})\,\nabla^{j}\phi(\hat{\bm{n}})\,+\mathcal{O}(\phi^{3})\,.$
(2.4)

Replacing the Fourier transform of the CMB temperature anisotropies and the
lensing potential as defined in Eq. (1.3), we find

$\widetilde{T}(\bm{\ell})=T(\bm{\ell})+\delta
T(\bm{\ell})+\delta^{2}T(\bm{\ell})+\mathcal{O}(\phi^{3})\,,$ (2.5)

where the $\mathcal{O}(\phi)$ and $\mathcal{O}(\phi^{2})$ terms are
respectively given by

$\displaystyle\delta T(\bm{\ell})$
$\displaystyle=-\int_{\bm{p}}\,\Big{[}\bm{p}\cdot(\bm{\ell}-\bm{p})\Big{]}\,\phi(\bm{\ell}-\bm{p})\,T(\bm{p})\,,$
(2.6) $\displaystyle\delta^{2}T(\bm{\ell})$
$\displaystyle=-\frac{1}{2}\int_{\bm{p}}\int_{\bm{q}}\,\Big{[}\bm{p}\\!\cdot\\!(\bm{p}+\bm{q}-\bm{\ell})\Big{]}\,(\bm{p}\\!\cdot\\!\bm{q})\,T(\bm{p})\,\phi(\bm{q})\,\phi(\bm{p}+\bm{q}-\bm{\ell})\,.$
(2.7)

Effectively, a fixed lens introduces a correlation between temperature modes
that would otherwise be independent. Averaging over the primary CMB
fluctuations for fixed lenses, we have at leading order

$\langle\widetilde{T}(\bm{\ell}_{1})\widetilde{T}(\bm{\ell}_{2})\rangle_{\text{CMB}}=(2\pi)^{2}\delta^{(2)}(\bm{\ell}_{1}+\bm{\ell}_{2})C_{\ell_{1}}^{TT}+f(\bm{\ell}_{1},\bm{\ell}_{2})\phi(\bm{\ell}_{1}+\bm{\ell}_{2})\,,$
(2.8)

where the lensing response function

$f(\bm{\ell}_{1},\bm{\ell}_{2})\equiv(\bm{\ell}_{1}+\bm{\ell}_{2})\cdot\left(\bm{\ell}_{1}C_{\ell_{1}}^{TT}+\bm{\ell}_{2}C_{\ell_{2}}^{TT}\right)\,,$
(2.9)

and $C_{\ell}^{TT}$ is the angular power spectrum of the unlensed CMB.

It is sometimes useful to consider a non-perturbative extension of the lensing
response function, defined by

$\left\langle\frac{\delta}{\delta\phi(\bm{L})}\left[\widetilde{T}(\bm{\ell}_{1})\widetilde{T}(\bm{\ell}_{2})\right]\right\rangle=\tilde{f}(\bm{\ell}_{1},\bm{\ell}_{2})\delta^{(2)}(\bm{L}-\bm{\ell}_{1}-\bm{\ell}_{2})\,,$
(2.10)

where the expectation value is now over both the unlensed CMB and the lenses.
It is shown in Ref. [42] that the non-perturbative response function has the
same form as Eq. (2.9) but with the unlensed CMB power spectrum replaced by
the spectrum $\tilde{C}_{\ell}^{T\nabla T}$, which describes the correlation
between the lensed temperature and the lensed gradient of the temperature. The
non-perturbative response function is useful as it captures (some) beyond-
leading-order terms due to lensing in the connected moments of the lensed CMB,
whose biases can easily be mitigated by an appropriate modification of the
normalisation of the quadratic estimator used for lensing reconstruction.
Throughout this work, we will adopt the non-perturbative approach and use Eq.
(2.10), which we rename $f(\bm{\ell}_{1},\bm{\ell}_{2})$ for convenience.

From the resulting off-diagonal terms in the 2-pt correlation function, we can
introduce an estimator to reconstruct the lensing potential from CMB maps that
is quadratic in the anisotropies:

$\hat{\phi}(\bm{L})=A_{L}\int_{\bm{\ell}_{1},\bm{\ell}_{2}}\,(2\pi)^{2}\delta^{(2)}(\bm{L}-\bm{\ell}_{1}-\bm{\ell}_{2})F(\bm{\ell}_{1},\bm{\ell}_{2})\,\widetilde{T}(\bm{\ell}_{1})\widetilde{T}(\bm{\ell}_{2})\,.$
(2.11)

Here, the weight function

$F(\bm{\ell}_{1},\bm{\ell}_{2})=\frac{f(\bm{\ell}_{1},\bm{\ell}_{2})}{2C^{TT,\,\mathrm{tot}}_{\ell_{1}}C^{TT,\,\mathrm{tot}}_{\ell_{2}}}$
(2.12)

is chosen to minimize the Gaussian variance of the estimator under the
constraint that the estimator is unbiased (see below).

The power spectrum $C^{TT,\,\mathrm{tot}}_{\ell}$ denotes the total
temperature power in the map, including CMB, noise, and foreground power.
Given that the goal of this work is to give a first estimate of the noise-
biases, we generally do not consider any particular CMB experimental
specifications and we assume that the total power spectrum is given by the
lensed power spectrum,
$C^{TT,\,\mathrm{tot}}_{\ell}=C^{\tilde{T}\tilde{T}}_{\ell}$.

The normalization of the quadratic estimator $A_{L}$ is derived assuming that
the resulting lensing map is unbiased, i.e.,
$\langle\delta\hat{\phi}(\bm{L})/\delta\phi(\bm{L}^{\prime})\rangle=\delta^{(2)}(\bm{L}-\bm{L}^{\prime})$
(cf. Eq. 2.10); it is given by

$A_{L}=\left[\int_{\bm{\ell}_{1}}\,F(\bm{\ell}_{1},\bm{L}-\bm{\ell}_{1})f(\bm{\ell}_{1},\bm{L}-\bm{\ell}_{1})\right]^{-1}\,.$
(2.13)

Finally, using Eq. (2.3) in Fourier space, we may write the quadratic
estimator for the lensing convergence field as

$\hat{\kappa}(\bm{L})=\frac{1}{2}L^{2}\hat{\phi}(\bm{L})=A_{L}^{\kappa}\int_{\bm{\ell}_{1},\bm{\ell}_{2}}\,(2\pi)^{2}\delta^{(2)}(\bm{L}-\bm{\ell}_{1}-\bm{\ell}_{2})F(\bm{\ell}_{1},\bm{\ell}_{2})\,\widetilde{T}(\bm{\ell}_{1})\widetilde{T}(\bm{\ell}_{2})\,,$
(2.14)

where $A_{L}^{\kappa}=L^{2}A_{L}/2$.

##  3 Noise-biases in the reconstructed lensing bispectrum

It is well known that statistical reconstruction noise, arising from chance
fluctuations in the CMB fields that mimic the effects of lensing, is present
in the OH quadratic estimator. This leads to “noise-biases” when using the
$n$-point correlation functions of the reconstructed $\hat{\kappa}$ to
estimate the correlation functions of the underlying lensing convergence (see,
e.g., [32] for the power spectrum case). In this paper, we calculate these
biases for the first time for the bispectrum of the reconstructed convergence.
Quantifying the biases carefully is crucial as their magnitude is typically
much larger than the signal we aim to detect. For power spectrum
reconstruction on data, the biases are now commonly removed with realization-
dependent methods that mix data and simulations [32, 33, 21]. These methods
reduce sensitivity to any mismatch between the statistical properties of the
simulations and reality (e.g., due to complex instrumental noise properties)
and also reduce covariance between multipoles of the reconstructed power
spectrum. Similar methods have been developed for a general $n$-pt function
and need to be tested for the bispectrum. However, as a first step here we
shall simply compute the dominant biases in the ideal case of an isotropic
(and, generally, noiseless) survey.

A generic $n$-pt function of the lensing convergence field estimator reads as

$\displaystyle\langle\hat{\kappa}(\bm{L}_{1})\hat{\kappa}(\bm{L}_{2})\ldots\hat{\kappa}(\bm{L}_{n})\rangle_{{\rm
CMB},\phi}$
$\displaystyle=\int^{\prime}_{\bm{\ell}_{1},\bm{\ell}_{2},\dots,\,\bm{\ell}_{n}}F(\bm{\ell}_{1},\bar{\bm{L}}_{11})\,F(\bm{\ell}_{2},\bar{\bm{L}}_{22})\ldots
F(\bm{\ell}_{n},\bar{\bm{L}}_{nn})$
$\displaystyle\qquad\qquad\qquad\times\langle\widetilde{T}_{1}\widetilde{T}_{11}\widetilde{T}_{2}\widetilde{T}_{22}\cdots\widetilde{T}_{n}\widetilde{T}_{nn}\rangle_{{\rm
CMB},\phi}\,,$ (3.1)

where the hat symbols refers to the reconstructed observable and

$\int^{\prime}_{\bm{\ell}_{1},\bm{\ell}_{2},\dots,\,\bm{\ell}_{n}}\equiv\left(\prod_{i=1}^{3}A_{L_{i}}^{\kappa}\right)\,\int_{\bm{\ell}_{1},\bm{\ell}_{2},\dots,\,\bm{\ell}_{n}}.$
(3.2)

The average is performed over many realizations of either the CMB or the
lensing potential, as indicated by the suffix ‘${\rm CMB},\phi$’. From now on,
we drop this notation the sake of brevity.

In order to compute Eq. (3), we need to evaluate the $2n$-pt function of the
lensed temperature fields

$\langle\widetilde{T}(\bm{\ell}_{1})\,\widetilde{T}(\bm{L}_{1}-\bm{\ell}_{1})\,\widetilde{T}(\bm{\ell}_{2})\,\widetilde{T}(\bm{L}_{2}-\bm{\ell}_{2})\,...\,\widetilde{T}(\bm{\ell}_{n})\,\widetilde{T}(\bm{L}_{n}-\bm{\ell}_{n})\rangle\,,$
(3.3)

which can be decomposed in connected $m$-pt functions ($m\leq 2n$). In the
case of the 3-pt function of $\hat{\kappa}$, the reconstructed lensing signal
receives contributions from the 6-pt function of lensed CMB modes

$\displaystyle\langle\widetilde{T}\,\widetilde{T}\,\widetilde{T}\,\widetilde{T}\,\widetilde{T}\,\widetilde{T}\rangle=\,$
$\displaystyle\langle\widetilde{T}\,\widetilde{T}\rangle_{C}\,\langle\widetilde{T}\,\widetilde{T}\rangle_{C}\,\langle\widetilde{T}\,\widetilde{T}\rangle_{C}+\langle\widetilde{T}\,\widetilde{T}\,\widetilde{T}\,\widetilde{T}\rangle_{C}\,\langle\widetilde{T}\,\widetilde{T}\rangle_{C}$
$\displaystyle+\langle\widetilde{T}\,\widetilde{T}\,\widetilde{T}\rangle_{C}\,\langle\widetilde{T}\,\widetilde{T}\,\widetilde{T}\rangle_{C}+\langle\widetilde{T}\,\widetilde{T}\,\widetilde{T}\,\widetilde{T}\,\widetilde{T}\,\widetilde{T}\rangle_{C}+(\mathrm{perm.})\,,$
(3.4)

where ‘perm.’ stands for permutations over the multipoles $\bm{\ell}_{i}$. In
the following, we shall ignore the 3-pt function of the lensed CMB,
$\langle\tilde{T}\tilde{T}\tilde{T}\rangle_{C}$. This is only sourced by the
correlation between lensing and the late-time integrated Sachs–Wolfe
contribution to the unlensed CMB, and should be negligible for the small-scale
CMB modes that dominate lensing reconstruction.

When expanding the connected moments in Eq. (3) perturbatively in $\phi$
(using Eq. (2.5)), we obtain the signal bispectrum and a number of noise-
biases:

$\begin{split}\langle\hat{\kappa}(\bm{L}_{1})\hat{\kappa}(\bm{L}_{2})\hat{\kappa}(\bm{L}_{3})\rangle=(2\pi)^{2}\delta(\bm{L}_{123})\bigg{[}B_{\kappa}(\bm{L}_{1},\bm{L}_{2},\bm{L}_{3})+N_{\mathrm{B}}^{(0)}+N_{\mathrm{B}}^{(1)}+N_{\mathrm{B}}^{(3/2)}+N_{\mathrm{B}}^{(2)}\bigg{]}\,.\end{split}$
(3.5)

The signal term $B_{\kappa}$ arises from the $\mathcal{O}(\phi^{3})$ part of
the connected 6-pt function of the lensed CMB in which the unlensed CMB fields
within each quadratic estimator are contracted. This involves terms like

$\langle\widetilde{T}_{1}\widetilde{T}_{11}\widetilde{T}_{2}\widetilde{T}_{22}\widetilde{T}_{3}\widetilde{T}_{33}\rangle_{C}\supset\langle\wick{\c{1}\delta
T_{1}\c{1}T_{11}\;\c{1}\delta T_{2}\c{1}T_{22}\;\c{1}\delta
T_{3}\c{1}T_{33}}\rangle_{\phi}\,,$ (3.6)

where the indicated (Gaussian) contractions are over the unlensed CMB and the
remaining expectation over $\phi$ returns the bispectrum. For the noise-biases
in Eq. (3.5), the superscripts denote the leading order at which they arise in
an expansion in $\phi$, with $N^{(j/2)}_{\text{B}}=\mathcal{O}(\phi^{j})$.
Similarly to lensing power spectrum reconstruction, $N_{\mathrm{B}}^{(0)}$
arises from the totally disconnected term in the decomposition of Eq. (3) and
would be present even in the absence of lensing. Here, we compute it non-
perturbatively using the lensed temperature power spectrum. The noise-bias
$N_{\mathrm{B}}^{(1)}$ comes from the connected 4-pt function up to order
$\mathcal{O}(\phi^{2})$, namely linear in the lensing potential power spectrum
(thus $j=2$). For instance, one such contributing term is

$\langle\widetilde{T}\widetilde{T}\widetilde{T}\widetilde{T}\rangle_{C}\langle\widetilde{T}\widetilde{T}\rangle_{C}\supset\langle\wick{\c{1}\delta
T\c{1}T\;\c{1}\delta
T\c{1}T}\rangle_{\phi}\langle\widetilde{T}\widetilde{T}\rangle_{C}\,,$ (3.7)

where, again, the contractions over the unlensed CMB are shown explicitly. The
terms of order $\mathcal{O}(\phi^{3})$ in the connected 4- and 6-pt functions
are captured by $N_{\mathrm{B}}^{(3/2)}$, which is linear in the lensing
potential bispectrum ($j=3$). Finally, $N_{\mathrm{B}}^{(2)}$ arises from the
part of the connected 6-pt function that is quadratic in the lensing potential
power spectrum ($j=4$) and from corrections at this order to the connected
4-pt function. These terms are therefore of order $\mathcal{O}(\phi^{4})$. See
Tab. 2 for a summary.

In the following we review two complementary approaches to evaluate
analytically these noise-biases, the Wick-expansion and a Feynman diagram
approach, and show the equivalence of the two methods by evaluating the
$N^{(0)}_{\mathrm{B}}$ and $N^{(1)}_{\mathrm{B}}$ noise-biases. We leave the
computation of $N_{\mathrm{B}}^{(2)}$ for future work. The largest terms that
are formally $\mathcal{O}(\phi^{4})$ should already be captured by our use of
non-perturbative response functions in the calculation of $N^{(1)}$ [42] and
the lensed CMB power spectra for
$\langle\widetilde{T}\widetilde{T}\rangle_{C}$ terms, and our numerical
results in Sec. 3 show that the residual $N^{(2)}_{\text{B}}$ noise-bias is
much smaller than the other bias terms .

| Noise-biases
---|---
| $\langle(\widetilde{T})^{2}\rangle_{C}\langle(\widetilde{T})^{2}\rangle_{C}\langle(\widetilde{T})^{2}\rangle_{C}$ | $\langle(\widetilde{T})^{4}\rangle_{C}\langle(\widetilde{T})^{2}\rangle_{C}$ | $\langle(\widetilde{T})^{6}\rangle_{C}$
all orders | $N_{\mathrm{B}}^{(0)}$ |  |
$\mathcal{O}(\phi^{2})$ |  | $N_{\mathrm{B}}^{(1)}$ |
$\mathcal{O}(\phi^{3})$ |  | $N_{\mathrm{B}}^{(3/2)}$
$\mathcal{O}(\phi^{4})$ |  | $N_{\mathrm{B}}^{(2)}$
Table 2: Noise-bias terms of
$\langle\hat{\kappa}\hat{\kappa}\hat{\kappa}\rangle$ and their origin in terms
of the decomposition (3).

###  3.1 Wick-expansion for correlation functions

In general, an $n$-pt function as in Eq. (3.3) can be computed with Wick’s
theorem in the case that the unlensed CMB and $\phi$ are Gaussian. In this
case, the only fundamental contraction is the 2-pt function defined in Eq.
(1.1):

$\contraction[1ex]{}{X}{(\bm{\ell}_{1})\,\,}{X}X(\bm{\ell}_{1})\,\,X(\bm{\ell}_{2})\,,\quad
X=\\{T,\phi\\}.$ (3.8)

On the other hand, if some level of non-Gaussianity arises in either the
temperature modes444We shall ignore this possibility in the rest of this
paper. or the lensing potential $\phi$, resulting in non-zero bispectra (1.2),
we have an additional fundamental contraction

$\contraction[1ex]{}{\phi}{(\bm{\ell}_{1})\,\,}{X}\contraction[1ex]{}{\phi}{(\bm{\ell}_{1})\,\,T(\bm{\ell}_{2})\,\,}{X}X(\bm{\ell}_{1})\,\,X(\bm{\ell}_{2})\,\,X(\bm{\ell}_{3})\,,\quad
X=\\{T,\phi\\}.$ (3.9)

Then, the lensed temperature $2n$-pt function (3.3) can be evaluated by
expanding the lensed temperature modes in terms of the lensing potential and
by taking all the terms allowed by the fundamental contractions Eqs.
(3.8–3.9).

####  $N^{(0)}_{\mathrm{B}}$ noise-bias: Wick approach.

The $N_{\mathrm{B}}^{(0)}$ noise-bias comes from the totally disconnected term
of the decomposition (3) and reflects the non-Gaussian nature of the
reconstruction noise in $\hat{\kappa}$. Here, we consider the full (non-
perturbative) term

$\begin{split}\langle\widetilde{T}(\bm{\ell}_{1})\widetilde{T}(\bm{L}_{1}-\bm{\ell}_{1})\rangle_{C}\,\langle\widetilde{T}(\bm{\ell}_{2})\widetilde{T}(\bm{L}_{2}-\bm{\ell}_{2})\rangle_{C}\,\langle\widetilde{T}(\bm{\ell}_{3})\widetilde{T}(\bm{L}_{3}-\bm{\ell}_{3})\rangle_{C}+\mathrm{(perms.)}.\,\end{split}$
(3.10)

In total, there are 15 permutations that need to be taken into account, given
that we must combine six momenta in groups of two. However, seven of these
permutations involve at least one contraction over fields in the same
quadratic estimator and are therefore proportional to $\delta(\bm{L}_{i})$.
For the $\bm{L}_{i}\neq 0$, we are left with eight terms, for example,

$\langle\widetilde{T}_{1}\widetilde{T}_{3}\rangle_{C}\,\langle\widetilde{T}_{2}\widetilde{T}_{11}\rangle_{C}\,\langle\widetilde{T}_{22}\widetilde{T}_{33}\rangle_{C}\propto\delta^{(2)}(\bm{\ell}_{1}+\bm{\ell}_{3})\delta^{(2)}(\bm{\ell}_{2}+\bar{\bm{L}}_{11})\delta^{(2)}(\bar{\bm{L}}_{22}+\bar{\bm{L}}_{33})C^{\widetilde{T}\widetilde{T}}_{\ell_{1}}\,C^{\widetilde{T}\widetilde{T}}_{\ell_{2}}\,C^{\widetilde{T}\widetilde{T}}_{|\bar{\bm{L}}_{22}|}\,.$
(3.11)

Using the fact that the integrand in Eq. (3) is symmetrical with respect to an
exchange of momenta $\bm{\ell}_{i}\leftrightarrow\bm{L}_{i}-\bm{\ell}_{i}$, it
can be shown that the remaining terms are equivalent and so the lowest-order
noise-bias is

$\displaystyle N^{(0)}_{\mathrm{B}}=8\int^{\prime}_{\bm{\ell}_{1}}\,$
$\displaystyle
F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{L}_{1}-\bm{\ell}_{1},\bm{L}_{3}+\bm{\ell}_{1})\,F(-\bm{\ell}_{1},\bm{L}_{3}+\bm{\ell}_{1})$
$\displaystyle\times
C^{\widetilde{T}\widetilde{T}}_{\ell_{1}}\,C^{\widetilde{T}\widetilde{T}}_{|\bm{L}_{1}-\bm{\ell}_{1}|}\,C^{\widetilde{T}\widetilde{T}}_{|\bm{L}_{3}+\bm{\ell}_{1}|}\,.$
(3.12)

Note that in the presence of instrumental noise, the lensed power spectra here
would be replaced by $C_{\ell}^{TT,\text{tot}}$.

####  $N^{(1)}_{\mathrm{B}}$ noise-bias: Wick approach.

The $N_{\mathrm{B}}^{(1)}$ noise-bias arises from the $4_{C}\times 2_{C}$ term
in Eq. (3) and its leading-order contribution is $\mathcal{O}(\phi^{2})$.Note
that terms of this order also appear in $N^{(0)}_{\mathrm{B}}$, but have
already been accounted for there by using the lensed temperature power spectra
in Eq. (3.1). Therefore, we are left with terms in $4_{C}\times 2_{C}$ such as

$\langle\widetilde{T}(\bm{\ell}_{1})\widetilde{T}(\bm{L}_{1}-\bm{\ell}_{1})\widetilde{T}(\bm{\ell}_{2})\widetilde{T}(\bm{L}_{2}-\bm{\ell}_{2})\rangle_{C}\,\langle\widetilde{T}(\bm{\ell}_{3})\widetilde{T}(\bm{L}_{3}-\bm{\ell}_{3})\rangle_{C}+\mathrm{(perms.)}\,.$
(3.13)

Let us first briefly show the computation of the connected 4-pt function
(details can be found in Ref. [31]). At leading order, we have

$\displaystyle\langle\widetilde{T}(\bm{\ell}_{1})\widetilde{T}(\bm{\ell}_{2})\widetilde{T}(\bm{\ell}_{3})\widetilde{T}(\bm{\ell}_{4})\rangle_{C}$
$\displaystyle=\frac{1}{2}\langle\wick{\c{1}T(\bm{\ell}_{1})\c{2}T(\bm{\ell}_{2})\c{1}\delta
T(\bm{\ell}_{3})\c{2}\delta T(\bm{\ell}_{4})}\rangle_{\phi}+(\text{23
perms.})$
$\displaystyle=\frac{1}{8}f(\bm{\ell}_{1},\bm{\ell}_{3})f(\bm{\ell}_{2},\bm{\ell}_{4})\langle\phi(\bm{\ell}_{1}+\bm{\ell}_{3})\phi(\bm{\ell}_{2}+\bm{\ell}_{4})\rangle+(\text{23
perms.})$
$\displaystyle=\frac{1}{8}(2\pi)^{2}\delta^{(2)}(\bm{\ell}_{1234})f(\bm{\ell}_{1},\bm{\ell}_{3})f(\bm{\ell}_{2},\bm{\ell}_{4})C^{\phi\phi}_{|\bm{\ell}_{1}+\bm{\ell}_{3}|}C^{\phi\phi}_{|\bm{\ell}_{1}+\bm{\ell}_{3}|}+(\text{23
perms.})\,,$ (3.14)

where the contractions indicated in the top line are over the unlensed CMB and
we have used555In practice, we use the non-perturbative response functions to
account for some higher-order terms in $\phi$. In particular, for a given pair
of lensed temperature fields, say
$\widetilde{T}(\bm{\ell}_{1})\widetilde{T}(\bm{\ell}_{3})$, whose unlensed
fields are contracted, the non-perturbative response functions account for
contractions over additional pairs of $\phi$ fields within
$\widetilde{T}(\bm{\ell}_{1})\widetilde{T}(\bm{\ell}_{3})$ leaving one $\phi$
field uncontracted (i.e., as an external leg).

$\wick{\c{1}\delta
T(\bm{\ell}_{1})\c{1}T(\bm{\ell}_{2})}+\wick{\c{1}T(\bm{\ell}_{1})\c{1}\delta
T(\bm{\ell}_{2})}=f(\bm{\ell}_{1},\bm{\ell}_{2})\phi(\bm{\ell}_{1}+\bm{\ell}_{2})\,.$
(3.15)

Gathering the permutations, we have

$\langle\widetilde{T}(\bm{\ell}_{1})\widetilde{T}(\bm{\ell}_{2})\widetilde{T}(\bm{\ell}_{3})\widetilde{T}(\bm{\ell}_{4})\rangle_{C}=(2\pi)^{2}\delta^{(2)}(\bm{\ell}_{1234})\mathcal{T}(\bm{\ell}_{1},\bm{\ell}_{2},\bm{\ell}_{3},\bm{\ell}_{4})\,,$
(3.16)

where the trispectrum is

$\mathcal{T}(\bm{\ell}_{1},\bm{\ell}_{2},\bm{\ell}_{3},\bm{\ell}_{4})=C^{\phi\phi}_{|\bm{\ell}_{1}+\bm{\ell}_{2}|}f(\bm{\ell}_{1},\bm{\ell}_{2})f(\bm{\ell}_{3},\bm{\ell}_{4})+C^{\phi\phi}_{|\bm{\ell}_{1}+\bm{\ell}_{3}|}f(\bm{\ell}_{1},\bm{\ell}_{3})f(\bm{\ell}_{2},\bm{\ell}_{4})\\\
+C^{\phi\phi}_{|\bm{\ell}_{1}+\bm{\ell}_{4}|}f(\bm{\ell}_{1},\bm{\ell}_{4})f(\bm{\ell}_{2},\bm{\ell}_{3})\,.$
(3.17)

Given this result, we may identify all the terms contributing to Eq. (3.13)
when summing over all the possible permutations of the momenta. In principle,
we have 15 terms, but three of them are proportional to the
$\delta(\bm{L}_{i})$ since the $2_{C}$ term involves fields within the same
quadratic estimator (as in the term shown explicitly in Eq. 3.13). Then,
thanks to the symmetry of the integrand in Eq. (3) with respect to the
exchange of $\bm{\ell}_{i}$ and $\bm{L}_{i}-\bm{\ell}_{i}$, only three of the
remaining 12 terms turn out to be independent (with a factor of 4) and are
given by

$\displaystyle N^{(1)}_{\mathrm{B}}$
$\displaystyle=4\int^{\prime}_{\bm{\ell}_{1},\,\bm{\ell}_{2}}\,F(\bm{\ell}_{1},\bar{\bm{L}}_{11})\,F(\bm{\ell}_{2},\bar{\bm{L}}_{22})\,F(-\bm{\ell}_{2},\bm{L}_{3}+\bm{\ell}_{2})\,C^{\widetilde{T}\widetilde{T}}_{\ell_{2}}\Big{[}C^{\phi\phi}_{L_{1}}\,f(\bm{\ell}_{1},\bar{\bm{L}}_{11})\,f(\bar{\bm{L}}_{22},\bm{L}_{3}+\bm{\ell}_{2})$
$\displaystyle\quad+\,C^{\phi\phi}_{|\bm{\ell}_{1}+\bar{\bm{L}}_{22}|}\,f(\bm{\ell}_{1},\bar{\bm{L}}_{22})\,f(\bar{\bm{L}}_{11},\bm{L}_{3}+\bm{\ell}_{2})\,+\,C^{\phi\phi}_{|\bm{\ell}_{1}+\bm{\ell}_{2}+\bm{L}_{3}|}\,f(\bm{\ell}_{1},\bm{L}_{3}+\bm{\ell}_{2})\,f(\bar{\bm{L}}_{11},\bar{\bm{L}}_{22})\Big{]}$
$\displaystyle\quad+\,(\bm{L}_{1}\leftrightarrow\bm{L}_{2})+(\bm{L}_{1}\leftrightarrow\bm{L}_{3}).$
(3.18)

This expression can be further simplified by a change of variables
$\bm{\ell}_{1}^{\prime}=\bm{L}_{1}-\bm{\ell}_{1}$ and
$\bm{\ell}_{2}^{\prime}=\bm{L}_{2}-\bm{\ell}_{2}$ (then renaming
$\bm{\ell}_{1}^{\prime}=\bm{\ell}_{1}$,
$\bm{\ell}_{2}^{\prime}=\bm{\ell}_{2}$) and applying the triangle condition
$\bm{L}_{1}+\bm{L}_{2}+\bm{L}_{3}=0$. Ultimately we obtain

$N^{(1)}_{\mathrm{B}}=N^{(1)}_{\mathrm{B}}(\bm{L}_{1},\bm{L}_{2},\bm{L}_{3})+N^{(1)}_{\mathrm{B}}(\bm{L}_{2},\bm{L}_{1},\bm{L}_{3})+N^{(1)}_{\mathrm{B}}(\bm{L}_{3},\bm{L}_{2},\bm{L}_{1})\,,$
(3.19)

with

$\displaystyle N^{(1)}_{\mathrm{B}}(\bm{L}_{1},\bm{L}_{2},\bm{L}_{3})=$
$\displaystyle
4\int^{\prime}_{\bm{\ell}_{1},\,\bm{\ell}_{2}}\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})F(\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{1}+\bm{\ell}_{2})\,$
$\displaystyle\qquad\qquad\qquad\times
C^{\widetilde{T}\widetilde{T}}_{|\bm{L}_{2}-\bm{\ell}_{2}|}\Big{[}C^{\phi\phi}_{L_{1}}\,f(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,f(\bm{\ell}_{2},-\bm{L}_{1}-\bm{\ell}_{2})\,$
$\displaystyle\qquad\qquad\qquad\;+\,2C^{\phi\phi}_{|\bm{\ell}_{1}+\bm{\ell}_{2}|}\,f(\bm{L}_{1}-\bm{\ell}_{1},-\bm{L}_{1}-\bm{\ell}_{2})\,f(\bm{\ell}_{1},\bm{\ell}_{2})\Big{]}\,\,.$
(3.20)

By making the change of variables
$\bm{\ell}^{\prime}_{1}=\bm{L}_{1}-\bm{\ell}_{1}$ and
$\bm{\ell}^{\prime}_{2}=-\bm{\ell}_{2}-\bm{L}_{1}$, it can be shown that
$N^{(1)}_{\mathrm{B}}(\bm{L}_{1},\bm{L}_{2},\bm{L}_{3})$ is symmetric under
exchange of $\bm{L}_{2}$ and $\bm{L}_{3}$ when the triangle condition
$\bm{L}_{1}+\bm{L}_{2}+\bm{L}_{3}=0$ is imposed. This symmetry ensures that
the right-hand side of Eq. (3.19) is symmetric under all permutations of
$\bm{L}_{1}$, $\bm{L}_{2}$ and $\bm{L}_{3}$ when the triangle condition holds.

###  3.2 Feynman rules

Similarly to the perturbative computation of scattering amplitudes that is
based on the wide application of the Wick’s theorem, one can introduce a set
of Feynman rules from which Feynman diagrams can be constructed, each
associated to a possible contribution to the $n$-pt functions in Eq. (3). We
will not give the explicit derivation of these rules, and we refer the reader
to the original publications [39, 40] for more details. In summary, this
method is applied as follows:

  * •

First, for each capital momentum $\bm{L}_{i}$ that appears in the $n$-pt
function (3) we draw a “rectangular” block, each including two vertexes that
correspond to the lensed temperatures. To it, we can associate the following
Feynman rule666Note that this Feynman rule is slightly different by what shown
in the original reference. This is because the lensing convergence field has a
different normalization convention than the lensing potential field.

$\bm{L_{i}-\ell}$$\bm{\ell}$$=A_{L_{i}}^{\kappa}\,F(\bm{\ell},\bm{L_{i}-\ell})$
(3.21)

where $F(\bm{\ell},\bm{L}_{i}-\bm{\ell})$ is the (optimal) filter in Eq.
(2.12). Therefore, when computing the 3-pt correlation function of
$\hat{\kappa}$, we have three external fundamental blocks, for a total of 6
circular vertexes. Each vertex corresponds to a lensed temperature mode and
the total number of vertexes coincides with the number of lensed temperature
modes in the expectation value of Eq. (3).

  * •

Coming out of each vertex, we can draw one straight line (that corresponds to
the unlensed temperature field) and an arbitrary number of wiggly lines
representing the order in the lensing potential field expansion. The
corresponding Feynman rule for the $n$-wiggly lines vertex is

$\bm{\ell}$$\bm{m}$$\bm{k_{1}}$$\bm{k_{2}}$$\bm{k_{n}}$$=\ \prod_{i}-(\bm{m\
\cdotp\ k_{i}})$ (3.22)

where $\bm{\ell}=\bm{m}+\sum_{i}\bm{k}_{i}$ due to the momentum conserving
delta-function at the vertex.

  * •

Next, we define the expression of the propagators, which correspond to
contractions between either straight or wiggly lines. Assuming Gaussian
distributed lensing potential and unlensed temperature fields, we get the
following Feynman rules

$=C^{TT}_{\ell}\,,$$=C^{\phi\phi}_{\ell}\,.$$\bm{\ell}$$\bm{\ell}$ (3.23)

  * •

In the presence of a non-zero bispectrum for the lensing potential, we get the
following additional Feynman rule

$\bm{\ell_{1}}$$\bm{\ell_{2}}$$\bm{\ell_{3}}$$=B_{\phi}(\ell_{1},\ell_{2},\ell_{3})$
(3.24)

where the momenta $\bm{\ell}_{i}$ are incoming.

  * •

The sum of the momenta coming inside each vertex must give zero due to
momentum conservation. As a result, not all the small momenta $\bm{\ell}_{j}$
are independent, but some of them turn out to be fixed in terms of the capital
momenta $\bm{L}_{i}$ and other small independent momenta $\bm{\ell}_{i}$. For
each of these one has to perform an integration
$\int\mathrm{d^{2}}\bm{\ell}_{i}/(2\pi)^{2}$. Outside the momenta integration,
a factor $(2\pi)^{2}\delta(\sum_{i}\bm{L}_{i})$ has to be added to the final
result.

With these fundamental Feynman rules, one can compute the $N^{(j/2)}$
contribution ($j\geq 0\in N$)777Notice that there is a $N^{(1/2)}$ term which
is proportional to the cross-correlation $C^{T\phi}_{\ell}$, but we ignore
such contributions given that both $C^{T\phi}_{\ell}/C^{TT}_{\ell}\ll 1$ and
$C^{T\phi}_{\ell}/C^{\phi\phi}_{\ell}\ll 1$. to the $n$-pt function of the
estimated $\hat{\kappa}$ by writing down all the possible diagrams with $j$
total wiggly lines coming off the lensed vertexes in the building blocks.
These rules are particularly convenient when dealing with the computation of
$j/2>1$ noise-biases, where one must include the $\mathcal{O}(\phi^{2})$ or
even $\mathcal{O}(\phi^{3})$ terms of the expansion (2.5), for which the Wick
expansion becomes challenging. Moreover, the Feynman diagrams allow us to
straightforwardly generalize the computations when including the CMB
polarization fields. We refer to Refs. [39, 40] for more details in this
regard.

Besides the fundamental Feynman rules, in Ref. [40] some more composite ones
have been derived. In the following, we report one which is useful for the
computation of $N_{\mathrm{B}}^{(3/2)}$ (see App. A). More specifically, block
(3.25) corresponds to all the possible ways in which we can connect two lensed
temperature vertexes with only 1 uncontracted wiggly line coming off one of
the two vertexes, diagrammatically

$\bm{\ell_{1}}$$\bm{\ell_{1}+\ell_{2}}$$\bm{\ell_{2}}$$=f(\bm{\ell_{1},\ell_{2}})$+$=$
(3.25)

(a) $2_{C}\times 2_{C}\times 2_{C}$

(b) $4_{C}\times 2_{C}$

(c) $3_{C}\times 3_{C}$

(d) $6_{C}$

Figure 1: Example of Feynman diagrams decomposed in the connected components.
Each block enclosed by a dotted line of a different color represents a
different connected component.

Finally, Feynman diagrams can be decomposed in connected contributions, which
are defined as subsets of the diagram where the connection between a set of
vertexes can be removed only by erasing one or more lines (either wiggly or
straight lines). As an example, in Fig. 1 we show some examples of Feynman
diagrams contributing to each of the terms appearing in the decomposition (3).

####  $N^{(0)}_{\mathrm{B}}$ noise-bias: Feynman approach.

The $N^{(0)}_{\mathrm{B}}$ noise-bias comes from the bispectrum of the
reconstructed $\hat{\kappa}$ field, thus the diagrams representing its
contribution comprise three building blocks (3.21). Moreover,
$N^{(0)}_{\mathrm{B}}$ arises from the connected 2-pt function of lensed
temperature anisotropies, therefore it features three propagators (3.23)
corresponding to the lensed temperature power spectra888Notice that the
propagator introduced in Sec. 3.2 corresponds to the unlensed temperature
power spectrum, however for our purposes of taking the non-perturbated term,
here we consider it as lensed.. Of the 15 diagrams, the ones that contain an
“interaction” of the type

$\bm{L_{i}-\ell_{i}}$$\bm{\ell_{i}}$ (3.26)

will not contribute, given that the sum of the momenta must be zero at each
vertex. The remaining 8 diagrams are equivalent because each block is
symmetrical for an exchange of momenta
$\bm{\ell}_{i}\leftrightarrow\bm{L}_{i}-\bm{\ell}_{i}$.

Figure 2: The only independent Feynman diagram contributing to the
$N_{\mathrm{B}}^{(0)}$ noise-bias of
$\langle\hat{\kappa}\hat{\kappa}\hat{\kappa}\rangle$.

Therefore, we are left with one diagram (for example, the one depicted in Fig.
2), which, according to the Feynman rules, corresponds to

$\displaystyle N^{(0)}_{\mathrm{B}}=8\int^{\prime}_{\bm{\ell}_{1}}\,$
$\displaystyle
F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{L}_{1}-\bm{\ell}_{1},\bm{L}_{3}+\bm{\ell}_{1})\,F(-\bm{\ell}_{1},\bm{L}_{3}+\bm{\ell}_{1})$
$\displaystyle\times
C^{\widetilde{T}\widetilde{T}}_{\ell_{1}}\,C^{\widetilde{T}\widetilde{T}}_{|\bm{L}_{1}-\bm{\ell}_{1}|}\,C^{\widetilde{T}\widetilde{T}}_{|\bm{L}_{3}+\bm{\ell}_{1}|}\,.$
(3.27)

where we used the relations $\bm{\ell}_{2}=\bm{L}_{1}-\bm{\ell}_{1}$ and
$\bm{\ell}_{3}=-\bm{\ell}_{1}$ coming from momentum conservation. This is
exactly what we found in Sec. 3.1.

####  $N^{(1)}_{\mathrm{B}}$ noise-bias: Feynman approach.

In order to compute $N^{(1)}_{\mathrm{B}}$ using the Feynman approach, we need
to add two wiggly lines coming off two vertexes corresponding to the $\phi$
fields involved.

Once we have defined the basic elements that constitute the diagrams, it is
possible to visualize the symmetries of the system. Indeed, as shown in Fig.
3, we have only 4 independent Feynman diagrams, as all other ones can be
obtained by exchanging the vertexes inside the building blocks or by
exchanging $\bm{L}_{i}$. The power of the Feynman approach lies in the fact
that we do not have to go through lengthy calculations to identify the
contributing terms.

(a) $N^{(1)}_{a}$

(b) $N^{(1)}_{b}$

(c) $N^{(1)}_{c}$

(d) $N^{(1)}_{d}$

Figure 3: Independent Feynman diagrams contributing to the
$N^{(1)}_{\mathrm{B}}$ noise-bias. The internal momenta are fixed by the
momenta-conserving delta function at each vertex, as explicitly shown in the
figure. Note that only two of the three $\ell_{i}$’s turn out to be
independent.

Using the Feynman rules outlined in Sec. 3.2, we can translate each of the
diagrams in Fig. 3 into the corresponding expressions (which account for
symmetry factors)

$\displaystyle N^{(1)}_{a}=8\,$
$\displaystyle\int^{\prime}_{\bm{\ell}_{1},\,\bm{\ell}_{2}}\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})F(\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{1}+\bm{\ell}_{2})$
$\displaystyle\qquad\times\,(\bm{\ell}_{1}\cdot\bm{L}_{1})\Big{[}\bm{L}_{1}\cdot(\bm{L}_{1}+\bm{\ell}_{2})\Big{]}C^{\phi\phi}_{L_{1}}\,C^{TT}_{|\bm{L}_{1}+\bm{\ell}_{2}|}C^{TT}_{|\bm{L}_{2}-\bm{\ell}_{2}|}C^{TT}_{\ell_{1}}+\mbox{perms($\bm{L}_{i}$)}\,,$
(3.28) $\displaystyle N^{(1)}_{b}=8\,$
$\displaystyle\int^{\prime}_{\bm{\ell}_{1},\,\bm{\ell}_{2}}\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})F(\bm{L}_{3}+\bm{\ell}_{2},-\bm{\ell}_{2})$
$\displaystyle\qquad\times\,\Big{[}(-\bm{L}_{2}+\bm{\ell}_{2})\cdot(\bm{L}_{3}+\bm{\ell}_{1}+\bm{\ell}_{2})\Big{]}\Big{[}(\bm{L}_{3}+\bm{\ell}_{1}+\bm{\ell}_{2})\cdot(\bm{L}_{3}+\bm{\ell}_{2})\Big{]}$
$\displaystyle\qquad\times\,C^{\phi\phi}_{|\bm{L}_{3}+\bm{\ell}_{1}+\bm{\ell}_{2}|}\,C^{TT}_{|\bm{L}_{3}+\bm{\ell}_{2}|}C^{TT}_{|\bm{L}_{2}-\bm{\ell}_{2}|}C^{TT}_{\ell_{2}}+(\bm{L}_{1}\leftrightarrow\bm{L}_{2})+(\bm{L}_{1}\leftrightarrow\bm{L}_{3})\,,$
(3.29) $\displaystyle N^{(1)}_{c}=-8\,$
$\displaystyle\int^{\prime}_{\bm{\ell}_{1},\,\bm{\ell}_{2}}\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})F(-\bm{\ell}_{1},\bm{L}_{3}+\bm{\ell}_{1})$
$\displaystyle\qquad\times\,\Big{[}(\bm{L}_{2}-\bm{\ell}_{2})\cdot(\bm{L}_{3}+\bm{\ell}_{1}+\bm{\ell}_{2})\Big{]}\Big{[}(\bm{L}_{3}+\bm{\ell}_{1}+\bm{\ell}_{2})\cdot(\bm{L}_{3}+\bm{\ell}_{1})\Big{]}$
$\displaystyle\qquad\times\,C^{\phi\phi}_{|\bm{\ell}_{1}+\bm{\ell}_{2}+\bm{L}_{3}|}\,C^{TT}_{|\bm{L}_{2}-\bm{\ell}_{2}|}C^{TT}_{|\bm{L}_{3}+\bm{\ell}_{1}|}C^{TT}_{\ell_{1}}+\mbox{perms($\bm{L}_{i}$)}\,,$
(3.30) $\displaystyle N^{(1)}_{d}=-8\,$
$\displaystyle\int^{\prime}_{\bm{\ell}_{1},\,\bm{\ell}_{3}}\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})F(\bm{\ell}_{3},\bm{L}_{3}-\bm{\ell}_{3})F(\bm{L}_{1}-\bm{\ell}_{1},\bm{L}_{3}+\bm{\ell}_{1})$
$\displaystyle\qquad\times\,\Big{[}(\bm{L}_{3}-\bm{\ell}_{3})\cdot(\bm{\ell}_{1}+\bm{\ell}_{3})\Big{]}\Big{[}(\bm{\ell}_{1}+\bm{\ell}_{3})\cdot\bm{\ell}_{3}\Big{]}$
$\displaystyle\qquad\times\,C^{\phi\phi}_{|\bm{\ell}_{1}+\bm{\ell}_{3}|}\,C^{TT}_{|\bm{L}_{3}-\bm{\ell}_{3}|}C^{TT}_{|\bm{L}_{1}-\bm{\ell}_{1}|}C^{TT}_{\ell_{3}}+(\bm{L}_{1}\leftrightarrow\bm{L}_{3})+(\bm{L}_{2}\leftrightarrow\bm{L}_{3})\,.$
(3.31)

The internal momenta, explicitly shown in Fig. 3, are fixed by momentum
conservation at the vertexes. We make a few change of variables to rearrange
the terms999In particular, $N^{(1)}_{a}$ stays the same. In $N^{(1)}_{b}$ we
made the change of variable $\bm{L}_{3}+\bm{\ell}_{2}=\bm{\ell}_{2}^{\prime}$.
In $N^{(1)}_{c}$ we made the double change of variable
$\bm{L}_{3}+\bm{\ell}_{1}=\bm{\ell}_{2}^{\prime}$,
$\,\bm{\ell}_{2}=\bm{\ell}_{1}^{\prime}$. In $N^{(1)}_{d}$ we made the double
change of variable $\bm{\ell}_{3}=\bm{\ell}_{1}^{\prime}$,
$\,\bm{\ell}_{1}=\bm{\ell}_{2}^{\prime}$., we can rewrite some of the
contributions as

$\displaystyle N^{(1)}_{b}=8\,$
$\displaystyle\int^{\prime}_{\bm{\ell}_{1},\,\bm{\ell}_{2}}\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})F(\bm{\ell}_{2},\bm{L}_{3}-\bm{\ell}_{2})F(\bm{L}_{3}-\bm{\ell}_{2},\bm{L}_{1}+\bm{\ell}_{2})$
$\displaystyle\times\,\Big{[}(\bm{L}_{1}+\bm{\ell}_{2})\cdot(\bm{\ell}_{1}+\bm{\ell}_{2})\Big{]}\Big{[}(\bm{\ell}_{1}+\bm{\ell}_{2})\cdot\bm{\ell}_{2}\Big{]}C^{\phi\phi}_{|\bm{\ell}_{1}+\bm{\ell}_{2}|}\,C^{TT}_{|\bm{L}_{1}+\bm{\ell}_{2}|}C^{TT}_{|\bm{L}_{3}-\bm{\ell}_{2}|}C^{TT}_{\ell_{2}}$
$\displaystyle\,+(\bm{L}_{1}\leftrightarrow\bm{L}_{2})+(\bm{L}_{1}\leftrightarrow\bm{L}_{3})\,,$
(3.32) $\displaystyle N^{(1)}_{c}=-8\,$
$\displaystyle\int^{\prime}_{\bm{\ell}_{1},\,\bm{\ell}_{2}}\,F(\bm{\ell}_{1},\bm{L}_{2}-\bm{\ell}_{1})F(\bm{\ell}_{2},\bm{L}_{3}-\bm{\ell}_{2})F(\bm{L}_{3}-\bm{\ell}_{2},\bm{L}_{2}+\bm{\ell}_{2})$
$\displaystyle\times\,\Big{[}(\bm{L}_{2}-\bm{\ell}_{1})\cdot(\bm{\ell}_{1}+\bm{\ell}_{2})\Big{]}\Big{[}(\bm{\ell}_{1}+\bm{\ell}_{2})\cdot\bm{\ell}_{2}\Big{]}\,C^{\phi\phi}_{|\bm{\ell}_{1}+\bm{\ell}_{2}|}\,C^{TT}_{|\bm{L}_{2}-\bm{\ell}_{1}|}C^{TT}_{|\bm{L}_{3}-\bm{\ell}_{2}|}C^{TT}_{\ell_{2}}$
$\displaystyle+\mbox{perms($\bm{L}_{i}$)}\,,$ (3.33) $\displaystyle
N^{(1)}_{d}=-8\,$
$\displaystyle\int^{\prime}_{\bm{\ell}_{1},\,\bm{\ell}_{2}}\,F(\bm{\ell}_{1},\bm{L}_{3}-\bm{\ell}_{1})F(\bm{\ell}_{2},\bm{L}_{1}-\bm{\ell}_{2})F(\bm{L}_{1}-\bm{\ell}_{2},\bm{L}_{3}+\bm{\ell}_{2})$
$\displaystyle\times\,\Big{[}(\bm{L}_{3}-\bm{\ell}_{1})\cdot(\bm{\ell}_{1}+\bm{\ell}_{2})\Big{]}\Big{[}(\bm{\ell}_{1}+\bm{\ell}_{2})\cdot\bm{\ell}_{2}\Big{]}C^{\phi\phi}_{|\bm{\ell}_{1}+\bm{\ell}_{2}|}\,C^{TT}_{|\bm{L}_{3}-\bm{\ell}_{1}|}C^{TT}_{|\bm{L}_{1}-\bm{\ell}_{2}|}C^{TT}_{\ell_{1}}$
$\displaystyle+(\bm{L}_{1}\leftrightarrow\bm{L}_{3})+(\bm{L}_{2}\leftrightarrow\bm{L}_{3})\,.$
(3.34)

It is now straightforward to verify that summing over all these contributions
leads to Eq. (3.1). This can be done by exploiting the definition in Eq. (2.9)
and expanding Eq. (3.1) in terms of scalar products between momenta and
temperature power spectra. Here we show this explicitly for one term as an
example. Consider the following contribution to Eq. (3.1)

$\begin{split}4\times\int_{\bm{\ell}_{1},\,\bm{\ell}_{2}}\,&F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})\,F(\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{1}+\bm{\ell}_{2})\\\
&\times\,C^{\phi\phi}_{L_{1}}\,f(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,f(\bm{\ell}_{2},-\bm{L}_{1}-\bm{\ell}_{2})\,C^{TT}_{|\bm{L}_{2}-\bm{\ell}_{2}|}\,.\end{split}$
(3.35)

By expanding the two response functions, and taking the product, we get

$\displaystyle
f(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,f(\bm{\ell}_{2},-\bm{L}_{1}-\bm{\ell}_{2})=$
$\displaystyle\Big{[}\bm{\ell}_{1}\cdot\bm{L}_{1}\Big{]}\Big{[}\bm{L}_{1}\cdot(\bm{L}_{1}+\bm{\ell}_{2})\Big{]}\,C^{TT}_{|\bm{L}_{1}+\bm{\ell}_{2}|}C^{TT}_{\ell_{1}}$
$\displaystyle+\Big{[}\bm{L}_{1}\cdot(\bm{L}_{1}-\bm{\ell}_{1})\Big{]}\Big{[}\bm{L}_{1}\cdot(\bm{L}_{1}+\bm{\ell}_{2})\Big{]}\,C^{TT}_{|\bm{L}_{1}+\bm{\ell}_{2}|}C^{TT}_{|\bm{L}_{1}-\bm{\ell}_{1}|}$
$\displaystyle-\Big{[}\bm{\ell}_{1}\cdot\bm{L}_{1}\Big{]}\Big{[}\bm{\ell}_{2}\cdot\bm{L}_{1}\Big{]}\,C^{TT}_{\ell_{2}}C^{TT}_{\ell_{1}}$
$\displaystyle-\Big{[}\bm{L}_{1}\cdot(\bm{L}_{1}-\bm{\ell}_{1})\Big{]}\Big{[}\bm{\ell}_{2}\cdot\bm{L}_{1}\Big{]}\,C^{TT}_{\ell_{2}}C^{TT}_{|\bm{L}_{1}-\bm{\ell}_{1}|}\,.$
(3.36)

If we replace Eq. (3.2) into Eq. (3.35), we get the following 4 terms

$\displaystyle\text{(A):}\quad
4\times\int_{\bm{\ell}_{1},\,\bm{\ell}_{2}}\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})\,F(\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{1}+\bm{\ell}_{2})$
$\displaystyle\qquad\;\quad\times\,\Big{[}\bm{\ell}_{1}\cdot\bm{L}_{1}\Big{]}\Big{[}\bm{L}_{1}\cdot(\bm{L}_{1}+\bm{\ell}_{2})\Big{]}\,C^{TT}_{|\bm{L}_{1}+\bm{\ell}_{2}|}\,C^{TT}_{\ell_{1}}\,C^{TT}_{|\bm{L}_{2}-\bm{\ell}_{2}|}\,C^{\phi\phi}_{L_{1}}\,,$
(3.37) $\displaystyle\text{(B):}\quad
4\times\int_{\bm{\ell}_{1},\,\bm{\ell}_{2}}\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})\,F(\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{1}+\bm{\ell}_{2})$
$\displaystyle\qquad\;\quad\times\,\Big{[}\bm{L}_{1}\cdot(\bm{L}_{1}-\bm{\ell}_{1})\Big{]}\Big{[}\bm{L}_{1}\cdot(\bm{L}_{1}+\bm{\ell}_{2})\Big{]}\,C^{TT}_{|\bm{L}_{1}+\bm{\ell}_{2}|}\,C^{TT}_{|\bm{L}_{1}-\bm{\ell}_{1}|}\,C^{TT}_{|\bm{L}_{2}-\bm{\ell}_{2}|}\,C^{\phi\phi}_{L_{1}}\,,$
(3.38)
$\displaystyle\text{(C):}\quad-4\times\int_{\bm{\ell}_{1},\,\bm{\ell}_{2}}\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})\,F(\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{1}+\bm{\ell}_{2})$
$\displaystyle\qquad\;\;\qquad\times\,\Big{[}\bm{\ell}_{1}\cdot\bm{L}_{1}\Big{]}\Big{[}\bm{\ell}_{2}\cdot\bm{L}_{1}\Big{]}\,C^{TT}_{\ell_{2}}\,C^{TT}_{\ell_{1}}\,C^{TT}_{|\bm{L}_{2}-\bm{\ell}_{2}|}\,C^{\phi\phi}_{L_{1}}\,,$
(3.39)
$\displaystyle\text{(D):}\quad-4\times\int_{\bm{\ell}_{1},\,\bm{\ell}_{2}}\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})\,F(\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{1}+\bm{\ell}_{2})$
$\displaystyle\qquad\;\;\qquad\times\,\Big{[}\bm{L}_{1}\cdot(\bm{L}_{1}-\bm{\ell}_{1})\Big{]}\Big{[}\bm{\ell}_{2}\cdot\bm{L}_{1}\Big{]}\,C^{TT}_{\ell_{2}}\,C^{TT}_{|\bm{L}_{1}-\bm{\ell}_{1}|}\,C^{TT}_{|\bm{L}_{2}-\bm{\ell}_{2}|}\,C^{\phi\phi}_{L_{1}}\,.$
(3.40)

By a change of variables $\bm{\ell}_{1}^{\prime}=\bm{L}_{1}-\bm{\ell}_{1}$ in
the integrals Eq. (3.38) and Eq. (3.40), they become equivalent to Eq. (3.37)
and Eq. (3.39), respectively. Therefore, in the end we have the two
independent contributions

$\begin{split}8\times\int_{\bm{\ell}_{1},\,\bm{\ell}_{2}}\,&F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})\,F(\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{1}+\bm{\ell}_{2})\\\
&\times\,\Big{[}\bm{\ell}_{1}\cdot\bm{L}_{1}\Big{]}\Big{[}\bm{L}_{1}\cdot(\bm{L}_{1}+\bm{\ell}_{2})\Big{]}\,C^{TT}_{|\bm{L}_{1}+\bm{\ell}_{2}|}\,C^{TT}_{\ell_{1}}\,C^{TT}_{|\bm{L}_{2}-\bm{\ell}_{2}|}\,C^{\phi\phi}_{L_{1}}\,,\end{split}$
(3.41)

and

$\begin{split}-8\times\int_{\bm{\ell}_{1},\,\bm{\ell}_{2}}\,&F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})\,F(\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{1}+\bm{\ell}_{2})\\\
&\times\,\Big{[}\bm{\ell}_{1}\cdot\bm{L}_{1}\Big{]}\Big{[}\bm{\ell}_{2}\cdot\bm{L}_{1}\Big{]}\,C^{TT}_{\ell_{2}}\,C^{TT}_{\ell_{1}}\,C^{TT}_{|\bm{L}_{2}-\bm{\ell}_{2}|}\,C^{\phi\phi}_{L_{1}}\,.\end{split}$
(3.42)

Now, it is evident that Eq. (3.41) gives exactly the first term of Eq. (3.2),
without accounting for permutations. When making an additional changes of
variables $\bm{\ell}_{2}^{\prime}=-\bm{\ell}_{2}-\bm{L}_{1}$, Eq. (3.42)
becomes

$\begin{split}8\times\int_{\bm{\ell}_{1},\,\bm{\ell}_{2}}\,&F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{\ell}_{2},\bm{L}_{3}-\bm{\ell}_{2})\,F(\bm{L}_{3}-\bm{\ell}_{2},\bm{L}_{1}+\bm{\ell}_{2})\\\
&\times\,\Big{[}\bm{\ell}_{1}\cdot\bm{L}_{1}\Big{]}\Big{[}\bm{L}_{1}\cdot(\bm{L}_{1}+\bm{\ell}_{2})\Big{]}\,C^{TT}_{|\bm{L}_{1}+\bm{\ell}_{2}|}\,C^{TT}_{\ell_{1}}\,C^{TT}_{|\bm{L}_{3}-\bm{\ell}_{2}|}\,C^{\phi\phi}_{L_{1}}\,,\end{split}$
(3.43)

which corresponds to the permutation $\bm{L}_{2}\leftrightarrow\bm{L}_{3}$ in
Eq. (3.2).

Once again, in order to obtain the non-perturbative noise-bias, we replace the
unlensed power spectra in the response function Eq. (2.9) and the ones
explicitly shown in the integrand respectively with the gradient lensed
$C_{\ell}^{\tilde{T}\nabla\tilde{T}}$ and the lensed
$C_{\ell}^{\widetilde{T}\widetilde{T}}$.

As seen in Sec. 3.1, in order to obtain $N^{(1)}_{\mathrm{B}}$ with the Wick
method, we need to compute all the possible contractions of the fields, which
involve a large number of permutations. On the other hand, here we show how
the relevant symmetries of the calculations become evident with the Feynman
approach, allowing us to write down all the relevant terms directly from the
diagrams. This makes the Feynman diagram method more efficient than the Wick
approach, particularly when the order of the $n$-pt function and of the noise-
bias $N^{(j/2)}$ increases.

##  4 Numerical results

In this section we present the results of evaluating the “Gaussian” noise-
biases $N_{\mathrm{B}}^{(0)}$ and $N_{\mathrm{B}}^{(1)}$ numerically. The
$N_{B}^{(3/2)}$ noise-bias appears if there is a non-zero matter bispectrum
inducing a lensing bispectrum; we calculate $N_{B}^{(3/2)}$ with the Feynman
formalism in App. A but leave its numerical evaluation for future work.

We make a number of simplifying assumptions as mentioned earlier: we consider
only CMB temperature modes, and we neglect the effect of primordial non-
Gaussianities and of the cross-correlations between temperature modes and the
lensing potential

$B_{T}(\bm{\ell}_{1},\bm{\ell}_{2},\bm{\ell}_{3})=0\,,\qquad
C^{T\phi}_{\ell}=0\,.$ (4.1)

We present the numerical evaluation of the two largest noise-biases,
$N_{\mathrm{B}}^{(0)}$ and $N_{\mathrm{B}}^{(1)}$ in two different slicings
(or configurations) of the bispectrum: the equilateral slicing, where
$L_{1}=L_{2}=L_{3}\equiv L$, and the folded one, given by $L_{1}\equiv L$ and
$L_{2}=L_{3}\equiv L/2$. We do not consider the squeezed limit for two main
reasons: the lensing bispectrum has low signal-to-noise in this configuration;
and it is significantly affected by multipole binning (see Ref. [43] for
details).

The numerical computation of these noise-biases, which include up to 4D
integrals over many multipoles, can be challenging to perform in a reasonable
amount of time. This problem was already highlighted in Ref. [34] for the
$N^{(3/2)}$ bias in lensing power spectrum reconstruction. Here, we tackle
this issue with a state-of-the-art numerical method, the VEGAS algorithm [44],
which gives Monte Carlo estimates of arbitrary multidimensional integrals by
using importance sampling. The algorithm locates those parts of the integrand
that contribute most to the final integral through a known or approximate
probability distribution function and then concentrates sampling in those
areas. This allows the computation of high-dimensional integrals with good
precision while keeping the computational cost low.

Similarly to the power spectrum reconstruction, $N_{\mathrm{B}}^{(0)}$ comes
from the $2_{C}\times 2_{C}\times 2_{C}$ of the connected 6-pt function (3),
and is given in Eq. (3.1). In Fig. 4, we show the magnitude of
$N^{(0)}_{\mathrm{B}}$ in the equilateral and folded slicings, and compare it
to the theoretically predicted lensing bispectrum induced by non-linear
evolution of the matter field and post-Born corrections [25]. For the former,
we employ the fitting formulae for the matter bispectrum given by Scoccimaro
and Couchman in Ref. [45], hereafter denoted by “SC”, and by Gil-Marín et al.
in Ref. [46], denoted below by “GM” 101010Recently, Ref. [47] provides a more
accurate fitting formula for the matter bispectrum calibrated by high-
resolution cosmological N-body simulations. Although such an accurate fitting
formula is desirable in practice, the difference between the new and SC/GM
fitting formula is not so large compared to the noise bias. Since we are
interested in the level of noise bias, this paper focuses on the SC and GM
formulas.. We find that the reconstructed convergence bispectrum is dominated
by $N_{\mathrm{B}}^{(0)}$ by several orders of magnitude, up to $10^{5}$ on
small scales for the equilateral slicing.

Figure 4: Left panel: Numerical evaluation of $N^{(0)}_{\mathrm{B}}$ in the
equilateral and folded slicings. The dashed line refers to negative values.
Right panel: Ratio between $N^{(0)}_{\mathrm{B}}$ and the theoretical
$\kappa$-bispectrum. Notice that the noise-bias in the folded slicing is
negative, while the signal bispectrum is positive over a large range of $L$
values.

The next-order noise-bias is $N^{(1)}_{\mathrm{B}}$. The expression in Eq.
(3.1) comprises two terms, one “separable”, where the bias term can be split
into a product of 2D integrals giving

$\displaystyle N^{(1)}_{\mathrm{B,\,sep}}=$ $\displaystyle
4C^{\phi\phi}_{L_{1}}A_{L_{1}}^{\kappa}\int_{\bm{\ell}_{1}}\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})f(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,$
$\displaystyle\qquad\times
A_{L_{2}}^{\kappa}A_{L_{3}}^{\kappa}\int_{\bm{\ell}_{2}}F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})F(\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{1}+\bm{\ell}_{2})f(\bm{\ell}_{2},-\bm{L}_{1}-\bm{\ell}_{2})\,C^{\widetilde{T}\widetilde{T}}_{|\bm{L}_{2}-\bm{\ell}_{2}|}$
$\displaystyle\qquad+(\bm{L}_{1}\leftrightarrow\bm{L}_{2})+(\bm{L}_{1}\leftrightarrow\bm{L}_{3})\,,$
(4.2)

and one “coupled”, where the full 4D integral has to be evaluated, giving

$\displaystyle N^{(1)}_{\mathrm{B,\,coupled}}=$ $\displaystyle
8\int^{\prime}_{\bm{\ell}_{1},\,\bm{\ell}_{2}}\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})\,F(\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{1}+\bm{\ell}_{2})$
$\displaystyle\quad\times
f(\bm{L}_{1}-\bm{\ell}_{1},-\bm{L}_{1}-\bm{\ell}_{2})f(\bm{\ell}_{1},\bm{\ell}_{2})\,C^{\phi\phi}_{|\bm{\ell}_{1}+\bm{\ell}_{2}|}C^{\widetilde{T}\widetilde{T}}_{|\bm{L}_{2}-\bm{\ell}_{2}|}+(\bm{L}_{1}\leftrightarrow\bm{L}_{2})+(\bm{L}_{1}\leftrightarrow\bm{L}_{3})\,.$
(4.3)

Note that the first integral in Eq. (4) is simply $1/A_{L_{1}}^{\kappa}$, by
the normalisation condition (2.13), and so cancels with the normalisation
$A_{L_{1}}^{\kappa}$.

Figure 5: Numerical evaluation of $N^{(1)}_{\mathrm{B,\,sep}}$ (dashed black)
and $N^{(1)}_{\mathrm{B,\,coupled}}$ (solid black) for the equilateral (left)
and folded (right) slicings. The total $N^{(1)}_{\mathrm{B}}$ is also shown in
each slicing (magenta for equilateral and green for folded). Note that the
noise-bias in the folded slicing is negative. Figure 6: As Fig. 4, but for
$N^{(1)}_{\mathrm{B}}$.

As shown in Fig. 5, for both slicings, the separable solution contributes the
most at low multipoles. Intuitively, at low values of $L$ the lensing
reconstruction weights couple tightly the small-scale CMB modes that dominate
the reconstruction and the separable term (4) allows a larger 4D integration
volume than the coupled term (4). On the other hand, at high $L$ the lensing
weights couple the CMB modes less tightly and the two terms become comparable.
Furthermore, in Fig. 6 we show that $N^{(1)}_{\mathrm{B}}$ is several orders
of magnitude larger than the $\kappa$-bispectrum and, like
$N^{(0)}_{\text{B}}$, needs to be carefully accounted for when reconstructing
the bispectrum.

It is interesting to note that both $N^{(0)}_{\mathrm{B}}$ and
$N^{(1)}_{\mathrm{B}}$ are negative in the folded slicing, while the signal
itself has positive values over most of the $L$ range (the negative ones are
due to Post-Born corrections). While in general, we must subtract any leading-
order noise-biases as they are much larger than the signal, in the folded
configuration this feature constitutes a consistency check of the
reconstruction process.

##  5 Validation with simulations

In this section, we use simulated temperature maps to validate our
calculations of the $N_{\mathrm{B}}^{(0)}$ and $N_{\mathrm{B}}^{(1)}$ noise-
biases. The simulation pipeline for each bias consists of:

  1. 1.

generating full-sky temperature maps;

  2. 2.

reconstructing $\hat{\kappa}$ from the above simulated temperature map with
the quadratic estimator; and

  3. 3.

computing the binned bispectrum of the resulting $\hat{\kappa}$.

We use the publicly available code cmblensplus [48] for the quadratic
estimator and binned bispectrum estimator, whose details are described in Ref.
[43].

In general, we bin the simulation results as it reduces the computational
costs, but we do not bin the analytical results as this would require costly
evaluation at many neighbouring triangle configurations. However, when binning
squeezed and folded bispectra, it is known that significant differences
between the binned and unbinned analytical results can occur. As shown in
[43], the binned bispectrum of a specific shape is less contaminated by
different triangle configurations if the number of multipole bins increases.
We discuss the binning dependency for the folded shape in Sec. 5.1 below.

Finally, in App. B we outline an alternative method for estimating these bias
terms that uses paired simulations to isolate the specific Wick contractions
that contribute to each term. This has the advantage of allowing each term in
the expansion in connected moments, Eq. (3) to be estimated individually.

###  5.1 Validation of $N^{(0)}_{\mathrm{B}}$

As discussed in Sec. 4, $N^{(0)}_{\mathrm{B}}$ arises from the $2_{C}\times
2_{C}\times 2_{C}$ term in Eq. (3). Therefore, in order to estimate such bias
from simulations, it suffices to generate Gaussian temperature maps using a
fiducial lensed temperature power spectrum computed with CAMB [49]. We use the
healpy package [50] for generating a full-sky temperature map with the lensed
temperature power spectrum. Consequently, the bispectrum of the reconstructed
convergence field, $\hat{\kappa}$, from each realization is proportional to
$N^{(0)}_{\mathrm{B}}$ only.

Figure 7: Numerical evaluation of $N^{(0)}_{\text{B}}$ (Eq. 3.1) for the
equilateral slicing (solid line), compared to binned results averaged over
$10^{4}$ Gaussian simulations. The lower panel shows the fractional residuals.
The error-bars are very small due to the large number of simulations. Figure
8: As Fig. 7, but for the folded slicing. Figure 9: Simulation results for
$N^{(0)}_{\text{B}}$ in the folded slicing computed with the binned bispectrum
estimator with increasing numbers of bins. In each case, 800 simulations are
used. The lower panel shows the residual plot. It is clear that finer binning
improves the match with the (unbinned) analytical results.

In Figs. 7 and 8 we show a comparison between simulations and the numerical
evaluation of the analytical result in Eq. (3.1) for two slicings –
equilateral and folded. We average over $10^{4}$ realizations to reduce the
simulation error in each data point, and we show the fractional residual
defined as

$\text{residual}=\frac{(\text{simulation}-\text{analytical})}{\text{analytical}}.$
(5.1)

We find good agreement between the analytical flat-sky result and the full-
sky, binned simulation results in the equilateral shape. We believe that the
small residuals, around $2\,\%$, are a consequence of the flat-sky
approximation employed in the analytical results. The folded slicing is
additionally affected by binning, as discussed above. In Fig. 9 we explicitly
show how increasing the number of bins in the bispectrum estimate improves the
match between the simulations and our analytical calculations. Finally, a
similar result is obtained using the paired-simulation approach (see App.
B.1).

###  5.2 Validation of $N^{(1)}_{\mathrm{B}}$

In order to test $N^{(1)}_{\mathrm{B}}$, which explicitly depends on the
lensing potential power spectrum, we generate unlensed CMB maps and
subsequently lens them with a Gaussian realization of $\phi$. For this step,
we employ the publicly available code
lenspyx111111https://github.com/carronj/lenspyx. In principle, the resulting
bispectrum contains all the “Gaussian” noise-biases121212The signal
$B_{\kappa}$ is zero given that we lens with a Gaussian $\phi$, as is
$N^{(3/2)}_{\text{B}}$.

$B_{\hat{\kappa}}=N^{(0)}_{\mathrm{B}}+N^{(1)}_{\mathrm{B}}+N^{(2)}_{\mathrm{B}}+\cdots\,.$
(5.2)

However, the results that we present below (see also App. B) suggest that
$N^{(2)}_{\text{B}}$ (and higher-order noise-biases) are small compared to
$N^{(1)}_{\text{B}}$. In this subsection, therefore, we show the residual,
$B_{\hat{\kappa}}-N^{(0)}_{\mathrm{B}}$, and denote it as
$N^{(1)}_{\mathrm{B}}$.

Figure 10: Numerical evaluation of $N^{(1)}_{\text{B}}$ from Eq. (3.1) for the
equilateral slicing compared to binned simulation results averaged over
$10^{4}$ realizations. In the lower panel, we show the fractional residuals.
Figure 11: As Fig. 10, but for the folded slicing.

Generally, there are many ways to subtract $N^{(0)}_{\mathrm{B}}$. For an
ideal survey, such as the one presented here, it would suffice to compute
$N^{(0)}_{\mathrm{B}}$ with a lensed power spectrum that agrees well with
reality (such as a smoothed version of the empirical power spectrum of our
sky). However, given the discrepancy between flat-sky analytical results and
the full-sky simulations as well as the binning effect discussed in Sec. 5.1,
we subtract a simulation-based131313Specifically, from each realization of the
lensed CMB, we generate Gaussian CMB maps drawn from the empirical angular
power spectrum measured from that simulation, and we then estimate
$N^{(0)}_{\mathrm{B}}$ from these Gaussian maps as done in Sec. 5.1. Note that
a small bias is introduced in this estimate of $N^{(0)}_{\mathrm{B}}$ since
its expectation value involves
$\langle\hat{C}^{TT}_{\ell_{1}}\hat{C}^{TT}_{\ell_{2}}\hat{C}^{TT}_{\ell_{3}}\rangle$,
the average of three products of the empirical power spectrum. This differs
from the true value,
$\widetilde{C}^{TT}_{\ell_{1}}\widetilde{C}^{TT}_{\ell_{2}}\widetilde{C}^{TT}_{\ell_{3}}$,
because of cosmic variance:
$\langle\hat{C}^{TT}_{\ell_{1}}\hat{C}^{TT}_{\ell_{2}}\hat{C}^{TT}_{\ell_{3}}\rangle=\widetilde{C}^{TT}_{\ell_{1}}\widetilde{C}^{TT}_{\ell_{2}}\widetilde{C}^{TT}_{\ell_{3}}\left[1+\left(\frac{2}{2\ell_{2}+1}\delta_{\ell_{2}\ell_{3}}+(\text{2
perms.})\right)+\frac{8}{(2\ell_{1}+1)^{2}}\delta_{\ell_{1}\ell_{2}\ell_{3}}\right]\,,$
where the last term enforces $\ell_{1}=\ell_{2}=\ell_{3}$. However, the
difference is suppressed by factors that are $\mathcal{O}(\ell_{i})$, which
are small for the CMB multipoles that dominate the lensing reconstruction, and
by the restricted integration volume from requiring two or three of the
$\ell_{i}$ to be equal. The benefit of subtracting $N^{(0)}_{\text{B}}$ with
this approach is that it is unaffected by any difference between the ensemble-
average of the angular power spectra of the lensed CMB maps and the fiducial
lensed power spectrum, which might arise from small inaccuracies in simulating
the lensing operation.  estimate of $N^{(0)}_{\mathrm{B}}$.

In Figs. 10 and 11, we show the comparison between the evaluation of Eq. (3.1)
and the binned simulation results (averaged over $10^{4}$ realizations) in the
equilateral and folded configurations, finding that the analytical curves
generally show good agreement with the simulations for both configurations.
The equilateral slicing has oscillatory features at high $L$ that cannot be
easily captured by the binned bispectrum, and we need a large number of
realizations, or an increase in the number of bins, to capture these features
accurately. Moreover, the error-bars on the binned simulation result are about
an order of magnitude larger than those of the folded configuration. This is
likely due to the fact that the equilateral $N^{(0)}_{\mathrm{B}}$ is larger
than the folded one by roughly an order of magnitude, hence when subtracting a
simulation-based $N^{(0)}_{\mathrm{B}}$ any variance in it from the finite
number of simulations used will leave additional variance in
$N^{(1)}_{\mathrm{B}}$. In order to reduce the error in the equilateral
slicing, we would need around 100 times more simulations.

In the folded configuration, we find a small mismatch, around $5\,\%$, which
likely mostly originates from the binning effect as we discussed in Sec. 5.1
and/or a small contribution of $N^{(2)}_{\mathrm{B}}$. The good agreement
shown in Figs. 10 and 11 suggests that $N^{(2)}_{\mathrm{B}}$ is much smaller
than $N^{(0)}_{\mathrm{B}}$ and $N^{(1)}_{\mathrm{B}}$, although further work
is needed to assess whether it is negligible compared to the bispectrum
signal. In App. B, we reach similar conclusions using the alternative, paired-
simulation approach. We show that the $4_{C}\times 2_{C}$ term is well
described by our analytical expression for $N^{(1)}_{\text{B}}$ and that the
$6_{C}$ term (which partly sources $N^{(2)}_{\text{B}}$) is significantly
smaller.

##  6 Conclusions and future perspectives

In the past decade, weak gravitational lensing of the CMB has been established
as a robust probe of late-time cosmological evolution. Data from upcoming CMB
experiments are expected to be sensitive to the nonlinear clustering of matter
at low redshift, which effectively induces non-Gaussianity in the lensing
convergence field. This signal, if detected, provides additional information
about structure growth, and has the potential to break degeneracies between
cosmological parameters.

In this paper, we calculated for the first time the noise-biases that arise
when estimating the lensing convergence bispectrum from CMB temperature
anisotropies via the Hu & Okamoto quadratic estimator. We provided a
theoretical model of the leading noise-biases using two methods: the standard
Wick contraction of correlation functions; and the Feynman diagram formalism
applied to CMB lensing. Regarding the latter approach, we adapted the results
of Refs. [39, 40] to the bispectrum case. The diagrammatic approach allows one
to organise very effectively the calculation of the noise-biases up to
$j$-orders of the lensing potential $\phi$ and is particularly helpful for
higher-order biases.

We further numerically evaluated the leading noise-biases,
$N^{(0)}_{\mathrm{B}}$ and $N^{(1)}_{\mathrm{B}}$, in the equilateral and
folded configurations of the bispectrum. We showed that they represent the
largest contribution to the reconstructed bispectrum and dominate over the
signal by several orders of magnitude at small scales ($L\gg 100$) for a
temperature-based reconstruction. We finally tested our theoretical results
against simulations of temperature maps, finding excellent agreement.

The general scope of this work was to give a first quantification of the
reconstruction noise-biases of the CMB lensing bispectrum in view of a
potential detection with ground-based experiments. Therefore, we leave a
number of generalizations and tests for future work:

  1. 1.

We have limited our analysis to temperature anisotropies only, however
polarization will be increasingly important for forthcoming surveys. We should
therefore extend our analysis to quantify the noise-biases for reconstructions
combining temperature and polarization anisotropies. A promising feature of
polarization-based bispectrum reconstruction is that for the $EB$ quadratic
estimator, which will dominate in the CMB-S4 era [23], the
$N^{(0)}_{\text{B}}$ noise-bias vanishes by parity.

  2. 2.

While our simulation results suggest that the higher-order
$N^{(2)}_{\text{B}}$ bias is much smaller than $N^{(0)}_{\text{B}}$ and
$N^{(1)}_{\text{B}}$, further work is needed to determine whether
$N^{(2)}_{\text{B}}$ is negligible compared to the signal bispectrum. More
converged simulation-based estimates and an analytic understanding of
$N^{(2)}_{\text{B}}$ will both be helpful in this regard.

  3. 3.

While we have provided a theoretical model of all the terms that contribute to
$N^{(3/2)}_{\mathrm{B}}$, a numerical evaluation of this noise-bias is needed
to assess its contribution to the reconstructed bispectrum. Given that it
contains the lensing convergence bispectrum, the analysis will depend on the
matter bispectrum modelling. Future work should also compare the numerical
results against estimates of the same non-Gaussian bias from ray-traced
$N$-body simulations.

  4. 4.

Our analysis ignores the effects of the curl-modes in the deflection angle
introduced by post-Born corrections, therefore it would be interesting to
assess how they affect our results.

  5. 5.

In the case of the reconstruction of the lensing power spectrum from data, the
biases are commonly removed with realization-dependent methods that mix data
and simulations. These make the subtraction more robust to mismodelling in the
simulations and reduce covariance between the reconstructed power at different
multipoles. The realization-dependent estimators for a general $n$-pt
correlation function has been developed in Ref. [51], but future work should
test the bispectrum estimator for specific experiments, such as SO. Given the
size of the $N^{(0)}_{\text{B}}$ and $N^{(1)}_{\text{B}}$ biases, achieving
highly converged bias corrections will be very important to avoid significant
additional scatter in the reconstructed bispectrum.

  6. 6.

It will also be important to quantify the impact on bispectrum reconstruction
of residuals of extragalactic foregrounds in maps of the temperature
anisotropies. As for the reconstructed power spectrum, polarization-based
estimators will be significantly less affected by these foreground signals.

Given the potential of the lensing bispectrum as a probe of cosmological
evolution, it will be critical to explore all of the work mentioned above.

##### Acknowledgments

We are grateful to Vanessa Böhm, Julien Carron, William Coulton, Mathew
Madhavacheril, Joel Meyers, and Alexander van Engelen for fruitful
discussions. A.K., G.O. and P.D.M acknowledge support from the Netherlands
organization for scientific research (NWO) VIDI grant (dossier 639.042.730).
A.K. also acknowledges the Simons Foundation for support in the latest stage
of this work. A.C. acknowledges support from the STFC (grant numbers
ST/N000927/1 and ST/S000623/1). T.N. acknowledges support from JSPS KAKENHI
Grant No. JP20H05859 and No. JP22K03682. The Kavli IPMU is supported by World
Premier International Research Center Initiative (WPI Initiative), MEXT,
Japan.

##  Appendix A Detailed computation of $N_{B}^{(3/2)}$

The noise-bias $N_{B}^{(3/2)}$ receives contribution from connected functions
both linear and non-linear in $\phi$. In the following we evaluate such terms
separately, using the Feynman diagrams.

#### Contributions from $4_{C}\times 2_{C}$ (linear order):

At linear order, there are 4 independent Feynman diagrams that correspond to

$\langle\delta T\delta T\delta
TT\rangle_{C}\,\langle\widetilde{T}\widetilde{T}\rangle_{C}\,,$ (A.1)

and are depicted in Fig. 12. These diagrams correspond to

(a) $N^{(3/2)}_{{\rm lin},a}$

(b) $N^{(3/2)}_{{\rm lin},b}$

(c) $N^{(3/2)}_{{\rm lin},c}$

(d) $N^{(3/2)}_{{\rm lin},d}$

Figure 12: Independent Feynman diagrams contributing to the
$N^{(3/2)}_{\mathrm{B}}$ noise-bias of
$\langle\hat{\kappa}\hat{\kappa}\hat{\kappa}\rangle$ with no lensed
temperature vertexes expanded beyond the linear order (2.6), $4_{C}\times
2_{C}$ terms.

$\begin{split}N_{{\rm
lin},a}^{(3/2)}=4\int^{\prime}_{\bm{\ell}_{1},\bm{\ell}_{2}}&\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})\,F(-\bm{\ell}_{2},\bm{L}_{3}+\bm{\ell}_{2})\,C^{\widetilde{T}\widetilde{T}}_{\ell_{2}}\\\
&\times\mathcal{F}(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1},\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{3}+\bm{\ell}_{2})+(\bm{L}_{1}\leftrightarrow\bm{L}_{2})+(\bm{L}_{1}\leftrightarrow\bm{L}_{3})\,,\end{split}$
(A.2) $\begin{split}N_{{\rm
lin},b}^{(3/2)}=8\int^{\prime}_{\bm{\ell}_{1},\bm{\ell}_{2}}&\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})\,F(-\bm{\ell}_{2},\bm{L}_{3}+\bm{\ell}_{2})\,C^{\widetilde{T}\widetilde{T}}_{\ell_{2}}\\\
&\times\mathcal{F}(\bm{\ell}_{1},\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{1}-\bm{\ell}_{1},\bm{L}_{3}+\bm{\ell}_{2})+(\bm{L}_{1}\leftrightarrow\bm{L}_{2})+(\bm{L}_{1}\leftrightarrow\bm{L}_{3})\,,\end{split}$
(A.3) $\begin{split}N_{{\rm
lin},c}^{(3/2)}=8\int^{\prime}_{\bm{\ell}_{1},\bm{\ell}_{2}}&\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})\,F(-\bm{\ell}_{2},\bm{L}_{3}+\bm{\ell}_{2})\,C^{\widetilde{T}\widetilde{T}}_{\ell_{2}}\\\
&\times\mathcal{F}(\bm{\ell}_{1},\bm{L}_{3}+\bm{\ell}_{2},\bm{L}_{1}-\bm{\ell}_{1},\bm{L}_{2}-\bm{\ell}_{2})+(\bm{L}_{1}\leftrightarrow\bm{L}_{2})+(\bm{L}_{1}\leftrightarrow\bm{L}_{3})\,,\end{split}$
(A.4) $\begin{split}N_{{\rm
lin},d}^{(3/2)}=4\int^{\prime}_{\bm{\ell}_{1},\bm{\ell}_{2}}&\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})\,F(-\bm{\ell}_{2},\bm{L}_{3}+\bm{\ell}_{2})\,C^{\widetilde{T}\widetilde{T}}_{\ell_{2}}\\\
&\times\mathcal{F}(\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{3}+\bm{\ell}_{2},\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})+(\bm{L}_{1}\leftrightarrow\bm{L}_{2})+(\bm{L}_{1}\leftrightarrow\bm{L}_{3})\,,\end{split}$
(A.5)

where we defined
$\mathcal{F}(\bm{\ell}_{1},\bm{\ell}_{2},\bm{\ell}_{3},\bm{\ell}_{4})=f(\bm{\ell}_{1},\bm{\ell}_{2})I(\bm{\ell}_{3},\bm{\ell}_{4})$,
with

$\begin{split}I(\bm{\ell}_{i},\bm{\ell}_{j})=\int_{\bm{p}}\left[\bm{p}\cdot(\bm{\ell}_{i}-\bm{p})\right]\,\left[\bm{p}\cdot(\bm{\ell}_{j}+\bm{p})\right]\,B_{\phi}(|\bm{\ell}_{i}-\bm{p}|,|\bm{\ell}_{j}+\bm{p}|,|\bm{\ell}_{i}+\bm{\ell}_{j}|)\,C^{\widetilde{T}\widetilde{T}}_{\bm{p}}\,.\end{split}$
(A.6)

Therefore, the final contribution from the $4_{C}\times 2_{C}$ terms reads as

$N^{(3/2)}_{B,\rm lin}=N^{(3/2)}_{B,\rm
lin}(\bm{L}_{1},\bm{L}_{2},\bm{L}_{3})+N^{(3/2)}_{B,\rm
lin}(\bm{L}_{2},\bm{L}_{1},\bm{L}_{3})+N^{(3/2)}_{B,\rm
lin}(\bm{L}_{3},\bm{L}_{1},\bm{L}_{2})\,,$ (A.7)

where

$\displaystyle N_{B,\rm lin}^{(3/2)}$
$\displaystyle(\bm{L}_{1},\bm{L}_{2},\bm{L}_{3})=4\int^{\prime}_{\bm{\ell}_{1},\bm{\ell}_{2}}\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})\,F(-\bm{\ell}_{2},\bm{L}_{3}+\bm{\ell}_{2})\,C^{\widetilde{T}\widetilde{T}}_{\ell_{2}}$
$\displaystyle\qquad\qquad\times\Big{[}\mathcal{F}(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1},\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{3}+\bm{\ell}_{2})+2\mathcal{F}(\bm{\ell}_{1},\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{1}-\bm{\ell}_{1},\bm{L}_{3}+\bm{\ell}_{2})$
$\displaystyle\qquad\qquad+2\mathcal{F}(\bm{\ell}_{1},\bm{L}_{3}+\bm{\ell}_{2},\bm{L}_{1}-\bm{\ell}_{1},\bm{L}_{2}-\bm{\ell}_{2})+\mathcal{F}(\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{3}+\bm{\ell}_{2},\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\Big{]}\,.$
(A.8)

#### Contributions from $6_{C}$ (linear order):

At linear order, the contribution to $N_{\mathrm{B}}^{3/2}$ comes from

$\langle\delta T\delta T\delta TTTT\rangle_{C}$ (A.9)

In terms of Feynman diagrams, depicted in Fig. 13, we have

$\displaystyle N_{{\rm lin},e}^{(3/2)}=$ $\displaystyle
2\int^{\prime}_{\bm{\ell}_{1},\bm{\ell}_{2},\bm{\ell}_{3}}F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})F(\bm{\ell}_{3},\bm{L}_{3}-\bm{\ell}_{3})$
$\displaystyle\times\,f(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,f(\bm{\ell}_{2},\bm{\ell}_{3})\,f(\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{3}-\bm{\ell}_{3})$
$\displaystyle\times
B_{\phi}(L_{1},|\bm{\ell}_{2}+\bm{\ell}_{3}|,|\bm{L}_{2}-\bm{\ell}_{2}+\bm{L}_{3}-\bm{\ell}_{3}|)+(\bm{L}_{1}\leftrightarrow\bm{L}_{2})+(\bm{L}_{1}\leftrightarrow\bm{L}_{3})\,,$
(A.10) $\displaystyle N_{{\rm lin},f}^{(3/2)}=$ $\displaystyle
8\int^{\prime}_{\bm{\ell}_{1},\bm{\ell}_{2},\bm{\ell}_{3}}F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})F(\bm{\ell}_{3},\bm{L}_{3}-\bm{\ell}_{3})\,f(\bm{\ell}_{1},\bm{\ell}_{2})\,f(\bm{L}_{1}-\bm{\ell}_{1},\bm{\ell}_{3})$
$\displaystyle\times\,f(\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{3}-\bm{\ell}_{3})B_{\phi}(|\bm{\ell}_{1}+\bm{\ell}_{2}|,|\bm{L}_{1}-\bm{\ell}_{1}+\bm{\ell}_{3}|,|\bm{L}_{2}-\bm{\ell}_{2}+\bm{L}_{3}-\bm{\ell}_{3}|).$
(A.11)

The expansion of $6_{C}$ function at linear order also includes the bispectrum
signal $B_{\kappa}$ that we aim to reconstruct, and it is as shown Fig. 14.

(a) $N^{(3/2)}_{{\rm lin},e}$

(b) $N^{(3/2)}_{{\rm lin},f}$

Figure 13: Independent Feynman diagrams contributing to the
$N^{(3/2)}_{\mathrm{B}}$ noise-bias of
$\langle\hat{\kappa}\hat{\kappa}\hat{\kappa}\rangle$ with no lensed
temperature vertexes expanded beyond the linear order (2.6), $6_{C}$ terms.
Figure 14: Feynman diagram corresponding to the bispectrum signal
$B_{\kappa}$.

#### Contributions from the second order expansion of the lensed temperature:

(a) $N^{(3/2)}_{a}$

(b) $N^{(3/2)}_{b}$

(c) $N^{(3/2)}_{c}$

(d) $N^{(3/2)}_{d}$

(e) $N^{(3/2)}_{e}$

(f) $N^{(3/2)}_{f}$

(g) $N^{(3/2)}_{g}$

(h) $N^{(3/2)}_{h}$

Figure 15: Independent Feynman diagrams contributing to the
$N_{\mathrm{B}}^{(3/2)}$ noise-bias of
$\langle\hat{\kappa}\hat{\kappa}\hat{\kappa}\rangle$ with at least one lensed
temperature vertex expanded beyond the linear order (2.6).

The estimation of $N_{\mathrm{B}}^{(3/2)}$ includes also non-zero
contributions from the $\mathcal{O}(\phi^{2})$ and $\mathcal{O}(\phi^{3})$
terms in (2.7). In particular, the independent Feynman diagrams providing
these additional contributions are depicted in Fig. 15 and correspond to

$\begin{split}\text{(a):}&\quad\langle\delta^{3}TT\rangle_{C}\,\langle\widetilde{T}\widetilde{T}\rangle_{C}\,\langle\widetilde{T}\widetilde{T}\rangle_{C}\,,\\\
\text{(d):}&\quad\langle\delta^{2}T\delta
T\rangle_{C}\,\langle\widetilde{T}\widetilde{T}\rangle_{C}\,\langle\widetilde{T}\widetilde{T}\rangle_{C}\,,\\\
\text{(b, c, e, f, g, h):}&\quad\langle\delta^{2}T\delta
TTT\rangle_{C}\,\langle\widetilde{T}\widetilde{T}\rangle_{C}\,.\end{split}$
(A.12)

In the following, by exploiting the Feynman rules in Sec. 3.2, we give their
analytical expression including symmetry factors and permutations over the
capital momenta $\bm{L}_{i}$’s.

$\displaystyle
N^{(3/2)}_{a}=8\int^{\prime}_{\bm{\ell}_{1},\bm{\ell}_{2},\bm{\ell}_{3}}$
$\displaystyle\,F(\bm{\ell}_{1}+\bm{L}_{1},-\bm{\ell}_{1})\,F(\bm{\ell}_{1},\bm{L}_{2}-\bm{\ell}_{1})\,F(\bm{\ell}_{1}-\bm{L}_{2},\bm{L}_{1}-\bm{\ell}_{1})$
$\displaystyle\times\Big{[}(\bm{\ell}_{2}+\bm{\ell}_{3})\cdot\bm{\ell}_{1}\Big{]}(\bm{\ell}_{1}\cdot\bm{\ell}_{2})\,(\bm{\ell}_{3}\cdot\bm{\ell}_{1})\,B_{\phi}(|\bm{\ell}_{2}+\bm{\ell}_{3}|,\ell_{2},\ell_{3})$
$\displaystyle\times
C^{\widetilde{T}\widetilde{T}}_{|\bm{L}_{1}+\bm{\ell}_{1}|}\,C^{\widetilde{T}\widetilde{T}}_{\ell_{1}}\,C^{\widetilde{T}\widetilde{T}}_{|\bm{L}_{2}-\bm{\ell}_{1}|}+\mbox{perms($\bm{L}_{i}$)}\,.$
(A.13)

By definition our noise-bias has to be symmetric under the exchange of the
$\bm{L}_{i}$’s. By making the change of variable
$\bm{\ell}_{1}=-\bm{\ell}_{1}^{\prime}$ in the previous term, we easily show
that

$N^{(3/2)}_{a}=-N^{(3/2)}_{a}=0\,.$ (A.14)

The result is that Feynman diagrams with three wiggly lines coming out from a
single vertex sum up to 0 and do not contribute to the $N_{B}^{(3/2)}$ noise-
bias. Therefore, only Feynman diagrams with two wiggly lines coming out from a
given vertex contribute to higher orders in the lensing expansion. Their
contribution is given in the following

$\displaystyle
N^{(3/2)}_{b}=8\int^{\prime}_{\bm{\ell}_{1},\bm{\ell}_{2},\bm{\ell}_{3}}$
$\displaystyle\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{\ell}_{2},\bm{L}_{2}-\bm{\ell}_{2})\,F(\bm{L}_{2}-\bm{\ell}_{2},\bm{L}_{1}+\bm{\ell}_{2})$
$\displaystyle\times\Big{[}(\bm{\ell}_{2}+\bm{L}_{1})\cdot(\bm{\ell}_{2}-\bm{\ell}_{1}+\bm{L}_{1})\Big{]}\Big{[}\bm{\ell}_{2}\cdot\bm{\ell}_{3}\Big{]}\Big{[}\bm{\ell}_{2}\cdot(\bm{\ell}_{2}-\bm{\ell}_{1}-\bm{\ell}_{3}+\bm{L}_{1})\Big{]}$
$\displaystyle\times
B_{\phi}(|\bm{\ell}_{2}-\bm{\ell}_{1}-\bm{\ell}_{3}+\bm{L}_{1}|,\ell_{3},|\bm{\ell}_{2}-\bm{\ell}_{1}+\bm{L}_{1}|)$
$\displaystyle\times
C^{\widetilde{T}\widetilde{T}}_{|\bm{L}_{1}+\bm{\ell}_{2}|}\,C^{\widetilde{T}\widetilde{T}}_{\ell_{2}}\,C^{\widetilde{T}\widetilde{T}}_{|\bm{L}_{2}-\bm{\ell}_{2}|}+\mbox{perms($\bm{L}_{i}$)}\,,$
(A.15) $\displaystyle
N^{(3/2)}_{d}=8\int^{\prime}_{\bm{\ell}_{1},\bm{\ell}_{2},\bm{\ell}_{3}}$
$\displaystyle\,F(-\bm{\ell}_{1}-\bm{L}_{3},\bm{\ell}_{1}-\bm{L}_{2})\,F(\bm{\ell}_{1},\bm{L}_{2}-\bm{\ell}_{1})\,F(-\bm{\ell}_{1},\bm{L}_{3}+\bm{\ell}_{1})$
$\displaystyle\times\Big{[}\bm{\ell}_{3}\cdot(\bm{L}_{2}-\bm{\ell}_{1}-\bm{\ell}_{3})\Big{]}\Big{[}\bm{\ell}_{3}\cdot\bm{\ell}_{2}\Big{]}\Big{[}\bm{\ell}_{3}\cdot(\bm{\ell}_{1}+\bm{\ell}_{3}-\bm{\ell}_{2}-\bm{L}_{2})\Big{]}$
$\displaystyle\times
B_{\phi}(|\bm{\ell}_{1}+\bm{\ell}_{3}-\bm{\ell}_{2}-\bm{L}_{2}|,\ell_{2},|\bm{L}_{2}-\bm{\ell}_{1}-\bm{\ell}_{3}|)$
$\displaystyle\times
C^{\widetilde{T}\widetilde{T}}_{|\bm{L}_{3}+\bm{\ell}_{1}|}\,C^{\widetilde{T}\widetilde{T}}_{\ell_{1}}\,C^{\widetilde{T}\widetilde{T}}_{\ell_{3}}+\mbox{perms($\bm{L}_{i}$)}\,,$
(A.16) $\displaystyle
N^{(3/2)}_{e}=8\int^{\prime}_{\bm{\ell}_{1},\bm{\ell}_{2},\bm{\ell}_{3}}$
$\displaystyle\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{L}_{3}-\bm{\ell}_{3},\bm{\ell}_{3}+\bm{L}_{1})\,F(\bm{\ell}_{3},\bm{L}_{3}-\bm{\ell}_{3})$
$\displaystyle\times\Big{[}\bm{\ell}_{1}\cdot(\bm{\ell}_{1}-\bm{\ell}_{3}-\bm{L}_{1})\Big{]}\Big{[}\bm{\ell}_{3}\cdot\bm{\ell}_{2}\Big{]}\Big{[}\bm{\ell}_{3}\cdot(\bm{L}_{1}-\bm{\ell}_{1}+\bm{\ell}_{3}-\bm{\ell}_{2})\Big{]}$
$\displaystyle\times
B_{\phi}(|\bm{L}_{1}-\bm{\ell}_{1}+\bm{\ell}_{3}-\bm{\ell}_{2}|,\ell_{2},|\bm{\ell}_{1}-\bm{\ell}_{3}-\bm{L}_{1}|)$
$\displaystyle\times
C^{\widetilde{T}\widetilde{T}}_{|\bm{L}_{3}-\bm{\ell}_{3}|}\,C^{\widetilde{T}\widetilde{T}}_{\ell_{1}}\,C^{\widetilde{T}\widetilde{T}}_{\ell_{3}}+\mbox{perms($\bm{L}_{i}$)}\,,$
(A.17) $\displaystyle
N^{(3/2)}_{f}=8\int^{\prime}_{\bm{\ell}_{1},\bm{\ell}_{2},\bm{\ell}_{3}}$
$\displaystyle\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(-\bm{\ell}_{1},\bm{\ell}_{1}-\bm{L}_{1}-\bm{L}_{3})\,F(\bm{\ell}_{3},\bm{L}_{3}-\bm{\ell}_{3})$
$\displaystyle\times\Big{[}(\bm{L}_{3}-\bm{\ell}_{3})\cdot(\bm{\ell}_{1}-\bm{\ell}_{3}-\bm{L}_{1})\Big{]}\Big{[}\bm{\ell}_{3}\cdot\bm{\ell}_{2}\Big{]}\Big{[}\bm{\ell}_{3}\cdot(\bm{L}_{1}-\bm{\ell}_{1}+\bm{\ell}_{3}-\bm{\ell}_{2})\Big{]}$
$\displaystyle\times
B_{\phi}(|\bm{L}_{1}-\bm{\ell}_{1}+\bm{\ell}_{3}-\bm{\ell}_{2}|,\ell_{2},|\bm{\ell}_{1}-\bm{\ell}_{3}-\bm{L}_{1}|)$
$\displaystyle\times
C^{\widetilde{T}\widetilde{T}}_{|\bm{L}_{3}-\bm{\ell}_{3}|}\,C^{\widetilde{T}\widetilde{T}}_{\ell_{1}}\,C^{\widetilde{T}\widetilde{T}}_{\ell_{3}}+\mbox{perms($\bm{L}_{i}$)}\,,$
(A.18) $\displaystyle
N^{(3/2)}_{g}=8\int^{\prime}_{\bm{\ell}_{1},\bm{\ell}_{2},\bm{\ell}_{3}}$
$\displaystyle\,F(\bm{\ell}_{1},\bm{L}_{1}-\bm{\ell}_{1})\,F(\bm{L}_{3}-\bm{\ell}_{3},\bm{L}_{1}+\bm{\ell}_{3})\,F(\bm{\ell}_{3},\bm{L}_{3}-\bm{\ell}_{3})$
$\displaystyle\times\Big{[}\bm{\ell}_{1}\cdot(\bm{\ell}_{2}-\bm{L}_{1})\Big{]}\Big{[}\bm{\ell}_{1}\cdot\bm{\ell}_{2}\Big{]}\Big{[}\bm{\ell}_{3}\cdot\bm{L}_{1}\Big{]}B_{\phi}(|\bm{L}_{1}-\bm{\ell}_{2}|,\ell_{2},|-\bm{L}_{1}|)$
$\displaystyle\times
C^{\widetilde{T}\widetilde{T}}_{|\bm{L}_{3}-\bm{\ell}_{3}|}\,C^{\widetilde{T}\widetilde{T}}_{\ell_{1}}\,C^{\widetilde{T}\widetilde{T}}_{\ell_{3}}+\mbox{perms($\bm{L}_{i}$)}\,,$
(A.19) $\displaystyle
N^{(3/2)}_{h}=8\int^{\prime}_{\bm{\ell}_{1},\bm{\ell}_{2},\bm{\ell}_{3}}$
$\displaystyle\,F(\bm{L}_{3}-\bm{\ell}_{3},\bm{L}_{2}+\bm{\ell}_{3})\,F(\bm{\ell}_{1},\bm{L}_{2}-\bm{\ell}_{1})\,F(\bm{\ell}_{3},\bm{L}_{3}-\bm{\ell}_{3})$
$\displaystyle\times\Big{[}\bm{\ell}_{3}\cdot\bm{\ell}_{2}\Big{]}\Big{[}\bm{\ell}_{3}\cdot(-\bm{L}_{2}-\bm{\ell}_{2})\Big{]}\Big{[}\bm{\ell}_{1}\cdot\bm{L}_{2}\Big{]}B_{\phi}(|-\bm{L}_{2}-\bm{\ell}_{2}|,\ell_{2},\ell_{2})$
$\displaystyle\times
C^{\widetilde{T}\widetilde{T}}_{|\bm{L}_{3}-\bm{\ell}_{3}|}\,C^{\widetilde{T}\widetilde{T}}_{\ell_{1}}\,C^{\widetilde{T}\widetilde{T}}_{\ell_{3}}+\mbox{perms($\bm{L}_{i}$)}\,.$
(A.20)

The contribution from $N^{(3/2)}_{c}$ turns out to be zero. By summing all
these terms to those in the previous subsections we get the full
$N^{(3/2)}_{B}$ noise-bias of Eq. (3.5).

##  Appendix B Alternative (paired) simulations

In this appendix, we describe an alternative form of simulations that allow us
to estimate the $N^{(0)}_{\text{B}}$ noise-bias and the $4_{C}\times 2_{C}$
bias term in the cumulant expansion (3) of the reconstructed lensing
bispectrum. These simulations allow us to target individual terms in the
expansion separately, and even particular Wick contractions of the unlensed
CMB that contribute to each term if required. The simulations here use
exclusively non-Gaussisan lensed CMB maps, which has the advantage over the
method used in Sec. 5 that there is no reliance on the power spectra of the
lensed simulations matching the theory lensed power spectrum to very high
accuracy on all scales.

As we shall need to identify explicitly the two temperature fields that enter
the quadratic estimator141414We shall refer to these as the two “legs” of the
quadratic estimator. for $\kappa(\bm{L})$, we adopt the notation
$\hat{\kappa}_{\bm{L}}(i,j)$ for the quadratic estimator applied to
temperature fields $i$ and $j$. We use a symmetrised form of the quadratic
estimator, as in Eq. (2.14), so that
$\hat{\kappa}_{\bm{L}}(i,j)=\hat{\kappa}_{\bm{L}}(j,i)$ for any temperature
fields $i$ and $j$.

In this appendix, the codes to generate the lensed temperature maps and to
compute $\hat{\kappa}$ and the binned bispectrum are the same as those
described in Sec. 4. Note that we include instrumental noise in the CMB maps
here. The noise model used is that relevant for the Simons Observatory large-
aperture telescope with a white noise level of $10\,\mu\text{K-arcmin}$ and a
beam size (full-width at half-maximum) of $1.4\,\text{arcmin}$ [22].

###  B.1 $N_{\mathrm{B}}^{(0)}$ bias

The $N_{\mathrm{B}}^{(0)}$ bias originates from the $2_{C}\times 2_{C}\times
2_{C}$ term in Eq. (3). For the $\bm{L}_{i}\neq 0$, the only contractions that
contribute couple fields across different quadratic estimators, giving eight
possible permutations. These are all equal because of the symmetry of the
quadratic estimator, so we have, symbolically,

$N^{(0)}_{\text{B}}\sim
8\langle\hat{\kappa}_{\bm{L}_{1}}(I,J)\hat{\kappa}_{\bm{L}_{2}}(K,I)\hat{\kappa}_{\bm{L}_{3}}(J,K)\rangle_{IJK}\,,$
(B.1)

where the expectation value is over independent realizations $I$, $J$, $K$ of
the noisy, lensed CMB temperature anisotropies. Note how the independence of
$I$, $J$ and $K$ ensures that only 2-pt correlations arise, despite the fields
being non-Gaussian.

We estimate Eq. (B.1) by averaging over sets of three independent simulations
of the noisy and lensed CMB. Specifically, we use 624 simulations and take
$J=I+1$ and $K=I+2$ (cyclically). Fig. 16 shows the comparison between this
method, using (B.1), and the numerical evaluation of $N_{\mathrm{B}}^{(0)}$
using Eq. (3.1), for the equilateral configuration. The two are in good
agreement. As seen in Sec. 5.1, the percent-level differences in Fig. 16 could
be due to the flat-sky approximation used in the analytical results.

Figure 16: Comparison of $N_{\mathrm{B}}^{(0)}$ for the equilateral
configuration calculated numerically from Eq. (3.1) and estimated from (B.1)
with 624 noisy and lensed simulations.

###  B.2 Bias from $4_{C}\times 2_{C}$ terms

We now estimate the bispectrum noise-bias from the $4_{C}\times 2_{C}$ term in
Eq. (3). This term contributes to $N^{(1)}_{\text{B}}$, but also to
$N^{(2)}_{\text{B}}$ (and higher-order noise-biases) for the case of Gaussian
lenses $\phi$. In Sec. 5.2, we compared our leading-order analytic prediction
for the $4_{C}\times 2_{C}$ bias (computed with the non-perturbative response
functions), Eq. (3.1), with the lensing bispectrum computed directly from
lensed CMB simulations with appropriate subtraction of $N^{(0)}_{\text{B}}$.
That comparison includes higher-order contributions to $4_{C}\times 2_{C}$,
along with those from the $6_{C}$ term. Here, we instead use suitably paired
simulations to isolate the $4_{C}\times 2_{C}$ term (at all orders in $\phi$),
removing contamination from the $6_{C}$ term in our comparison with the
leading-order anayltic result.

There are three independent permutations contributing to $4_{C}\times 2_{C}$,
which we label $A$, $B$ and $C$, respectively. These correspond to whether the
$2_{C}$ contraction is between fields with one appearing in
$\hat{\kappa}_{\bm{L}_{1}}$ and the other in $\hat{\kappa}_{\bm{L}_{2}}$, or
the $\bm{L}_{1}$–$\bm{L}_{3}$ or $\bm{L}_{2}$–$\bm{L}_{3}$ pairings.
Symbolically, we have

$4_{C}\times 2_{C}\sim
4\underbrace{\langle\widetilde{T}_{11}\widetilde{T}_{22}\widetilde{T}_{3}\widetilde{T}_{33}\rangle_{C}\langle\widetilde{T}_{1}\widetilde{T}_{2}\rangle_{C}}_{\text{perm.
$A$}}+4\underbrace{\langle\widetilde{T}_{11}\widetilde{T}_{22}\widetilde{T}_{2}\widetilde{T}_{33}\rangle_{C}\langle\widetilde{T}_{1}\widetilde{T}_{3}\rangle}_{\text{perm.
$B$}}+4\underbrace{\langle\widetilde{T}_{1}\widetilde{T}_{11}\widetilde{T}_{22}\widetilde{T}_{33}\rangle_{C}\langle\widetilde{T}_{2}\widetilde{T}_{3}\rangle_{C}}_{\text{perm.
$C$}}\,.$ (B.2)

The factors of $4$ account for whether the first and first, first and second,
second and first, or second and second legs of each quadratic estimator are
contracted across (they are all equal by symmetry of the quadratic estimator).
For each of the permutations $A$, $B$ and $C$, there are three possible Wick
contractions over the unlensed CMB in the connected 4-pt function. For
example, for $C$ we have

$\langle\widetilde{T}_{1}\widetilde{T}_{11}\widetilde{T}_{22}\widetilde{T}_{33}\rangle_{C}=\wick[offset=1.5em]{\langle\c{1}{\widetilde{T}}_{1}\c{1}{\widetilde{T}}_{11}\c{1}{\widetilde{T}}_{22}\c{1}{\widetilde{T}}_{33}\rangle_{\phi,C}}+\wick[offset=1.5em]{\langle\c{1}{\widetilde{T}}_{1}\c{2}{\widetilde{T}}_{11}\c{1}{\widetilde{T}}_{22}\c{2}{\widetilde{T}}_{33}\rangle_{\phi,C}}+\wick[offset=1.5em]{\langle\c{1}{\widetilde{T}}_{1}\c{2}{\widetilde{T}}_{11}\c{2}{\widetilde{T}}_{22}\c{1}{\widetilde{T}}_{33}\rangle_{\phi,C}}\,,$
(B.3)

where the contractions over the unlensed CMB fields that appear (linearly) in
each $\tilde{T}$ are shown explicitly. The second and third terms are equal
because of the symmetry of $\hat{\kappa}_{\bm{L}_{1}}$.

To extract the different contractions of the unlensed CMB fields in Eq. (B.3),
we need simulations with independent realizations of the unlensed temperature
field but the same realization of the lensing potential. This ensures that the
unlensed temperature fields are correctly paired, but the same $\phi$ fields
are being averaged in all lensed CMB fields. We must, however, remove any
contractions of $\phi$ that would give only disconnected contributions to the
4-pt function. We therefore estimate, for example, the contribution of the
first term in Eq. (B.3) to permutation $C$ as

$\int^{\prime}_{\bm{\ell}_{1},\bm{\ell}_{2},\bm{\ell}_{3}}F(\bm{\ell}_{1},\bar{\bm{L}}_{11})F(\bm{\ell}_{2},\bar{\bm{L}}_{22})F(\bm{\ell}_{3},\bar{\bm{L}}_{33})\wick[offset=1.5em]{\langle\c{1}{\widetilde{T}}_{1}\c{1}{\widetilde{T}}_{11}\c{1}{\widetilde{T}}_{22}\c{1}{\widetilde{T}}_{33}\rangle_{\phi,C}}\langle\widetilde{T}_{2}\widetilde{T}_{3}\rangle_{C}\\\
=\langle\hat{\kappa}_{\bm{L}_{1}}(ii,ii)\hat{\kappa}_{\bm{L}_{2}}(Jj,ki)\hat{\kappa}_{\bm{L}_{3}}(Jj,ki)\rangle_{Jj,ii,ki}-\langle\hat{\kappa}_{\bm{L}_{1}}(ii,ii)\hat{\kappa}_{\bm{L}_{2}}(Jj,kk)\hat{\kappa}_{\bm{L}_{3}}(Jj,kk)\rangle_{Jj,ii,kk}\,.$
(B.4)

Here, the first index in a field pair labels the realization of the unlensed
temperature field entering $\widetilde{T}$ and the second the realization of
the lensing potential. Capital indices indicate that instrumental noise is
added to the lensed field. Hence, the first leg of the quadratic estimator
$\hat{\kappa}_{\bm{L}_{2}}(Jj,ki)$ is a noisy lensed temperature field and the
second leg a noiseless lensed temperature field constructed from an
independent unlensed temperature realization and lensed by an independent
realization of the lensing potential. Since the noise model is Gaussian, the
simulations within the connected 4-pt correlation function can be noiseless.
Note that in the case where there is no lensing, the right-hand side of Eq.
(B.4) vanishes in every realization (not just in the mean), reducing the
scatter in our simulation-based estimates. The final expression for the noise-
bias to the lensing convergence bispectrum originating from the $4_{C}\times
2_{C}$ term is, symbolically,

$\displaystyle 4_{C}\times 2_{C}$ $\displaystyle\sim
4\langle\hat{\kappa}_{\bm{L}_{1}}(ii,ii)\hat{\kappa}_{\bm{L}_{2}}(Jj,ki)\hat{\kappa}_{\bm{L}_{3}}(Jj,ki)\rangle_{Jj,ii,ki}$
$\displaystyle\mbox{}-4\langle\hat{\kappa}_{\bm{L}_{1}}(ii,ii)\hat{\kappa}_{\bm{L}_{2}}(Jj,kk)\hat{\kappa}_{\bm{L}_{3}}(Jj,kk)\rangle_{Jj,ii,kk}$
$\displaystyle\mbox{}+8\hat{\kappa}_{\bm{L}_{1}}(ii,ki)\hat{\kappa}_{\bm{L}_{2}}(Jj,ii)\hat{\kappa}_{\bm{L}_{3}}(Jj,ki)\rangle_{Jj,ii,ki}$
$\displaystyle\mbox{}-8\langle\hat{\kappa}_{\bm{L}_{1}}(ii,kk)\hat{\kappa}_{\bm{L}_{2}}(Jj,ii)\hat{\kappa}_{\bm{L}_{3}}(Jj,kk)\rangle_{Jj,ii,kk}+(\bm{L}_{1}\leftrightarrow\bm{L}_{2})+(\bm{L}_{1}\leftrightarrow\bm{L}_{3})\,.$
(B.5)

Figure 17: Contribution from $4_{C}\times 2_{C}$ terms to the noise-bias on
the convergence bispectrum in the equilateral slicing, estimated from Eq.
(B.5) with 624 paired simulations (data points). The solid line shows the
flat-sky analytical result for $N^{(1)}_{\text{B}}$, Eq. (3.1), which agrees
well with the data points to within the simulation errors (estimated from the
scatter across simulations). The bottom panel shows the fractional residuals.

In Figs. 17 and 18 we show results of using 624 paired simulations to estimate
the contribution from $4_{C}\times 2_{C}$ terms to the noise-bias on the
bispectrum in the equilateral and folded configurations. These estimates are
compared to the flat-sky analytical expression for $N^{(1)}_{\text{B}}$, Eq.
(3.1). The results are in good agreement and support the findings from the
alternative simulation approach described in Sec. 5.2 that the bispectrum
noise-bias that originates from $4_{C}\times 2_{C}$ terms is well described by
the leading-order approximation, $N^{(1)}_{\text{B}}$. Note that the folded
slicing suffers from a binning effect as discussed in Sec. 5.1, which is
expected to become less significant as the number of bins used in the
bispectrum estimator is increased.

Figure 18: As Fig. 17, but for the folded configuration of the bispectrum.

###  B.3 Bias from $6_{C}$ term

Finally, we estimate the bias from the connected 6-pt function in the case of
Gaussian $\phi$. This term contributes to $N^{(2)}_{\text{B}}$ (and higher-
order biases), along with the $4_{C}\times 2_{C}$ terms. We estimate the
$6_{C}$ noise-bias to the bispectrum with151515An alternative estimator built
from the possible Wick contractions of the unlensed CMB fields, analogous to
Eq. (B.5), can easily be constructed and would likely have less scatter than
the simpler form adopted here. However, the latter is adequate for our initial
estimate of the relative importance of the bias from $6_{C}$.

$6_{C}\sim\langle\hat{\kappa}_{\bm{L}_{1}}(I,I)\hat{\kappa}_{\bm{L}_{2}}(I,I)\hat{\kappa}_{\bm{L}_{3}}(I,I)\rangle_{I}-4_{C}\times
2_{C}-N^{(0)}_{\text{B}}\,.$ (B.6)

We evaluate the $4_{C}\times 2_{C}$ term with Eq. (B.5) and
$N^{(0)}_{\text{B}}$ from Eq. (B.1). Gaussian instrumental noise does not
affect the connected 6-pt function of the noisy CMB so Eq. (B.6) could also be
evaluated without including noise in any of the simulations used for any of
the terms. Since the lensing potential $\phi$ is Gaussian there is no
$N_{\mathrm{B}}^{(3/2)}$ bias.

Fig. 19 compares the estimated $6_{C}$ term from Eq. (B.6) to the size of the
$N^{(1)}_{\text{B}}$ bias in the equilateral slicing. Although our estimate is
very noisy, it suggests that $6_{C}$ is around an order of magnitude smaller
than $N^{(1)}_{\text{B}}$ for $L<1000$. However, given the size of
$N^{(1)}_{\text{B}}$ compared to the bispectrum signal, we cannot rule out the
$6_{C}$ term being comparable to the signal.

Figure 19: Comparison of the $N_{\mathrm{B}}^{(1)}$ bias calculated
numerically from Eq. (3.1) and the bias from $6_{C}$ estimated with Eq. (B.6)
in the equilateral slicing using 20 multipole bins. The bias from $6_{C}$ is
less significant by around an order of magnitude than the
$N_{\mathrm{B}}^{(1)}$ bias for $L<1000$.

## References

  * [1] A. Lewis and A. Challinor, “Weak gravitational lensing of the CMB”, Physics Reports 429 no. 1, (Jun, 2006) 1–65, arXiv:astro-ph/0601594.
  * [2] M. Kamionkowski and E. D. Kovetz, “The Quest for B Modes from Inflationary Gravitational Waves”, Annual Review of Astronomy and Astrophysics 54 no. 1, (Sep, 2016) 227–269, arXiv:1510.06042.
  * [3] L. Knox and Y.-S. Song, “Limit on the Detectability of the Energy Scale of Inflation”, Physical Review Letters 89 no. 1, (Jun, 2002) , arXiv:astro-ph/0202286.
  * [4] U. Seljak and C. M. Hirata, “Gravitational lensing as a contaminant of the gravity wave signal in the CMB”, Physical Review D 69 no. 4, (Feb, 2004) , arXiv:astro-ph/0310163.
  * [5] K. M. Smith, O. Zahn, and O. Doré, “Detection of gravitational lensing in the cosmic microwave background”, Physical Review D 76 no. 4, (Aug, 2007) .
  * [6] S. Das et al., “Detection of the Power Spectrum of Cosmic Microwave Background Lensing by the Atacama Cosmology Telescope”, Physical Review Letters 107 no. 2, (Jul, 2011) , arXiv:1103.2124.
  * [7] S. Das et al., “The Atacama Cosmology Telescope: temperature and gravitational lensing power spectrum measurements from three seasons of data”, Journal of Cosmology and Astroparticle Physics 2014 no. 04, (Apr, 2014) 014–014, arXiv:1301.1037.
  * [8] P. A. R. Ade et al., “Planck2013 results. XVII. Gravitational lensing by large-scale structure”, Astronomy & Astrophysics 571 (Oct, 2014) A17, arXiv:1303.5077.
  * [9] P. A. R. Ade et al., “Measurement of the Cosmic Microwave Background Polarization Lensing Power Spectrum with the POLARBEAR Experiment”, Physical Review Letters 113 no. 2, (Jul, 2014) , arXiv:1312.6646.
  * [10] A. van Engelen et al., “A measurement of gravitational lensing of the Microwave Background using the South Pole Telescope data”, The Astrophysical Journal 756 no. 2, (Aug, 2012) 142, arXiv:1202.0546.
  * [11] K. T. Story et al., “A Measurement of the Cosmic Microwave Background Gravitational Lensing Potential from 100 Square Degrees of SPTpol Data”, The Astrophysical Journal 810 no. 1, (Aug, 2015) 50, arXiv:1412.4760.
  * [12] P. A. R. Ade et al., “BICEP2 / Keck Array VIII: Measurement of gravitational lensing from large-scale B-mode polarization”, Astrophys. J. 833 no. 2, (2016) 228, arXiv:1606.01968 [astro-ph.CO].
  * [13] N. Aghanim et al., “Planck2018 results. VIII. Gravitational lensing”, Astronomy & Astrophysics 641 (Sep, 2020) A8, arXiv:1807.06210.
  * [14] B. D. Sherwin et al., “Evidence for Dark Energy from the Cosmic Microwave Background Alone Using the Atacama Cosmology Telescope Lensing Measurements”, Physical Review Letters 107 no. 2, (Jul, 2011) , arXiv:1105.0419.
  * [15] M. Kaplinghat, L. Knox, and Y.-S. Song, “Determining Neutrino Mass from the Cosmic Microwave Background Alone”, Physical Review Letters 91 no. 24, (Dec, 2003) , arXiv:astro-ph/0303344.
  * [16] R. A. Battye and A. Moss, “Evidence for Massive Neutrinos from Cosmic Microwave Background and Lensing Observations”, Physical Review Letters 112 no. 5, (Feb, 2014) , arXiv:1308.5870.
  * [17] U. Seljak, “Extracting Primordial Non-Gaussianity without Cosmic Variance”, Physical Review Letters 102 no. 2, (Jan, 2009) , arXiv:0807.1770.
  * [18] T. Giannantonio and W. J. Percival, “Using correlations between cosmic microwave background lensing and large-scale structure to measure primordial non-Gaussianity”, Monthly Notices of the Royal Astronomical Society: Letters 441 no. 1, (Mar, 2014) L16–L20, arXiv:1312.5154.
  * [19] C. Modi, M. White, and Z. Vlah, “Modeling CMB lensing cross correlations with CLEFT”, Journal of Cosmology and Astroparticle Physics 2017 no. 08, (Aug, 2017) 009–009, arXiv:1706.03173.
  * [20] A. Banerjee, B. Jain, N. Dalal, and J. Shelton, “Tests of neutrino and dark radiation models from galaxy and CMB surveys”, Journal of Cosmology and Astroparticle Physics 2018 no. 01, (Jan, 2018) 022–022, arXiv:1612.07126.
  * [21] M. Schmittfull and U. Seljak, “Parameter constraints from cross-correlation of CMB lensing with galaxy clustering”, Physical Review D 97 no. 12, (Jun, 2018) , arXiv:1710.09465.
  * [22] P. Ade et al., “The Simons Observatory: science goals and forecasts”, Journal of Cosmology and Astroparticle Physics 2019 no. 02, (Feb, 2019) 056–056, arXiv:1808.07445.
  * [23] K. N. Abazajian et al., “CMB-S4 Science Book, First Edition”, arXiv:1610.02743.
  * [24] T. Namikawa, “CMB lensing bispectrum from nonlinear growth of the large scale structure”, Physical Review D 93 no. 12, (Jun, 2016) , arXiv:1604.08578.
  * [25] G. Pratten and A. Lewis, “Impact of post-Born lensing on the CMB”, Journal of Cosmology and Astroparticle Physics 2016 no. 08, (Aug, 2016) 047–047.
  * [26] J. Liu, J. C. Hill, B. D. Sherwin, A. Petri, V. Böhm, and Z. Haiman, “CMB lensing beyond the power spectrum: Cosmological constraints from the one-point probability distribution function and peak counts”, Phys. Rev. D 94 no. 10, (2016) 103501, arXiv:1608.03169 [astro-ph.CO].
  * [27] E. Calabrese et al., “Precision epoch of reionization studies with next-generation CMB experiments”, Journal of Cosmology and Astroparticle Physics 2014 no. 08, (Aug, 2014) 010–010, arXiv:1406.4794.
  * [28] A. Suzuki et al., “The Polarbear-2 and the Simons Array Experiments”, Journal of Low Temperature Physics 184 no. 3–4, (Jan, 2016) 805–810, arXiv:1512.07299.
  * [29] W. Hu and T. Okamoto, “Mass Reconstruction with Cosmic Microwave Background Polarization”, The Astrophysical Journal 574 no. 2, (Aug, 2002) 566–574.
  * [30] T. Okamoto and W. Hu, “Cosmic microwave background lensing reconstruction on the full sky”, Physical Review D 67 no. 8, (Apr, 2003) , arXiv:astro-ph/0301031.
  * [31] M. Kesden, A. Cooray, and M. Kamionkowski, “Lensing reconstruction with CMB temperature and polarization”, Physical Review D 67 no. 12, (Jun, 2003) .
  * [32] D. Hanson, A. Challinor, G. Efstathiou, and P. Bielewicz, “CMB temperature lensing power reconstruction”, Physical Review D 83 no. 4, (Feb, 2011) , arXiv:1008.4403.
  * [33] T. Namikawa, D. Hanson, and R. Takahashi, “Bias-hardened CMB lensing”, Monthly Notices of the Royal Astronomical Society 431 no. 1, (Feb, 2013) 609–620, arXiv:1209.0091.
  * [34] V. Böhm, M. Schmittfull, and B. D. Sherwin, “Bias to CMB lensing measurements from the bispectrum of large-scale structure”, Physical Review D 94 no. 4, (Aug, 2016) , arXiv:1605.01392.
  * [35] V. Böhm, B. D. Sherwin, J. Liu, J. C. Hill, M. Schmittfull, and T. Namikawa, “Effect of non-Gaussian lensing deflections on CMB lensing measurements”, Phys. Rev. D 98 no. 12, (2018) 123510, arXiv:1806.01157 [astro-ph.CO].
  * [36] G. Fabbian, A. Lewis, and D. Beck, “CMB lensing reconstruction biases in cross-correlation with large-scale structure probes”, JCAP 10 (2019) 057, arXiv:1906.08760 [astro-ph.CO].
  * [37] R. K. Sachs and A. M. Wolfe, “Perturbations of a Cosmological Model and Angular Variations of the Microwave Background”, Astrophysical Journal 147 (Jan, 1967) 73.
  * [38] D. M. Goldberg and D. N. Spergel, “Microwave background bispectrum. II. A probe of the low redshift universe”, Physical Review D 59 no. 10, (Apr, 1999) .
  * [39] E. E. Jenkins, A. V. Manohar, W. J. Waalewijn, and A. P. S. Yadav, “Gravitational Lensing of the CMB: a Feynman Diagram Approach”, Phys. Lett. B 736 (2014) 6–10, arXiv:1403.2386.
  * [40] E. E. Jenkins, A. V. Manohar, W. J. Waalewijn, and A. P. S. Yadav, “Higher-Order Gravitational Lensing Reconstruction using Feynman Diagrams”, JCAP 09 (2014) 024, arXiv:1403.4607.
  * [41] The Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters”, arXiv:1807.06209.
  * [42] A. Lewis, A. Challinor, and D. Hanson, “The shape of the CMB lensing bispectrum”, Journal of Cosmology and Astroparticle Physics 2011 no. 03, (Mar, 2011) 018–018.
  * [43] T. Namikawa, B. Bose, F. m. c. R. Bouchet, R. Takahashi, and A. Taruya, “CMB lensing bispectrum: Assessing analytical predictions against full-sky lensing simulations”, Phys. Rev. D 99 (Mar, 2019) 063511, arXiv:1812.10635.
  * [44] G. P. Lepage, “Adaptive multidimensional integration: vegas enhanced”, Journal of Computational Physics 439 (Aug, 2021) 110386, arXiv:2009.05112.
  * [45] R. Scoccimarro and H. M. P. Couchman, “A fitting formula for the non-linear evolution of the bispectrum”, Monthly Notices of the Royal Astronomical Society 325 no. 4, (Aug, 2001) 1312–1316, arXiv:astro-ph/0009427.
  * [46] H. Gil-Marín, C. Wagner, F. Fragkoudi, R. Jimenez, and L. Verde, “An improved fitting formula for the dark matter bispectrum”, Journal of Cosmology and Astroparticle Physics 2012 no. 02, (Feb, 2012) 047–047.
  * [47] R. Takahashi, T. Nishimichi, T. Namikawa, A. Taruya, I. Kayo, K. Osato, Y. Kobayashi, and M. Shirasaki, “Fitting the nonlinear matter bispectrum by the Halofit approach”, Astrophys. J. 895 no. 2, (2020) 113, arXiv:1911.07886 [astro-ph.CO].
  * [48] T. Namikawa, “cmblensplus: Cosmic microwave background tools.” Astrophysics source code library, record ascl:2104.021, Apr., 2021.
  * [49] A. Lewis, A. Challinor, and A. Lasenby, “Efficient Computation of Cosmic Microwave Background Anisotropies in Closed Friedmann-Robertson-Walker Models”, The Astrophysical Journal 538 no. 2, (Aug, 2000) 473–476, arXiv:astro-ph/9911177.
  * [50] A. Zonca, L. Singer, D. Lenz, M. Reinecke, C. Rosset, E. Hivon, and K. Gorski, “healpy: equal area pixelization and spherical harmonics transforms for data on the sphere in Python”, Journal of Open Source Software 4 no. 35, (Mar, 2019) 1298.
  * [51] T. Namikawa, “CMB internal delensing with general optimal estimator for higher-order correlations”, Physical Review D 95 no. 10, (May, 2017) .