# Constraints on ${\tau_{\text{NL}}}$ from Planck temperature and polarization

${}^{\text{}}$Kareem Marzouk  0000-0003-2060-8956 Department of Physics &
Astronomy, University of Sussex, Brighton BN1 9QH, United Kingdom
${}^{\text{}}$Antony Lewis  0000-0001-5927-6667 Department of Physics &
Astronomy, University of Sussex, Brighton BN1 9QH, United Kingdom
${}^{\text{}}$Julien Carron  0000-0002-5751-1392 Université de Genève,
Département de Physique Théorique et CAP, 24 Quai Ansermet, CH-1211 Genève 4,
Switzerland  Department of Physics & Astronomy, University of Sussex, Brighton
BN1 9QH, United Kingdom

###### Abstract

We update constraints on the amplitude of the primordial trispectrum, using
the final Planck mission temperature and polarization data. In the squeezed
limit, a cosmological local trispectrum would be observed as a spatial
modulation of small-scale power on the CMB sky. We reconstruct this signal as
a source of statistical anisotropy via quadratic estimator techniques. We
systematically demonstrate how the estimated power spectrum of a reconstructed
modulation field can be translated into a constraint on ${\tau_{\text{NL}}}$
via likelihood methods, demonstrating the procedures effectiveness by
inferring known ${\tau_{\text{NL}}}$ signal(s) from simulations. Our baseline
results constrain ${\tau_{\text{NL}}}<1700$ at the 95% confidence level,
providing the most stringent constraints to date.

##  I Introduction

A wide range of observations appear to be broadly consistent with a Universe
that started the hot Big Bang phase with purely Gaussian adiabatic super-
horizon scalar perturbations in an otherwise Friedmann-Robertson-Walker
background. Gaussian adiabatic perturbations have only one locally varying
degree of freedom (‘single clock’), which is not locally observable on super-
horizon scales since it corresponds to a local re-definition of the time
coordinate or scale factor [1, 2, 3]. Gaussian adiabatic models are simple and
highly predictive, for example, the small-scale CMB perturbations in each
independent Hubble patch at recombination are expected to be statistically
identical. This would be consistent with inflation, where perturbation
generation is dominated by the effect of a single effective scalar field. Any
deviation from this simple prediction would be powerful evidence for a more
complex multi-field inflationary evolution, or other new early universe
physics involving more than one local degree of freedom. In this paper we
constrain large-scale modulation of small-scale CMB fluctuations by
reconstructing the modulation field using the latest Planck data.

In the language of primordial non-Gaussianity, a large-scale modulation of
small-scale power corresponds to ‘squeezed’-shape non-Gaussianity (see Ref.
[4, 5, 6] for reviews). We model the initial non-Gaussian curvature
perturbation $\zeta$ as a locally-modulated version of an underlying Gaussian
field $\zeta_{0}$, so that in real space

$\zeta({\mbox{\boldmath{$x$}}})\approx\zeta_{0}({\mbox{\boldmath{$x$}}})(1+\phi({\mbox{\boldmath{$x$}}})),$
(1)

where $\phi$ is some modulating field, and we assume
$|\phi({\mbox{\boldmath{$x$}}})|\ll 1$. If $\phi$ is correlated to
$\zeta_{0}$, there will be a non-zero local (squeezed) curvature bispectrum,
as commonly parameterized by $f_{\text{NL}}$, corresponding in the CMB to the
amplitude of the small-scale CMB power spectrum varying in a way that’s
correlated to the large-scale temperature. However, whether or not $\phi$ is
correlated to $\zeta_{0}$, there will always be a non-zero squeezed
trispectrum, corresponding to the small-scale $\zeta$ power varying spatially
with $\phi$. In the inflationary context this is usually parameterized by
${\tau_{\text{NL}}}$, where
${\tau_{\text{NL}}}\geq\left(\frac{5}{6}f_{\text{NL}}\right)^{2}$ depending on
the degree of correlation of $\phi$ and $\zeta_{0}$ [7]. The local bispectrum
constraint $f_{\text{NL}}=-0.9\pm 5.1$ from Planck [8] means that any
${\tau_{\text{NL}}}$ from a nearly scale-invariant $\phi\propto\zeta_{0}$ must
be very small, however there remains the possibility of an uncorrelated
modulation giving rise to a more substantial ${\tau_{\text{NL}}}$.

The best current empirical constraint comes from the Planck 2013 nominal-
mission temperature data, ${\tau_{\text{NL}}}<2800\,(95\%{\rm CL})$ [9], where
${\tau_{\text{NL}}}\sim 500$ would correspond to a modulation at the
$|\phi|\sim 10^{-3}$ level111We quote constraints assuming no other non-
Gaussian signals apart from CMB lensing (and other effects modelled in the
simulations) are present. Ref. [10] find a weaker joint constraint with
$g_{\text{NL}}$, but also used a more suboptimal estimator.. This result was a
conservative upper limit, accounting for evidence of frequency-dependent
foregrounds or systematic modulation in power (particularly in the octopole).
With full-mission Planck data [11] the temperature noise is substantially
reduced and the CMB polarization data give additional constraining power;
combined with better cleaning of foregrounds and modelling of systematics, it
should therefore be possible to substantially improve on the first Planck
constraint. The updated NPIPE analysis pipeline [12] is also now available
using slightly more data, lowering the noise further, as well as offering the
prospect of a better control of systematics and more reliable simulations.

In practice, only large-scale modes of the modulation field can be constrained
well, since small-scale variations in power are impossible to distinguish from
random fluctuations in local realizations of the modes [13, 14]. We focus on
modulation modes corresponding to multipoles $L\lesssim 100$, where at
recombination the modulation field remains largely super-horizon and
approximately constant through the last-scattering surface. In this case, a
primordial curvature modulation translates directly into a large-scale
modulation of the small-scale CMB [14, 9]

$X(\hat{{\mbox{\boldmath{$n$}}}})\approx
X_{g}(\hat{{\mbox{\boldmath{$n$}}}})[1+\phi(\hat{{\mbox{\boldmath{$n$}}}},r_{*})]\equiv
X_{g}(\hat{{\mbox{\boldmath{$n$}}}})[1+f(\hat{{\mbox{\boldmath{$n$}}}})],$ (2)

where $X$ is the small-scale CMB temperature perturbation or polarization
field, $X_{g}$ is the corresponding standard unmodulated linear Gaussian
field, and $r_{*}$ is the radial distance to last scattering. The two-
dimensional modulation field
$f(\hat{{\mbox{\boldmath{$n$}}}})\equiv\phi(\hat{{\mbox{\boldmath{$n$}}}},r_{*})$
parameterizes the modulation directly in terms of the CMB observables, and can
be reconstructed from the data using quadratic estimator techniques [15]. For
measuring ${\tau_{\text{NL}}}$, where most of the signal is in modulation
modes $L\ll 100$, the super-horizon approximation is accurate at the percent
level, allowing a greatly simplified analysis [14] that is nearly optimal.
With no additional computational cost we can also constrain
$f({\hat{{\mbox{\boldmath{$n$}}}}})$ more generally, without assuming
${\tau_{\text{NL}}}$ form, allowing for more general scale dependence of the
modulation field. The analysis is then equivalent to looking for a power
modulation, as possibly suggested by hints for statistical anisotropy [16, 15,
17, 18]. We shall however assume that all the small-scale modes are modulated
as implied by the model of Eq. (2), rather than allowing for $X_{g}$ to be
modulated in a different way depending on the CMB harmonic scale. This allows
us to place good constraints using the full power of high-resolution Planck
data, and avoids issues regarding possible a posteriori choices for more
general parameterizations (for a discussion of scale-dependent models, see
e.g., Ref. [19, 20]).

The signal modulation of Eq. (2) is of exactly the same structure as having an
optical depth to last scattering ($\Delta\tau$) that varies with line of
sight, with

$X(\hat{{\mbox{\boldmath{$n$}}}})=e^{-\Delta\tau(\hat{{\mbox{\boldmath{$n$}}}})}X_{g}(\hat{{\mbox{\boldmath{$n$}}}})\approx
X_{g}(\hat{{\mbox{\boldmath{$n$}}}})[1-\Delta\tau(\hat{{\mbox{\boldmath{$n$}}}})].$
(3)

The modulation reconstruction can therefore also be used to constrain patchy
reionization [21, 22, 23, 24, 25, 26]. The main difference is that the
standard ${\tau_{\text{NL}}}$ model imposes a specific scale-dependence to the
modulation power spectrum, where the signal of interest is concentrated on
very large scales. Optical-depth perturbations by contrast have signal out to
small scales, since the modulation more closely follows the density spectrum
than the curvature/potential spectrum. Interpretation of the patchy
reionization constraint is also more complicated as it depends on the details
of the reionization model. There is also a significant correlation with CMB
lensing [27, 28], and polarization can also be generated by scattering at
reionization. In this work we focus on constraining ${\tau_{\text{NL}}}$, but
our full modulation reconstruction power spectrum (which is dominated by the
temperature modulations) demonstrates a non-detection which could also be
translated into a constraint on specific reionization models. For Planck-level
sensitivity, the patchy reionization signal is not expected to significantly
contaminate ${\tau_{\text{NL}}}$ at the level that could be detected.

##  II Theory and Methodology

###  II.1 A squeezed ${\tau_{\text{NL}}}$ trispectrum

We construct the general trispectrum from four wavevectors
$\textbf{l}_{1},\textbf{l}_{2},\textbf{l}_{3},\textbf{l}_{4}$ that define a
closed quadrilateral in harmonic space. We will consider $X=X(\textbf{l})$ to
be a random field, which may correspond to, for example, the temperature
fluctuations of the CMB at last scattering. The four-point function $\langle
X(\textbf{l}_{1})X(\textbf{l}_{2})X(\textbf{l}_{3})X(\textbf{l}_{4})\rangle$
will receive two contributions, which are referred to as “connected” and
“disconnected” pieces. Regardless of the statistical properties of the field
$X$, the disconnected component can be factorized into the products of two-
point function $\langle X(\textbf{l}_{i})X(\textbf{l}_{j})\rangle$, i.e.,
power spectra, which are typically non-zero even if $X$ is a Gaussian random
field. By comparison, the connected part of the four-point function is exactly
zero when $X$ is Gaussian. We can therefore probe the non-Gaussianity of the
field $X$ by analysing the connected trispectrum.

In this work, we consider a diagonally squeezed form of the trispectrum. To
understand this, let us first consider a triangle in harmonic space. If we
take this triangle, and squeeze one of its sides such that
$|\mathbf{L}|\ll|\textbf{l}_{1}|,|\textbf{l}_{2}|$ and
$|\mathbf{L}|\approx|\textbf{l}_{1}+\textbf{l}_{2}|$, then we necessarily have
one large-scale mode and two small-scale modes. If the bispectrum positive, in
regions of space with a positive large-scale mode, the small-scale power over
those regions must be enhanced; correspondingly, it must be suppressed if the
bispectrum is negative; and vice versa for a negative large-scale mode. The
small-scale power therefore appears to be modulated by the large-scale mode,
with a correlation between the small-scale power and the large-scale mode.

By extension to this triangular picture, we could append two further modes
($\textbf{l}_{3},\textbf{l}_{4}$) to the squeezed side of the triangle to form
a closed quadrilateral. So long as these additional modes are much longer than
the common diagonal length,
$|\mathbf{L}|\ll|\textbf{l}_{3}|,|\textbf{l}_{4}|$, then we will form a
diagonal-squeezed trispectrum. We display this geometry in Figure 1. In this
case, the field with large-scale mode $L$ is not one of the observed modes
defining the measured trispectrum, and does not have to be directly observable
itself. If there is a correlation between the large-scale mode and the
curvature perturbation, there will potentially be an observable bispectrum.
However, the large-scale mode may be uncorrelated, or weakly correlated, in
which case there may be a substantial trispectrum without a corresponding
observable bispectrum.

If the large-scale mode is associated with a scalar field, as expected for
example in simple multi-field inflation models, the squeezed-diagonal
trispectrum is independent of the relative orientation of the small and large-
scale modes. The small-scale power is therefore expected to remain locally
isotropic, but to have a spatial modulation determined by the large-scale
mode. Since the large-scale mode may not be observable, the connected
trispectrum is what quantifies this observable modulation in small-scale
power.

Figure 1: Left: Example trispectrum. Wavevectors form a closed quadrilateral
in harmonic space. The quadrilateral can be described equivalently by two
triangles related by a common wavevector $L$. Right: Schematic of a squeezed
trispectrum, corresponding to the ${\tau_{\text{NL}}}$ signal peak; $L\ll
l_{1},l_{2},l_{3},l_{4}$ and
$L\approx|\mathbf{l}_{1}+\mathbf{l}_{2}|\approx|\mathbf{l}_{3}+\mathbf{l}_{4}|$.

This squeezed diagonal trispectrum is what is known as the
${\tau_{\text{NL}}}$-shape. Its definition is usually accompanied with another
shape constructed via squeezing one outer length of the quadrilateral such
that $|\textbf{l}_{1}|\ll,|\textbf{l}_{2}|,|\textbf{l}_{3}|,|\textbf{l}_{4}|$,
rather than across its diagonal. This “$g_{\text{NL}}$-form” trispectrum
produces a large-scale modulation of a small-scale bispectrum [4]. Together,
we conventionally write the connected trispectrum as the sum

$\langle
X(\textbf{l}_{1})X(\textbf{l}_{2})X(\textbf{l}_{3})X(\textbf{l}_{4})\rangle_{c}=\langle
X(\textbf{l}_{1})X(\textbf{l}_{2})X(\textbf{l}_{3})X(\textbf{l}_{4})\rangle^{g_{\text{NL}}}_{c}+\langle
X(\textbf{l}_{1})X(\textbf{l}_{2})X(\textbf{l}_{3})X(\textbf{l}_{4})\rangle^{{\tau_{\text{NL}}}}_{c}.$
(4)

The ${\tau_{\text{NL}}}$ and $g_{\text{NL}}$ trispectra can be produced by
local quadratic and cubic corrections to an otherwise Gaussian initial
perturbation field. It is corrections of this type that are computed in the
$\delta N$ formalism for inflationary perturbations [29, 19, 5, 30], and these
are often used as key predictors for models of inflation described by multiple
fields. The fact that the peak signals for ${\tau_{\text{NL}}}$ and
$g_{\text{NL}}$ are attained in different momentum configurations allow us to
form constraints almost independently of one another. In the remainder of this
paper, we will only consider the ${\tau_{\text{NL}}}$ trispectrum, and refer
the reader to the 2018 constraints on $g_{\text{NL}}$ by the Planck
collaboration [8]. It is usually also assumed that the large-scale modes are
approximately scale invariant (have a power spectrum proportional to
$P_{\zeta_{0}}$, as would be the case for local $\zeta_{0}^{2}$ and
$\zeta_{0}^{3}$ corrections), which then determines the overall wavenumber-
dependence of the signals.

The explicit form of the ${\tau_{\text{NL}}}$ trispectrum in Eq. (4) connects
the primordial fluctuations to those on the last scattering surface. We can
write the trispectrum in its reduced form (i.e. with redundant symmetries
eliminated) as [31]

$p^{l_{1}l_{2}}_{l_{3}l_{4}}(L)={\tau_{\text{NL}}}\int\text{d}r_{1}\text{d}r_{2}r_{1}^{2}r_{2}^{2}F_{L}(r_{1},r_{2})\Xi(r_{1},l_{1},l_{2})\Xi(r_{2},l_{3},l_{4}),$
(5)

where $F_{L}(r_{1},r_{2})=4\pi\int\text{d}\ln
k\mathcal{P}_{\zeta}(k)j_{L}(kr_{1})j_{L}(kr_{2})$ and $\Xi$ contains the
radiative transfer functions and further spherical Bessel functions $j_{L}$
required for projecting onto the CMB surface. Exact results can be obtained
numerically via Boltzmann codes, but there are analytic simplifications. In
particular at recombination the kernel
$F(r_{1}=r_{*},r_{2}=r_{*})=C_{L}^{\zeta_{*}}$ and in the squeezed limit the
(approximate) constancy of the large-scale mode $L$ relative to the short
scale fluctuations means that
$\int\text{d}rr^{2}\Xi(r,l_{i},l_{j})\approx(C_{l_{i}}+C_{l_{j}})$ [14]. The
squeezed ${\tau_{\text{NL}}}$ trispectrum at recombination is therefore
approximated as

$p^{l_{1}l_{2}}_{l_{3}l_{4}}(L)\approx{\tau_{\text{NL}}}C_{L}^{\zeta_{*}}(C_{l_{1}}+C_{l_{2}})(C_{l_{3}}+C_{l_{4}}).$
(6)

This is exactly the form of squeezed trispectrum expected from the modulated
Gaussian field model of Eq. (2), where the modulation field is proportional to
the curvature perturbation at recombination, $f\propto\zeta_{*}$.

###  II.2 The modulation field estimator

We denote the observed CMB temperature and polarization fields with the block
vector

$\mathbf{\hat{X}}=({}_{0}\mathbf{\hat{X}},\ {}_{+2}\mathbf{\hat{X}},\
{}_{-2}\mathbf{\hat{X}})^{T}=(\mathbf{\hat{T}},\ {}_{+2}\mathbf{\hat{P}},\
{}_{-2}\mathbf{\hat{P}})^{\mathrm{T}},$ (7)

where we have adopted the use of “spin indices” to indicate the spin-0
temperature field $\mathbf{\hat{T}}$ and spin-2 polarization fields ${}_{\pm
2}\mathbf{\hat{P}}=\mathbf{\hat{Q}}\pm i\mathbf{\hat{U}}$ in order to write
the fields in a unified notation. CMB fields that appear in bold typeface
refer to the vector of pixel-data, and if non-bold refer to their scalar value
when evaluated at a specific location on the sky, e.g.
$\mathrm{\hat{T}}(\hat{{\mbox{\boldmath{$n$}}}})\in\mathbf{T}$.
$\mathbf{\hat{Q}}$ and $\mathbf{\hat{U}}$ are the Stokes polarization maps
which are obtained at the detector level. Each data component in Eq. (7) is
modelled as the linear combination of the cosmological signal and noise
$\mathbf{\hat{X}}=\mathbf{X}+\mathbf{N}$, each of which have the same block
vector structure: $\mathbf{X}=({}_{0}\mathbf{X},\ {}_{+2}\mathbf{X},\
{}_{-2}\mathbf{X})^{\mathrm{T}}$, $\mathbf{N}=({}_{0}\mathbf{N},\
{}_{+2}\mathbf{N},\ {}_{-2}\mathbf{N})^{\mathrm{T}}$. A large-scale modulation
that is approximately constant through the last-scattering light cone will
alter each component of the signal vector in the same way, which we define as
a linear real-space operation
${}_{s}\mathrm{X}(\hat{{\mbox{\boldmath{$n$}}}})=[1+f(\hat{{\mbox{\boldmath{$n$}}}})]{}_{s}\mathrm{X}_{0}(\hat{{\mbox{\boldmath{$n$}}}})$
(where the modulation acting on the otherwise unmodulated signal
$\mathbf{X}_{0}=({}_{0}\mathbf{X}_{0},\ {}_{+2}\mathbf{X}_{0},\
{}_{-2}\mathbf{X}_{0})^{\mathrm{T}}$). We then model the observed data as

$\mathbf{\hat{X}}=\mathbfcal{B}\mathbf{X}+\mathbf{N},$ (8)

where $\mathbfcal{B}$ denotes the operation of the instrumental beam (a
convolution in real space) and pixel window function.

Assuming that the unmodulated signal and noise components are reliably
approximated as Gaussian random fields, we can use maximum likelihood
estimates (MLE) to build estimates of $f$. We follow the methods developed by
Hirata and Seljak [32] (for CMB lensing), and later by Hanson and Lewis [15]
(for generalized weak statistical anisotropy), defining the (negative) log-
likelihood function

$\displaystyle\begin{split}\mathcal{L}[\mathbf{f}]&\equiv-\ln
P(\mathbf{\hat{X}}|\mathbf{f})\\\
&=\frac{1}{2}\mathbf{\hat{X}}^{\dagger}\left(\mathbf{C^{\hat{X}\hat{X}}}\right)^{-1}\mathbf{\hat{X}}+\frac{1}{2}\ln\mathrm{det}\mathbf{C^{\hat{X}\hat{X}}}+\text{const}.\end{split}$
(9)

The conditional probability $P(\mathbf{\hat{X}}|\mathbf{f})$ appearing in the
log-likelihood is the standard multivariate Gaussian density with data
covariance $\mathbf{C^{\hat{X}\hat{X}}}$, which is not isotropic as it depends
on the realization of the statistical anisotropy field $\mathbf{f}$. Note that
we have represented the modulation $f$ as the vector $\mathbf{f}$ (here
containing only one item) for notational compatibility with the formalism.

We can find $\mathbf{f}$ corresponding to the maximum-likelihood point by
identifying the minimum of Eq. (9). This can be computed directly by finding
where the functional derivative of the likelihood function is zero,
$\delta\mathcal{L}/\delta\mathbf{f^{\dagger}}=0$. As explicitly given in [15],
the derivative may be conveniently expressed in terms of the combination

$\mathbfcal{H}=\frac{1}{2}\left[\mathbfcal{B}^{\dagger}\left(\mathbf{C^{\hat{X}\hat{X}}}\right)^{-1}\mathbf{\hat{X}}\right]^{\dagger}\frac{\mathbf{\delta
C^{\hat{X}\hat{X}}}}{\delta\mathbf{f^{\dagger}}}\left[\mathbfcal{B}^{\dagger}\left(\mathbf{C^{\hat{X}\hat{X}}}\right)^{-1}\mathbf{\hat{X}}\right],$
(10)

such that

$\frac{\delta\mathcal{L}[\mathbf{f}]}{\delta\mathbf{f^{\dagger}}}=\langle\mathbfcal{H}\rangle-\mathbfcal{H}=0.$
(11)

The expectation value $\langle\mathbfcal{H}\rangle$ corresponds to a “mean
field”, which comes from the derivative of the determinant term. It can
conveniently be thought of as a data simulation average.222From the identity
$\mathrm{Tr}\mathbf{A}=\langle\mathbf{x}^{\dagger}\mathbf{A}\mathbf{C}^{-1}\mathbf{x}\rangle$,
which holds for a random vector $\mathbf{x}$ with
$\langle\mathbf{x}\mathbf{x}^{\dagger}\rangle=\mathbf{C}$ and any (square)
matrix $\mathbf{A}$, we have
$2\langle\mathbfcal{H}\rangle=\left\langle\mathbf{\hat{X}}^{\dagger}\Big{(}\mathbf{C^{\hat{X}\hat{X}}}\Big{)}^{-1}\frac{\mathbf{\delta
C^{\hat{X}\hat{X}}}}{\delta\mathbf{f^{\dagger}}}\Big{(}\mathbf{C^{\hat{X}\hat{X}}}\Big{)}^{-1}\mathbf{\hat{X}}\right\rangle=\mathrm{Tr}\left[\Big{(}\mathbf{C^{\hat{X}\hat{X}}}\Big{)}^{-1}\frac{\mathbf{\delta
C^{\hat{X}\hat{X}}}}{\delta\mathbf{f^{\dagger}}}\right].$ Solutions for
$\mathbf{f}$ that satisfy Eq. (11) can be found iteratively, e.g. with a
Newton-Raphson approach. This is particularly efficient for weak sources of
anisotropy since $\mathbf{f}\approx\mathbf{0}$, meaning that a single
iteration from $\mathbf{f}=\mathbf{0}$ will yield a good approximate solution
to the anisotropy. Explicitly, the first iteration produces the estimator

$\mathbf{\hat{f}}=-\left[\frac{\delta}{\delta\mathbf{f^{\dagger}}}(\langle\mathbfcal{H}\rangle-\mathbfcal{H})^{\dagger}\Big{|}_{0}^{\dagger}\right]^{-1}(\langle\mathbfcal{H}_{0}\rangle-\mathbfcal{H}_{0}),$
(12)

where the lower index “0” indicates evaluation with $\mathbf{f}=\mathbf{0}$,
corresponding to the initial step. The functional derivative in Eq. (12) can
be approximated with the Fisher information

$\mathbfcal{F}_{0}=\left\langle\frac{\delta}{\delta\mathbf{f^{\dagger}}}(\langle\mathbfcal{H}_{0}\rangle-\mathbfcal{H}_{0})\right\rangle,$
(13)

which provides a natural form of normalization for the likelihood based
estimator. The normalized quadratic estimator that we use then takes the form

$\mathbf{\hat{f}}=\mathbfcal{F}^{-1}_{0}(\mathbfcal{H}_{0}-\langle\mathbfcal{H}_{0}\rangle),$
(14)

where

$\mathbfcal{H}_{0}=\frac{1}{2}\mathbf{\mathbf{\bar{X}}}^{\dagger}\left[\frac{\delta\mathbf{C^{XX}}}{\delta\mathbf{f}^{\dagger}}\right]_{0}\mathbf{\mathbf{\bar{X}}},$
(15)

and we have introduced the inverse (co-)variance filtered fields, defined as
$\mathbf{\mathbf{\bar{X}}}\equiv\mathbfcal{B}^{\dagger}(\mathbf{C^{\hat{X}\hat{X}}})^{-1}\big{|}_{0}\mathbf{\hat{X}}$.
An analytic expression for the Fisher normalization is in general non-trivial,
so in practice we consider the isotropic full-sky limit of the reconstruction
to obtain a diagonal simplification that we then correct for with a Monte
Carlo correction. We implement a normalization that is consistent333The
normalization differs in the case of the joint estimators as defined in
Eq.(24), but is identical otherwise. with that given in Refs. [15, 14].

The covariance $\mathbf{C^{{X}{X}}}$ is (in the general case) a $3\times 3$
block matrix, whose elements
$\mathbf{C}^{s_{1}s_{2}}=\langle{}_{s_{1}}{\mbox{\boldmath{$X$}}}\
{}_{s_{2}}{\mbox{\boldmath{$X$}}}^{\dagger}\rangle$ represent the covariance
between the spin fields with fixed anisotropy. We can rewrite the product of
matrices as the sum over spin-indices to obtain

$\mathbfcal{H}_{0}=\frac{1}{2}\sum_{{s_{1}},{s_{2}}}{}_{s_{1}}\mathbf{\bar{X}}^{\dagger}\left[\frac{\delta\mathbf{C}^{{s_{1}}{s_{2}}}}{\delta\mathbf{f}^{\dagger}}\right]_{0}{}_{s_{2}}\mathbf{\bar{X}}.$
(16)

We calculate (up to linear order in modulation) the derivative of the
covariance to be

$\left[\frac{\delta\mathrm{C}^{{s_{1}}{s_{2}}}(\hat{{\mbox{\boldmath{$n$}}}}_{1},\hat{{\mbox{\boldmath{$n$}}}}_{2})}{\delta
f(\hat{{\mbox{\boldmath{$n$}}}})}\right]=\langle{}_{s_{1}}\mathrm{X}(\hat{{\mbox{\boldmath{$n$}}}}_{1})\
{}_{s_{2}}\mathrm{X}^{*}(\hat{{\mbox{\boldmath{$n$}}}}_{2})\rangle(\delta(\hat{{\mbox{\boldmath{$n$}}}}-\hat{{\mbox{\boldmath{$n$}}}}_{1})+\delta(\hat{{\mbox{\boldmath{$n$}}}}-\hat{{\mbox{\boldmath{$n$}}}}_{2})),$
(17)

where
$\delta(\hat{{\mbox{\boldmath{$n$}}}}-\hat{{\mbox{\boldmath{$n$}}}}_{i})$
appearing on the right-hand side result from locality of the modulation. By
expanding ${}_{s}\mathrm{X}(\hat{{\mbox{\boldmath{$n$}}}})$ into spherical
harmonics Eq. (17) can be written in terms of the fiducial theoretical power
spectra associated with the spin fields

$\langle{}_{s_{1}}\mathrm{X}(\hat{{\mbox{\boldmath{$n$}}}}_{1})\
{}_{s_{2}}\mathrm{X}^{*}(\hat{{\mbox{\boldmath{$n$}}}}_{2})\rangle=\sum_{\ell_{1},m_{1}}\sum_{\ell_{2},m_{2}}\mathcal{B}_{\ell_{1}}\mathcal{B}_{\ell_{2}}\
{}_{s_{1}}Y_{\ell_{1}m_{1}}(\hat{{\mbox{\boldmath{$n$}}}}_{1})\
{}_{s_{2}}Y_{\ell_{2}m_{2}}^{*}(\hat{{\mbox{\boldmath{$n$}}}}_{2})\langle{}_{s_{1}}X_{\ell_{1}m_{1}}(\hat{{\mbox{\boldmath{$n$}}}}_{1})\
{}_{s_{2}}X_{\ell_{2}m_{2}}^{*}(\hat{{\mbox{\boldmath{$n$}}}}_{2})\rangle,$
(18)

where $\mathcal{B}_{\ell}$ is the harmonic space representation of the
instrumental beam and pixel window $\mathbfcal{B}$. The quadratic building
block for the modulation estimator is then

${}_{\mathcal{S}}\mathcal{H}_{0}^{XX}(\hat{{\mbox{\boldmath{$n$}}}})=\sum_{s\in\mathcal{S}}{}_{s}\bar{X}(\hat{{\mbox{\boldmath{$n$}}}})^{*}{}_{s}X^{\mathrm{WF}}(\hat{{\mbox{\boldmath{$n$}}}}),$
(19)

where $\mathcal{S}$ is the set of CMB field spin indices to sum over. In
practice, we implement a version of conjugate-gradient descent to obtain an
approximate iterative solution for the Wiener-filtered (WF) fields
${{\mbox{\boldmath{$X$}}}}^{\mathrm{WF}}=\bm{C}\mathbf{\mathbf{\bar{X}}}\equiv\bm{C}\mathbfcal{B}^{\dagger}(\mathbf{C^{\hat{X}\hat{X}}})^{-1}\big{|}_{0}\mathbf{\hat{X}}$,
where $\mathbf{C}$ is the covariance of fiducial power spectra, for both data
and the simulations [33, 34]. The inverse-variance filtered (IVF) fields are
simply related to the WF fields by dividing the matrix of fiducial spectra
$\mathbf{C}$. The quadratic estimator remains unbiased, though slightly
suboptimal [35], if an approximate version of $\mathbf{C^{\hat{X}\hat{X}}}$ is
used, as long as the normalization is computed consistently. For simplicity,
we take $C^{TE}=0$ in the theory part of the covariance
$\mathbf{C^{\hat{X}\hat{X}}}$ to allow separate filtering of temperature and
polarization maps, and approximate the noise as isotropic over the unmasked
area. The inverse noise is taken to be zero in the masked regions, so that
masked pixels do not contribute to the WF maps. The modulation estimator does
respond to $B$ modes, but for Planck they contribute negligibly to the
${\tau_{\text{NL}}}$ constraint.444A scalar primordial non-Gaussianity signal
should not produce $B$ modes, but a modulation does. In the ultra-squeezed
limit in which the modulation approximation is exact, the $B$-mode
contribution also vanishes. Including the $B$ modes is correct for modulation
reconstruction, but for futuristic measurements where the $B$-mode
contribution may become non-negligible due to low delensed cosmic variance,
the $B$ mode contribution should strictly be dropped when using the modulation
model as an approximation to constrain ${\tau_{\text{NL}}}$.  The equivalent
harmonic space representation of (19) can be directly obtained via the spin-
zero harmonic transformation

${}_{\mathcal{S}}\mathcal{H}_{0,LM}^{XX}=\int\mathrm{d}\hat{{\mbox{\boldmath{$n$}}}}\
{}_{\mathcal{S}}\mathcal{H}_{0}^{XX}(\hat{{\mbox{\boldmath{$n$}}}})_{0}Y_{LM}^{*}(\hat{{\mbox{\boldmath{$n$}}}}).$
(20)

By defining the set of spin-indices, $\mathcal{S}$ appearing in Eq. (19) we
can construct different estimators. For example, we obtain the (unnormalized)
temperature estimator given in [15] by setting
$\mathcal{S}_{\mathrm{TT}}=\\{0\\}$

${}_{\mathcal{S}_{\mathrm{TT}}}\mathcal{H}_{0}^{XX}(\hat{{\mbox{\boldmath{$n$}}}})=\bar{{T}}(\hat{{\mbox{\boldmath{$n$}}}})T^{\mathrm{WF}}(\hat{{\mbox{\boldmath{$n$}}}}),$
(21)

whereas the polarization and minimum variance estimators are defined by
$\mathcal{S}_{\mathrm{PP}}=\\{+2,-2\\}$ and
$\mathcal{S}_{\mathrm{MV}}=\\{0,+2,-2\\}$ respectively. These latter
estimators have simple expressions analogous to Eq. (21). Each estimator can
equivalently be expressed in terms of the CMB harmonics
$\\{T_{lm},E_{lm},B_{lm}\\}$, though appear more cumbersome.

The quadratic estimators can reconstruct the modulation field up to the limits
of instrumental noise and the cosmic variance of the CMB fields themselves.
These Gaussian sources produce a reconstruction noise, $N_{L}^{(0)}$, which
corresponds directly to the disconnected trispectrum contributions discussed
in II.1, and should be subtracted from any modulation field power spectrum
estimate. This can be seen by computing the covariance between two quadratic
estimators, which will result in Gaussian terms, composed of products of noise
and (theoretical) CMB power spectra. Under idealized conditions, this
reconstruction noise is statistically isotropic and equivalent to the
isotropic Fisher normalization, such that it is straightforward to account
for. We normalize using the full-sky analytic result, but correct our final
spectra by a Monte Carlo normalization to account for inaccuracies in the
analytic result. The reconstruction noise bias and mean field are also
calculated from simulations, as described in the following section. Using
harmonic representation of the full-sky analytic Fisher normalization gives us
the normalized quadratic estimator

${}_{\mathcal{S}}\hat{f}_{LM}^{XX}=({}_{\mathcal{S}}\mathcal{F}^{XX}_{0,L})^{-1}{}_{\mathcal{S}}\mathcal{H}_{0,LM}^{XX},$
(22)

where we have appended indices to the normalization factor in order to keep
track of the estimator properties, which will alter the required normalization
in general.

We can generalize the quadratic estimator in the Eq. (19) to a “joint
estimator”, which combines information from pairs of maps which may correspond
to e.g. distinct observations of the CMB. For map pairs
$X(\hat{{\mbox{\boldmath{$n$}}}}),Z(\hat{{\mbox{\boldmath{$n$}}}})$ we define
the unnormalized estimator analogous to (19)

${}_{\mathcal{S}}\mathcal{H}_{0}^{XZ}(\hat{{\mbox{\boldmath{$n$}}}})=\frac{1}{2}\sum_{s\in\mathcal{S}}\left[{}_{s}\bar{Z}^{*}(\hat{{\mbox{\boldmath{$n$}}}}){}_{s}X^{\mathrm{WF}}(\hat{{\mbox{\boldmath{$n$}}}})+{}_{s}\bar{X}^{*}(\hat{{\mbox{\boldmath{$n$}}}}){}_{s}Z^{\mathrm{WF}}(\hat{{\mbox{\boldmath{$n$}}}})\right],$
(23)

where the symmetrized mean accounts for the fact that in general the quadratic
building block is not symmetric under exchange $X\leftrightarrow Z$.
Similarly, we write the symmetrized normalization as

${}_{\mathcal{S}}\mathcal{A}^{XZ}_{L}=\frac{1}{2}\Big{(}{}_{\mathcal{S}}\mathcal{F}^{XZ}_{0,L}+{}_{\mathcal{S}}\mathcal{F}^{ZX}_{0,L}\Big{)},$
(24)

to obtain

${}_{\mathcal{S}}\hat{f}_{LM}^{XZ}=({}_{\mathcal{S}}\mathcal{A}^{XZ}_{L})^{-1}{}_{\mathcal{S}}\mathcal{H}_{0,LM}^{XZ}.$
(25)

Note that if $X=Z$ these definitions automatically recover the former result
of Eq. (22). In the remainder of this paper we will use power spectrum
estimates of the modulation field computed from the normalized quadratic
estimator(s) described in this section.

###  II.3 Power spectrum estimator

Power spectrum estimates of the modulation field are computed from the cross-
spectrum of pairs of quadratic estimators
${}_{\mathcal{S}_{1}}\hat{f}_{LM}^{X_{1}Z_{1}}$ and
${}_{\mathcal{S}_{2}}\hat{f}_{LM}^{X_{2}Z_{2}}$. Each of these components
require mean field subtraction as given in Eq. (11), which in practice is
calculated as the average quadratic estimate over simulations with $f=0$. The
mean field serves as a map-level correction to the data that accounts for the
main sources of statistical anisotropy that are not primordial. These include
Doppler aberration and modulation, anisotropic noise and beam anisotropy, as
well as estimator response to masked regions of the sky. The accuracy to which
the mean field can be estimated is limited by the number and fidelity of the
simulations. The Doppler-induced effects are important for
${\tau_{\text{NL}}}$ since most of their signal is at very low $L$, but have a
fixed known direction such that they are expected to be accounted for
accurately by the simulations.

Monte Carlo noise on the mean field limits accuracy of the reconstruction, but
can on average be nulled at the power spectrum level by ensuring the mean
field subtracted from the data reconstruction on each leg has independent
Monte Carlo noise. Since realizations of CMB and instrumental noise are
uncorrelated from simulation-to-simulation, this can be achieved
straightforwardly by constructing two separate mean field estimates from
independent sets of simulations. The mean field (on each component
$i\in\\{1,2\\}$) is

${}_{\mathcal{S}_{i}}f_{LM}^{X_{i}Z_{i},\mathrm{MF}}=\big{\langle}{}_{\mathcal{S}_{i}}\hat{f}^{X_{i}Z_{i}}_{LM}\big{\rangle}_{\mathrm{sims}_{i}}$
(26)

such that the leading contribution to the modulation power spectrum is given
by555Our use of “$\supset$” in this section implies that whatever appears on
the right-hand side is a partial contribution to the full quantity given on
the left; which is added to linearly.

$\bar{C}_{L}^{\mathcal{S}_{1}X_{1}Z_{1},\
\mathcal{S}_{2}X_{2}Z_{2}}\supset\frac{1}{(2L+1)f_{\rm
sky}}\sum_{M}\left({}_{\mathcal{S}_{1}}\hat{f}_{LM}^{X_{1}Z_{1}}-{}_{\mathcal{S}_{1}}f_{LM}^{X_{1}Z_{1},\mathrm{MF}}\right)\left({}_{\mathcal{S}_{2}}\hat{f}_{LM}^{X_{2}Z_{2}}-{}_{\mathcal{S}_{2}}f_{LM}^{X_{2}Z_{2},\mathrm{MF}}\right)^{*}.$
(27)

If we assume that the simulations reliably capture the expected lensing
signal666The simulations have independent random lensing realizations, so
lensing produces no mean field on average, but does have a trispectrum. In
this paper we do not attempt a joint reconstruction of the lensing and the
modulation. the estimator can then be made unbiased by simply subtracting the
mean over simulations with zero modulation
$\bar{C}_{L}^{\mathcal{S}_{1}X_{1}Z_{1},\
\mathcal{S}_{2}X_{2}Z_{2}}-\langle\bar{C}_{L}^{\mathcal{S}_{1}X_{1}Z_{1},\
\mathcal{S}_{2}X_{2}Z_{2}}\rangle_{f=0}$. However, we can instead use the
optimal perturbative trispectrum estimator that does better than this, by
correcting for the realization of the Gaussian power in the data map compared
to the mean in the simulations [36, 37, 38, 39]. This realization-dependent
correction, $N_{L}^{{\rm RD},(0)}$, reduces the variance of the estimator, and
also makes it perturbatively self-correcting to misestimates of the Gaussian
power in the simulations. In the idealized full sky case this can be done by
adding the linear correction

$\bar{C}_{L}^{\mathcal{S}_{1}X_{1}Z_{1},\
\mathcal{S}_{2}X_{2}Z_{2}}\supset\langle\bar{C}_{L}^{\mathcal{S}_{1}X_{1}Z_{1},\
\mathcal{S}_{2}X_{2}Z_{2}}\rangle_{f=0}-\sum_{ab,l}\frac{\partial\langle\bar{C}_{L}^{\mathcal{S}_{1}X_{1}Z_{1},\
\mathcal{S}_{2}X_{2}Z_{2}}\rangle_{f=0}}{\partial\bar{C}_{\ell,\rm{expt}}^{ab}}\left(\bar{C}_{\ell,\rm{expt}}^{ab}-\bar{C}_{\ell,\rm{expt,fid}}^{ab}\right),$
(28)

where $\bar{C}_{\ell,\rm{expt}}^{ab}$ is the CMB power spectrum measured on
the data and $\bar{C}_{\ell,\rm{expt},fid}^{ab}$ the fiducial model of the
data power spectrum (i.e. lensed power spectrum + isotropic noise). The
dominant part of the derivative correction comes from the response of
${N_{L}^{(0)}}$ to the power, so we can approximate
$\langle\bar{C}_{L}^{\mathcal{S}_{1}X_{1}Z_{1},\
\mathcal{S}_{2}X_{2}Z_{2}}\rangle_{f=0}\approx{N_{L}^{(0)}}$ in that term; the
first term depends on ${N_{L}^{(0)}}$ and the lensing 4-point contribution to
the modulation power.

Accounting for these bias corrections, our estimator for the modulation power
spectrum is given by

$\bar{C}_{L}^{\mathcal{S}_{1}X_{1}Z_{1},\
\mathcal{S}_{2}X_{2}Z_{2}}=\frac{1}{(2L+1)f_{\rm
sky}}\sum_{M}\left({}_{\mathcal{S}_{1}}\hat{f}_{LM}^{X_{1}Z_{1}}-{}_{\mathcal{S}_{1}}f_{LM}^{X_{1}Z_{1},\mathrm{MF}}\right)\left({}_{\mathcal{S}_{2}}\hat{f}_{LM}^{X_{2}Z_{2}}-{}_{\mathcal{S}_{2}}f_{LM}^{X_{2}Z_{2},\mathrm{MF}}\right)^{*}-\Delta\bar{C}_{L}^{(0)}-N_{L}^{{\rm
RD},(0)},$ (29)

where $\Delta\bar{C}_{L}^{(0)}=\langle\bar{C}_{L}^{\mathcal{S}_{1}X_{1}Z_{1},\
\mathcal{S}_{2}X_{2}Z_{2}}\rangle_{f=0}-{N_{L}^{(0)}}$ is the connected
lensing contribution. In the general non-full-sky case of power spectra based
on four input maps $\\{X_{1},Z_{1},X_{2},Z_{2}\\}$ the explicit expressions
for these terms are:

$\begin{split}\Delta\bar{C}_{L}^{(0)}=\Big{\langle}\bar{C}_{L}^{\mathcal{S}_{1}{X_{1}}(i){Z_{1}}(i),\
\mathcal{S}_{2}{X_{2}}(i){Z_{2}}(i)}-\Big{(}\bar{C}_{L}^{\mathcal{S}_{1}{X_{1}}(i){Z_{1}}(j),\
\mathcal{S}_{2}{X_{2}}(j){Z_{2}}(i)}+\bar{C}_{L}^{\mathcal{S}_{1}{X_{1}}(i){Z_{1}}(j),\
\mathcal{S}_{2}{X_{2}}(i){Z_{2}}(j)}\Big{)}\Big{\rangle}_{i\neq
j,f=0}\end{split}$ (30)

and

$\begin{split}\bar{N}_{L}^{{\rm
RD},(0)}=\Big{\langle}&\bar{C}_{L}^{\mathcal{S}_{1}{X_{1}}(i){Z_{1}}(d),\
\mathcal{S}_{2}{X_{2}}(i){Z_{2}}(d)}+\bar{C}_{L}^{\mathcal{S}_{1}{X_{1}}(i){Z_{1}}(d),\
\mathcal{S}_{2}{X_{2}}(d){Z_{2}}(i)}+\bar{C}_{L}^{\mathcal{S}_{1}{X_{1}}(d){Z_{1}}(i),\
\mathcal{S}_{2}{X_{2}}(i){Z_{2}}(d)}\\\
+&\bar{C}_{L}^{\mathcal{S}_{1}{X_{1}}(d){Z_{1}}(i),\
\mathcal{S}_{2}{X_{2}}(d){Z_{2}}(i)}-\Big{(}\bar{C}_{L}^{\mathcal{S}_{1}{X_{1}}(i){Z_{1}}(j),\
\mathcal{S}_{2}{X_{2}}(j){Z_{2}}(i)}+\bar{C}_{L}^{\mathcal{S}_{1}{X_{1}}(i){Z_{1}}(j),\
\mathcal{S}_{2}{X_{2}}(i){Z_{2}}(j)}\Big{)}\Big{\rangle}_{i\neq
j,f=0}\end{split}$ (31)

where the indices $(i)$ and $(j)$ imply the maps are Monte Carlo simulations,
and $(d)$ the data. In practice we implement the convention $j=i+1$ to obtain
uncorrelated Monte Carlo pairs, though in principle further permutations of
indices could lower the Monte Carlo noise even further. The power spectra have
the same definition as in Eq. (27) other than the mean field term, which is
set to zero. The exception to this is the first term appearing in Eq. (30)
when each map carries the same index. This term includes lensing 4-point
information and will play a further role in our likelihood analysis, it is
computed from individual Monte Carlo simulations as

$\bar{N}_{L}^{{\rm
MC},(0)}=\big{\langle}\bar{C}_{L}^{\mathcal{S}_{1}{X_{1}}(i){Z_{1}}(i),\
\mathcal{S}_{2}{X_{2}}(i){Z_{2}}(i)}\big{\rangle}_{f=0}.$ (32)

Finally we introduce a calibration factor $k_{L}$ to form the final estimate

$\hat{C}_{L}^{\mathcal{S}_{1}X_{1}Z_{1},\
\mathcal{S}_{2}X_{2}Z_{2}}=k_{L}\bar{C}_{L}^{\mathcal{S}_{1}X_{1}Z_{1},\
\mathcal{S}_{2}X_{2}Z_{2}}.$ (33)

The calibration factor is derived empirically from with-signal simulations
such that the ensemble average reconstruction spectra match the theoretical
input
$C_{L}^{f,\mathrm{Theory}}=k_{L}\langle\bar{C}_{L}^{\mathcal{S}_{1}X_{1}Z_{1},\
\mathcal{S}_{2}X_{2}Z_{2}}\rangle_{\rm sims}$. Note that when the calibration
factor multiplies the reconstruction noise spectra we apply the same
convention; e.g. $\hat{N}_{L}^{{\rm MC},(0)}=k_{L}\bar{N}_{L}^{{\rm MC},(0)}$.

The calibration corresponds to an $\mathcal{O}(1)$ multiplicative correction
after appropriate binning and smoothing is applied, though is slightly scale-
dependent. The dipole receives the largest correction corresponding to the
greatest power loss after masking the sky, and towards smaller scales the
power is typically reduced by up to $15\%$. We assess the calibration
procedure in more detail in Section IV, and have confirmed that it leads to
more accurate ${\tau_{\text{NL}}}$ estimates, and does not lead to a false
detection when applied to without-signal reconstructions. In the remainder of
this paper, spectra that are uncalibrated or calibrated will be indicated with
a bar or hat, i.e. $\bar{C}_{L}$ or $\hat{C}_{L}$ respectively.

Figure 2: Power spectrum of the modulation mean field, reconstruction noise
biases, and lensing bias, for the NPIPE (A, B) analysis. The mean field power
is computed as the cross-spectrum of the independent mean field components
given in Eq. (26). For comparison, the dotted line shows the modulation power
expected if ${\tau_{\text{NL}}}=1500$; i.e. $C_{L}^{f}=1500C_{L}^{\zeta_{*}}$.

###  II.4 ${\tau_{\text{NL}}}$ estimation and likelihood analysis

The simplest estimate that we can construct for ${\tau_{\text{NL}}}$
corresponds to the modulation power spectrum at a fixed scale $L$, normalized
by the curvature power spectrum at recombination

$\hat{\tau}_{\mathrm{NL}}(L)=\frac{\hat{C}_{L}^{f}}{C_{L}^{\zeta_{*}}},$ (34)

where we assume that we know $C_{L}^{\zeta_{*}}$ theoretically; and can be
readily computed with e.g. CAMB. In general
$\hat{C}_{L}^{f}=\hat{C}_{L}^{\mathcal{S}_{1}X_{1}Z_{1},\
\mathcal{S}_{2}X_{2}Z_{2}}$ (i.e. is the full modulation power spectrum
estimate) but we have omitted indices for convenience. The calibrated
estimates are related to the uncalibrated estimates multiplicatively
$\hat{C}_{L}^{f}=k_{L}\bar{C}_{L}^{f}$, such that
$\hat{\tau}_{\mathrm{NL}}(L)=k_{L}\bar{\tau}_{\mathrm{NL}}(L)$. We adopt the
same bar/hat convention for all quantities in our likelihood analysis to make
this distinction clear; including the covariance and Gaussian noise bias.
Estimates at each $L$ can be combined using an inverse variance weighting. For
small signals, an approximately scale-invariant primordial spectrum, and
approximating the reconstruction noise as white over the modulation multipole
range $L_{\mathrm{min}}\leq L\leq L_{\mathrm{max}}$ that is used, this gives
the approximately-optimal combined estimator [14]

$\hat{\tau}_{\mathrm{NL}}=\frac{1}{R}\sum_{L=L_{\mathrm{min}}}^{L_{\mathrm{max}}}\frac{2L+1}{L^{2}(L+1)^{2}}\frac{\hat{C}_{L}^{f}}{C_{L}^{\zeta_{*}}},$
(35)

where the normalization $R$ accounts for the finite multipole range,
$R={L_{\mathrm{min}}}^{-2}-(L_{\mathrm{max}}+1)^{-2}$. Regardless of whether
or not there is in fact a primordial modulation field, the reconstruction is
expected to produce an approximately Gaussian random field on large scales. In
the null case, the reconstruction should correspond to a realization of the
disconnected Gaussian noise bias, $N_{L}^{(0)}$, with an approximately white
power spectrum. If there is a primordial modulation, this too is likely to
correspond to an approximately Gaussian random field (e.g., if it is generated
by quantum fluctuations early in inflation). If it is of ${\tau_{\text{NL}}}$
form, the spectrum would peak on large scales if the signal-to-noise of the
modulation is sufficiently large to be distinguishable from Gaussian
background. We therefore expect that the estimators formed as a sum over
$\hat{\tau}_{\mathrm{NL}}(L)$ are approximately $\chi^{2}$-distributed. Since
the signal peaks at low multipoles, this is a very skewed distribution, with
$95\%$ of the signal in ${\hat{\tau}_{\rm NL}}$ at $L\leq 4$ [14].

To constrain the true value of ${\tau_{\text{NL}}}$ given the high level of
realization dependence, we use an approximate model for the posterior
distribution of the set of $\hat{\tau}_{\mathrm{NL}}$ estimates. As pointed
out in Ref. [40], the statistical properties of a ${\tau_{\text{NL}}}$ signal
are analogous to those assumed when estimating the statistics of the CMB
temperature field. We largely borrow from this “single-field” procedure and
construct an approximate log-likelihood function on the cut-sky as described
in Ref. [41]

$-2\ln P\left(\\{{\tau_{\text{NL}}}(L)\\}|\\{{\hat{\tau}_{\rm
NL}}(L)\\}\right)\approx\sum_{LL^{\prime}}\Big{[}g(\hat{x}(L))\hat{N}^{\mathrm{MC},(0)}_{{\tau_{\text{NL}}}}(L)\Big{]}\Big{[}\hat{M}^{-1}\Big{]}_{LL^{\prime}}\Big{[}\hat{N}^{\mathrm{MC},(0)}_{{\tau_{\text{NL}}}}(L^{\prime})g(\hat{x}(L^{\prime}))\Big{]}+\text{const,}$
(36)

using a restricted multipole range $L_{\mathrm{min}}\leq L\leq
L_{\mathrm{max}}$. In the above ${\hat{\tau}_{\rm NL}}(L)$ is defined in Eq.
(34), and $\hat{N}_{\tau_{\text{NL}}}^{\mathrm{MC},(0)}(L)$ is of the same
form but using the reconstruction noise power spectrum in place of the
modulation power spectrum (and therefore defines the Gaussian noise bias on
${\tau_{\text{NL}}}$). $\hat{M}$ is the covariance matrix of ${\hat{\tau}_{\rm
NL}}(L)$ estimates inferred from without-signal simulation reconstructions,
and is truncated to the prescribed multipole interval before inversion. To
better account for the error incurred by the inverse-covariance due to the
finite number of simulations, we rescale the inverse by a constant (Hartlap
factor [42]) which leads to a slight broadening on the final posterior
distribution. Similarly, we account for potential Monte Carlo errors (i.e.
error on the mean field and reconstruction noise bias) by including the
additive corrections to the diagonal of the covariance [34]

$\sigma_{\mathrm{MC},L}^{2}=\left(\frac{2}{N_{\mathrm{MF}}}+\frac{9}{N_{\mathrm{bias}}}\right)\frac{2(N_{\tau_{\text{NL}}}^{(0)})^{2}}{(2L+1)f_{\rm
sky}},$ (37)

where $N_{\mathrm{MF}}$ and $N_{\mathrm{bias}}$ are the number of simulations
used to compute the mean field and realization dependent noise bias
respectively. The likelihood further relies on the functions

$g(x)=\mathrm{sign}(x-1)\sqrt{2(x-\ln(x)-1)},\quad\mathrm{where}\quad\hat{x}(L)=\frac{{\hat{\tau}_{\rm
NL}}(L)+\hat{N}_{\tau_{\text{NL}}}^{\mathrm{RD},(0)}(L)}{{\tau_{\text{NL}}}(L)+\hat{N}_{\tau_{\text{NL}}}^{\mathrm{RD},(0)}(L)},$
(38)

where $\hat{N}_{\tau_{\text{NL}}}^{\mathrm{RD},(0)}(L)=k_{L}\bar{N}_{L}^{{\rm
RD},(0)}/C_{L}^{\zeta_{*}}$ accounts for realization dependent Gaussian noise
fluctuations in the data.777An exception to this may be “TTxPP” estimators. In
this case, the resulting $N_{L}^{{\rm RD},(0)}$ spectra have a small
amplitude, which may be negative for some multipoles. For small
${\tau_{\text{NL}}}$ the logarithm in Eq. (38) will then become undefined.
Alternative approximations to the reconstruction noise should then instead be
used, e.g. $N^{(0)}_{\tau_{\text{NL}}}(L)\simeq\sqrt{(2L+1)f_{\rm
sky}M_{LL}/2}$ which is non-negative by definition. Although we reconstruct
the modulation field with general $L$ dependence, for constraining
${\tau_{\text{NL}}}$ we follow the standard definition by assuming no scale
dependence in ${\tau_{\text{NL}}}$ itself;
${\tau_{\text{NL}}}(L)={\tau_{\text{NL}}}$.

Once we have obtained our posterior distribution of ${\tau_{\text{NL}}}$
estimates, it is straightforward to construct constraints on
${\tau_{\text{NL}}}$ by numerically integrating the distribution from
${\tau_{\text{NL}}}=0$ up until some desired limit. Recall that
${\tau_{\text{NL}}}$ is positive definite by definition in the context of
squeezed primordial non-Gaussianity.

##  III Data and simulations

The nominal Planck 2013 results placed constraints on
${\tau_{\text{NL}}}$using the 143 and 217 GHz intensity maps to reconstruct
the modulation field. The foreground signal was partly mitigated by projecting
a dust template from higher frequencies, but the Planck 2013 analysis still
showed evidence of important residuals. There is also a large spatial
variation of the noise, which can mimic the effect of spatially varying power
from ${\tau_{\text{NL}}}$ if it cannot be simulated and corrected for
accurately. This noise modulation can be mitigated by correlating different
maps with independent noise; in Planck 2013 a large quadrupole was
substantially reduced by performing the joint (143, 217) GHz reconstruction.
The joint reconstruction is obtained by direct application of the estimators
discussed in section II.2, where one copy of the inverse filtered harmonics
(e.g., as given in Eq. (21)) is replaced with those corresponding to a
different map. More generally, the noise from observations at different times
and different detectors is also approximately uncorrelated, and we can also
use joint reconstruction between different data subsets to reduce sensitivity
to the spatially-varying noise power.

Advancing from the 2013 to 2018, the data products are derived from a longer
observation period via repeated scans of the sky, leading to a lower effective
noise in all frequency channels. In addition to this, the release of
polarization data enables us to expand the range of estimators at our
disposal; from temperature only (TT), to polarization (PP) and minimum
variance (MV) combinations, as defined in Section II.2. PP estimators
typically offer weak constraining power in comparison to TT and MV, since the
small-scale signal is dominated by noise for Planck. However, the MV estimate
(which combines information from polarization with temperature), by definition
reduces the reconstruction noise when compared to TT and PP (if implemented
optimally), and we use it for our best constraints. Other cross-power spectrum
estimators can be built by considering the cross between estimators of
different types, such as “TT$\times$PP”. Whilst this estimator has low
reconstruction noise due to the small noise cross-power spectrum between the
temperature and polarization fields, this comes at the expense of a larger
estimator variance such that the resulting spectra provide relatively weak
constraints.

To see the impact of the changes since the 2013 result, we first produce
constraints from updated 143 GHz and 217 GHz temperature-only data, to
illustrate improvements associated with lowering the effective noise, whilst
remaining subject to the same frequency-dependent signals as in 2013.
Secondly, we then use the component-separated SMICA maps in temperature and
polarization to reduce foreground contamination, using the joint
reconstruction between observations in first and second half of the
observational period (HM1, HM2) to mitigate noise modelling uncertainty.
Comparison between temperature reconstruction and the minimum-variance joint
constraint with polarization (MV) shows the impact of the polarization signal,
while also have different dependence on foreground residuals (which are
expected to be lower in polarization on small scales). Finally, we use the
latest NPIPE data and simulation set, which have (in most respects) higher-
fidelity noise-modelling and signal processing. We use the joint MV
reconstruction between NPIPE time-period data splits A and B to obtain our
baseline results.

In the remainder of this section, we provide further technical details
regarding the processing of the data and simulations that propagate into our
estimators, and finally the likelihood analysis.

###  III.1 Simulations and Corrections

The Full Focal Plane (FFP) simulations [43, 44] aim to replicate the
statistics of the observed CMB sky. The end products from the map-making
process are a set of lensed scalar realizations of the cosmological CMB along
with simulated noise + residual systematics. The CMB simulations rely on
fiducial power spectra from which scalar unlensed CMB realizations are
constructed. The unlensed CMBs are then lensed by Gaussian realizations of a
lensing potential that are independent between simulations. The input
cosmological parameters used in the process of generating and subsequently
lensing the scalar maps are based on Planck public release (PR) data. For
example, in the case of FFP10, which are the set of CMB and noise simulations
that accompany the Planck PR3 data, the CMB simulations assume a base
$\Lambda$CDM cosmology with parameters inferred from previous PR2 [45] data.
Note that the PR3 data are also referred to as the “DX12” data - the latter of
which usually appears in technical documentation - and we use these
interchangeably throughout this work. The simulated CMB signal contains
additional frequency-dependent sources. These include Rayleigh scattering, a
second order dipole-induced quadrupole [46], and a Doppler modulation and
(frequency-independent) aberration [47]. Once the lensed scalar maps are
generated, they are convolved with FEBeCoP [48], which accounts for realistic
beaming and instrumental pointing.

In order to generate the accompanying FFP10 noise maps a fiducial realization
of the scalar CMB is first constructed. The process further relies on
sophisticated foreground modelling, via the Planck Sky Model [49], which
contains galactic thermal dust, spinning dust, synchrotron, free-free
scattering and CO line emission based on data-based models, as well as random
realizations of CIB and SZ that are (unrealistically) uncorrelated to the
fiducial CMB lensing potential realization. The nominal input sky map in
conjunction with the fiducial CMB realization is then processed to include
systematic effects such as detector pointing, and known point sources [50].
From this point, independent realizations of noise are added to mimic the
observational data stream [51]. Finally, by subtracting out the fiducial CMB
and foregrounds, the noise (+ systematics) simulations are obtained.

In addition to the DX12 data and FFP10 simulations, in this work we also make
use of the PR4 NPIPE data [12] and its corresponding simulations. At the level
of the data, there are improvements in the polarization calibration procedure,
and additional data ($\sim 8\%$) from the (previously neglected) satellite
repointing periods, leading to a slightly reduced noise level. Enhanced data
processing further accounts for systematic effects by building templates in
the time-ordered data stream. These template amplitudes can be inferred via
marginalization over the sky, and lead to iteratively “cleaned” maps as the
template amplitudes converge to zero. With respect to the simulations, the
lensed scalar CMB remains unchanged from FFP10, but the end-to-end simulations
emulate many features from the data processing pipeline, such as the removal
of systematics in time ordered information previously discussed. Finally, the
NPIPE processing allows for the construction of maximally independent
detector-set split maps, which are referred to as the A and B splits. These
should in principle provide improved noise decorrelation with respect to one
another than that between the half mission splits.

Despite the sophisticated instrumentation model used in FFP10 and NPIPE, small
discrepancies between the power spectra of the simulations and those observed
in the equivalent data will persist. Small spatially varying differences in
power could bias our inference of primordial power modulation. Fortunately,
the realization-dependent disconnected noise correction ($N_{L}^{{\rm
RD},(0)}$) that we apply automatically corrects for the difference between the
true and simulation covariance to leading order, making the estimators
relatively robust. However, for higher-order differences to be small, we do
still need to ensure that the power in the simulations is at least
perturbatively close to the power in the data.

We therefore adjust the simulations so that the power spectra of the filtered
data and filtered simulations match, targeting consistency to the sub-percent
level. In order to do this, we compute the power difference from the WF data
and the ensemble average spectrum from WF simulations: $\Delta C_{\ell}\equiv
C_{\ell}^{\text{data}}-\langle C_{\ell}^{\text{sims}}\rangle$. The corrections
are Gaussian and isotropic realizations of power that are added to the
simulations and/or data. We defer explicit details of this prescription to
Appendix. A where we explain how we simultaneously correct the auto and cross-
spectra between datasets. We illustrate the effectiveness of this corrective
procedure for the DX12 143/217 GHz data and simulations in Figure 11.

##  IV Pipeline validation and calibration

We now validate our pipeline by simulating CMB skies modulated with a known
${\tau_{\text{NL}}}$ modulation field before adding noise. We check that we
recover the modulation spectrum correctly, and that the corresponding
posterior is consistent with the input value. Since the reconstruction on a
single realization is noisy, we check for biases over a set of simulations
with the same input ${\tau_{\text{NL}}}$. The mean of the reconstructed
modulation power spectra should match the input spectrum, and the product of
the posteriors should peak close to the input ${\tau_{\text{NL}}}$.

We perform these tests using 300 FFP10 HM1/HM2 and NPIPE A/B simulations on
the cut-sky, using the same masks as in the final data analysis. In each
simulation set, we use common scalar CMB simulations (from FFP10) modulated
with realizations of the modulation field $f(\hat{{\mbox{\boldmath{$n$}}}})$
with theoretical spectrum
$C_{L}^{f}={\tau_{\text{NL}}}C_{L}^{\zeta^{*}}=1500C_{L}^{\zeta^{*}}$. Key
differences in the simulation sets are found in their respective noise maps,
as described in Section III. After building the modulated simulations, we
compute the WF and IVF fields, and then construct ensembles of modulation
field estimators from map pairs as given in Eq. (25). In order to compute the
mean field and $N_{L}^{(0)}$ corrections we repeat this procedure for the same
CMB + noise simulations without the modulation signal. Each modulation
reconstruction is then mean field and $N_{L}^{{\rm MC},(0)}$ subtracted. We
use the same realizations of the noise and unmodulated lensed CMB for the
with-signal and without-signal simulations, so that the Monte Carlo and cosmic
variance noise largely cancels. This choice enables us to measure any required
calibration factors out to small scales, even where the signal is well below
the reconstruction noise.

Figure 3: Modulation power spectra from with-signal ${\tau_{\text{NL}}}$
simulation minimum variance reconstructions: Left: NPIPE (A, B), Right: FFP10
(HM1, HM2). Dotted line indicates input theory spectrum, and solid blue the
mean reconstructed spectrum; with shaded contours indicating the 99, 95 and
68% empirical C.L. The calibration factor is shown within the inset for each
simulation set.

With the ensemble of reconstructed spectra we compute the calibration factor
by comparing their mean against the theoretical input spectrum:
$k_{L}=C_{L}^{f,\mathrm{Theory}}/\langle\bar{C}_{L}^{f}\rangle_{\mathrm{sims}}$.
Since the calibration is empirical (and therefore depends on the details of
the simulations, mask and processing) it is not universally defined, and
requires expensive calculation. We therefore only compute the calibration for
the simulation sets discussed in this section; which include those relevant
for our baseline results. In Figure 3 we plot the mean over reconstructed
power spectra after $k_{L}$-correction, and in Figure 4 plot the product of
the corresponding posterior distributions (normalized by their maximum value).

The $k_{L}$ factors are plotted within the inset panels of Figure 4 for MV
reconstructions using the FFP10 and NPIPE sets. We further find that the TT
calibration factors are approximately equal to those given for MV. In order to
compute $k_{L}$ we first bin the raw empirical result. For the largest scales
of interest $L=1,2,3,4$ we retain unit multipole bins, before incrementally
expanding the bin width in multiples of three, weighting by an $L(L+1)$
kernel. The binned solutions are then interpolated to build smooth functions
over the smaller scales. This strategy is designed to avoid over-fitting
fluctuations in the Monte Carlo noise, whilst modelling the large-scales
accurately where the bulk of the modulation signal resides. At first pass we
found that the large-scale calibration leads to over-correction of the
posterior estimates, and we therefore decrease the amplitude of $k_{L}$ by 5%
on these scales.

In Figure 4 we plot the calibrated, uncalibrated and “idealized” posterior
distributions. In the idealized case, we use the (known) realizations of input
full-sky modulation power spectra in place of the reconstructions appearing in
the estimates of (38). This serves as a second reference point in addition to
the input ${\tau_{\text{NL}}}=1500$ marker, with which the product of
posteriors should be consistent. Notably this posterior does not peak at the
precise fiducial value, which is driven by the cosmic variance and the limited
sample size. After calibration, the product of reconstructed simulation
posteriors are highly compatible with the product of the idealized posteriors,
both of which agree with the theoretical value (${\tau_{\text{NL}}}=1500$) to
within $1\sigma$. We have further checked that the calibrated posteriors of
without-signal simulations are consistent with ${\tau_{\text{NL}}}=0$ to the
same degree.

Figure 4: Product of posterior distributions (normalized to one) from 300
with-signal ${\tau_{\text{NL}}}$ simulation minimum variance reconstructions:
Left: NPIPE (A, B); Right: FFP10 (HM1, HM2). Vertical dash-dotted lines
indicates the input simulation value. Solid black posterior assumes perfect
full-sky reconstruction of the input modulation power. Red and green lines
indicate the uncalibrated and calibrated reconstructions respectively. All
posterior distributions are constructed over the multipole range $L\in[1,50]$.

##  V Results

We first report our results for the TT modulation reconstruction obtained from
the 143 and 217 GHz channels, which are derived from DX12 data and FFP10
simulations. In Figure 5 we plot the uncalibrated power spectra for the auto-
and joint-reconstructions between these data. In red we show the modulation
field power spectrum after mean field and $N_{L}^{{\rm RD},(0)}$ subtraction,
the latter of which is given in green for reference. In orange bands, we plot
the 68%, 95% and 99% confidence intervals obtained directly from the
equivalent modulation power spectra of without-signal simulations. This figure
is directly comparable to that given in Figure 20 of the Planck 2013 non-
Gaussianity constraints [9]. Beyond obvious differences in the data and
simulations, we note that the methodology presented in this work is slightly
different. In Planck 2013 the disconnected noise bias was computed from
simulations alone ($N_{L}^{{\rm MC},(0)}$), which is not as robust as the
realization-dependent estimate ($N_{L}^{{\rm RD},(0)}$) used in this work.
Moreover in this work we adjust the Gaussian simulation power measured against
the data, which makes our ${N_{L}^{(0)}}$ and mean field subtraction more
optimal. When comparing Figure 5 to its equivalent in 2013, it should be
notated that “(143, 217) GHz” in this work equates to “143$\times$217 GHz” in
the other. In Planck 2013 the 857 GHz intensity map was projected from the 143
and 217 GHz intensity maps whilst obtaining the IVF maps. This projection
acted as a template for dusty foreground removal. We project 857 GHz in the
same way but now using the DX12 data.

Figure 5: Uncalibrated modulation reconstruction power spectra from DX12 143
and 217 GHz data. From left-to-right we plot the auto 143 GHz, joint (143,
217) GHz and auto 217 GHz reconstructions. Red and green lines indicate the
modulation and $\bar{N}_{L}^{{\rm RD},(0)}$ power spectra, and orange bands
indicate the 68%, 95% and 99% C.L. of the simulation distribution with no
signal.

With this in mind we note that there are differences between the reconstructed
spectra, most vividly between the amplitude and sign of the largest scale
reconstruction modes. In Planck 2013 there was an $\sim 4\times$ dipolar power
difference between the auto-reconstructions, with the power of the joint
reconstruction roughly matching the 143 GHz auto. We find qualitatively
similar results; though the difference is $\sim 2\times$ and the sign of our
143 GHz dipole is now negative. In both analyses, it is likely that the 217
GHz auto signal is strongly contaminated by non-cosmological signal which is
not present in the 143 GHz observations. The differences in the dipole sign
are partially explained by the realization dependent noise modelling, which
slightly boosts the amplitude of ${N_{L}^{(0)}}$ in this work resulting in a
more heavily subtracted signal. The octopole signal similarly has a frequency-
dependent amplitude, which in the joint reconstruction approximately match the
143 GHz auto.

The most striking difference is in the quadrupole. Whilst we report that there
is very little reconstruction signal in this mode across the auto- and joint-
reconstructions, in Planck 2013 this dominated the autos before being heavily
mitigated in the joint. The strong frequency-dependent amplitude seen in the
Planck 2013 analysis already indicated that the bulk of the signal was not
from primordial modulation. The difference we observe in this work is most
likely explained by better quality (and power-corrected) simulations, which
enable us to subtract the frequency-dependent components of the mean field
more reliably. Though varying by a lesser extent, we observe a reduced
amplitude in the 143 GHz auto, which could be explained in the same way.

We now consider the foreground cleaned SMICA maps, for which we compare DX12
and NPIPE data. It is helpful to consider both of these data (and
corresponding simulation) products since it enables us to identify possible
systematic effects arising from the map-making procedure that could lead to a
false detection. When comparing these data we perform modulation
reconstruction for a variety of different sky fractions, the masks of which
are shown in Figure 6. By analysing the response of the modulation
reconstruction to lesser sky coverage, we check for possible galactic
foreground contaminants. The baseline mask maintains $f_{\mathrm{sky}}\sim
67\%$, and is almost identical to that used in the Planck lensing constraints.
Additional masks are then obtained by taking a threshold cut in the 353 GHz
intensity emission (after smoothing the data to avoid removing regions of
high-intensity CMB). We vary the threshold in order to obtain masks with an
incremental 10% reduction in $f_{\mathrm{sky}}$.

Figure 6: Masks used to test foreground contamination. Black corresponds to
the fiducial mask used in the Planck lensing analysis, coloured masks layered
below are generated from a threshold cut in the smoothed intensity of emission
in 353 GHz. By construction, the larger masks are a superset of the smaller
ones, so that
$\mathrm{M37}\subseteq\mathrm{M47}\subseteq\mathrm{M57}\subseteq\mathrm{M67}$.

For NPIPE analyses we mask a small number of additional pixels, which
correspond to the pre-processing mask, applied during our implementation of
the SMICA map making process; this essentially leaves the sky fraction
unchanged. For DX12 analyses (when we use the HM splits) the data contains
missing pixels, therefore we apply further masking which results in a $1-3\%$
reduction in $f_{\mathrm{sky}}$. In the following text, these masks will be
referred to as M67, M57, M47 and M37, where any additional masking should be
interpreted from the data it is applied to and the number convention indicates
the sky fraction of the baseline masks in Figure 6.

Although NPIPE has more end-to-end simulations available than FFP10, to keep
these initial comparisons on equal footing we construct mean field and
$\bar{N}_{L}^{{\rm RD},(0)}$ corrections using the same number of simulations
(80 for mean field, 220 for $\bar{N}_{L}^{{\rm RD},(0)}$). In Figure 7 we
present the results for $\bar{C}_{L}^{f}$ computed from minimum variance
modulation reconstruction. In the top and bottom rows we plot the power
spectra of joint reconstructions NPIPE (A, B) and DX12 (HM1, HM2)
respectively, and left-to-right the reduced $f_{\rm sky}$. For a given
baseline mask, NPIPE and DX12 have similar spectra, with the amplitude of the
large-scale reconstruction modes of NPIPE being generically lower. All of the
spectra have a shape consistent with the (143, 217) GHz reconstruction, in
which most of the power is found in the octopole; which has decreasing power
as a function of $f_{\rm sky}$. This systematic decrease in the octopole power
suggests that there may be a galactic foreground contaminant in the modulation
signal, which is partially removed as the size of the galactic mask is
increased. However, as we decrease the available sky fraction the statistical
significance of the reconstructed signal decreases; as can be seen by the
error bars expanding left-to-right. We do not present the TT reconstruction
spectra since these are essentially the same as the MV, though the spectra and
errors have a slightly larger amplitude.

Figure 7: Uncalibrated modulation power spectra from NPIPE (A, B) and DX12
(HM1, HM2) joint-reconstructions. In each panel we plot the MV reconstruction,
and from left-to-right reduce the sky fraction; corresponding to M67, M57, M47
and M37 as shown in Figure 6. The 68%, 95% and 99% confidence intervals are
represented by the orange contours.

In the upper part of Table 1 we provide estimates for the uncalibrated
modulation power spectrum reconstructions. For each power spectrum we compute
an uncalibrated $\bar{\tau}_{\mathrm{NL}}$ estimate via direct application of
Eq. (35), and provide a rough estimate of the error
$\sigma(\bar{\tau}_{\mathrm{NL}})$ from applying the same estimator to
reconstructed modulation spectra of without-signal simulations and computing
the resulting standard deviation. For each spectrum we form an estimate with
$L_{\mathrm{min}}=1$ and consider $L_{\mathrm{max}}=2,10$ and $50$. We provide
results for both TT and MV reconstructions.

Uncalibrated inverse-variance weighted estimates:
$\bar{\tau}_{\mathrm{NL}}\pm\sigma(\bar{\tau}_{\mathrm{NL}})$ for $1\leq L\leq
L_{\mathrm{max}}$
---
| Temperature | Minimum variance
Data | Mask | $L_{\mathrm{max}}=2$ | $L_{\mathrm{max}}=10$ | $L_{\mathrm{max}}=50$ | $L_{\mathrm{max}}=2$ | $L_{\mathrm{max}}=10$ | $L_{\mathrm{max}}=50$
DX12 (143, 217) GHz | M67 | $-340\pm 267$ | $-56\pm 272$ | $-94\pm 274$ | N/A | N/A | N/A
DX12 (HM1, HM2) | M67 | $-260\pm 273$ | $95\pm 274$ | $50\pm 277$ | $-230\pm 269$ | $113\pm 267$ | $49\pm 269$
DX12 (HM1, HM2) | M57 | $-459\pm 301$ | $-154\pm 302$ | $-185\pm 306$ | $-430\pm 302$ | $-145\pm 294$ | $-200\pm 297$
DX12 (HM1, HM2) | M47 | $-492\pm 319$ | $-177\pm 330$ | $-208\pm 332$ | $-500\pm 323$ | $-241\pm 322$ | $-291\pm 323$
DX12 (HM1, HM2) | M37 | $-384\pm 346$ | $-166\pm 354$ | $-207\pm 357$ | $-425\pm 357$ | $-302\pm 357$ | $-369\pm 360$
NPIPE (A, B) | M67 | $-46\pm 260$ | $129\pm 259$ | $79\pm 261$ | $1\pm 248$ | $167\pm 248$ | $107\pm 248$
NPIPE (A, B) | M57 | $-306\pm 276$ | $-145\pm 274$ | $-185\pm 275$ | $-288\pm 267$ | $-134\pm 266$ | $-186\pm 266$
NPIPE (A, B) | M47 | $-313\pm 282$ | $-107\pm 286$ | $-134\pm 287$ | $-290\pm 273$ | $-125\pm 278$ | $-165\pm 279$
NPIPE (A, B) | M37 | $-170\pm 296$ | $-66\pm 299$ | $-120\pm 301$ | $-130\pm 285$ | $-103\pm 285$ | $-174\pm 285$
Calibrated inverse-variance weighted estimates: $\hat{\tau}_{\mathrm{NL}},\
\sigma(\hat{\tau}_{\mathrm{NL}})$ for $1\leq L\leq L_{\mathrm{max}}$
| Temperature | Minimum variance
Data | Mask | $L_{\mathrm{max}}=2$ | $L_{\mathrm{max}}=10$ | $L_{\mathrm{max}}=50$ | $L_{\mathrm{max}}=2$ | $L_{\mathrm{max}}=10$ | $L_{\mathrm{max}}=50$
DX12 HM (HM1, HM2) | M67 | $-321\pm 326$ | $12\pm 317$ | $-30\pm 319$ | $-287\pm 321$ | $39\pm 309$ | $-21\pm 310$
NPIPE AB (A, B) | M67 | $-62\pm 290$ | $86\pm 279$ | $36\pm 280$ | $-63\pm 293$ | $86\pm 283$ | $36\pm 283$
Table 1: Summary of the inverse-variance weighted ${\tau_{\text{NL}}}$
estimates. For each combination of data, mask and modulation estimator, we
compute $\bar{\tau}_{\mathrm{NL}}$ estimates from $L=1$ to $L_{\mathrm{max}}$
using Eq. (35). Each entry in the table(s) corresponds to the data estimate
and standard deviation of equivalent estimates from without-signal
simulations. Upper: Uncalibrated estimates from DX12 and NPIPE data. Lower:
Calibrated estimates.

When only $L\in[1,2]$ define the estimator we sometimes have large negative
estimates, which are driven directly by the negative mode(s) in the
reconstructed power spectra; these estimates are un-physical. For the
remaining estimators when $L_{\mathrm{max}}\in[10,50]$ we see that all of the
data estimates are consistent with ${\tau_{\text{NL}}}=0$. It should be noted
that using the standard deviation of $\bar{\tau}_{\mathrm{NL}}$ from
simulations is not the most rigorous approach to estimating the error, since
the distribution is slightly skewed by the $\chi^{2}$ properties of the
reconstruction modes. Our estimated error is however still useful as an
approximate means of assessing the data consistency with
${\tau_{\text{NL}}}=0$.

The error bars for the $\bar{\tau}_{\mathrm{NL}}$ constraints in Table 1
nearly saturate by $L_{\rm max}=10$, since the bulk of the signal is in the
lowest multipoles. However, the measured values mostly shift downward moving
from $L_{\rm max}=10$ to $L_{\rm max}=50$, by an amount that look anomalous
given the little information that should be available in that range (hence the
error bars being almost unchanged). This shift is dominated by the multipole
range $10\lesssim L\lesssim 40$, corresponding to the run of mostly negative
modulation power multipole estimates clearly seen by eye in e.g. Fig. 7.
However, we find consistency with the simulations to $\lesssim 3\sigma$ by
applying Eq. (35) to each of the sub-intervals $5+\delta L\leq L\leq 15+\delta
L$ for $\delta L\in[0,10,20,30]$. We have also found that NPIPE (A, B) is just
over $3\sigma$ low for $9\leq L\leq 35$. While visually striking, negative
signals cannot be produced by a physical modulation signal, and this multipole
range is chosen a posteriori, so we do not regard this is significant evidence
of a problem.

In addition to the data provided in the table, we report
$\bar{\tau}_{\mathrm{NL}}\pm\sigma(\bar{\tau}_{\mathrm{NL}})=$ $-150\pm 579$,
$173\pm 577$ and $156\pm 578$ for $L_{\mathrm{min}}=1$ with
$L_{\mathrm{max}}=2,10$ and $50$, from the TTxPP NPIPE (A, B) reconstruction
spectrum using the M67 mask. As discussed in Section III, estimators of this
type have a large variance, and this can be seen directly by comparing the
error on this result to the analogous case of MV NPIPE (A, B) given in Table
1, which is more than $2\times$ smaller.

In Figure 8 we show the corresponding posterior distribution of estimates for
a range of $L_{\mathrm{max}}$ with fixed $L_{\mathrm{min}}=1$. All of the
posterior distributions are consistent with ${\tau_{\text{NL}}}=0$ with very
long tails. In a few instances, the posterior peaks slightly away from zero,
but in most cases peaks at zero. The dashed vertical lines in each panel
indicate the upper bound on ${\tau_{\text{NL}}}$ inferred from the 95% C.L. on
the $L_{\mathrm{max}}=50$ posterior distribution. This bound shifts as a
function of $f_{\rm sky}$ but is relatively stable for the larger-sky
fractions. Note that these posteriors have no calibration factors of $k_{L}$,
but do contain the Monte Carlo and Hartlap corrections to their covariances as
described in Section II.4. Without applying the full calibration these
estimates are slightly too small.

Figure 8: Uncalibrated posterior distribution of ${\tau_{\text{NL}}}$
estimates, from MV reconstructed modulation power spectra from NPIPE (A, B)
and DX12 (HM1, HM2) data. The underlying data power spectra correspond to
those given in Figure 7. All posterior distributions are computed from $L=1$
to an $L_{\mathrm{max}}$, as indicated by the colour of the line. The single
dashed line in each panel indicates the corresponding 95% C.L.
${\tau_{\text{NL}}}$ constraint obtained by integrating the
$L_{\mathrm{max}}=50$ posterior distribution.

.

All of the uncalibrated posterior constraints are provided in the upper
portion of Table 2. We provide the constraints from using data up to
$L_{\mathrm{max}}=2,10,50$ in both temperature and polarization (except for
DX12 (143, 217) GHz). In general, the upper-bound decreases as more multipoles
are included in the likelihood analysis; with the tightest constraints arising
when we preserve the larger sky-fractions. This latter trend is associated
with the increased (co-)variance of modulation estimators as the sky fraction
is decreased, leading to broader posteriors. For M67 and M57, we find that the
95% upper bound lies approximately between 1500-2000.

Uncalibrated $\tau_{\mathrm{NL}}$ constraints: 95% C.L. on
$P(\tau_{\mathrm{NL}}|\bar{\tau}_{\mathrm{NL}})$ for $1\leq L\leq
L_{\mathrm{max}}$
---
| Temperature | Minimum variance
Data | Mask | $L_{\mathrm{max}}=2$ | $L_{\mathrm{max}}=10$ | $L_{\mathrm{max}}=50$ | $L_{\mathrm{max}}=2$ | $L_{\mathrm{max}}=10$ | $L_{\mathrm{max}}=50$
DX12 (143, 217) GHz | M67 | 3300 | 2600 | 2000 | N/A | N/A | N/A
DX12 (HM1, HM2) | M67 | 3400 | 3100 | 2200 | 3200 | 2600 | 1800
DX12 (HM1, HM2) | M57 | 4100 | 2300 | 1800 | 3500 | 1800 | 1300
DX12 (HM1, HM2) | M47 | 5200 | 4100 | 3100 | 5000 | 3000 | 1900
DX12 (HM1, HM2) | M37 | 8700 | 3900 | 3600 | 9800 | 3000 | 2500
NPIPE (A, B) | M67 | 3600 | 2200 | 2000 | 3200 | 2300 | 1600
NPIPE (A, B) | M57 | 4000 | 2200 | 1600 | 3400 | 2200 | 1600
NPIPE (A, B) | M47 | 3500 | 3300 | 2600 | 3400 | 2700 | 2100
NPIPE (A, B) | M37 | 6200 | 2600 | 2200 | 8600 | 2100 | 1500
Calibrated $\tau_{\mathrm{NL}}$ constraints: 95% C.L. on
$P(\tau_{\mathrm{NL}}|\hat{\tau}_{\mathrm{NL}})$ for $1\leq L\leq
L_{\mathrm{max}}$
| Temperature | Minimum variance
Data | Mask | $L_{\mathrm{max}}=2$ | $L_{\mathrm{max}}=10$ | $L_{\mathrm{max}}=50$ | $L_{\mathrm{max}}=2$ | $L_{\mathrm{max}}=10$ | $L_{\mathrm{max}}=50$
DX12 (HM1, HM2) | M67 | 16000 | 3900 | 2800 | 14000 | 3300 | 2300
NPIPE (A, B) | M67 | 5900 | 2200 | 1900 | 5800 | 2200 | 1700
Table 2: Summary of upper limits on ${\tau_{\text{NL}}}$ from posterior
probabilities. Upper: Uncalibrated data constraints from DX12 and NPIPE data;
Sep. refers to the separate frequency constraint from the (143, 217) GHz
cross-estimator. Lower: 95% C.L. ${\tau_{\text{NL}}}$ constraints from
calibrated analyses. The NPIPE (A, B) constraint in the lower table is our
baseline result.

We finally produce results for the calibrated analyses of TT and MV modulation
reconstruction using DX12 (HM1, HM2) and NPIPE (A, B) map pairs. For the
calibrated NPIPE analysis, we have now used all available simulations; 120 and
480 simulations for the mean field and noise bias estimates respectively. In
the lower part of Table 1 we provide the inverse-variance weighted estimates
for these data. For the $L_{\mathrm{max}}=2$ constraint, we see that the
estimate is pulled more negative compared to the uncalibrated result. This is
due to the dipole and quadrupole receiving multiplicative corrections with
$k_{L}>1$ (see inset on Figure 3) and the fact that the reconstructed spectra
are negative for these two modes. Once further multipoles are included, the
estimates become more consistent with zero. For all of the calibrated analyses
the error estimate becomes slightly larger; which reflects the high weighting
that the low multipoles receive for this particular form of the estimator and,
as previously noted, that $k_{L}>1$ on these scales.

In Table 2 we present the calibrated ${\tau_{\text{NL}}}$ constraints on the
posterior probability, which exhibit relatively large changes from their
respective uncalibrated constraints. This is especially notable from for DX12
95% C.L. constraints for $L_{\mathrm{max}}=2$. For these calibrated results,
the posterior peaks at zero, but has a long and shallow tail, which leads to
the relatively weak constraining power. As in the uncalibrated analyses, all
of the posterior distributions remain consistent with zero. After calibrating,
however, the posterior peaks are pulled slightly positive, such that the
normalized distributions exhibit greater shifts as a function of
$L_{\mathrm{max}}$. For our baseline result we advocate the calibrated MV
NPIPE (A, B) analysis, for which we report ${\tau_{\text{NL}}}<1700$ 95% C.L.
The full-calibrated result for the modulation power spectrum and the
corresponding posterior distributions are shown in Figure 9.

Figure 9: Left: Calibrated MV modulation power spectrum reconstruction for
NPIPE (A, B) analysis. Right: Corresponding posterior distribution of
${\tau_{\text{NL}}}$ estimates. This figure represents the key data for our
baseline constraints.

Finally, we attempt to gauge the impact improved data quality has on
$\sigma({\tau_{\text{NL}}})$ estimates using a semi-analytic approximation of
the reconstruction noise bias ${N_{L}^{(0)}}$. These ${N_{L}^{(0)}}$ can be
computed cheaply from the power spectra of the IVF maps of different data, and
we do so for the NPIPE, DX12 and Planck 2013 (PR1) data. Using these spectra
we simulate $10^{6}$ realizations of the Gaussian noise bias following a
$\chi^{2}$ distribution, and propagate the spectra into Eq. (35) to build
$\bar{\tau}_{\mathrm{NL}}$ estimates. The standard deviation from these
ensembles is then used to build approximate errors for each dataset, and by
further splitting the estimates into ($10^{3}$) subsets before computing
$\sigma({\tau_{\text{NL}}})$ we estimate the sample-variance from the scatter
on the errors. The results from this exercise are given in Figure 10 for
$L_{\mathrm{min}}=1$ and $L_{\mathrm{max}}=10$, and we see that these are
highly consistent with those computed in Table 1 for the full end-to-end
simulations.

Moving from PR1 to DX12 for (143, 217) GHz analyses, there is a large
improvement from the overall reduction in noise. However when moving from TT
to MV (relevant for NPIPE (A, B) and DX12 (HM1, HM2)) there is a relatively
small improvement, which highlights that at Planck polarization noise levels
the reconstruction signal is largely accounted for by the temperature modes.
The forecast suggests that there is $\sim 20\%$ reduction on the error
(assuming a null hypothesis) when using NPIPE rather than PR1 data.
Interestingly the errors for DX12 (143, 217) GHz seem to outperform those from
DX12 (HM1, HM2). In the analytic model, this is explained by the competing
factors of noise levels and sky fraction: Although the HM maps have slightly
lower noise properties, they require additional masking due to missing pixel
data, and this few percent difference in $f_{\rm sky}$ propagates directly
into the error estimate. It should further be noted that the component-
separated maps have a better control over galactic and astrophysical
contaminants, the details of which are not folded into this analytic
calculation. The component-separated maps are therefore expected to produce
cleaner estimates of primordial modulation despite the seemingly competitive
reconstruction noise bias of the separate frequency maps.

Figure 10: Semi-analytic error forecasts for $\bar{\tau}_{\mathrm{NL}}$
estimates with different data choices: Diamond markers correspond to MV and
the remaining circular markers to TT estimators.

##  VI Conclusions

Using the latest Planck data, we found the significantly tightened constraint
${\tau_{\text{NL}}}<1700$ 95% C.L., reducing the upper limit by a factor of
$\sim 1.7$ compared to previous results. A Universe with purely Gaussian
primordial fluctuations remains consistent with the data, though the upper
limit remains substantially larger than would be predicted in most multi-field
inflation models.

In addition to making use of new data scans and polarization, and using the
improved NPIPE processing, we also made several significant improvements to
the analysis methodology, including a more optimal and robust model for the
connected biases and use of foreground-separated maps. The quadratic estimator
pipeline is very similar to that for CMB lensing, but with the additional
requirement to consistently model the expected lensing four-point signal.
Several lessons from this work have already been directly useful to inform the
analogous updated NPIPE CMB lensing analysis (in prep.). The constraint is
still dominated by information from the temperature due to the relatively low
signal to noise in the Planck polarization, but could be significantly
improved in future with small-scale low-noise polarization data. More dramatic
improvements may be possible in conjunction with low shot noise large-scale
structure surveys and small-scale kinetic SZ observations [52].

The current analysis also could be further optimized by more optimal filtering
and handling of the noise anisotropy (e.g. see Ref. [53] in the context of
lensing), but given the signal is temperature dominated, the constraint is not
expected to improve dramatically. Future work could also easily extend the
analysis to scale-dependent ${\tau_{\text{NL}}}$-like signals.

## Acknowledgements

KM acknowledges support from STFC training grant ST/S505766/1. AL and JC
acknowledge support from the European Research Council under the European
Union’s Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement No.
[616170], and AL support by the UK STFC grants ST/P000525/1 and ST/T000473/1.
JC acknowledges support from a SNSF Eccellenza Professorial Fellowship (No.
186879). This research used resources of the National Energy Research
Scientific Computing Center (NERSC), a U.S. Department of Energy Office of
Science User Facility located at Lawrence Berkeley National Laboratory. KM is
grateful to Reijo Keskitalo for the provision of - and information pertaining
to - Planck simulation products available on NERSC. KM further thanks Maude Le
Jeune and Jean-Francois Cardoso for helpful discussion regarding SMICA
processed simulations. Some of the key results from this paper were derived
using the plancklens and healpy (healpix) [54, 55] software packages.

##  Appendix A Spectral corrections

Figure 11: Temperature spectral improvements $\Delta C_{\ell}^{TT}$, as
measured from the DX12 143/217 GHz data and simulations. The power spectra are
computed directly from the WF maps, and have been binned in multipole
intervals of $\Delta\ell=50$ and weighted by a $L(L+1)$ kernel. In each panel
the fractional difference of the fiducial and corrected maps is given in red
and green respectively. For reference, two grey shaded bands indicate the 1%
and 2% difference contours. Left: Auto-power spectra from 143 GHz . Centre:
Cross-power spectra from 143 GHz and 217 GHz. Right: Auto-power spectra from
217 GHz.

Accurate computation of the $N_{L}^{{\rm RD},(0)}$spectral correction
described in Section II.3 rely on the Gaussian power between the data and
simulations being perturbatively close. In order to ensure this, we add a
small amount of Gaussian power to the simulations and/or data. We quantify
this difference by measuring the power spectra of the WF data and the ensemble
averaged power spectra of WF simulation maps; $\Delta
C_{\ell}=C_{\ell,\mathrm{expt}}-\langle C_{\ell}\rangle_{\mathrm{sims}}\approx
0$. The WF maps (rather than the raw maps) provide a better benchmark for
assessing the power differences since it is the filtered maps that enter the
quadratic estimator. Moreover, in the case of e.g. the (143, 217) GHz analysis
where we perform a template projection of the 857 GHz intensity map during the
filtering stage, there is a reduction in power that would not otherwise be
measured on the raw (unfiltered) maps.

Since we consider the joint modulation reconstruction between map pairs (e.g.
see discussion towards the bottom of Section II.2) it is important that we not
only correct the Gaussian auto power spectrum within a given map set, but also
the cross power between different map sets entering the quadratic estimator.
By map sets, we are collectively referring to a given data and its
corresponding set of simulations. A concrete example is ensuring that $\Delta
C_{\ell}\approx 0$ for 143x143 GHz, 217x217 GHz and 143x217 GHz
simultaneously. We achieve this by solving the system of equations

$\displaystyle C_{\ell}[P_{\ell m}^{\mathrm{a,dat}},P_{\ell
m}^{\mathrm{b,dat}}]-\langle C_{\ell}[P_{\ell m}^{\mathrm{a,sim}},P_{\ell
m}^{\mathrm{b,sim}}]\rangle$ $\displaystyle=\Delta C_{\ell}^{\mathrm{ab}},$
(39) $\displaystyle C_{\ell}[P_{\ell m}^{\mathrm{a,dat}},P_{\ell
m}^{\mathrm{a,dat}}]-\langle C_{\ell}[P_{\ell m}^{\mathrm{a,sim}},P_{\ell
m}^{\mathrm{a,sim}}]\rangle$ $\displaystyle=\Delta C_{\ell}^{\mathrm{aa}},$
(40) $\displaystyle C_{\ell}[P_{\ell m}^{\mathrm{b,dat}},P_{\ell
m}^{\mathrm{b,dat}}]-\langle C_{\ell}[P_{\ell m}^{\mathrm{b,sim}},P_{\ell
m}^{\mathrm{b,sim}}]\rangle$ $\displaystyle=\Delta C_{\ell}^{\mathrm{bb}},$
(41)

where $P_{\ell m}\in\\{T_{\ell m}^{\mathrm{WF}},E_{\ell
m}^{\mathrm{WF}},B_{\ell m}^{\mathrm{WF}}\\}$ and ‘a’ and ‘b’ refer to
distinct map sets. We want to construct solutions such that $\Delta
C_{\ell}^{ab}=\Delta C_{\ell}^{aa}=\Delta C_{\ell}^{bb}=0$. To this end, we
assume the following form for the corrections in harmonic space

Map set a: $\displaystyle\begin{cases}P_{\ell
m}^{\mathrm{a,dat}}\rightarrow&\hat{P}_{\ell m}^{\mathrm{a,dat}}=P_{\ell
m}^{\mathrm{a,dat}}+\theta_{\ell m}^{\mathrm{a}}+\phi_{\ell m}\\\ P_{\ell
m}^{\mathrm{a,sim}}\rightarrow&\hat{P}_{\ell m}^{\mathrm{sim,a}}=P_{\ell
m}^{\mathrm{a,sim}}+\delta_{\ell m}^{\mathrm{a}}+\epsilon_{\ell m}\end{cases}$
(42) Map set b: $\displaystyle\begin{cases}P_{\ell
m}^{\mathrm{b,dat}}\rightarrow&\hat{P}_{\ell m}^{\mathrm{b,dat}}=P_{\ell
m}^{\mathrm{b,dat}}+\theta_{\ell m}^{\mathrm{b}}+\phi_{\ell m}\\\ P_{\ell
m}^{\mathrm{b,sim}}\rightarrow&\hat{P}_{\ell m}^{\mathrm{sim,b}}=P_{\ell
m}^{\mathrm{b,sim}}+\delta_{\ell m}^{\mathrm{b}}+\epsilon_{\ell m}\end{cases}$
(43)

where all of the additive corrections are only correlated with themselves.
This means that (by construction) $\phi_{\ell}$ and $\epsilon_{\ell}$ fix the
cross map set power spectra difference, and the remaining the auto map set
power difference.

Substituting the corrected harmonics into our ansatz Eq. (39), the power
difference can be written in terms of the isotropic power spectrum of the
harmonic corrections

$\displaystyle\Delta C_{\ell}^{\mathrm{ab}}$
$\displaystyle=(C_{\ell}^{\epsilon}-C_{\ell}^{\phi}),$ (44)
$\displaystyle\Delta C_{\ell}^{\mathrm{aa}}$
$\displaystyle=(C_{\ell}^{\epsilon}-C_{\ell}^{\phi})+(C_{\ell}^{\delta,\mathrm{a}}-C_{\ell}^{\theta,\mathrm{a}}),$
(45) $\displaystyle\Delta C_{\ell}^{\mathrm{bb}}$
$\displaystyle=(C_{\ell}^{\epsilon}-C_{\ell}^{\phi})+(C_{\ell}^{\delta,\mathrm{b}}-C_{\ell}^{\theta,\mathrm{b}}).$
(46)

We can obtain solutions in terms of the combinations given in parentheses (…)
by treating them as single variables. The solutions can then be unpacked
further by requiring power spectra to be non-negative via a piecewise
definition, e.g.

$\displaystyle C_{\ell}^{\phi}$
$\displaystyle=\begin{dcases}|(C_{\ell}^{\epsilon}-C_{\ell}^{\phi})|\quad&\text{if}\quad(C_{\ell}^{\epsilon}-C_{\ell}^{\phi})<0\\\
\quad 0\quad&\text{if}\quad(C_{\ell}^{\epsilon}-C_{\ell}^{\phi})\geq 0\\\
\end{dcases},$ (47) $\displaystyle C_{\ell}^{\epsilon}$
$\displaystyle=\begin{dcases}\phantom{|}(C_{\ell}^{\epsilon}-C_{\ell}^{\phi})\phantom{|}\quad&\text{if}\quad(C_{\ell}^{\epsilon}-C_{\ell}^{\phi})\geq
0\\\ \quad 0\quad&\text{if}\quad(C_{\ell}^{\epsilon}-C_{\ell}^{\phi})<0\\\
\end{dcases}.$ (48)

The remaining components $C_{\ell}^{\delta,\mathrm{a}}$,
$C_{\ell}^{\theta,\mathrm{a}}$, $C_{\ell}^{\delta,\mathrm{b}}$ and
$C_{\ell}^{\theta,\mathrm{b}}$ are obtained in the same way. Since the
solutions for the corrective power spectra depend on the empirical power
difference, they are in general not smooth functions of $\ell$. We therefore
bin the solutions and interpolate over $\ell$ in order to construct smooth
continuous functions. 888A bin width of $\Delta\ell=51$ is used throughout.
Note that to minimize boundary effects incurred during binning and
interpolation, we use the multipole range $30\leq\ell\leq 2048$, then retain
multipoles $100\leq\ell\leq 2018$ for final analyses.

Recall that these spectral corrections are built from the power spectra of WF
maps. We want the equivalent result for unfiltered maps on the full sky, such
that we can simulate and apply corrections to the raw data and simulations
before reprocessing. The required transformation is (approximately) given by

$C_{\ell}^{\mathrm{Raw}}\simeq
C_{\ell}^{\mathrm{WF}}\left(1+\frac{N_{\ell}}{C_{\ell}^{\mathrm{Fid.}}}\right)^{2}\frac{\mathcal{B}_{\ell}^{2}}{f_{\rm
sky}},$ (49)

where $\mathcal{B}_{\ell}$ is the beam transfer function and
$N_{\ell}=\Delta^{2}\mathcal{B}_{\ell}^{-2}$ is the approximate isotropic
noise power spectrum; for a fixed noise level $\Delta$. For uncorrelated
corrections, we can simulate the corresponding isotropic map corrections
directly from Eq. (49) using healpix. For correlated components, we should
instead simulate harmonic coefficients from the WF spectra, then transform
these via an appropriate definition of the instrumental beam and noise (if
they differ) for each map set;

$a_{\ell m}^{\mathrm{Raw}}\simeq a_{\ell
m}^{\mathrm{WF}}\left(1+\frac{N_{\ell}}{C_{\ell}^{\mathrm{Fid.}}}\right)\frac{\mathcal{B}_{\ell}}{\sqrt{f\mathrm{sky}}},$
(50)

before finally transforming the harmonics into real space maps.

## References

  * Maldacena [2003] J. M. Maldacena, Non-Gaussian features of primordial fluctuations in single field inflationary models, JHEP 05, 013, arXiv:astro-ph/0210603 .
  * Creminelli and Zaldarriaga [2004] P. Creminelli and M. Zaldarriaga, CMB 3-point functions generated by non-linearities at recombination, Phys. Rev. D 70, 083532 (2004), arXiv:astro-ph/0405428 .
  * Senatore and Zaldarriaga [2013] L. Senatore and M. Zaldarriaga, On Loops in Inflation II: IR Effects in Single Clock Inflation, JHEP 01, 109, arXiv:1203.6354 [hep-th] .
  * Lewis [2011] A. Lewis, The real shape of non-Gaussianities, JCAP 1110, 026 (2011), arXiv:1107.5431 [astro-ph.CO] .
  * Byrnes and Choi [2010] C. T. Byrnes and K.-Y. Choi, Review of local non-Gaussianity from multi-field inflation, Adv. Astron. 2010, 724525 (2010), arXiv:1002.3110 [astro-ph.CO] .
  * Bartolo  _et al._ [2004] N. Bartolo, E. Komatsu, S. Matarrese, and A. Riotto, Non-Gaussianity from inflation: Theory and observations, Phys. Rept. 402, 103 (2004), arXiv:astro-ph/0406398 .
  * Suyama and Yamaguchi [2008] T. Suyama and M. Yamaguchi, Non-Gaussianity in the modulated reheating scenario, Phys. Rev. D 77, 023505 (2008), arXiv:0709.2545 [astro-ph] .
  * Akrami  _et al._ [2020a] Y. Akrami _et al._ (Planck), Planck 2018 results. IX. Constraints on primordial non-Gaussianity, Astron. Astrophys. 641, A9 (2020a), arXiv:1905.05697 [astro-ph.CO] .
  * Ade  _et al._ [2014a] P. A. R. Ade _et al._ (Planck), Planck 2013 Results. XXIV. Constraints on primordial non-Gaussianity, Astron. Astrophys. 571, A24 (2014a), arXiv:1303.5084 [astro-ph.CO] .
  * Feng  _et al._ [2015] C. Feng, A. Cooray, J. Smidt, J. O’Bryan, B. Keating, and D. Regan, Planck Trispectrum Constraints on Primordial Non-Gaussianity at Cubic Order, Phys. Rev. D92, 043509 (2015), arXiv:1502.00585 [astro-ph.CO] .
  * Aghanim  _et al._ [2020a] N. Aghanim _et al._ (Planck), Planck 2018 results. I. Overview and the cosmological legacy of Planck, Astron. Astrophys. 641, A1 (2020a), arXiv:1807.06205 [astro-ph.CO] .
  * Akrami  _et al._ [2020b] Y. Akrami _et al._ (Planck), $Planck$ intermediate results. LVII. Joint Planck LFI and HFI data processing, Astron. Astrophys. 643, A42 (2020b), arXiv:2007.04997 [astro-ph.CO] .
  * Kogo and Komatsu [2006] N. Kogo and E. Komatsu, Angular Trispectrum of CMB Temperature Anisotropy from Primordial Non-Gaussianity with the Full Radiation Transfer Function, Phys. Rev. D 73, 083007 (2006), arXiv:astro-ph/0602099 .
  * Pearson  _et al._ [2012] R. Pearson, A. Lewis, and D. Regan, CMB lensing and primordial squeezed non-Gaussianity, JCAP 1203, 011 (2012), arXiv:1201.1010 [astro-ph.CO] .
  * Hanson and Lewis [2009] D. Hanson and A. Lewis, Estimators for CMB Statistical Anisotropy, Phys. Rev. D 80, 063004 (2009), arXiv:0908.0963 [astro-ph.CO] .
  * Hansen  _et al._ [2004] F. K. Hansen, A. J. Banday, and K. M. Gorski, Testing the cosmological principle of isotropy: Local power spectrum estimates of the WMAP data, Mon. Not. Roy. Astron. Soc. 354, 641 (2004), arXiv:astro-ph/0404206 .
  * Contreras  _et al._ [2018] D. Contreras, J. Hutchinson, A. Moss, D. Scott, and J. P. Zibin, Closing in on the large-scale CMB power asymmetry, Phys. Rev. D 97, 063504 (2018), arXiv:1709.10134 [astro-ph.CO] .
  * Akrami  _et al._ [2020c] Y. Akrami _et al._ (Planck), Planck 2018 results. VII. Isotropy and Statistics of the CMB, Astron. Astrophys. 641, A7 (2020c), arXiv:1906.02552 [astro-ph.CO] .
  * Byrnes  _et al._ [2010] C. T. Byrnes, M. Gerstenlauer, S. Nurmi, G. Tasinato, and D. Wands, Scale-dependent non-Gaussianity probes inflationary physics, JCAP 1010, 004 (2010), arXiv:1007.4277 [astro-ph.CO] .
  * Adhikari  _et al._ [2018] S. Adhikari, A.-S. Deutsch, and S. Shandera, Statistical anisotropies in temperature and polarization fluctuations from a scale-dependent trispectrum, Phys. Rev. D98, 023520 (2018), arXiv:1805.00037 [astro-ph.CO] .
  * Dvorkin and Smith [2009] C. Dvorkin and K. M. Smith, Reconstructing Patchy Reionization from the Cosmic Microwave Background, Phys. Rev. D 79, 043003 (2009), arXiv:0812.1566 [astro-ph] .
  * Dvorkin  _et al._ [2009] C. Dvorkin, W. Hu, and K. M. Smith, B-mode CMB Polarization from Patchy Screening during Reionization, Phys. Rev. D 79, 107302 (2009), arXiv:0902.4413 [astro-ph.CO] .
  * Gluscevic  _et al._ [2013] V. Gluscevic, M. Kamionkowski, and D. Hanson, Patchy Screening of the Cosmic Microwave Background by Inhomogeneous Reionization, Phys. Rev. D 87, 047303 (2013), arXiv:1210.5507 [astro-ph.CO] .
  * Namikawa [2018] T. Namikawa, Constraints on Patchy Reionization from Planck CMB Temperature Trispectrum, Phys. Rev. D 97, 063505 (2018), arXiv:1711.00058 [astro-ph.CO] .
  * Namikawa  _et al._ [2021] T. Namikawa, A. Roy, B. D. Sherwin, N. Battaglia, and D. N. Spergel, Constraining reionization with the first measurement of the cross-correlation between the CMB optical-depth fluctuations and the Compton y-map, Phys. Rev. D 104, 063514 (2021), arXiv:2102.00975 [astro-ph.CO] .
  * Roy  _et al._ [2022] A. Roy, A. van Engelen, V. Gluscevic, and N. Battaglia, Probing the circumgalactic medium with CMB polarization statistical anisotropy, (2022), arXiv:2201.05076 [astro-ph.CO] .
  * Feng and Holder [2019] C. Feng and G. Holder, Searching for patchy reionization from cosmic microwave background with hybrid quadratic estimators, Phys. Rev. D 99, 123502 (2019), arXiv:1808.01592 [astro-ph.CO] .
  * Fidler  _et al._ [2019] C. Fidler, J. Lesgourgues, and C. Ringeval, Lensing anomalies from the epoch of reionisation, JCAP 10, 042, arXiv:1906.05042 [astro-ph.CO] .
  * Byrnes  _et al._ [2006] C. T. Byrnes, M. Sasaki, and D. Wands, The primordial trispectrum from inflation, Phys. Rev. D 74, 123519 (2006), arXiv:astro-ph/0611075 .
  * Smidt  _et al._ [2010] J. Smidt, A. Amblard, C. T. Byrnes, A. Cooray, A. Heavens, _et al._ , CMB Constraints on Primordial non-Gaussianity from the Bispectrum ($f_{NL}$) and Trispectrum ($g_{NL}$ and $\tau_{NL}$) and a New Consistency Test of Single-Field Inflation, Phys. Rev. D 81, 123007 (2010), arXiv:1004.1409 [astro-ph.CO] .
  * Okamoto and Hu [2002] T. Okamoto and W. Hu, The angular trispectra of CMB temperature and polarization, Phys. Rev. D66, 063008 (2002), arXiv:astro-ph/0206155 [astro-ph] .
  * Hirata and Seljak [2003] C. M. Hirata and U. Seljak, Analyzing weak lensing of the cosmic microwave background using the likelihood function, Phys. Rev. D 67, 043001 (2003), arXiv:astro-ph/0209489 [astro-ph] .
  * Smith  _et al._ [2007] K. M. Smith, O. Zahn, and O. Dore, Detection of Gravitational Lensing in the Cosmic Microwave Background, Phys. Rev. D 76, 043510 (2007), arXiv:0705.3980 [astro-ph] .
  * Aghanim  _et al._ [2020b] N. Aghanim _et al._ (Planck), Planck 2018 results. VIII. Gravitational lensing, Astron. Astrophys. 641, A8 (2020b), arXiv:1807.06210 [astro-ph.CO] .
  * Maniyar  _et al._ [2021] A. S. Maniyar, Y. Ali-Haïmoud, J. Carron, A. Lewis, and M. S. Madhavacheril, Quadratic estimators for CMB weak lensing, Phys. Rev. D 103, 083524 (2021), arXiv:2101.12193 [astro-ph.CO] .
  * Regan  _et al._ [2010] D. M. Regan, E. P. S. Shellard, and J. R. Fergusson, General CMB and Primordial Trispectrum Estimation, Phys. Rev. D 82, 023520 (2010), arXiv:1004.2915 [astro-ph.CO] .
  * Hanson  _et al._ [2011] D. Hanson, A. Challinor, G. Efstathiou, and P. Bielewicz, CMB temperature lensing power reconstruction, Phys. Rev. D 83, 043005 (2011), arXiv:1008.4403 [astro-ph.CO] .
  * Planck Collaboration XV [2016] Planck Collaboration XV (Planck), Planck 2015 results. XV. Gravitational lensing, A&A 594, A15 (2016), arXiv:1502.01591 [astro-ph.CO] .
  * Namikawa  _et al._ [2013] T. Namikawa, D. Hanson, and R. Takahashi, Bias-hardened CMB lensing, MNRAS 431, 609 (2013), arXiv:1209.0091 [astro-ph.CO] .
  * Ade  _et al._ [2014b] P. A. R. Ade _et al._ (Planck), Planck 2013 Results. XXIV. Constraints on primordial non-Gaussianity, Astron. Astrophys. 571, A24 (2014b), arXiv:1303.5084 [astro-ph.CO] .
  * Hamimeche and Lewis [2008] S. Hamimeche and A. Lewis, Likelihood Analysis of CMB Temperature and Polarization Power Spectra, Phys. Rev. D 77, 103013 (2008), arXiv:0801.0554 [astro-ph] .
  * Hartlap  _et al._ [2007] J. Hartlap, P. Simon, and P. Schneider, Why your model parameter confidences might be too optimistic. Unbiased estimation of the inverse covariance matrix, A&A 464, 399 (2007), astro-ph/0608064 .
  * Ade  _et al._ [2016] P. A. R. Ade _et al._ (Planck), Planck 2015 results. XII. Full Focal Plane simulations, Astron. Astrophys. 594, A12 (2016), arXiv:1509.06348 [astro-ph.CO] .
  * Aghanim  _et al._ [2020c] N. Aghanim _et al._ (Planck), Planck 2018 results. III. High Frequency Instrument data processing and frequency maps, Astron. Astrophys. 641, A3 (2020c), arXiv:1807.06207 [astro-ph.CO] .
  * Planck Collaboration XIII [2016] Planck Collaboration XIII (Planck), Planck 2015 results. XIII. Cosmological parameters, A&A 594, A13 (2016), arXiv:1502.01589 [astro-ph.CO] .
  * Kamionkowski and Knox [2003] M. Kamionkowski and L. Knox, Aspects of the cosmic microwave background dipole, Phys. Rev. D 67, 063001 (2003), arXiv:astro-ph/0210165 .
  * Aghanim  _et al._ [2014] N. Aghanim _et al._ (Planck), Planck 2013 results. XXVII. Doppler boosting of the CMB: Eppur si muove, Astron. Astrophys. 571, A27 (2014), arXiv:1303.5087 [astro-ph.CO] .
  * Mitra  _et al._ [2011] S. Mitra, G. Rocha, K. M. Gorski, K. M. Huffenberger, H. K. Eriksen, M. A. J. Ashdown, and C. R. Lawrence, Fast Pixel Space Convolution for CMB Surveys with Asymmetric Beams and Complex Scan Strategies: FEBeCoP, Astrophys. J. Suppl. 193, 5 (2011), arXiv:1005.1929 [astro-ph.CO] .
  * Delabrouille  _et al._ [2013] J. Delabrouille, M. Betoule, J.-B. Melin, M.-A. Miville-Deschênes, J. Gonzalez-Nuevo, M. Le Jeune, G. Castex, G. de Zotti, S. Basak, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, J.-P. Bernard, F. R. Bouchet, D. L. Clements, A. da Silva, C. Dickinson, F. Dodu, K. Dolag, F. Elsner, L. Fauvet, G. Faÿ, G. Giardino, S. Leach, J. Lesgourgues, M. Liguori, J. F. Macías-Pérez, M. Massardi, S. Matarrese, P. Mazzotta, L. Montier, S. Mottet, R. Paladini, B. Partridge, R. Piffaretti, G. Prezeau, S. Prunet, S. Ricciardi, M. Roman, B. Schaefer, and L. Toffolatti, The pre-launchplancksky model: a model of sky emission at submillimetre to centimetre wavelengths, Astronomy & Astrophysics 553, A96 (2013).
  * Reinecke  _et al._ [2006] M. Reinecke, K. Dolag, R. Hell, M. Bartelmann, and T. Ensslin, A simulation pipeline for the planck mission, Astron. Astrophys. 445, 373 (2006), arXiv:astro-ph/0508522 .
  * Aghanim  _et al._ [2020d] N. Aghanim _et al._ (Planck), Planck 2018 results. III. High Frequency Instrument data processing and frequency maps, Astron. Astrophys. 641, A3 (2020d), arXiv:1807.06207 [astro-ph.CO] .
  * Kumar  _et al._ [2022] N. A. Kumar, G. Sato-Polito, M. Kamionkowski, and S. C. Hotinli, Primordial trispectrum from kSZ tomography, (2022), arXiv:2205.03423 [astro-ph.CO] .
  * Mirmelstein  _et al._ [2019] M. Mirmelstein, J. Carron, and A. Lewis, Optimal filtering for CMB lensing reconstruction, Phys. Rev. D100, 123509 (2019), arXiv:1909.02653 [astro-ph.CO] .
  * Zonca  _et al._ [2019] 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, 1298 (2019).
  * Górski  _et al._ [2005] K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann, HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere, Astrophys. J. 622, 759 (2005), arXiv:astro-ph/0409513 .